NOISE IN NONLINEAR MICRO-RESONATORS
By
Nicholas Miller

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
DOCTOR OF PHILOSOPHY
Mechanical Engineering
Physics
2012

ABSTRACT
NOISE IN NONLINEAR MICRO-RESONATORS
By
Nicholas Miller
In this work we examine several noise-induced phenomena that are found in nanoand micro-electromechanical systems. Our aim is to exploit these phenomena, where
possible, for system identification and sensing applications. The exception being our
discussion of nonlinear oscillators, where our goal is to remove the effect of noise
as much as possible. In general, however, our viewpoint is one of embracing the
noise and nonlinearity that are an inevitable part of miniaturization and seek new
modes of micro- and nano-system operation applicable where noise and nonlinearity
are fundamentally unavoidable.
In this work we discuss three different, but related, topics. First, we discuss the
dynamics of a system subject to a parameter sweep through a bifurcation in the presence of noise. We develop approximate distributions for the time at which a system
swept through a subcritical pitchfork or a saddle-node bifurcation will leave the local
region of phase space, leading to a sudden jump in response amplitude. These distributions are developed in the limit of “fast” sweeping, where the sweep rate is large
compared to the noise strength, which is relevant to most MEMS applications. We
then present the results of sweeping experiments performed on a MEMS resonator
and show the value of our analytic predictions. Our primary conclusion is that delayed bifurcation is prominent in MEMS devices, and we provide predictive tools for
predicting the resulting distribution of jump events.
Next, we consider a novel dynamical bridge sensing paradigm. This sensing
methodology employs a vibrating bistable system arranged, by tuning of parameters,
such that the rate of noise-activated escape out of the two stable basins of attraction
are identical. The ratio of occupation probabilities can become extremely sensitive at

this point when the noise is weak. We calculate the sensitivity and detection time of
the balanced dynamical bridge for use as a general use detector and we develop the
conditions required and demonstrate its application as a detector of non-Gaussian
noise. We illustrate the measurement of the parameters of a shot-noise process using
a one-dimensional bridge model.
Finally, we discuss phase noise in oscillators employing nonlinear frequency selective elements. We employ the method of averaging to develop an expression for the
spectrum of fractional frequency fluctuations in an oscillator. The expression for the
spectrum of frequency fluctuations is quite general, although it neglects dynamics in
the feedback elements of the oscillator. The expression contains elements of the noise
model and the nominal limit cycle shape of the oscillator. We demonstrate the utility
of the results by exploring the parameter space of a prototypical oscillator modeled
by a biased Duffing equation. We find that one can reduce the phase diffusion constant of the oscillator by tuning the system to a zero-dispersion point of the resonator
element, thus eliminating the action-angle coupling at zero-dispersion points of the
resonator frequency can lower the phase diffusion constant.

ACKNOWLEDGMENT
First and foremost, this work could not have been accomplished without the efforts of
my advisors: Professor Steven W. Shaw and Professor Mark I. Dykman. Beyond the
technical guidance and instruction Professors Shaw and Dykman have been mentors
and exemplars of scientific integrity. I am honored to have had the opportunity to
work with them both. The example set by these gentlemen will continue to guide
me for the remainder of my career. Thanks also to the committee who oversaw this
work; Professor Brian Feeny, Professor Norman Birge, and Professor Jongeun Choi.
I have had the good fortune to collaborate with Professor Kimberly L. Turner and
Chris Burgner of UCSB. I would like to thank Chris for his creativity, grounding
practicality, and hard work. Collaborating with him has been very educational and
enjoyable. I am indebted to my fellow lab-mates and former lab-mates for many
meaningful discussions about various interesting academic problems. These include
Jeff Rhoads, Umar Farooq, Ryan Monroe, Brendan Vidmar, Venkat Ramakrishnan,
and Pavel Polunin. I would especially like to thank Scott Strachen who is always
ready to discuss some sticky point of condensed matter physics or philosophy or
anything else. Our conversations have been particularly useful and enjoyable. Finally,
I would like to thank the funding agencies that made this work possible. The National
Science Foundation under grants CMMI-0758419 and CMMI-0900666 and the Defense
Advanced Research Agency under the DARPA-MTO-DEFYS project.

iv

TABLE OF CONTENTS

List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1 Introduction . . . . . . . .
1.1 MEMS and NEMS . . . .
1.2 Noise and Nonlinearity . .
1.3 Noise-Induced Phenomena
1.4 Summary of Results . . .

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

.
.
.
.
.

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

.
.
.
.
.

.
.
.
.
.

1
2
4
6
16

2 Escape Statistics for Parameter Sweeps Through Bifurcations
2.1 Stochastic Normal Forms . . . . . . . . . . . . . . . . . . . . . . .
2.2 Kramers Escape Rate . . . . . . . . . . . . . . . . . . . . . . . . .
2.3 Escape Near a Saddle-Node Bifurcation . . . . . . . . . . . . . . .
2.4 Escape Near a Subcritical Pitchfork Bifurcation . . . . . . . . . .
2.5 Subcritical Pitchfork Bifurcation in a MEMS Resonator . . . . . .
2.6 Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.

.
.
.
.
.
.
.

19
21
26
30
39
47
54

3 The
3.1
3.2
3.3
3.4
3.5
3.6

56
59
61
64
67
71
76

Dynamical Bridge . . . . . .
Activated Escape . . . . . . . . .
Dynamical Bridge Sensitivity . .
Detection of Non-Gaussian Noise
The Duffing Resonator Bridge . .
Shot Noise Measurement . . . . .
Outlook . . . . . . . . . . . . . .

.
.
.
.
.

.
.
.
.
.
.
.

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

.
.
.
.
.
.

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

4 Phase Noise in Nonlinear Oscillators .
4.1 Action-Angle Coordinates . . . . . . .
4.2 Microscopic Resonator Model . . . . .
4.3 Oscillator Frequency Fluctuations . . .
4.4 The Biased Duffing Oscillator . . . . .
4.5 Outloook . . . . . . . . . . . . . . . . .

.
.
.
.
.
.

.
.
.
.
.

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

.
.
.
.
.

.
.
.
.
.

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

.
.
.
.
.

.
.
.
.
.

vii

.
.
.
.
.
.
.

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

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.

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

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.

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

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.

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

.
.
.
.
.
.

.
.
.
.
.
.

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

.
.
.
.
.
.

.
.
.
.
.
.

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

.
.
.
.
.
.

. 78
. 81
. 84
. 92
. 100
. 109

5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
A Stochastic Processes . . . . . . . . . . . . .
A.1 The Kinetic Equation for a Stochastic Process .
A.2 The Langevin Equation . . . . . . . . . . . . . .
A.3 Ito and Stratonovich Integrals . . . . . . . . . .
v

.
.
.
.

.
.
.
.

. .
. . .
. . .
. . .

.
.
.
.

.
.
.
.

. .
. . .
. . .
. . .

.
.
.
.

.
.
.
.

114
114
118
121

B Stochastic Averaging . . . . . . . . . . . .
B.1 A Near-Identity Coordinate Transformation . .
B.2 Stochastic Averaging Part I: O( 1/2 ) . . . . . .
B.3 Stochastic Averaging Part II: O( ) . . . . . . .
B.4 Time Scaling . . . . . . . . . . . . . . . . . . .
B.5 Approximation by a Markov Process . . . . . .
C Derivation of Equation (2.65)
Bibliography

.
.
.
.
.
.

.
.
.
.
.
.

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

.
.
.
.
.
.

.
.
.
.
.
.

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

.
.
.
.
.
.

.
.
.
.
.
.

126
128
129
132
134
136

. . . . . . . . . . . . . . . . . . 140

. . . . . . . . . . . . . . . . . . . . . . . . . . 150

vi

LIST OF FIGURES

1.1

Activated escape in one dimension. . . . . . . . . . . . . . . . . . . .

6

1.2

Most probable paths (dashed lines) of activated escape in a two dimensional system. Invariant stable and unstable manifolds (solid lines)
represent underlying noise free dynamics. . . . . . . . . . . . . . . . .

10

2.1

Example potential governing normal form dynamics.

. . . . . . . . .

26

2.2

Illustration of the evolution of the potential (a), and the standard
bifurcation diagram (b), for a saddle-node bifurcation. . . . . . . . . .

30

c∞ boundary for linear sweeping, µ(t) = µ0 + rt, with µ0 < 0. The
first two branches are shown along with the variance (shaded region)
of c for 2 = 0.05. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

35

2.3

2.4

Mean (a), variance (b), and sample distributions (c) of the escape time
away from a swept saddle-node bifurcation. Numeric solution (solid
line) and approximations (dashed) for fast and slow sweeping are shown. 38

2.5

Illustration of the evolution of the potential (a), and the standard
bifurcation diagram (b), for a subcritical pitchfork bifurcation. . . . .

39

Mean (a), variance (b), skewness (c), and example distributions (d) of
the distribution of escape times away from a swept subcritical pitchfork bifurcation. Numeric solution (solid lines) and approximations
(dashed) for fast and slow sweeping are shown. Dots represent data
points from Monte-Carlo simulations of the nonlinear Mathieu equation, equations (C.17) and (C.17), swept through the subcritical pitchfork bifurcation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

46

2.6

vii

2.7

(a) SEM of MEMS gyroscope [89]. (b) Amplitude and phase of drive
axis vibration during an example escape event under parameter sweep.
For interpretation of the references to color in this and all other figures,
the reader is referred to the electronic version of this dissertation. . .

48

(Top) Frequency response for the MEMS resonator. Blue data points
track the sweep up and red data points track the sweep down. The
green data points represent the mean jump frequency for various sweep
rates. Black lines represent fitted stable (solid) and unstable (dashed)
frequency response branches. (Bottom) Real and imaginary parts of
the eigenvalues belonging to the fixed point at zero. Dashed line represents the linear sweep approximation. . . . . . . . . . . . . . . . . .

51

Experimentally measured distributions of escape events for ambient
noise and sweep rates ranging from 0.005 to 0.3 Hz/sec. The inset
shows one such distribution, artifically shifted, and compared against
the theoretical distribution. The purpose of the shift is discussed in the
text. The three distributions in the inset correspond each to different
sweep rate and noise intensities yet all three with the same ratio. This
illustrates that, while the normalized statistics depend only on the
ratio, the true statistics depend on both values. . . . . . . . . . . . .

52

2.10 Normalized mean (a) and variance (b) of the bifurcation parameter at
escape. The solid line represents the mean of the distribution for a linear sweep predicted from the normal form, as given in equation (2.62).
The data points were calculated from the measured distributions shown
in figure 2.9. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

54

2.8

2.9

3.1

3.2

The blue and magenta lines depict heteroclinic solutions connecting
the operation points to the saddle point for system (3.43), as projected
onto the (q1 , q2 ) plane. Solid lines represent the p = 0 solutions
and the dashed lines are the noise-free (p = 0) stable and unstable
manifolds of the saddle point. Parameter values: Ω = 4.706182 and
β = 8.501802 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

70

Heteroclinic solutions for system (3.43) projected onto the (p1 , p2 )
plane. For the solution with Ω = 4.706182 and β = 8.501802. Note
that these heteroclinics start and end at noise-free points where p1 =
p2 = 0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

71

viii

3.3

Two dimensional parameter space for the Duffing resonator. Solid
lines are saddle-node branches which enclose the bistable region. The
dashed line is the balance curve on which the escape actions for two
stable states are equal, S1 = S2 . In the weak noise limit this line is
very close to the line where the escape rates are equal. . . . . . . . .

72

3.4

Double-well potential of equation (3.45); q1 and q2 denote positions of
stable fixed points, while qS refers to the unstable fixed point. Poisson
pulses with positive area g change the switching rates such that r1 >
r2 , resulting in a shift of the occupation probability from q1 towards q2 . 73

3.5

Change in probability ratio due to the addition of zero mean shot noise,
equations (3.53). Data are results from Monte-Carlo simulations with
D = 0.04 and ν = 0.5. The inset shows the same quantities over a
larger range. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

75

Noise roadmap (a) illustrating the different routes by which fluctuations may affect the angle. Prototypical oscillator block diagram (b).

79

Illustration depicting how noise is filtered and down-converted. The
noise is selected near the oscillator harmonics, that is, in the shaded
regions, weighted by the Fourier coefficients in equation (4.51), and
modulated down to baseband by the oscillator motion. . . . . . . . .

97

4.1

4.2

4.3

Frequency dispersion for the biased Duffing resonator for several values
of the bias. Red dashed line is the locus of critical (zero dispersion)
points. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101

4.4

Several closed orbits belonging to the potential (4.65) with (a) λ = 0
and (b) λ = 4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102

4.5

The logarithm of the normalized long term phase diffusion coefficient,
log10 Sy (0)/Denv , as a function of operating action and bias force.
Dashed white line is the locus of zero dispersion (∂I ω(I) = 0) points.
For interpretation of the references to color in this and all other figures,
the reader is referred to the electronic version of this dissertation. . . 108

4.6

Mean oscillator frequency as a function of operating action and bias
force. Dashed white line is the locus of zero dispersion (∂I ω(I) = 0)
points. For interpretation of the references to color in this and all other
figures, the reader is referred to the electronic version of this dissertation.109

ix

Chapter 1
Introduction
It is well known that small fluctuations in nonlinear systems can have large repercussions in terms of the system response. One means by which this occurs is noise
activated escape, where a series of pulses from a weak noise source may conspire in
such a way as to push a system out of a potential well. For small noise, this process
is a rare event which can be described by an average rate of occurrence. For classical
systems the rate, W , of activated escape is given by the Arrhenius law [100]

W ∝ exp(−R/D),

(1.1)

where R is the activation energy and D is the noise strength. From this expression, we
can see, for small noise (small D), that a small change in the activation energy, ∆R,
corresponding to a small change in the system or its environment, can produce a large
change in W if the ratio ∆R/D is not small. This change in escape rate is described as
a logarithmic susceptibility [100], a nomenclature which implies that the logarithm of
W changes linearly with respect to a small perturbation acting on the system. Thus
highly sensitive measurements can be, and have been [2, 91, 116, 77], made using
the principle of activated escape. However, the full utility of such measurements in
1

the field of MEMS and NEMS has yet to be explored. The aim of this work is to
investigate noise phenomena in micro- and nano- electromechanical systems (MEMS
and NEMS). This includes the investigation of noise activated escape in N/MEMS
resonators for sensing, analysis of noise-dispersed delayed bifurcations for bifurcation
point detection, and phase diffusion in N/MEMS based nonlinear oscillators.
In this chapter we provide motivation and a brief overview of the main contributions described in this dissertation. In section 1.1 MEMS and NEMS are introduced
and demonstrated to be a rich and interesting area of research. In section 1.2 we discuss the intrinsic nature of noise and nonlinear dynamics in MEMS and NEMS that
makes them an appropriate area in which to exploit noise induced phenomena in nonlinear dynamic systems. In section 1.3 we give introductions to the phenomena that
are the primary foci of this work: noise activated escape, delayed bifurcations, noise
in slow-fast systems, etc. In section 1.4 we give a summary of the major contributions
of this work.

1.1

MEMS and NEMS

The field of microelectromechanical systems (MEMS) emerged during the years between the 1960s to 1980s, arising from the integrated circuit industry [75]. It has since
expanded rapidly, driven by commercial successes and the development of new enabling technologies. The ability to fabricate even smaller and smaller devices has led
to the related field of nanoelectromechanical systems (NEMS) [26, 41]. The ability to
batch fabricate such small mechanical structures and integrate them with electronics
has had a large impact on both commercial electronics and research in fundamental
physics alike. Devices such as pressure sensors, accelerometers, gyroscopes, ink jet
printer nozzles, and DLP projectors have found a great deal of success in the commercial arena [75]. Devices such as cell sorters, retinal and cochlear implants, neural
2

probes, and other in vivo sensors have had significant impact on biological research
[51]. Atomic force microscopy [13] and the use of nanoresonators to probe mesoscopic
electrical systems [24, 116, 6] have had significant impact on fundamental and applied
physics research. The list goes on, but this is sufficient to make the point that the
field of MEMS and NEMS is rich and promising.
MEMS and NEMS can generally be split up into various subfields including sensors, actuators, and signal processing devices such as filters and oscillators. In this
work, we focus on sensors and oscillators, specifically, systems that rely on vibration
for operation.
The oscillators we consider employ nonlinear resonators as frequency selective
devices. These resonators are models of N/MEMS devices. Our goal is to optimize
the oscillator performance, in terms of phase noise, by tuning the feedback parameters
and designing the nonlinearity in the N/MEMS resonators.
The sensors we consider are parametric sensors. A parametric sensor is a general
type of sensor in which the signal to be detected influences the dynamic behavior of
the sensing device. The system is then interrogated to retrieve the signal. Examples
include gyroscopes [71] and resonant chemical mass sensors [118]. Gyroscopes encode
the rotation rate in a coupling parameter that governs the interaction between two
resonant modes of vibration. These two modes are designed to have identical natural
frequencies and to be uncoupled save for the coupling arising from Coriolis effects.
One mode is driven, and the other is monitored, or sensed. The amplitude of vibration
of the sensed mode is then correlated with the rotation rate of the gyro. Resonant
chemical mass sensors encode the mass of the substance being detected in the natural
frequency, specifically by shifts in the natural frequency, of a resonator. This is
typically monitored using a self-excited oscillator scheme based on a phase locked loop
[40, 49]. In such a setup, the resonator is made to oscillate at its natural frequency
via feedback, and the frequency of oscillation is monitored and correlated with the
3

measured mass.
While the above examples of dynamical sensors employ linear dynamic systems,
nonlinear dynamics systems are often exploited for properties that might either avoid
the problems inherent in the linear systems, attain more sensitivity, or simply provide
functionality that is not possible with linear systems. For example, Oropeza-Ramos
et al. constructed a gyroscope that employs nonlinear dynamics to circumvent the
difficulty of matching the modal frequencies of its drive and sense modes [89]. Buks
and Yurke explored the idea of using nonlinear dynamics in a chemical mass sensor to
increase sensitivity by operating at a point of infinite slope in the frequency response
curve [15]. A frequency response with infinite slope can only be produced by a
nonlinear system. It was found that the increase in sensitivity came at the price
of a slower response time. For a more complete review of applications of nonlinear
dynamics in MEMS and NEMS, see [73, 92].

1.2

Noise and Nonlinearity

The use of nonlinear dynamics in MEMS and NEMS is in no way artificial. Such small
devices are subject to significant thermal and electrical noise forces. Accordingly, it
is often necessary to drive the devices hard enough to elicit a nonlinear response in
order to obtain a response with good signal to noise ratio. It is by virtue of the
smallness of these devices that the linear dynamic range is made small and one must
account for the noise and nonlinear effects. To illustrate this, consider the following
simple example. The equation describing a simple harmonic oscillator coupled to a
thermal bath is [26]
mω0
2
x+
˙
m¨ + mω0 x = −
x
Q
4

2mω0 kB T
ξ(t),
Q

(1.2)

where x is the displacement, m is the mass, ω0 is the natural frequency, Q is the
quality factor, kB is Boltzmann’s constant, T is the temperature, and ξ(t) is a white
Gaussian noise force with ξ(t)ξ(t ) = δ(t − t ). The noise force, ξ(t), and the dissimω
pative force, − Q 0 x, result from coupling of the resonator to the environment and
˙
are responsible for keeping the oscillator in a state of thermal equilibrium. Dividing
equation (1.2) through by m gives
ω
2
x + ω0 x = − 0 x +
¨
˙
Q

2ω0 kB T
ξ(t).
mQ

(1.3)

As devices are made smaller the mass, m, clearly must decrease if we continue to
use the same materials. The quality factor, Q, also tends to decrease [81, 74, 43].
Meanwhile, the natural frequency, ω0 , tends to increase. Thus, we can see from (1.3)
that as devices decrease in size, the effective noise strength increases. Meanwhile, for
a fixed-fixed beam, for example, the amplitude of vibration corresponding with the
onset nonlinearity goes as tQ−1/2 , where t is the thickness of the beam [41]. This
upper limit of linear dynamic range will decrease with device size if t decreases faster
than Q1/2 . While it is unclear what the result of this competition is in general, we
can concede that for many small devices the linear dynamic range is not large. Thus,
in order to take advantage of the many attractive properties of small devices, one
must consider models that go beyond the linear dynamic range and account for both
nonlinear effects and random forcing. In this simple example, we have only considered
thermal noise acting on the oscillator. However, MEMS and NEMS devices are often
subject to additional noise sources. The reader is referred to [41] for a more complete
discussion of linear dynamic range in NEMS, and to [25, 82] for a review of noise in
MEMS and NEMS.
The host of phenomena created by the interplay of random fluctuations and nonlinearity in dynamic systems is quite large. In this work we will only touch on a
5

U (x)

B
∆U
x
A
Figure 1.1: Activated escape in one dimension.
few examples, all of which can be found in nonlinear MEMS or NEMS resonators
subject to noise, but they are also relevant to other physical systems of interest. In
the following section we will introduce the phenomena of interest and cite some of
the other interesting contexts in which they arise.

1.3

Noise-Induced Phenomena

Perhaps the most basic noise-induced phenomenon is diffusion. Diffusion arises from
the sum of a large number of random variables. The value taken by this sum is
itself random with a variance that grows with the number of variables in the sum.
Brownian motion is one example [39]. Brownian motion can be illustrated by a
particle suspended in a fluid. The random force acting on the particle resulting from
many collisions with the fluid molecules is integrated in time to produce a wandering
trajectory. A similar situation occurs in clocks [72]. A clock cannot have knowledge
of a reference point in true time and so random perturbations on the clock accumulate
causing the phase to diffuse. This results in the clock’s estimate of true time drifting.
Quartz and MEMS clocks are ubiquitous in electronic devices and the reduction
of phase noise in these clocks can have significant impact on communication and
navigation for example.
Another very important noise-induced phenomenon is activated escape. As men6

tioned above, activated escape is a rare event where a series of noise pulses drives a
system out of a confining potential well, or metastable state. To understand this in
more detail, but still from a qualitative perspective, it is helpful to consider a simple
one dimensional case. Suppose we have a particle trapped in a potential well, as
illustrated in figure 1.1, and subject to a random force. In the overdamped limit, the
particle behaves as if it has no inertia and its dynamics are captured by the Langevin
equation

x = −U (x) + ξ(t)
˙

(1.4)

where U (x) is the potential and ξ(t) is the random force. If the noise, ξ, is small, and
the particle starts inside the potential well, then it will very likely relax to near the
bottom of the well at point A. The particle will then rattle around the vicinity of
point A for some time as the noise pushes it back and forth. On occasion a rare large
outburst of noise may overcome the potential force, −U (x), and push the particle
up and over the hill, past point B. It will then roll down the hill and will have thus
escaped the potential well. If the probability of such an outburst is nonzero, then it
will inevitably occur, if one waits long enough. The question is: what are the statistics
of the escape times for a given noise process? We can answer this question for a few
cases. The simplest of which is where ξ(t) is a white Gaussian noise. Considering this
case, we take ξ(t) = 0 and ξ(t)ξ(t ) = 2Dδ(t − t ). This is just restating that ξ is
white and adding that it has zero mean and the strength of this noise is 2D, which
we assume to be small. The escape rate, W , in this case is given by the celebrated
Kramers result [64, 93]

W =

1
2π

∆U
U (xA )|U (xB )| exp −
.
D

7

(1.5)

A rate of this form is often referred to as an Arrhenius equation, or Arrhenius law,
which gives the rate constant, k, of chemical reactions in terms of the temperature,
T , and activation energy, Ea ,

k = A exp −

Ea
RT

(1.6)

where A is the prefactor and R is the gas constant.
Near a simple bifurcation point with a zero eigenvalue the dynamics of a system,
e.g., a noisy nonlinear resonator, squeeze down onto a one dimensional slow manifold,
under certain conditions [62, 19]. If a resonator is operated just before a saddle-node
bifurcation, for example, the dynamics are governed by a one dimensional equation
with a potential similar to that shown in figure 1.1. Such systems are used in the
readout of superconducting quantum bits, or qubits. Readout of qubits is difficult with
conventional amplifiers because the amplifier must minimally disturb the system while
maintaining high sensitivity [98]. However, a method proven to be successful is to
employ what is called a Josephson junction bifurcation amplifier [98, 97, 116, 78, 76].
The method is to operate a Josephson junction in a bistable regime near a saddle-node
bifurcation. The Josephson junction is used because “the superconducting tunnel
junction is the only electronic dipolar circuit element whose nonlinearity remains
unchanged at arbitrary low temperatures” [98]. The resonator is coupled to the qubit
so that the actual position of the bifurcation point depends on the qubit state. In
other words, the height of the energy barrier preventing escape depends on the qubit
state. The resonator is operated in one of its states near the bifurcation point, held
there for some time. It is then backed away from the bifurcation point, after which the
state of the resonator is measured. If activated escape occurs in the resonator, that
is, the response switches from one metastable state to the other, the amplitude and
phase of the resonator will be found to have been changed. Because the probability
8

of such a state switch in the resonator by activated escape is related to the state of
the measured qubit through the energy barrier, the state of the qubit can be inferred
with some degree of confidence by determining whether or not the resonator switched
during the measurement time.
Josephson junctions are commonly used in many practical uses of activated escape, because this switching phenomenon is intrinsic in their behavior. In the early
study of the Josephson junction, it was found that the low current features of the
current-voltage relationship across the junction, which are not predicted by a noisefree theory, are the result of activated escape events in a potential shaped like a tilted
washboard [4, 102]. Thus, Josephson junctions make excellent test-beds for experimental observation of escape phenomena, including so-called resonant activation [28],
and switching between dynamic states over a dissipation barrier [117]. The washboard
potential of a current biased Josephson junction is of the form

U (x) = −ax − cos x,

(1.7)

which is analogous to a classical pendulum subject to a constant torque. When
operated in the overdamped regime, these devices exhibit activated escape in one
dimension, similar to the discussion above. The use of these junctions was proposed
for the experimental measurement of the full counting statistics of electrical current
[109]. However, the details of such an experiment prevent the use of the overdamped
regime, requiring underdamped junctions to be used instead [105, 88]. Nevertheless,
the third cumulant of the current fluctuations can be measured from one dimensional
escape event statistics [56, 70] using the asymmetry of the switching rate with respect
to the bias current. In the study of current fluctuations using activated escape, the
random force is no longer white and Gaussian, but has an arbitrary set of cumulants
that we are trying to measure. The problem of escape due to an arbitrary noise force
9

Figure 1.2: Most probable paths (dashed lines) of activated escape in a two dimensional system. Invariant stable and unstable manifolds (solid lines) represent underlying noise free dynamics.
is still open. However, the problem of poisson noise induced escape in one dimension
has been recently solved [31].
Activated escape out of a multidimensional basin of attraction, as in the case of
the underdamped Josephson junction, is more complex than escape in one dimension.
For the underdamped Josephson junction the basin is two-dimensional, where the
dimensions are the system coordinate and conjugate momentum. The coordinate
in this case is the phase difference of the superconducting order parameters across
the junction. The added complexity in the multi-dimensional case comes from the
necessity to consider the path of escape in the higher dimensional phase space. In the
one dimensional case we know the system must travel from point A to point B, from
the local potential minimum to the local maximum, and it can only do so in a straight
line. In the multidimensional case the escape event occurs along some path through
phase space [37]. In the case of small noise, these paths are tightly packed around
an optimal path which can be found by solving the variational problem for the most
probable path. Such paths have been experimentally observed in MEMS resonators
[21, 20]. Nevertheless, there are details that prevent a complete solution of the escape
problem. The escape rate can only be found up to an approximately constant prefactor. Escape in distributed systems amounts to escape in infinite dimensions. This
10

phenomenon arises in nucleation [100] at the onset of changes in thermodynamic
phase, for example, in boiling. Figure 1.2 shows an example of optimal paths of
escape for a bistable two dimensional system (described in more detail in Chapter 3).
The figure shows the invariant manifolds (solid lines) of the saddle point in a system
containing one saddle and two stable foci. The invariant manifolds give a picture of
the underlying noise free dynamics. The stable manifold of the saddle point divides
the phase space into the basins of attraction of the two stable points. Escape occurs
when noise drives the system across this line. For weak noise, the escape events
traverse paths very close to the optimal escape paths, the first segments of which
are shown as dashed lines in the figure. Once the neighborhood of the saddle point
is reached along such a path, the system will ride a trajectory close to one of the
branches of the unstable manifold of the saddle point, and either escape or return to
the region from which it started.
An interesting situation, related to the above, is escape out of a time varying
potential well. Examples include periodically modulated systems [32] and systems
swept through a bifurcation, where the potential well may vanish, as in the saddlenode or subcritical pitchfork bifurcations, or split into two wells, as in the supercritical
pitchfork. Stochastic-dynamic bifurcations are particularly interesting candidates for
application in MEMS and NEMS devices as they can easily be employed for device
interrogation by setting up simple parameter sweeping experiments. The concept of
a bifurcation in a stochastic system must be treated somewhat delicately, however.
The reason is that, unlike a deterministic system, the solution of a stochastic system,
described by a probability distribution over phase space, does not undergo any sudden
change in form or structure as a function of the system parameters. Therefore, in
some sense, the concept of a bifurcation point in a stochastic system is meaningless
and the notion must be replaced with a bifurcation region [79]. Still, the underlying
deterministic drift in a stochastic system still retains the structure of a deterministic
11

system, and so we can speak of the bifurcation points of a corresponding noise free
system. In this work this is the perspective we take. When we speak of an escape
event occurring before or after the bifurcation point, we mean that it occurs at a time
during the sweep before or after the bifurcation of the deterministic drift term. For
the present investigation, a suitable partition of stochastic-dynamic bifurcations is to
distinguish between the concepts of escape before (activated escape) or after (delayed
bifurcation) the noise-free bifurcation point is reached.

Escape before the bifurcation point in a swept system occurs with very high probability when the system is swept in a quasi-static fashion. That is, when the time-rate
of change of the potential well is much slower than the escape attempt frequency. For
small noise intensity a system will form a quasi-stationary distribution at the bottom
of a potential well. Only rarely does a large noise pulse come along that pushes the
system farther away from the minimum than the width of this quasi-stationary distribution. The escape attempt frequency is the rate of occurrence of these large noise
pulses. Only a few of these large pulses will be large enough to cause escape. In this
limit of slow change with respect to the escape attempt frequency, the escape rate in
a time-varying potential can be well approximated by the quasi-static rate. For example, if we have a time-varying quasi-static escape rate, W (t), then the probability
to escape at time t, Pe (t), is governed by the differential equation
˙
Pe (t) = W (t)(1 − Pe (t)).

(1.8)

If we take initial conditions Pe (t0 ) = 0, i.e. assume that the system starts inside the
potential well, then the solution to equation (1.8) is
t
Pe (t) = 1 − exp −
12

W (t )dt
t0

.

(1.9)

Measuring the distribution Pe (t) by many repeated sweeping experiments can be
a useful tool for gaining information about a system in the vicinity of a bifurcation
point. For example, this experimental technique has often been employed in the study
of switching fields in magnetic material, which has applications in magnetic recording media [115], high speed magnetoelectronics [63], and magnetic RAM cells [106]
among other things. In the late 1980s a fundamental question was determining the
power law dependence of a coercive field strength, as measured from the critical field,
on the switching rate for a single particle magnet. The original model, due to N´el
e
[115], suggested a power law dependence of (∆H)2 . However, the refined N´el-Brown
e
model, suggested a dependence of (∆H)3/2 [120]. Another relevant issue is whether
switching dynamics actually follow a simple Arrhenius law, or is it more complicated?
Experimental confirmation of the 3/2 dependence and agreement with Arrhenius law
was obtain using sweeping experiments [120] following the theory of Kurkij¨rvi [67].
a
However, in submicron-sized magnetic thin films a different behavior is observed. Activated escape experiments, where the system is placed close to the bifurcation point,
and then the time until switching is recorded, showed that the behavior in magnetic
thin films can be explained as an “energy ladder” wherein the system climbs through
repeated thermally activated escape events up a ladder of states of increasing energy
[63]. In [106] the authors use sweeping experiments to characterize the switching
dynamics of submicron magnetic tunneling junctions, used in magnetic RAM cells,
and found the Arrhenius law to hold. In [84, 65] the authors use thermally activated
escape to experimentally investigate magnetic reversal due to spin-polarized currents.
In [119] the authors use sweeping experiments based on the theory of Kurkij¨rvi [67]
a
to determine the magnetic properties of single-molecule magnets: the energy barrier, magnetic anisotropy constant, spin, and crossover temperature from classical to
quantum regime.
Similar sweeping experiments have also been use to study superconducting inter13

ference devices (SQUIDs), atomic friction, and rare molecular events such as protein
unfolding and ligand dissociation. The study of escape is important in determining
limits of sensitivity of weak link superconducting devices [47], such as SQUIDs, which
are very sensitive magnetometers. In a superconducting ring with a weak link junction, however, the fluctuations in the internal flux can be very small. Thus, it proves
convenient to investigate such systems through dynamic bifurcation or sweeping experiments [67, 58]. In experimental studies of friction between two surfaces, a friction
microscope is used where a tip is dragged at constant slow velocity across a surface.
The tip exhibits stick-slip behavior which can be described by activated escape out of
a diminishing potential well [95, 83], similar to the magnetic and superconducting ring
experiments mentioned above. This rate-state model for atomic friction is, of course,
just one of many models for friction [113]. In the study of rare molecular events,
pulling experiments are often used [57]. The technique is also referred to as dynamic
force spectroscopy [29]. These experiments are carried our by using an atomic force
microscope or laser tweezers to pull on an anchored molecule via a pulling spring,
possibly a linker molecule, until some molecular transition occurs, e.g. unfolding,
unwrapping, or dissociation. These experiments are also analogous to the sweeping
experiments discussed above.
Escape after the bifurcation point occurs with high probability when the sweep
rate is sufficiently fast such that the time it takes the system to respond in a significant
way is greater than the time it takes to sweep past the bifurcation point. This
occurs with certainty in situations where the bifurcating fixed point is stable and not
metastable before the bifurcation point [87]. The concept of escape, when discussing
escape after the bifurcation point, is, however, somewhat poorly defined because
after the bifurcation point there is typically no longer a well defined region from
which the system can be said to have escaped. This phenomenon is also called
delayed bifurcation, though this term is usually used in the literature to refer to
14

the noise free situation. The term delayed bifurcation is also ambiguous as there
is no well defined counter-part to the static bifurcation point in the swept (time
varying) system. Thus multiple different definitions for the bifurcation time have
been suggested. These include first threshold crossing time [104, 110], last threshold
crossing time [59], time at which the variance of escape trajectories reaches a threshold
[122, 14], and time of escape to infinity. The latter definition is employed in the
present investigation, and to our knowledge, it is not used by anyone else. One
reason for this might be that the bulk of the literature addresses delayed bifurcation
in noisy systems that are globally bounded, since the systems considered exhibit
supercritical pitchfork or Hopf bifurcations. Both of these bifurcations have one
stable equilibrium before the bifurcation, and after the bifurcation this equilibrium
is unstable, and there exist stable states that emerge from this equilibrium at the
bifurcation. In the case of the supercritical pitchfork, the stable post-bifurcation
states are a symmetric pair of stable equilibria, as in bucking of a symmetric elastic
beam. In contrast, the supercritical Hopf bifurcation gives rise to a stable periodic
orbit. The literature focuses on these supercritical bifurcations because they describe
the important problem of laser turn-on [107, 103, 104, 122, 110, 111, 14]. These
bifurcations can also be used to describe the firing of nuerons [59]. In contrast,
the saddle-node and subcritical pitchfork bifurcations we consider in this work both
exhibit metastable states before the bifurcation and escape from the region of this
state after the bifurcation. The time of escape to infinity is finite for the generic local
models of these bifurcations. Since these bifurcation models describe a local picture
of the dynamics of the full system model in the region of interest (in both the phase
space and parameter space), the (finite) time of escape to infinity for the local models
provide good estimates of the times of escape from the local region.

15

1.4

Summary of Results

In this work we discuss three different, but related, topics. The first topic is the
dynamics of a system subject to a parameter sweep through a bifurcation in the
presence of noise. This topic is examined in chapter 2. In this chapter we review
the normal form models of dynamic systems near bifurcation points. We then review
noise-activated escape in one-dimensional systems subject to white Gaussian noise.
Our primary contributions follow this review material. We first develop approximate
distributions for the time at which a system swept through a subcritical pitchfork
or a saddle-node bifurcation will leave the local region of phase space. The general
problem must be solved numerically, and we do so for comparison with the approximate solution. Our approximate expressions for these distributions are developed
in the limit of “fast” sweeping. We determine that the ratio of sweep rate to noise
strength is the important parameter condition that separates fast and slow sweeping.
These results are compared with the numerical solution of the general case. We then
present the results of sweeping experiments performed on a MEMS resonator, carried
out in collaboration with colleagues at the University of California-Santa Barbara.
Our primary conclusion is that delayed bifurcation is prominent in MEMS devices.
While small, the noise levels we expect to see in MEMS devices is not large enough
to cause an activated escape event before the bifurcation is reached, unless the sweep
rate is extremely slow, to the point of being infeasible for some devices. Therefore,
understanding this phenomenon is important if one wishes to locate the bifurcation
point for system identification, sensing, or for use as a bifurcation amplifier. We however also conclude from our study that performing many sweeping experiments is an
overly time-consuming way to retrieve information from a MEMS resonator. Ongoing work makes use of feedback to track the bifurcation point for bifurcation-based
sensing strategies [17].
16

Our second topic is the novel dynamical bridge sensing paradigm, which we discuss
in chapter 3. This sensing methodology employs a bistable system arranged such that
the rate of noise-activated escape out of each stable basin of attraction is identical.
The ratio of occupation probabilities can become extremely sensitive at this point
if the noise is weak. Our contributions to this sensing concept are divided into two
parts. First, we calculate the sensitivity and detection time of the dynamical bridge
for use as a general use detector. Second, we develop the conditions required and
demonstrate the application of the dynamical bridge as a detector of non-Gaussian
noise. We illustrate the measurement of the parameters of a shot-noise process with
a one-dimensional bridge model. In terms of the bridge as a general use detector, we
conclude that it is a time-consuming sensing method. The reason is that the bridge
sensitivity scales with the inverse of the noise strength while the measurement time
is exponential in this same quantity. As a detector of the full-counting statistics of a
shot-noise process, the bridge should perform comparable to other methods with the
added convenience that one does not need to reset the system.
Our third topic is phase noise in oscillators employing nonlinear frequency selective
elements. This is discussed in chapter 4. We review the theory of a nonlinear oscillator
coupled to a medium [36], and employ this microscopic resonator model to gain
insight into the origin of phase noise. Our primary contribution is the use of the
method of averaging to develop an expression for the spectrum of fractional frequency
fluctuations in an oscillator. The method of averaging applied to a noisy system is
discussed in appendix B. The expression for the spectrum of frequency fluctuations
is quite general, although it does ignore dynamics in the feedback elements of the
oscillator. The expression contains elements of the noise model and nominal limit
cycle shape of the oscillator. We demonstrate its utility by exploring the parameter
space of a prototypical oscillator. It is shown that the action-angle decoupling at
zero-dispersion points of the resonator frequency provide an operating point that
17

reduces the phase diffusion constant. Our conclusion is that these zero dispersion
points are near the optimal operating points for strongly nonlinear oscillators. The
reduction in phase noise gleaned from oscillator tuning should be complimented by
the reduction of noise influence at the source. These results are relevant to the design
of high performance N/MEMS frequency sources.

18

Chapter 2
Escape Statistics for Parameter
Sweeps Through Bifurcations
In this chapter we discuss the dynamics of nonlinear N/MEMS resonators near codimension one bifurcation points that result in large changes in system response. Of
the phenomena discussed in this dissertation, we choose to begin with this because
the one-dimensional models we use provide ease of visualization and yield more complete solutions than the more complicated systems considered in later chapters. Our
primary purpose is to examine the initial transient that occurs when the parameters of a system are swept through a saddle-node or subcritical pitchfork bifurcation.
Since these bifurcations exhibit a loss of local stability, and have no stable states
emerging from the bifurcation, we can qualitatively say that the system will undergo
a transient, a “jump,” or an “escape”, to another region in phase space. These escape events are highly visible transient processes and therefore useful for detection.
However, they can also be strongly influenced by noise and therefore must be understood as random events. In this chapter we formulate the general problem, determine
conditions which separate fast and slow sweeping, and then develop the distribution
of escape times in the limit of fast sweeping through a saddle-node or a subcritical
19

pitchfork bifurcation. These are new results, relevant for many N/MEMS resonators,
and possibly for other systems as well.
This work is motivated by the applications of escape statistics to sensing, parameter fitting, and system identification. For example, circuits containing Josephson junctions often exhibit saddle-node bifurcations. Interrogating such systems by
sweeping through the bifurcation point has been used to fit the critical current parameter of the junction [67, 47, 58]. In the study of nanomagnetic systems, the measured
distribution of escape events was used to determine the type of bifurcation occurring
in the the system, leading to support for certain system models [115, 120]. More recently escape due to a parameter sweep trajectory approaching a bifurcation has been
used as a method for reading the state of a qubit coupled to a nonlinear resonator
[116]. Another important example is the supercritical pitchfork or Hopf bifurcation
encountered during parameter sweeps during laser turn-on [14, 110, 122, 8]. This
example, however, does not exhibit loss of local stability and so differs qualitatively
from the others.
As one might expect from these example applications, noise-activated escape near
bifurcations has been the subject of many previous investigations. Many of these
focus on systems in which parameters remain fixed [85, 108] or vary slowly [67, 48,
44, 116, 44, 54]. However, systems with comparatively weak noise typically exhibit
delayed bifurcation during parameter sweeps [16]. It is this regime in which we
are primarily interested. Berglund and Gentz [9] developed scaling laws which are
applicable to the distribution of delayed bifurcation trajectories. However, more
quantitative results are required for bifurcation-based sensing schemes, for example,
in dynamic M/NEMS [16, 66].
We will begin, in section 2.1, by discussing the modeling of a nonlinear N/MEMS
resonator by a reduced model, the normal form, near the bifurcation point. This
section concludes by recalling the normal forms for the saddle-node and subcritical
20

pitchfork bifurcations. We then discuss, in section 2.2, how noise can cause the
system to escape from the region of the bifurcating equilibria by reviewing Kramers
[64, 93] escape formulation, which is valid in the limiting case of slow (adiabatic)
sweeping. We then proceed with our solution to the problem of delayed saddle-node
and subcritical pitchfork bifurcations. The saddle-node is discussed in section 2.3
and the pitchfork in section 2.4. Experimental demonstration of a delayed subcritical
pitchfork bifurcation is presented in section 2.5 and future research directions related
to this subject are discussed in section 2.6.

2.1

Stochastic Normal Forms

In this section we sketch the transformation of a weakly nonlinear resonator to a
reduced one-dimensional stochastic differential equation that describes the resonator
dynamics close to a bifurcation point. It is well known that the resonantly forced
Duffing and nonlinear Mathieu resonators exhibit saddle-node and pitchfork bifurcations of equilibria, respectively, in an averaged rotating frame in which equilibria
represent nearly-harmonic periodic system responses. In appendix C we develop the
pitchfork normal form arising in the nonlinear Mathieu type resonator. In this section, however, we leave our discussion more general. For more background material,
the reader is directed to [52] for the noise-free case. For additional discussion of
stochastic normal forms, see [62].
A resonator with weak damping and weak nonlinearity can be described by
2
q + ω0 q = F (q, q, t),
¨
˙

(2.1)

where q is the resonator coordinate and the function F represents small forces acting
on the resonator including damping, nonlinear effects, periodic forcing, stochastic
21

processes, etc. From equation (2.1) we proceed to make a series of approximations
based on the relative size of various terms. The first step is to transform equation
(2.1) into a set of coupled first order equations with small right hand side. This
can be done with the well-known Van der Pol transformation. A polar form of the
transformation is given by
q = a(t) cos(ωt + φ(t)),
(2.2)
q = −ωa(t) sin(ωt + φ(t)), .
˙
where ω is the frequency of the harmonic excitation in F . This transformation is
analogous to the method of variation of parameters. Applying equations (2.2) to
(2.1) gives
F (a, φ, t)
,
ω
F (a, φ, t)
˙
φ = − cos(ωt + φ) σ cos(ωt + φ) +
.
ωa
a = − sin(ωt + φ) σa cos(ωt + φ) +
˙

(2.3)

2
where σ = (ω 2 − ω0 )/ω is the detuning, which is assumed small, and F (a, φ, t) is
understood to imply F (q(a, φ, t), q(a, φ, t), t). Owing to the smallness of the right
˙
hand side, equations (2.3) imply that a and φ change on a much longer time scale
than the period of oscillation 2π/ω. Accordingly, fast oscillating periodic terms in
(2.3) integrate to approximately zero and can be ignored. A more complete discussion
of this approximation, called the method of averaging, is given in appendix B, where
particular attention is given to dealing with stochastic processes. The result is an
approximate set of envelop equations describing the slow variation of the resonator’s
amplitude and phase. We may assume that the equations (2.3) can be approximated
by (see appendix B)

˙
x = f (x, µ) + G(x, µ)ξ(t),
22

(2.4)

where we have lumped a and φ into a state vector, x, µ is some resonator parameter, ξ(t) is an n-dimensional vector of zero mean stochastic processes, and G is a
2-by-n dimensional matrix. Now, suppose there exists a point, x∗ and a corresponding parameter value, µ∗ such that f (x∗ , µ∗ ) = 0 and the Jacobian at this point,
Dx f (x∗ , µ∗ ), has one zero eigenvalue and one negative eigenvalue, −λ. The bifurcation point is understood to be the pair (x∗ , µ∗ ). The bifurcation value is µ∗ alone,
and the term bifurcation describes a change in the sign of the real part of an eigenvalue
of an equilibrium point resulting from a change of the parameter µ near µ∗ .

Let u be the deviation from x∗ along the stable eigenvector, i.e., the eigenvector associated with eigenvalue −λ, and let v be the deviation from x∗ along the
marginally stable eigenvector, i.e., the eigenvector associate with the zero eigenvalue.
We express equation (2.4) in the (u, v) coordinate system and expand in both the
coordinates and the parameter about (u, v, µ) = (0, 0, µ∗ ). Keeping leading order
terms, the equations for u and v take the form

u = −λu + a1 ξ1 (t) + · · · ,
˙

(2.5)

v = b1 v n1 + b2 uv n2 + b3 u2 v n3 + cv n4 δµ + a2 ξ2 (t) + · · · ,
˙

(2.6)

where ai ξi are the weighted projections of ξ onto the two eigenvectors and δµ = µ−µ∗ .
For sufficiently small noise and small u, v, and δµ the u dynamics evolve much faster
than the v dynamics owing to the nonzero eigenvalue −λ. Thus we assume that u is
quite small. Accordingly, we keep the mean contribution from those terms containing
powers of u in equation (2.6), that is the terms that arise in the deterministic normal
form [52], but ignore the fluctuating parts of u, as it is dominated by the fluctuating
term a2 ξ2 (t). This amounts to simply averaging equation (2.6) over u, and it results
23

in the stochastic normal form
v ≈ b1 v n1 + b3 u2 v n3 + cv n4 δµ + a2 ξ2 (t).
˙

(2.7)

Expanding equation (2.4) in powers of u and v is called the tangent space approximation since u and v describe a deviation in the eigenspace of the bifurcation point. A
more detailed reduction can be performed by seeking the deterministic center manifold u = h(v) + such that u = h (v)v for ξ = 0. h(v) can be found by expressing it
˙
˙
as a series in v and then using the condition u = h (v)v to solve for the coefficients by
˙
˙
matching like powers. This provides a more complete picture of the slow manifold in
phase space. We take this approach when we consider the nonlinear Mathieu equation
in appendix C.
The exponents ni in equation (2.7) define the type of bifurcation. The coefficients
bi and c can be given standard values by a scaling of the coordinate v and time
and re-normalizing the bifurcation parameter, yielding a normal form with minimal
parameters.
For the saddle-node bifurcation n1 = 2 and n4 = 0. For the pitchfork bifurcation
n1 = 3 and n4 = 1. In the nonlinear Mathieu resonator we also have n3 = 1, which
simply shifts the bifurcation point, an effect that can be accounted for by re-defining
the bifurcation parameter, δµ [85, 86]. This effect can also be neglected if the noise
is sufficiently weak. Thus, in the saddle-node and pitchfork bifurcations, only the
value of a2 and the relationship between δµ and the original resonator parameters
are specific to any given system undergoing such a bifurcation. For this reason the
model is very general and valid for a wide class of systems undergoing the bifurcations
of interest.
There are two types of pitchfork bifurcations, the supercritical, for which b1 > 0,
and the subcritical, for which b1 < 0. As mentioned above the supercritical pitchfork
24

arises in laser dynamics, but also in parametric resonance of the nonlinear Mathieu
resonator, when the system restoring force is hardening. The subcritical pitchfork
likewise can be found in parametric resonance, for systems with softening restoring
force. Again, we are interested in the subcritical bifurcation because it results in a
loss of local stability and an escape event.
Since the normal form, equation (2.7), is one-dimensional, we can write

v = −U (v) +
˙

√
Dξ(t),

(2.8)

where U is an effective potential, D is the effective noise strength, and ξ(t) is the
effective noise on the slow manifold arising from a combination of the components
of ξ. As we mentioned above, we have neglected any parametric noise terms by
appealing to an appropriate scaling, c.f. appendix C. Coordinate v can be understood
as describing the position of an overdamped particle in the potential U . We are
particularly interested in bifurcations for which U presents a finite barrier for the
system to escape to v → ±∞. The finite barrier allows the system to escape from
the local region of phase space and exhibit a highly visible transient that can be
used for sensing or system identification, as discussed above. Two examples are the
saddle-node and the subcritical pitchfork. The supercritical pitchfork does not have
this property, since responses stay in the neighborhood of the bifurcation point.
Now, we wish to consider the parameter sweeping experiments described above.
Thus, we will let δµ vary in time according to a prescribed parameter sweep trajectory.
In order for the reduced model to remain valid, this variation, as well as other effects,
including the noise, should not push the system far from the slow manifold. This
should be the case so long as the time rate of change of the eigenvectors is much
smaller than λ and the noise is sufficiently weak. For the remainder of the analysis,
we assume these assumptions hold. As we will see, when sweeping the bifurcation
25

U (v)

v1

vmax

vmin v2

a

v

Figure 2.1: Example potential governing normal form dynamics.
parameter linearly in time, only the ratio of the sweep rate to the noise strength is
important. Therefore we speak of fast and slow sweep rates according to whether
this ratio is large or small, but always assume the sweep rate is sufficiently slow (i.e.,
noise sufficiently small) so that the normal form model holds.

2.2

Kramers Escape Rate

Consider equation (2.8) and suppose that the potential U presents a finite barrier
beyond which the system response will become unbounded. For example, see figure
2.1. Since the barrier height is finite, there is some chance that the noise will push
the system out of the potential well and it will escape to ∞. Of course, equation (2.8)
is only a local model of the full system dynamics, obtained from the reduction to the
slow mode. As v grows, this model will eventually become invalid, and rather than
escape to ∞, the resonator state will move to another region of phase space that is
outside the domain of this local picture.
The probability per unit time of escape out of the well was first derived in the
classic 1940 work of Kramers [64]. This result applies to a time-independent potential, but can be applied to a slowly varying potential via an adiabatic approximation.
26

This important result is the counterpoint to the new results we develop in sections 2.3
and 2.4. These new fast sweeping solutions complement the slow sweeping (adiabatic)
solutions since they are the two asymptotic cases that bracket the complete picture.
For this reason we review the celebrated Kramers’ result and the corresponding adiabatic approximation in this section. Our development of Kramers’ escape rate comes
from Risken [93]. We assume that U is time independent, D is constant, and ξ is a
zero mean white Gaussian noise with ξ(t)ξ(t ) = 2δ(t − t ). Under the white noise
assumption, the probability distribution, ρ, for v is governed by the Fokker-Planck
(FP) equation [93]
∂S
∂ρ
= − ,
∂t
∂v
∂ U/D
S = −De−U/D
e
ρ ,
∂v

(2.9)
(2.10)

where S is the probability current. If the potential barrier is high, it is unlikely
that the system will escape and ρ will settle into a quasi steady state distribution.
Accordingly, ∂t ρ is small and S is approximately constant. Making this assumption,
we integrate equation (2.10) over the interval [vmin , a] (see figure 2.1)
a
S
vmin

eU/D dv ≈ D eU (vmin )/D ρ(vmin ) − eU (a)/D ρ(a) ,

(2.11)

where a an arbitrary point sufficiently beyond the barrier peak. We can also assume
that the probability of finding the system at point a is small. We therefore drop the
second term on the right hand side of equation (2.11). Thus, the probability current
is approximately given by
S ≈ DeU (vmin )/D ρ(vmin )

a

eU/D dv

−1
.

(2.12)

vmin

Returning to the initial assumption that ∂t ρ is small, we approximate ρ by its steady
27

state distribution

ρ ≈ ρ(vmin ) exp −[U (v) − U (vmin )]/D .

(2.13)

Thus, the probability for the system to be near vmin is
p=

v2
v1

ρdv ≈ ρ(vmin )eU (vmin )/D

v2

e−U/D dv.

(2.14)

v1

The inverse escape rate, W , is given by the ratio of the probability to be near the
bottom of the well to the probability current out of the well.
a
1
p
1 v2 −U/D
eU/D dv.
e
dv
≡ ≈
W
S
D v
vmin
1

(2.15)

The integrals can be approximated by the method of steepest decent. In this method
U is expanded to second order about the point that maximizes the integrand and
the bounds of the integration are taken to ±∞. For the first integral U is expanded
about vmin , and for the second integral U is expanded about vmax . The resulting
escape rate is thus found to be

W =

1
2π

U (vmin )|U (vmax )| exp −(U (vmax ) − U (vmin ))/D .

(2.16)

If the potential allows escape in both directions, the total escape rate is the sum of
the rates for escape in the two directions.

Now, suppose that the potential is changing slowly in time. We can model the
probability to escape the potential well per unit time with equation (2.16). Thus, the
probability for the system to remain in the well until time t, ρ, is governed by

ρne = −W (t)ρne ,
˙

ρne (t0 ) = 1,
28

(2.17)

where the subscript “ne” stands for “not escaped.” This equation admits the solution
t
ρne = exp

W dt .

(2.18)

t0

The probability to escape from the well by time t is ρe = 1 − ρne . The probability
density, P , of escape at time t, is ρe . This gives
˙
t
P = W exp

W dt .

(2.19)

t0
This result is valid in the slow sweeping limit. In this limit noise-activated escape
occurs with probability 1 by the time the the sweep has reached the bifurcation value.
If the sweep is fast, however, there must be nonzero probability to pass through
the bifurcation value before escape occurs. Results analogous to equation (2.19) for
the saddle-node and subcritical pitchfork bifurcations in the fast sweeping limit are
developed in the following sections.

Before proceeding to discuss escape in the non-adiabatic, or fast, sweeping condition, a few words must be said about what it means to “escape” in this context. In the
slow sweeping limit, discussed above, the duration of the escape event is small compared to the time scale of the sweep. On the long time scale of the sweep it appears
that escape happens instantaneously. Accordingly, the escape rate is independent of
the point a at which escape is defined to occur. However, when the sweep is fast, the
transient escape process may occur over a duration in which the bifurcation parameter changes significantly. Escape, therefore, does not appear instantaneous and the
probability density to escape per unit time depends on the threshold at which “escape” is said to occur. To deal with this complication, we observe that both normal
form models considered herein exhibit escape to ∞ in finite time. Accordingly, the
escape trajectories become extremely steep at (locally) large values of v. We thus
29

x
µ
µ<0

µ=0

µ>0

(a)

(b)

Figure 2.2: Illustration of the evolution of the potential (a), and the standard bifurcation diagram (b), for a saddle-node bifurcation.

define escape to occur when v → ±∞ and note that the arrival time at some large
value of v, taken to be the limit of the domain of validity of the normal form model,
is well approximated by the arrival time at ±∞.

2.3

Escape Near a Saddle-Node Bifurcation

When sweeping a parameter near a saddle-node bifurcation, the dynamics along the
slow manifold are described by a one-dimensional Langevin equation representing the
stochastic normal form. Switching notation from the previous section, the saddlenode normal form is given by [9],
x = µ(t) + x2 + ξ(t),
˙

(2.20)

where µ(t) is the bifurcation parameter and ξ(t) is zero mean, delta-correlated white
noise with autocorrelation ξ(t)ξ(t ) = 2δ(t − t ). For simplicity, we assume that µ
is non-decreasing. In this case, the sweep trajectory can be classified by the charac30

teristic sweep rate, r.

r=

µc − µ0
,
tc − t0

(2.21)

where µc is the bifurcation value (for equation (2.20) µc = 0), tc is defined by
µ(tc ) = µc , µ0 is the initial value of the bifurcation parameter, and t0 is the time
at which the sweep is initiated. The dynamics of equation (2.20) are conveniently
visualized by an overdamped particle moving in a time-varying potential, as illustrated
in figure 2.2a. The bifurcation diagram is shown in figure 2.2b. Our goal is to
determine the distribution of escape times T , i.e., those values of T for which x → ∞
as t → T . A numerical solution for this problem can be obtained by solving the FP
equation associated with equation (2.20),
∂
∂ 2ρ
∂ρ
= − [(µ + x2 )ρ] + 2
,
∂t
∂x
∂x2

(2.22)

with appropriate initial and boundary conditions at large, but finite, x, the results of
which are used to illustrate the connection between the fast and slow sweeping limits.
It can be shown by the following simple scaling argument that these limits are given
by 2

r and 2

r, respectively: Scaling time by r−1/3 and the coordinate by

r1/3 gives rise to a renormalized characteristic sweep rate of unity and a renormalized
noise strength of 2 /r.

Now, suppose that the system is initially near the local minimum of the potential
that exists for µ < 0 and the sweep is slow (or the noise is large). If µ is increased
slowly (for the given level of noise intensity), one can approximate the rate of escape
by the well-known adiabatic approximation due to Kramers [64, 54] as discussed in
31

the previous section. For the saddle-node, this rate is

W (t) =

−µ(t)
4(−µ(t))3/2
exp −
.
π
32

(2.23)

The probability to remain in the well diminishes according to the product of this
rate and the probability to remain in the well. Accordingly, the probability density
function for escape at time T in the limit of adiabatic sweeping is given by [67, 54],
T
P (T ) = W (T ) exp −

W (t)dt .

(2.24)

t0

Of particular interest to us is the case where the noise strength is small compared
to the characteristic rate of change of the bifurcation parameter, that is, 0 ≤ 2

r.

In this case, the potential well will disappear before noise activated escape is likely
to occur, and the system thus exhibits a delayed bifurcation, escaping to x → ∞ at
a finite time T , for which µ(T ) > 0. To approximate the distribution of escape times
for equation (2.20) in this situation, we employ a direct perturbation method. To
this end it is convenient to transform the equation using a variation of parameters
approach, in which the = 0 solution, x0 , plays the role of the homogeneous solution.
The transformation x0 = −u/u in equation (2.20) with = 0 yields,
˙
u + µ(t)u = 0,
¨

(2.25)

which has linearly independent solutions u1 and u2 , so that x0 can be written as,
u + cu2
˙
˙
x0 = − 1
,
u1 + cu2

(2.26)

where c is a constant of integration. To account for the presence of noise, we let c vary
in time. Utilizing expression (2.26) with c(t) in equation (2.20) gives the equation for
32

c(t) as,
c = − a(u1 + cu2 )2 ξ,
˙

(2.27)

where a−1 = u1 u2 − u1 u2 is the Wronskian. The variable c is a means of parameter˙
˙
izing the solutions of the noise-free problem, and noise causes the system to diffuse
among the noise-free trajectories according to equation (2.27).
The escape time, T , is given by the time at which x0 becomes unbounded, which,
according to equation (2.26), corresponds to c(t) reaching the condition,
u (t)
c(t) = c∞ (t) = − 1 .
u2 (t)

(2.28)

Here the first passage time problem reduces to a transformation of the stochastic
variable c(t) to the times T for which c(T ) = c∞ (T ), for some initial distribution
of initial conditions. To determine the distribution of T , it is useful to consider the
features of c∞ (t). The function c∞ (t) has branches, separated by vertical asymptotes
at the zeros of u2 (t), and c(t) evolves between these branches, as depicted in figure 2.3.
Also, note that since
1
dc∞
=
dt
au2
2

(2.29)

does not change sign, c∞ (t) is monotonic between branches. In addition, by considering equations (2.27) and (2.28), it is seen that c = 0 when c = c∞ . Thus, trajectories
˙
of c(t) cross each branch of c∞ (t) at a single time. Consequently, in considering the
first passage time of c(t) across a particular branch of c∞ (t), we do not need to account for trajectories that cross the branch multiple times. We can therefore do away
with the adsorbing wall boundary condition usually employed in first passage time
problems. Without an adsorbing wall boundary condition, a crossing corresponds to
33

reinjection at x = −∞. Hence, subsequent crossings of c∞ occur on immediately
subsequent branches only. These subsequent crossings occur when the system has
once again reached x = ∞.
The change of coordinates from x(t) to c(t), expressed in equation (2.26), makes
the problem amenable to direct perturbation analysis. To develop a first order asymptotic approximation of the escape current, we begin with the expansion,

c = c0 + c1 + · · · .

(2.30)

Substitution of this into equation (2.27) and expanding in powers of

yields,

c0 = 0,
˙

(2.31)

c1 = −a(u1 + c0 u2 )2 ξ.
˙

(2.32)

Note that the first order perturbation removes the state-dependent diffusion, and
thus Ito and Stratonovich calculus give the same result at this order. Solving equations (2.31-2.32) gives a time-dependent Gaussian probability density function (PDF)
for the noisy response of c, as follows,

Pc (c, t) =

4πa2 2

t
t0

−1/2
(u1 + c0 u2 )4 dt

(2.33)




2
(c − c0 )
.
× exp −
t (u + c u )4 dt
2 2
4a
0 2
t0 1
Using this result, the PDF of escape times can be expressed as,

P∞ (T ) = Pc (c∞ , T )

dc∞ (T )
.
dT

(2.34)

Combining equations (2.29), (2.33), and (2.34), we arrive at an expression for the
34

6 c
4

(2)
c∞

2
-2

-1

1
-2
-4

(1)
c∞

2

3

t

( 2 = 0.05)

-6
Figure 2.3: c∞ boundary for linear sweeping, µ(t) = µ0 + rt, with µ0 < 0. The first
two branches are shown along with the variance (shaded region) of c for 2 = 0.05.
PDF of escape times,
−1/2
T
4 2 u4 (T )
4 dt
P∞ (T ) =
4πa
(u1 + c0 u2 )
2
t0


2
(u1 (T ) + c0 u2 (T ))
.
× exp −
2 2 u2 (T ) T (u + c u )4 dt
4a
0 2
t0 1
2

(2.35)

Note that this result is quite general, since the details of the parameter sweep µ(t),
which dictates u1,2 (t), are not yet specified.

The simplest method for sweeping

the bifurcation parameter is linear in time, that is, µ = µ0 + rt, with µ0 < 0 and
r > 0. The affine relationship between the bifurcation parameter and time allows one
to identify the bifurcation parameter as a renormalized time variable. Thus, without
loss of generality, we take µ = t. With this form of sweeping, equation (2.25) becomes
Airy’s equation in backwards time, that is, u + tu = 0, with independent solutions,
¨

u1 (t) = Ai (−t),

(2.36)

u2 (t) = Bi (−t),

(2.37)

where Ai and Bi are the standard Airy functions, for which the standard parameter
a = −π [1]. In principle, u1 and u2 can be taken as linear combinations of Ai and
35

Bi , however, the present choice is most convenient here. The first two branches of
(1,2)
the c∞ boundary, c∞ , are shown in figure 2.3. In the weak noise case, only the
(1)
(1)
second boundary can be crossed, since c∞ is bounded above by zero (c∞ < 0) and
(1)
(1)
monotonically decreases (c∞ < 0), while c0 > c∞ must hold for initial conditions
˙
restricted to be near the minimum of the potential well. Thus, the only way to
(1)
cross c∞ is for c(t) to diffuse to +∞ and be reinjected at −∞, and this cannot
happen under the weak noise assumption. Therefore, computing escape times we
(2)
(2)
consider only those trajectories that reach c∞ . Note that c∞ spans the time interval
between the first and second zeros of Bi (−t), so that our approximation of the escape
distribution times is necessarily limited to this time interval, which is given below. In
principle, this time interval can be adjusted by taking a different linear combination
of Ai and Bi for u1 and u2 , since the time interval is specified by the zeros of u2 .
However, as the perturbation solution requires small noise, the probability of escape
outside this time interval is also small.
Equation (2.35) captures the escape distribution with the system initially situated
on the noise-free trajectory specified by c0 at time t0 . Generally, the dependence on
c0 and t0 becomes weak as |t0 | becomes large since
2
c ∼ exp − |t0 |3/2 ,
3

(2.38)

assuming the system begins the sweep near the bottom of the potential well. This
weak dependence on this class of initial conditions is a property of the deterministic
system and results from the annihilation of the fixed point after the bifurcation.
Diffusion only increases this effect as the system settles into a quasi-steady state
distribution early in the sweep. Thus, when the sweep is started well before the
bifurcation point, the initial conditions are forgotten.
With this understanding, the initial conditions for the linear sweep are taken to
36

be c0 = 0 and t0 = −∞. The resulting escape distribution then applies to all initial
−1. The corresponding escape current is

conditions where c0 is near zero and t0

computed using equations (2.35) and (2.37), resulting in,
−1/2
T
4 (−t)dt
5 2 B 4 (−T )
Ai
P∞ (T ) =
4π
i
−∞


2 (−T )
−Ai
,
× exp 
T
2
4π 2 2 Bi (−T ) −∞ A4 (−t)dt
i

(2.39)

with T restricted to the range between the first two zeros of Bi (−t), that is, 1.17 <
T < 3.27.
The first two moments of the escape distribution given by equation (2.39) are
illustrated in figure 2.4 along with two sample distributions, one representative of the
slow sweeping limit and one representative of the fast sweeping limit. The moments
are plotted over a wide range of noise strengths, spanning adiabatic (Kramers) escape
and delayed bifurcation. The solid lines represent the solution obtained by numerically
solving the FP equation over a finite domain. The equation was solved using the
Crank-Nicolson method with Dirichlet boundary conditions at large coordinate values
where drift dominates diffusion. The moments for the adiabatic approximation and
the approximation developed here are shown as dashed lines. Note that the two
limiting cases bracket an intermediate region near 2 = 1, where one must rely on the
numerical solution to obtain the distribution. Thus, the present analysis, in addition
to providing convenient approximations for the fast sweep rate escape distributions,
also offers information about the parameter limitations for the two limiting cases, and
a method for determining the distributions when neither approximation holds.
Note that the variance is dramatically larger in the adiabatic case, 2
compared to the non-adiabatic case, 2

1, when

1. This is the result of the different escape

mechanisms in these two limiting cases. In the adiabatic case, escape is a rare event
37

2

T

0
-2
-4
-6
-8
10−4 10−3 0.01 0.1
1
10
Noise to Sweep Ratio 2
(a)
103

100

Var(T )

10
0.1
10−3
10−5
10−4

0.01
1
100
2
Noise to Sweep Ratio
(b)
2 = 0.01
2 = 53

0.06

4

0.05
3

0.03

P∞

0.04

2

0.02
1

0.01
-60

-40 -20
0 1.8
T
(c)

2.2
T

2.6

Figure 2.4: Mean (a), variance (b), and sample distributions (c) of the escape time
away from a swept saddle-node bifurcation. Numeric solution (solid line) and approximations (dashed) for fast and slow sweeping are shown.

38

x

µ

µ<0

µ=0

µ>0

(a)

(b)

Figure 2.5: Illustration of the evolution of the potential (a), and the standard bifurcation diagram (b), for a subcritical pitchfork bifurcation.
in which the noise overcomes a potential barrier. Thus, the probability of escape
over a small time interval is small, but the sweep takes a long time, so escape is
virtually inevitable. In this way, the escape events have a wide distribution. In
the non-adiabatic case, escape is inevitable, since the potential well disappears, and
the solution closely follows the deterministic trajectory that would result from an
initial condition at the bottom of the potential well. The escape times are randomly
distributed about this trajectory due to diffusion. These observations also indicate
why the mean escape time is insensitive to the noise strength for small noise intensity,
but highly sensitive to the noise strength for large noise intensity.

2.4

Escape Near a Subcritical Pitchfork Bifurcation

Near a subcritical pitchfork bifurcation, the dynamics with noise along the slow manifold are described by the normal form [62],
x = 2µ(t)x + 4x3 + ξ(t),
˙
39

(2.40)

where, as before, µ is the bifurcation parameter, 2 is the noise strength, and ξ is zero
mean, delta-correlated white noise with autocorrelation ξ(t)ξ(t ) = 2δ(t − t ). It is
recognized that the presence of noise shifts the effective bifurcation parameter [85], so
µ here includes that shift. The illustration of this system as a particle in an evolving
potential is shown in figure 2.5a. The bifurcation diagram is shown in figure 2.5b.
The FP equation corresponding with equation (2.40) is given by,
∂
∂ 2ρ
∂ρ
= − [(2µx + 4x3 )ρ] + 2
,
∂t
∂x
∂x2

(2.41)

where ρ = ρ(x, t) is the probability distribution for the system state. As in the saddlenode case, we compare approximate results for fast and slow sweeping with numerical
solutions of this FP equation.
Suppose the system is initially near the local minimum x = 0 for µ < 0. If µ is
increased slowly, the rate of escape can be approximated by the well-known adiabatic
Kramers result [64, 54], given by,
√
8
µ2 (t)
.
|µ(t)| exp −
W (t) =
π
42

(2.42)

The probability to remain in the well diminishes according to the product of this
rate and the probability to remain in the well. Accordingly, the probability density
function for escape at time T in the limit of adiabatic sweeping is given by [67, 54],
T
P (T ) = W (T ) exp −

If µ is swept quickly, that is 0 ≤ 2

−t0

W (t)dt .

(2.43)

µ, the system response is dominated by delayed
˙

bifurcations. We develop an approximate expression for the escape current in this
case using asymptotic expansions.
Note that the method employed for the saddle-node bifurcation is not convenient
40

for the pitchfork, since the approach is not capable of capturing noise induced independence of initial conditions. This results from the fact that, for the pitchfork,
the noise-free system maintains an equilibrium state, albeit unstable, beyond the bifurcation point. For the noise-free system, the initial conditions enter the solution
multiplicatively and cannot be neglected. Nevertheless, diffusion removes the dependence on initial conditions when the parameter sweep is started sufficiently far from
the bifurcation point, such that the system settles into a quasi-steady-state distribution [122, 14, 104]. To capture this effect we must employ an alternative method.

In this case the solution of the corresponding FP equation is approximated by the
method of matched asymptotic expansions. Following Suzuki [107], we separate the
sweeping time into a first interval, over which nonlinear drift is ignored, and a second
interval, over which diffusion is ignored. We show that, under some assumptions,
there exists a time interval over which both approximations are valid, allowing one
to match the solutions. A change of coordinates in equation (2.41) makes clear the
two phases of response. Let

y =

x
σ

where

(2.44)


σ2 =

t
2 e4 t0 µ dt σ 2 + 2
0



t −4 t µ dt
t0
e
dt  .
t0

(2.45)

and where µ is assumed to grow monotonically with t. This particular σ 2 is the
variance of the solution to the linearized problem in which the nonlinear drift term
is neglected. By changing coordinates into one in which the linearized problem has
a constant solution, we more explicitly expose the roles of nonlinear drift and its
relationship to diffusion. σ has the following properties: it is positive, takes on the
value σ0 ≥ 0 at time t = t0 , it decreases until the time at which 2µσ 2 / 2 = −1, which
must occur before the bifurcation point is reached, and it increases monotonically for
41

σ0 < / −2µ0 . Taking y as the new coordinate, equation (2.41) becomes,
2 ∂ρ
2 ∂ 2ρ
∂ρ
∂
−
= −2µρ − 4σ 2
y3ρ +
y
.
∂t σ 2 ∂y
∂y
σ 2 ∂y 2

(2.46)

The second term on the LHS of equation (2.46) arises since y changes with time, t.
The three terms on the RHS correspond to, in order, linear drift, nonlinear drift, and
diffusion. In equation (2.46), the two regions of time become apparent by comparing
the coefficients of the nonlinear drift and diffusion terms, namely 4σ 2 and 2 /σ 2 . The
ratio of the two coefficients is,
4σ 4
.
2

(2.47)

In the initial phase of the response, this ratio is O( 2 ) (assuming σ0 is O(1)), implying
that diffusion dominates the nonlinear drift. At later times, the ratio is O( −2 ),
owing to the exponential growth of σ, and the nonlinear drift dominates diffusion. In
the intermediate region, the ratio is O(1), and thus σ 2 ∼ O( ). Delayed bifurcation
dominates the system response when this overlap region occurs beyond the bifurcation
point. Now, an upper bound for σ 2 at the bifurcation point, i.e., when µ = 0, is
2π 2 / min µ. Thus, the condition 2
˙

µ ensures that the overlap region occurs after
˙

the bifurcation point and the system will exhibit delayed bifurcation.

In the initial time region we ignore the nonlinear drift with respect to diffusion
and obtain an approximate solution, ρl , of equation (2.46) by solving,
2
∂ρl
= −2µρl +
∂τ
σ2

y
42

∂ρl ∂ 2 ρl
+
∂y
∂y 2

.

(2.48)

2
Taking a Gaussian initial condition with zero mean and variance σ0 , ρl is given by

ρl (y, τ ) =

1
2πσ 2

exp −

y2
.
2

(2.49)

At the other end of the process, in the second time region, which leads to escape,
diffusion can be ignored with respect to drift. An approximate solution, ρnl , of
equation (2.46) for this case is obtained by solving,
∂ρnl
∂ρ
= −2µρnl − 12σ 2 y 2 ρnl − 4σ 2 y 3 nl .
∂τ
∂y

(2.50)

This equation is solved using the method of characteristics, which yields,
t
y1 3
exp −2
µ dt
y
t1
−1/2
t
2 dt
2
,
8σ
y1 (y, t|t1 ) = y 1 + y
t1
(1)
ρnl (y, t) = ρ (y1 , t1 )
nl

,

(2.51)

(2.52)

(1)
(y , t ) captures the constant of integration along each characteristic. This
nl 1 1
corresponds to the initial conditions for equation (2.50). Now, time t1 need not be
where ρ

the starting time of the sweep, but can be taken to be any time before escape. To
match the linear and nonlinear solutions, we take t1 to be some time in the common
region where the quantity in equation (2.47) is O(1). To do the matching, we choose
(1)
ρ
such that ρnl and ρl are identical to leading order in the common region.
nl

We begin by rewriting equation (2.49) as,

ρl = √

t
t 2
1
y2
exp − − 2
µ dt −
dt
2
2πσ(t1 )
t1
t1 σ 2
43

.

(2.53)

In the common region, the last integral is O( ), and it can be ignored, resulting in,
t
1
y2
ρl ≈ √
µ dt
exp − − 2
2
2πσ(t1 )
t1

.

(2.54)

Note that this is valid because t is not too far from t1 as both are in the common
region. By the same argument, we drop the integral in equation (2.52). Thus, y1 ≈ y
and the nonlinear solution is approximated by,
t
(1)
µ dt
ρnl ≈ ρ (y, τ1 ) exp −2
nl
t1

.

(2.55)

2
y1
1
(1)
ρ (y1 , t1 ) = √
exp −
.
nl
2
2πσ(t1 )

(2.56)

The two distributions are matched by taking,

It only remains to calculate the probability current at ±∞, the escape current.
The escape current comes entirely from the nonlinear drift which is responsible for
escape to ∞ in finite time. Thus, the current is given by P∞ (T ) = limx→∞ 8x3 ρ.
This gives
8σ 2 (T )
√
P∞ =
2π

T
8
t1

−3/2
σ 2 dt


exp 

T 2
t1

σ2

T
dt −

16
t1

−1
σ 2 dt


 .(2.57)

Equation (2.57) depends weakly on the choice of t1 , and to leading order, this dependence can be removed, as follows. First, since 2 /σ 2 is small for times after t1 ,
we drop the first integral in the exponent. Second, since σ 2 is small for times before
t1 , we extend the lower bound of integration for the remaining integrals back to the
bifurcation point, which is a well-defined time and makes a convenient choice. The
44

resulting approximation is given by,
8σ 2 (T )
√
P∞ (T ) ≈
2π

T
8

−3/2
σ 2 dt


exp − 16

tb

T

−1
σ 2 dt

tb


 , (2.58)

where tb is the time at which the bifurcation point is reached.
To solve a specific example, we again consider the linear sweep, for which µ = t.
In this case,

σ2 = 2

π 2t2
e
2

4
√
√
2σ0 −2t2
0 + erf( 2t) − erf( 2t0 ) .
e
π

(2.59)

For large |t0 |, that is, starting far from the bifurcation point, σ becomes approximately
independent of σ0 and t0 . Under this assumption, we take the approximation,
σ2 ≈ 2

√
π 2t2
e
1 + erf( 2t) .
2

(2.60)

For µ = t, the bifurcation point is reached when t = tb = 0. Thus, the integration of
σ 2 gives,
T
8

σ 2 dt = 2π 2

0

√
4T 2
erfi( 2T) +
F (2T 2 ) ,
π 2 2

(2.61)

where erfi is the imaginary error function, and 2 F2 is the hypergeometric function
2
2 F2 (1, 1; 3/2, 2; 2T ) [1]. Thus, the escape distribution for a linear sweep is approximated by,
√
2(1 + erf( 2T ))
P∞ (T ) ≈ √
√
2π (π erfi( 2T ) + 4T 2 2 F2 (2T 2 ))3/2
× exp −

(2.62)

1
√
.
2 (π erfi( 2T ) + 4T 2 F (2T 2 ))
4
2 2

Figure 2.6 illustrates the first three moments and two example distributions of
45

100
•
•

T

0

•

Var(T )

2

-2

•

•

10

•
•

1
•

0.1

•

•

-4
1
10−4 0.001 0.01 0.1
Noise to Sweep Ratio, 2
(a)

10

0.08

-1
-2
10
10−410−3 0.01 0.1 1
Noise to Sweep Ratio, 2
(c)

102

0.06

0.6
0.4

0.02

0

0.8

0.04

1

P∞

Skew(T )

2

0.01
10 102
10−410−3 0.01 0.1 1
2
Noise to Sweep Ratio,
(b)
2 = 50
2 = 0.01
1.0
0.10

0.2

-60 -50 -40 -30 -20 0
T
(d)

1

2
T

3

4

Figure 2.6: Mean (a), variance (b), skewness (c), and example distributions (d) of
the distribution of escape times away from a swept subcritical pitchfork bifurcation.
Numeric solution (solid lines) and approximations (dashed) for fast and slow sweeping
are shown. Dots represent data points from Monte-Carlo simulations of the nonlinear Mathieu equation, equations (C.17) and (C.17), swept through the subcritical
pitchfork bifurcation.

46

the escape current for a linearly swept subcritical pitchfork bifurcation over a wide
range of noise intensities computed using numerical and approximate solutions of the
FP equation. Solid lines indicate the solution obtained by numerically solving the
FP equation over a finite domain. The moments for the adiabatic approximation
(large noise) and the delayed bifurcation (small noise) approximation developed in
this paper, given in equation (2.62), are shown as dashed lines. The dots represent
data points obtained by Monte-Carlo simulation of a nonlinear Mathieu equation
near the subcritical pitchfork bifurcation point. In appendix C we illustrate how this
system reduces to the one-dimensional model employed in this analysis. An interesting
feature of the distributions is that the skewness is positive for delayed bifurcations and
negative for activated escape, indicating that the sign of the skewness is an indicator
of the sweep rate relative to the noise strength.

2.5

Subcritical Pitchfork Bifurcation in a MEMS
Resonator

In an experimental setting, it is common to interrogate a nonlinear resonator by
slowly sweeping its drive frequency through resonance. This practice captures a portion of the frequency response curve from which device parameters can be fit. Given
the discussion in the previous sections of this chapter, we must acknowledge that
the branch end points, i.e., the bifurcation points, are not faithfully captured by this
experimental technique. To examine this further and compare equation (2.62) with
experimental data, in this section we consider a set of sweeping experiments carried
out by Chris Burgner at the University of California, Santa Barbara, under the supervision of Prof. Kimberly Turner [16]. These experiments where done on a MEMS
gyroscope which exhibits parametric resonance. A scanning electron microscope image of the device is shown in figure 2.7 (a). Figure 2.7 (b) shows the amplitude and
47

phase, (rad) amplitude, (µm)

2
1.5
1
0.5
0
0

5

10

15

5

10

15

2
0
-2
0

(a)

time, (s)
(b)

Figure 2.7: (a) SEM of MEMS gyroscope [89]. (b) Amplitude and phase of drive axis
vibration during an example escape event under parameter sweep. For interpretation
of the references to color in this and all other figures, the reader is referred to the
electronic version of this dissertation.

phase of the drive mode of the device during an escape event while the forcing frequency is swept. The device, described in more detail in [89], was employed for this
work because it was conveniently available and exhibits a subcritical pitchfork bifurcation in its subharmonic resonant response. The gyroscope possesses two nominally
orthogonal vibration modes, the drive mode and the sense mode. The wide frequency
range of the resonant response of this nonlinear gyro allows these modes to be significantly mistuned and so in our experiment, which focuses on the drive mode, the
sense mode can be ignored. We model the drive mode by
2
z + 2Γz + ωo (1 + δ1 + δ1 cos 2ωt)z + γ(1 + δ3 + δ3 cos 2ωt)z 3 = ξ(t) (2.63)
¨
˙
ξ(t) = 0

ξ(t)ξ(t ) =

4ΓkB T
m δ(t − t )

(2.64)

where Γ is the damping rate, ω0 is the natural frequency, ω is the drive frequency,
δ1 is the normalized variation of the linear stiffness due to parametric forcing, γ
is the nonlinear stiffness coefficient, δ3 is the normalized variation of the nonlinear
48

stiffness due to parametric forcing, m is the effective mass of the resonator, T is the
temperature, and ξ is a white Gaussian thermal noise. The reduction of this model
near the subcritical pitchfork is carried out in appendix C following the method
presented in [62]. The reduced model is
√
dq
= 2µq + 4q 3 + Dξ(τ ),
dτ
ξ = 0,

(2.65)

ξ(τ )ξ(τ ) = 2δ(τ − τ ),

where
τ = Γt
µ =
r =
D =

α =
β =

(2.66)

1
(−1 + β 2 − α2 ),
2
dδ
1 dµ
1
dω
=
(2α + β)ω0 1 − 4α
,
2
Γ dt
dt
dt
8Γ
kB T γ αβ 2 (3 + 2δ3 ) + δ3 β(β + α)2
,
32ω 3 Γ
(β 2 − α2 )3/2
2
ω 2 − ω0 (1 + δ1 )
ω δ
σ
≈− 0 1 + ,
2ωΓ
2Γ
Γ
2
δ1 ω0
ω δ
≈ 0 1.
4ωΓ
4Γ

(2.67)
(2.68)
(2.69)

(2.70)
(2.71)

In the reduced model, the dynamic variable q is a linear combination of averaged
quadrature components of z and z. The details of the transformation are left to ap˙
pendix C. τ is a scaled time variable, µ is the bifurcation parameter which depends on
the device parameters, r is the non-dimensional sweep rate of the bifurcation parameter assuming that the swept parameters are δ1 and/or ω, D is the non-dimensional
effective strength of the noise on the slow manifold, α is the effective detuning between
the drive and natural frequencies, and β is the effective parametric forcing strength.
Note that the effective noise strength is also a function of the detuning and forcing
strength. To capture the leading order term, these parameters should be evaluated at
2
2
the bifurcation point where βb = αb + 1. The non-dimensional noise to sweep ratio
49

is given by
2
αb βb
2 = D = 3ΓkB T γ
.
3
˙
r
˙
4ω0 (2αb + βb )ω0 δ1 − 4αb ω

(2.72)

In this expression we have dropped terms proportional to δ3 for convenience since
δ3
1 these terms don’t contribute much to the value of 2 . Note the overdots in
equation (2.72) denote derivative with respect to the lab frame time, not the nondimensionalized scaled time. The non-dimensional parameters αb and βb appear up
in equation (2.72) because the dynamics near the bifurcation point depend on the
system parameters. These values have been given the subscript b to denote that they
should be evaluated at the bifurcation point encountered during the parameter sweep.
2
2
Again, this is the point along the parameter sweep trajectory where βb = αb + 1.
˙
When the frequency is swept along, i.e. δ1 = 0, equation (2.72) simplifies to
2
2 = D = 3kB T |γ|δ1 ,
r
256ω0 Γ|ω|
˙

(2.73)

which indicates how 2 depends on physical parameters.

A frequency response for the device is shown in figure 2.8. The device was placed
in a vacuum chamber at 1mTorr and monitored via a laser vibrometer. At this vacuum
level the quality factor of the resonator is Q ≈ 3000 and the natural frequency is near
8.485kHz. An electrostatic potential of 15V was placed across the non-interdigitated
combs seen in figure 2.7 (a) and the frequency of the drive voltage was swept downward
from well above the subcritical pitchfork bifurcation point, which is seen in figure 2.8
to be near 8485Hz. Sweep rates between 0.005 and 0.3 Hz/sec were used and 1000
escape measurements were performed for each sweep rate to generate the family of
escape distributions shown in figure 2.9. A single example escape event is shown
in figure 2.7 (b).

This set of escape distributions was taken with ambient noise.
50

Displacement, µm
Normalized
Eigenvalues

1
0.5
0
8478

8480

8482 8484 8486
Frequency Hz

8488

5

8490

Re
Im

0
-5
8478

8480

8482 8484 8486
Frequency Hz

8488

8490

Figure 2.8: (Top) Frequency response for the MEMS resonator. Blue data points
track the sweep up and red data points track the sweep down. The green data points
represent the mean jump frequency for various sweep rates. Black lines represent
fitted stable (solid) and unstable (dashed) frequency response branches. (Bottom)
Real and imaginary parts of the eigenvalues belonging to the fixed point at zero.
Dashed line represents the linear sweep approximation.
Additional noise was added with base excitation using a piezo element, in order to
reach smaller values of the sweep rate to noise ratio. Recall that the normalized mean
and variance of escape are functions of this ratio when the sweep is linear in time.
The inset of figure 2.9 shows three different distributions that share the same value
of r/D in order to illustrate that, while the normalized statistics depend only on the
ratio, the true statistics depend on both values. The experimental distribution shown
in the inset has had its mean shifted for comparison. The reason this was required will
be discussed presently. The primary attribute we notice about the measured escape
distributions is that they all appear to be described by delayed bifurcations. Despite
sweeping so slowly that a single escape data point took several minutes to acquire,
activated escape over the barrier was not observed. The skew of the distributions was
expected to confirm this, but the data was inconclusive. This may be because the
data may have been somewhat contaminated by the drifting parameter values of the
51

35
20

Shifted
experimental
distribution.

30
15

0.005 Hz/sec

25
10
20
5
15
0
10
5

0

0.1 0.2
0.3
Bifn Parameter, µ

0.4

0.3 Hz/sec

0
8481 8481.5 8482 8482.5 8483 8483.5 8484 8484.5
Frequency (Hz)
Figure 2.9: Experimentally measured distributions of escape events for ambient noise
and sweep rates ranging from 0.005 to 0.3 Hz/sec. The inset shows one such distribution, artifically shifted, and compared against the theoretical distribution. The
purpose of the shift is discussed in the text. The three distributions in the inset
correspond each to different sweep rate and noise intensities yet all three with the
same ratio. This illustrates that, while the normalized statistics depend only on the
ratio, the true statistics depend on both values.

52

resonator. The drift observed was noted to be periodic and correlated to the HVAC
cycles of the lab environment. In the data presented here the drift was compensated
for by measuring the natural frequency prior to each sweeping experiment.

Figure 2.10 shows the normalized mean and variance of the escape distributions
plotted against the sweep rate to noise ratio. The results show qualitative agreement
with the result in equation (2.62). We believe the reason for the difference is a
deviation from the assumptions employed in the development of equation (2.62).
The assumption of a linear sweep requires that the eigenvalue of the bifurcating
fixed point change linearly in time. Figure 2.8 shows the quasistatic values of the
eigenvalues of the fixed point at zero amplitude during the sweep. The green dots
show the mean frequency of escape. It is clear that the linear approximation, shown
as the dashed line, is a poor representation of the actual eigenvalue. The eigenvalue
is smaller than anticipated and so the response does not grow as fast as expected.
Hence the mean occurs later than expected, as seen in figure 2.10(a). This is the
reason the example curve in the inset of figure 2.9 was shifted for comparison with
the theoretical distribution. The variance is also larger than predicted for most data
points. This may also arise from the failure of the linear sweep assumption, since
the eigenvalue was smaller than expected it took more time for the system to escape.
During this added time the system continued to diffuse, resulting in a larger spread
of escape times. In addition, over this wide range of parameter variation the effective
noise strength on the center manifold will depend on the bifurcation parameter. Thus,
while equation (2.58) could be employed to capture the nonlinear relationship between
the eigenvalue and the forcing frequency, it does not account for the change in effective
noise intensity.
53

×104
Normalized Variance

Normalized Mean

6
5
4
3
2
1
0

100
105
Sweep rate to Noise Ratio, D/r
(a)

106
104
102
100
100 102 104 106 108
Sweep rate to Noise Ratio, D/r
(b)

Figure 2.10: Normalized mean (a) and variance (b) of the bifurcation parameter at
escape. The solid line represents the mean of the distribution for a linear sweep
predicted from the normal form, as given in equation (2.62). The data points were
calculated from the measured distributions shown in figure 2.9.

2.6

Outlook

As mentioned at the beginning of this chapter, this work was motivated by applications to sensing, parameter fitting, and system identification. Escape experiments,
used for these applications, employ highly visible transients which can be useful in
very small systems that are difficult to measure [58]. However, when compared to
more conventional methods, measurement by collecting escape event data is quite
slow [17]. The reason is that there are three time scales acting in the system and
activated escape near a bifurcation point operates on the slowest time scale. Moreover, many escape events must be observed (with the exception of single-shot qubit
readout [78]) to obtain an estimate of the escape moments. This brings the added
complication of parameter drift over such long time scales. If bifurcation location
or tracking is to be used for sensing or system identification in a device that can
be monitored with reasonable precision, then feedback control [17] or experimental
continuation [99] may be more appropriate techniques.

54

Chapter 3
The Dynamical Bridge
In this chapter we discuss the use of noise activated escape in a parametric sensor. The past decade has seen significant advances in the employment of micro- and
nano-electromechanical (M/NEMS) resonators as parametric sensors. The small size
and high quality factors of these devices enables very sensitive coupling to and/or
detection of adsorbed masses on the resonator surface [40, 42, 121, 33, 92, 43, 66],
superconducting qubit states [23, 68, 112], and electron transport statistics [6, 5], to
cite a few examples. Many of the implementations of NEMS resonators for parametric sensing mirror those of Josephson junction circuits for the same applications
[55, 76, 88, 114, 70, 105, 78, 116, 96]. Accordingly, the sensing paradigms are very
similar to the methods developed in the 1970s to measure the Josephson junction
critical current [67, 58, 47].
The current biased Josephson junction is analogous to the simple pendulum where
the critical current is likened to the coefficient of the restoring torque (mg in a
pendulum) and the phase difference, δ, between the superconductor order parameters
across the junction is likened to the pendulum angle. For small δ the junction is
modeled, like many NEMS resonators, by Duffing’s equation. Accordingly, it is not
surprising to find many common experimental methods associated with both systems.
55

The measurement paradigm we highlight here is, roughly speaking, the use of
such a Duffing resonator as a threshold detector. While there are several variations
on this concept, the basic principle is as follows. The system is placed in the bistable regime and brought close to a bifurcation (the threshold). Noise in the system
blurs the threshold and so the system is placed sufficiently far from the threshold
so that noise only occasionally pushes the system over the edge. This occurs when
a rare large outburst of noise pushes the system out of the domain of attraction
belonging to the occupied attractor. When this happens the system is said to have
escaped. The probability per unit time, or rate, of escape is exponentially sensitive
to the distance from the operating point to the threshold and, by extension, to the
system parameters. Accordingly, monitoring the time until escape, and performing
many such escape experiments, allows one to calculate the mean escape rate and thus
determine parameter variation in the resonator with high precision. The threshold
position can be thought of as being detected through the measurement of the escape
rate. The escape events monitored in this method are often highly visible transient
dynamics of the resonator. This can enable a very precise measurement to be made
even when the detailed motion of the resonator is difficult to resolve, c.f. [58].
A single-shot variation of this method is employed in qubit readout, relevant for
quantum computing [116, 78]. If the state of the qubit, which is coupled to the resonator, displaces the bifurcation sufficiently, then an escape event can be strongly
correlated with the qubit state. The qubit state can then be inferred from the occurrence, or lack thereof, an escape event. This detection scheme has become known as
the Josephson bifurcation amplifier [116].
The phenomenon of noise activated escape, on which these measurement methodologies are based, has been studied in great detail. The escape rate has been calculated
for Gaussian noise [64], Poisson noise [31] and non-Gaussian perturbations superimposed on Gaussian noise [105, 11, 10, 12]. These and similar results have application
56

in the measurement of electron transport statistics [7, 55, 53, 88, 114]. It is also
known that the escape events occur along trajectories forming narrow tubes in phase
space [21, 20]. This gives experimental credence to the method of calculating the
escape rate by considering the most probable path.

In this chapter we discuss a variant on the detection schemes described above,
called the balanced dynamical bridge [30]. In the dynamical bridge paradigm a
bistable system subject to weak noise is arranged such that the probability to be
in either state is equal. A perturbing parameter variation or change in the noise
statistics can then be measured by the resulting shift in relative state populations,
or by the force required to rebalance the system. The balanced dynamical bridge
is so named because it is analogous to a bridge circuit in which the voltage across
the bridge connection is tuned to zero in order to measure, for example, a resistance.
However, here we consider a dynamical system driven by noise, hence the “dynamical”
designation.

Our discussion of the balanced dynamical bridge is organized in the following
way. In section 3.1 we review the theory of noise activated escape, treating both
Gaussian noise and a non-Gaussian perturbing process. In section 3.2 we discuss the
dynamical bridge as a general use detector, formulating expressions for the sensitivity
and measurement time of the device. In section 3.3 we discuss the dynamical bridge
as a detector of non-Gaussian noise. In section 3.4 we discuss the Duffing resonator
as an example of a prototypical bistable system. We present the locus of operating
points in a two-dimensional parameter space at which the bridge is balanced. In
section 3.5 we give an example of measuring non-Gaussian noise by considering a
one-dimensional dynamic bridge subject to weak shot noise. Concluding remarks and
a discussion of future work are given in section 3.6.
57

3.1

Activated Escape

In section 2.2 we reviewed activated escape by white Gaussian noise in one dimension.
In this section we review the general theory of noise activate escape due to additive
Gaussian noise perturbed by weak non-Gaussian noise; see [11]. This theory lays
the groundwork for our subsequent discussions of bridge sensitivity and detection of
non-Gaussian noise.

We begin with a general nonlinear system

˙
q = K(q) +

√
Df (t) + G(q)ξ(t),

(3.1)

where q ∈ Rn is the system state vector, K is the nonlinear vector field, f is a
stationary Gaussian noise vector with characteristic strength D, G ∈ Rn×m is the
non-Gaussian noise matrix coefficient, and ξ ∈ Rm is the non-Gaussian noise vector,
assumed to be independent of f . Suppose that the noise free system possesses two
stable fixed points, q a and q b , and a saddle point, q s , on the boundary between
the basins of attraction for q a and q b . If the noise is weak, then the system will
likely settle down to the vicinity of one stable fixed point, say q a . It may happen
subsequently that a large rare burst of noise pushes the system out of this basin and
into that belonging to q b . Suppose we know the path by which the system exits
the basin and the accompanying realization of ξ, then the Gaussian noise realization
required to make this transition follows from equation (3.1),
˙
f (t) = D−1/2 q − K(q) − G(q)ξ(t) .
58

(3.2)

The probability for such a noise outburst to occur is [46]
1
P [f (t)] ∝ exp − 2

ˆ
f T (t )F(t − t)f (t) dt dt ,

ˆ
f (t)f T (s) F(s) ds = Iδ(t),

(3.3)
(3.4)

where F is the inverse pair correlator of ξ(t). We assume that f possesses time
ˆ
ˆ
ˆ
reversal symmetry so that F(t) = F(−t) = F T (t). From equation (3.2) we see that
the integrand of (3.3) is inversely proportional to the characteristic noise strength.
Thus, for weak noise, the probability for a transition to occur is sharply peaked
about the most probable switching path, i.e., the path that minimizes P [f (t)] in
equation (3.3) where f is given by (3.2). The mean switching rate is approximately
proportional to the probability for the realization of the most probable switching
path. If ξ = 0, the rate is

r0 ∝ exp [−S0 /D]
1
S0 = min 2

(3.5)

ˆ
˙
˙
(q − K)T F(t − t)(q − K)t dt dt
t

(3.6)

˙
˙
where (q − K)t denotes q(t) − K(q(t)) and the minimum is taken over all paths that
take the system from q a to q b .
ˆ
When f (t)f T (s) = Iδ(t − s), then F(t) = Iδ(t) and S0 simplifies to
S0 = min

1
2

˙
|q − K|2 dt ,

(3.7)

which is the action belonging to an auxiliary Hamiltonian system
1
H = |p|2 + pT K(q),
2

˙
p ≡ q − K(q).

(3.8)

The costate, p, in this case, is proportional to the optimal noise realization, and the
59

constant of proportionality is

√

D. The path minimizing S0 in equation (3.7) is the

made up of two parts. The first part is the heteroclinic trajectory going from q a to
q s along which p = 0, since this transition requires noise. The second part is the
heteroclinic trajectory going from q s to q b along which p = 0, since this response is
noise free. This second part of the switching trajectory does not contribute to the
calculation of the switching probability and so our focus is on the first (noisy) part.
If f and ξ are independent noise processes and ξ is small compared to f , then
the mean switching rate can be expressed as [11].

r ∝ exp [−S[ξ]/D] ξ

(3.9)

where, to leading order, S[ξ] is
χ(t)T ξ(t) dt,

(3.10)

ˆ
˙
q(t ) − K(q(t )) T F(t − t)G(q(t))dt .

(3.11)

S[ξ] = S0 +
χT (t) = −

The quantity χ can be thought of as the susceptibility of the activation energy to a
perturbing force. The switching rate in the presence of this perturbing noise is

r ∝ Ar0 ,
1
A = exp − D

χT ξdt

(3.12)
= Φξ [iχ/D],

(3.13)

where Φξ is the characteristic functional of ξ.

3.2

Dynamical Bridge Sensitivity

We now turn to the dynamical bridge. In this section we set ξ = 0 and discuss the
bridge as a general use detector. We calculate the sensitivity and measurement time
60

assuming the rates developed in the previous section are known.

Equation (3.1) serves as the microscopic model for the dynamic bridge. The bridge
itself functions on a much longer time scale. For times scales exceeding that of the
switching dynamics and correlations of f , the coarsened dynamics may be modeled
by the master equation
˙
P1 = −r1 P1 + r2 P2 ,

(3.14)

˙
P2 = r1 P1 − r2 P2 ,
where P1 and P2 are the probabilities to find the system in the first or the second
states respectively, and ri is the switching rate out of the ith state, i = 1, 2. Normalization of the probability distribution requires P2 = 1 − P1 and so we will work with
P1 only for this analysis. The solution to equations (3.14) is
P1 (t) =

r2
1 − e−(r1 +r2 )t + P1 (0)e−(r1 +r2 )t ,
r1 + r2

(3.15)

from which we see that the characteristic relaxation time of the system is tr ∼
(r1 + r2 )−1 .

To make a measurement with the dynamic bridge we must measure the probability
P1 . This may be done by measuring some quantity x that takes on the value X1
when the system is in state 1 and X2 when the system is in state 2. The probability
to be in state 1 can be estimated from N measurements of x, where the system is
given time to relax between each measurement. Thus, the total measurement time is
on the order N tr and gives the probability estimate
1
ˆ
P1 =
N (X1 − X2 )

N
(xi − X2 ).
i=1

61

(3.16)

ˆ
The mean and variance of the estimator P1 are
ˆ
P1
ˆ
¯
(P1 − P1 )2

¯
= P1 ,
¯
¯2
P1 − P1
=
,
N

(3.17)
(3.18)

¯
where P1 = r2 /(r1 + r2 ) is the steady state value of P1 . Equations (3.17) and (3.18)
ˆ
¯
show that P1 is an unbiased estimator of P1 with uncertainty proportional to N −1/2 .

Now, suppose that some parameter of the bridge system, λ, changes by an amount
¯
∆λ. The sensitivity of the dynamic bridge is given by ∂λ P1 . In calculating this
quantity, we ignore the dependence of the rate prefactor on λ under the assumption
that D is small. Thus the dominant change in the transition rates comes from the
change in activation energy with respect to λ. Under this assumption, the bridge
sensitivity is given by
¯
r1 r2
∂ P1
=
∂λ
D(r1 + r2 )2

∂(S1 − S2 ) S1 − S2 ∂D
+
∂λ
D
∂λ

,

(3.19)

where Si is the activation energy (equation (3.6)) to escape from the ith state, i = 1, 2.
Equation (3.19) applies for any r1 and r2 . The most sensitive operating point is at
the balance point where r1 = r2 and S1 ≈ S2 . In this case equation (3.19) simplifies
¯
to ∂λ P1 ≈ D−1 ∂λ (S1 − S2 ), which quantifies how the dynamical bridge becomes
increasingly sensitive with decreasing noise.

¯
In order to detect a change ∆λ, we must measure P1 to sufficient accuracy. In
¯
order to see a shift in the occupation probability, the shift in P1 should be on the
order of the error in our estimate. Thus, we require
¯
∂ P1
ˆ
¯
( P1 − P1 ) 2 ∼
∆λ
∂λ
62

(3.20)

ˆ
The variance of P1 , equation (3.18), is bounded by 1/4N , so we use this value to
obtain an upper bound for the number of measurements that must be made on the
bridge. Solving for N and multiplying by the relaxation time gives the measurement
time

N tr =

(r1 + r2 )3 D2
2 2
4r1 r2 ∆λ2

≈ r

∂(S1 − S2 ) S1 − S2 ∂D −2
+
,
∂λ
D
∂λ

2
D
,
∆λ∂λ (S1 − S2 )

at the balance point.

(3.21)
(3.22)

From expressions (3.21) and (3.22), we find that the measurement time increases
exponentially with decreasing noise, owing to the exponential dependence of r1,2 on
D.
We conclude this section by again considering the case where f (t)f T (s) =
ˆ
Iδ(t − s). This gives F(t) = Iδ(t) and ∂λ Si can be written as a simple integral along
the most probable switching path.
∂S1,2
∂λ

=

∂K(q)
˙
q − K(q) T
dt,
∂λ

(3.23)

where q, for S1,2 , is evaluated along the heteroclinic trajectory taking the system
from q a,b to q s .

3.3

Detection of Non-Gaussian Noise

The dynamical bridge can also be used as a sensitive detector of non-Gaussian noise.
Suppose we have the system (3.1) and we have the ability to turn ξ on and off.
Then we can use the the system as a dynamical bridge to determine a degree of
non-Gaussian character in the noise ξ. The ability to turn ξ on and off allows us to
remove the mean by examining the mean fluctuation of the coordinate q about the
63

fixed points q a or q b . A bias can then be applied to compensate for the mean of ξ,
which is convenient for the development, as described below.
In general, the addition of the perturbing noise ξ will change the transition rates,
see section 3.1, and this will shift the relative occupancy ratio in the bridge. The
ratio of probabilities to be in state 1 and 2 is
(0)
P1
r2
A2 r2
=
=
,
P2
r1
A1 (0)
r1

(3.24)

(0)
where ri is the switching rate out of the ith state when ξ = 0. Balancing the bridge
(0)
(0)
before exposure to the perturbation sets r1 = r2 . We now introduce the random
variable
zi =

dtχT (t)ξ(t),
i

(3.25)

where χi (t) is given in equation (3.11) and the subscript denotes escape from the ith
(i)
state. Let κn be the nth cumulant of zi . From equation (3.13), the pre factors, Ai ,
are the moment generating functions for zi and the logarithm, ln Ai , is the cumulant
generating function. Thus,
∞
ln Ai =
n=1

(i)
κn
.
n!(−D)n

(3.26)

Note, if ξ is Gaussian, then zi is also Gaussian and κn = 0 for n > 2. Of course,
here we assume ξ is non-Gaussian, so we keep all cumulants. From equation (3.26),
the probability ratio becomes
P1
= exp
P2

∞
n=1

(2)
(1)
κn − κn
.
n!(−D)n

(3.27)

We have mentioned above how the mean of ξ can be removed by an applied bias. If
64

(1)
(2)
the second cumulants κ2 and κ2 can be balanced, they too can be removed from
equation (3.27). In this case, the population ratio will depend only on higher order
cumulants, i.e., on the non-Gaussian character of ξ. The second cumulant of zi can
be expressed as follows
(i)
κ2 =

dtdt χT (t) ξ(t)ξ T (t ) χi (t ).
i

(3.28)

(1)
(2)
From equations (3.11) and (3.28) it can be shown that κ2 and κ2 will be equal if
ˆ
ˆ
ˆ
dsds F (s − t)G(q i (s)) ξ(s)ξ T (s ) GT (q i (s ))F (s − t ) = cF (t − t ),

(3.29)

where q i is evaluated along the heteroclinic escape trajectory out of state i and c is
some constant, equal for i = 1, 2. This condition is satisfied if, for example, f is a
vector of uncorrelated Gaussian white noise processes, ξ is a vector of uncorrelated
Poisson noise processes, and G ∝ I. If the mean has been removed, and if the second
cumulants cancel, then the population ratio of the bridge will reflect a measure of
the non-Gaussian character of ξ, which is captured in the series of cumulants of z1,2 ,
that is, by the projection of the noise realizations onto the switching susceptibility.
This will be illustrated by an example in section 3.5. Note that while the second
cumulants drops out of equation (3.27), they remain in equation (3.26) and thus
change the switching rate. In this case, the second cumulant can be lumped together
with that of the Gaussian noise, f , and be understood as a renormalization of the
Gaussian noise strength.

65

3.4

The Duffing Resonator Bridge

In this section we discuss the lightly damped Duffing resonator and its use as a
dynamic bridge. We map out the bistable region of parameter space and numerically
determine the locus of points in this space where the bridge is balanced, i.e., r1 = r2 .
The resonator, subject to a periodic driving force and additive noise, is described by
2
ˆ
y + 2Γy + ω0 y + γy 3 = h cos ωt + f (t),
¨
˙

(3.30)

where, y is the oscillator coordinate, ω0 is its natural frequency, Γ is the friction
ˆ
coefficient, (Γ
ω0 ), γ is the nonlinearity parameter, and f (t) is Gaussian white
ˆ
ˆ ˆ
ˆ
noise with f (t) = 0 and f (t)f (t ) = 2Dδ(t − t ). We assume near resonant forcing
2
of amplitude h and frequency ω with ω 2 − ω0

ω 2 . It is convenient to switch

from y to a slowly-varying complex amplitude, u, using the well-known Van der Pol
transformation

y(t) =

2ωΓ iωt
ue
+ c.c.,
3|γ|

(3.31a)

ueiωt + c.c. = 0.
˙

(3.31b)

We apply this transformation to equation (3.30), neglect fast oscillating terms, and
move to a nondimensional time scale, τ , to obtain
u = −(1 + iΩ)u + i sgn(γ)|u|2 u − i
˙

β − if (τ ),

(3.32)

where the dot, (˙), is now interpreted as the derivative with respect to τ and the following nondimensional parameters have been defined, along with the nondimensional
66

noise f (τ ),

τ = Γt,
2
ω 2 − ω0
Ω =
,
2ωΓ
3|γ|h2
β =
,
32ω 3 Γ3
3|γ| −iωτ /Γ ˆ
f (τ /Γ),
e
8ω 3 Γ3

f (τ ) =

(3.33)
(3.34)
(3.35)
(3.36)

where f (τ ) refers to f (t) in (3.1). By neglecting the fast oscillating moments, we
approximate the real (fr ) and imaginary (fi ) parts of f as two independent white
Gaussian noise processes with [35]
fr (τ ) = fi (τ ) = fr (τ )fi∗ (τ ) = 0,

(3.37a)

∗
fr (τ )fr (τ ) = fi (τ )fi∗ (τ ) = Dδ(t − t ).

(3.37b)

where D is the nondimensional noise strength

D=

ˆ
3|γ|D
.
8ω 3 Γ2

(3.38)

In the absence of noise (f (τ ) = 0), the fixed point are found by setting equation
(3.32) to zero. Thus, the fixed points must satisfy
√

u=

β
.
sgn(γ)|u|2 − Ω + i

(3.39)

Taking the norm of (3.39) gives a cubic equation in |u|2 ,
|u|6 − 2 sgn(γ)Ω|u|4 + (Ω2 + 1)|u|2 − β = 0.
67

(3.40)

The fixed points are associated with the the positive real roots of this equation. The
bistable region is bounded by two saddle-node bifurcations. The amplitudes of u at
which these bifurcations occur can be found by taking the derivative of (3.40) with
respect to |u|2 and solving for the bifurcation value of |u|2
2
|u|2 = sgn(γ)Ω ±
B 3

Ω2 − 3,

(3.41)

where the subscript B denotes a saddle-node bifurcation value. This solution is
substituted back into (3.40) to determine the boundaries of the bistable region in
parameter space, which can be expressed as

βB =

2
sgn(γ)Ω(Ω2 + 9) ± (Ω2 − 3)3/2 ,
27

Ω2 ≥ 3.

(3.42)

This bistable region is illustrated in figure 3.3 as the region between the solid lines.
Plotting with axes Ω−2 and βΩ−3 allows us to visualize the bistable region as a
closed region of parameter space. Inside the bistable region the system posses two
stable fixed points and one saddle point. Returning to the notation of section 3.1, we
set the real and imaginary components of u to be the elements of the two component
vector q ∈ R2 . The costate, p ∈ R2 , is proportional to the optimal realization of the
noise quadratures; the real and imaginary components of f (τ ). The most probable or
optimal switching trajectories, as described in section 3.1, are heteroclinic solutions
of the Hamiltonian system, given here by

H=

|p|2
+ (|q|2 − Ω)(q ∧ p) − q T p − p2
2

β,

(3.43)

where ∧ is the wedge product defined as q ∧ p = q1 p2 − q2 p1 , and p2 is understood to
be the second component of the two dimensional vector p. The pair of heteroclinics
for a sample set of parameters is depicted in figures 3.1 and 3.2. In figure 3.1 the
68

1
0
q2 -1
-2
-3

-2

-1

0
q1

1

2

Figure 3.1: The blue and magenta lines depict heteroclinic solutions connecting the
operation points to the saddle point for system (3.43), as projected onto the (q1 , q2 )
plane. Solid lines represent the p = 0 solutions and the dashed lines are the noisefree (p = 0) stable and unstable manifolds of the saddle point. Parameter values:
Ω = 4.706182 and β = 8.501802

heteroclinic is projected onto the (q1 , q2 ) plane. The dashed lines are noise-free stable
and unstable manifolds of the saddle point. The solid lines show the most probable
path by which noise will drive the system from the stable fixed points to the saddle
point. Upon reaching the saddle point, the system will relax to one of the stable
fixed points, most likely following quite closely the corresponding noise-free unstable
manifold.

Figure 3.2 shows the heteroclinic projected onto the (p1 , p2 ) plane. The

trajectories represent the most probable realization of the noise quadratures that gives
rise to a switching event. These trajectories begin and end at the noise-free equilibria,
hence at p = 0.
The heteroclinics shown in figures 3.1 and 3.2 were obtained by the method of
shooting. The parameter values were chosen because the action along these two heteroclinics is equal. Thus, these heteroclinics represent a balanced bridge configuration
(at least up to the difference in prefactors). The locus of points in parameter space
for which the action is equal for both switching directions is plotted as the dashed line
69

1.0
0.5
p2 0.0
-0.5
-1.0
-0.4

-0.2 0.0

0.2
p1

0.4

0.6

0.8

Figure 3.2: Heteroclinic solutions for system (3.43) projected onto the (p1 , p2 ) plane.
For the solution with Ω = 4.706182 and β = 8.501802. Note that these heteroclinics
start and end at noise-free points where p1 = p2 = 0.

in figure 3.3. The balance curve in figure 3.3 was determined numerically. However,
the end points were obtained from results given in [36], where it was shown that (for
γ > 0) on this locus lim −2
βΩ−3 = 0.13 and near the cusp at Ω−2 = 1/3 the
Ω →0
balance curve lies approximately halfway between the two saddle-node branches, and
is given by βbl Ω−3 = 2(1 + 9Ω−2 )/27. These results set the stage for using the
Duffing system with additive Gaussian noise as a balanced bridge.

3.5

Shot Noise Measurement

Near the cusp (Ω−2 = 1/3) the system discussed above collapses onto a one-dimensional
slow manifold. On the balance curve, the system reduces to the classic overdamped
symmetric quartic double well potential. In this section we discuss the measurement
of shot noise statistics by such a system.
Shot noise is characterized by the pulse area, g, and pulse rate, ν. Thus, we
70

0.30
0.25
βΩ−3

0.20
0.15
0.10
0.05
0.00
0.00

0.05

0.10

0.15 0.20
Ω−2

0.25

0.30

Figure 3.3: Two dimensional parameter space for the Duffing resonator. Solid lines
are saddle-node branches which enclose the bistable region. The dashed line is the
balance curve on which the escape actions for two stable states are equal, S1 = S2 .
In the weak noise limit this line is very close to the line where the escape rates are
equal.

suppose our system is described by the equation

q = −U (q) +
˙

√
Df (t) + ξ(t) − νg,

(3.44)

where q is a scalar, f (t) is a white Gaussian noise with f (t)f (t ) = δ(t − t ), and
ξ(t) is a Poisson process. We assume that the shot noise acting on the resonator
is modulated by the resonator driving signal so that it appears as a simple Poisson
process on the slow manifold. Our system is assumed to move in a overdamped
double-well potential, as shown in Fig. 3.4, that can be described by
1
1
U (q) = − q 2 + q 4 .
2
4

(3.45)

The mean ξ = νg has been explicitly removed in order that the bridge, as discussed in section 3.3, will be sensitive to only the higher moments. It is assumed that
this mean can be measured and compensated for. This measurement also amounts
71

g>0
qs

q1

q2

Figure 3.4: Double-well potential of equation (3.45); q1 and q2 denote positions of
stable fixed points, while qS refers to the unstable fixed point. Poisson pulses with
positive area g change the switching rates such that r1 > r2 , resulting in a shift of
the occupation probability from q1 towards q2 .

to the first of the two data points required to calculate ν and g.

The cumulant generating function is known for the Poisson process. For the zero
mean Poisson process under consideration, ξ − νg, it is given by

ln Ai = ν

dt e−gχi (t)/D + gχi (t)/D − 1 .

(3.46)

Equating this term by term to equation (3.26), (repeated here)
∞
ln Ai =
n=1

(i)
κn
,
n!(−D)n

gives
(i)
κ1 = 0,
(i)
κn = νg n

(3.47)
dtχn (t),
i

n > 1.

(3.48)

As discussed in section 3.1, the susceptibility, χi , is a heteroclinic solution of the
72

auxiliary Hamiltonian system. For this system, the Hamiltonian is

H=

p2
− pU (q),
2

p = q + U (q).
˙

(3.49)

The heteroclinic trajectories of interest connect the fixed points (q, p) = (±1, 0) with
the saddle point (q, p) = (0, 0). By inspection, we see that H = 0 along these
trajectories. Accordingly, the heteroclinics are solutions to

q = U (q),
˙

(3.50)

that is, the reversed-time deterministic trajectories. For the potential given in equation (3.45), the solutions are
−1/2
(h)
q± (t) = ± 1 + e2t
,

(3.51)

and the attendant susceptibilities are given by
(h)
χ1,2 = −2q± (t) = ±2e2t 1 + e2t −3/2 .
˙

(3.52)

Using equation (3.50), the integrals in equation (3.48) can be evaluated explicitly.
Applying the result to equation (3.27) gives

P1
= exp −ν
P2

n=1


2n+1
2g
Γ(n + 1/2)
,
D
(2n + 1)Γ(3n + 3/2)

(3.53)

where Γ(n) is the Euler Gamma function. Figure 3.5 shows log(P2 /P1 )/ν as a function of g/D along with Monte-Carlo simulation results. The inset in the figure shows
the same quantities as the primary axes over a larger range. The Monte-Carlo simulations were carried out with parameter values D = 0.04 and ν = 0.5. The weak noise
73

2.0

Log [P2 /P1 ] /ν

10
8
• •
1.5 6
•
4
•
1.0 2
•
0 •
0 1 2 3 4 5 6 7
0.5
•
•
•
0.0
•

•
•
•

•

-0.5
0.0

0.5

1.0
g/D

1.5

2.0

Figure 3.5: Change in probability ratio due to the addition of zero mean shot noise,
equations (3.53). Data are results from Monte-Carlo simulations with D = 0.04 and
ν = 0.5. The inset shows the same quantities over a larger range.
limit assumed by the analysis requires D/2∆U

1, where ∆U is the barrier height.

We also assumed that the non-Gaussian noise was weak compared to the Gaussian
noise; this condition is νg 2
require g/D

D. The parameters used give D/2∆U = 0.08 and

7. Figure 3.5 shows reasonable agreement in this region.

For g/D not too large, the series can be truncated at n = 2. This accounts for up
to the fifth order moments of ξ. The error due to this truncation grows monotonically
with g/D, reaching approximately 2% error for g/D = 6. The truncation, and use of
the relation ξ = νg, give a quadratic equation for g 2 . Solving for g 2 and taking the
positive root gives

g2 = −

715 2
D +
32

715 2 2 225225D5 P1
D
−
ln
,
32
256 ξ
P2

(3.54)

where the label “1” applies to the stable fixed point at q = −1 and the label “2”
applies to the stable fixed point at q = 1. By considering the two cases of positive
and negative g we now show that equation (3.54) gives a positive real result for g 2 .
For positive g, ξ is also positive. Moreover, since the system is overdamped, pulses
74

pushing in the positive direction can only decrease P1 while increasing P2 . Thus,
for positive pulses, g > 0, ln P1 /P2 < 0 and ξ > 0, and so equation (3.54) gives a
positive real result for g 2 from which we take the positive root for g. On the other
hand, when g < 0, we have ln P1 /P2 > 0 and ξ < 0, and equation (3.54) still gives
a positive real result for g 2 from which we take the negative root for g. Generally,
the sign of g is found from the sign of ln(P1 /P2 ).
Equation (3.54) is our final result, from which g is found. Use of this equation
requires measurement of the quantities D, ξ , and P1 /P2 . The Gaussian noise
strength, D, and the non-Gaussian noise mean, ξ , can be measured from the quasi
steady state distribution around a fixed point of the drift vector. Specifically, the
Gaussian noise strength can be found from the width of the distribution about the
fixed point, and the mean can be found from the force required to rebalance the
distribution, as described above, and reiterated here: With ξ turned off, the mean of
this quasi steady state distribution is measured. With ξ turned on, this mean may
shift. A bias force is applied until the mean has been restored to the value found
with ξ turned off. The force to rebalance gives ξ . The final quantity required to
estimate g is the population ratio P1 /P2 . This is found by observing the system for
a period of time over which many escape events occur. The system can be sampled
with sufficiently long sampling period and the probabilities can be estimated by use
of equation (3.16). From these three measurements, the full counting statistics are
found.

3.6

Outlook

In this chapter we discussed a parametric sensing strategy based on noise-induced
escape. Parametric sensing is applicable to a wide range of problems. It only requires
that one can couple a quantity of interest to the parameters, or effective parameters,
75

of the device while preventing disturbances from doing the same. While this may
not be a trivial task, parametric sensors have found success in force and mass [121,
40, 42, 33, 43, 66] spectroscopy, and measurement of electron transport [6, 5], among
other applications. Likewise, the dynamic bridge may be employed for force or mass
spectroscopy, as it can be made extremely sensitive to changes in the natural frequency
of the resonator. However, this comes at the cost of long measurement times.
Alternatively, the bridge can be used for measurements of counting statistics, as
demonstrated in sections 3.3 and 3.5. The processes that may be counted with a
NEMS bridge depend on the mean rate and how it compares to the relaxation time
of the resonator. Pulse rates much faster than the relaxation time of the resonator will
appear approximately Gaussian, and therefore the higher moments will be difficult to
measure. Intriguing possible applications include measuring the transport of electrons
through a single electron transistor, or the transport of cells through a micro-channel,
both of which can be viewed as shot noise processes that can affect the dynamics of
a NEMS resonator.

76

Chapter 4
Phase Noise in Nonlinear
Oscillators
An oscillator is a device which produces a periodic signal. Such devices are central
to many areas of science and technology including timing, spectroscopy, communications, radar, GPS, etc. An oscillator is constructed from a frequency selective device,
which is inevitably lossy, and a feedback loop designed to compensate for this energy
loss in order to maintain a constant amplitude periodic signal. See, for example,
Figure 4.1(b). The frequency selective element may be, for example, a mechanical
pendulum as in a grandfather clock, a quartz crystal as in many electronic oscillator
circuits, or a microwave cavity as in a hydrogen MASER. In all cases there are random perturbations acting on the elements of the oscillator system stemming from the
environment. Each perturbation may knock the oscillator phase a little forward or
backward. Over time many small perturbations accumulate and the oscillator phase
diffuses. Figure 4.1(a) illustrates the different routes by which fluctuations affect the
oscillator phase. Fundamentally, this diffusion problem cannot be completely eliminated because an oscillator must be autonomous and thus have time-translational
symmetry. In other words, if we build an oscillator and set it up in a given initial
77

Environment
Resonator Noise Feedback Noise
Resonator
Action
Limiter

Amplifier

Phase Shifter

Angle
(a)

(b)

Figure 4.1: Noise roadmap (a) illustrating the different routes by which fluctuations
may affect the angle. Prototypical oscillator block diagram (b).

configuration, then we can observe its behavior. If we and set the oscillator in the
same initial configuration at a different time, we expect it to behave the same. This
property implies that the oscillator cannot know about a fixed reference point in time.
Likewise, it has no way to know if its perception of time has slipped a little bit from
“true” time. Thus, the only thing we can do to make a better oscillator is to reduce
the action of perturbations on the phase. There are two approaches one can take to
address this problem. The first approach, which is hardware specific, is to combat the
disturbances at their source. For example, we can put the oscillator in a temperature
controlled environment to diminish perturbations due to temperature fluctuations.
The second approach is to tune the oscillator parameters so that it operates optimally given the perturbations that are present. Returning to the noise roadmap in
figure 4.1(a), the first approach is analogous to reducing the quantity of noise descending along the four lines from the environment block. The second approach may
be used to, for example, cut the noise route connecting the action to the angle or reduce the noise strength relative to the oscillator’s inertia. In general both approaches
must be employed. In this chapter, however, we focus on some fundamental issues
related to the second approach.
78

Our goal is to analyze oscillators constructed with nano-eletro-mechanical systems (NEMS) as frequency selective elements. During the past decade NEMS have
demonstrated operation at high frequencies with high quality factors and low power
consumption [94]. These qualities put NEMS in a position to be a key component in
the development of compact, robust, low phase-noise frequency sources [45]. However,
due to their size, NEMS exhibit strong nonlinearity and relatively large noise levels,
which limit the linear dynamic range [90]. The realization of NEMS microwave sources
therefore requires engineering beyond the linear dynamic range of the frequency selective element. Thus, a fully nonlinear treatment of the problem is required. We
approach the problem by examining the slow dynamics in the action-angle coordinate
representation of an oscillator. For alternative analyses of phase noise, see [27, 22, 80].
Our treatment of this problem is arranged as follows. Our fundamental assumptions include neglecting the dynamics of the feedback loop and treating the resonator
as nearly conservative. We therefore model the oscillator with the resonator’s equations of motion. The most appropriate coordinate system in which to perform the
analysis is in terms of the action-angle coordinates of the resonator. Accordingly,
we review action-angle coordinates in section 4.1. We then discuss a microscopic
model for the resonator in section 4.2. This section provides insight into the noise
forces acting on an open-loop resonator. In section 4.3, we consider the closed loop
dynamics of a general nonlinear oscillator. Our primary result is an expression for
the spectrum of fractional frequency fluctuations in the oscillator. This expression is
useful for oscillator tuning and optimization. We illustrate this by example in section
4.4, in which we examine an oscillator employing a biased Duffing resonator. Using
noise and damping effects developed from the microscopic model presented in section
4.2, and a simple feedback model, we discuss oscillator tuning over a reduced twodimensional parameter space. The region of parameter space we consider is restricted
to reasonable values. In this restricted region, the minimum phase diffusion rate is
79

found when the resonator is maximally biased and operated at a particular amplitude where the frequency of the resonator is locally independent of the amplitude.
Following this example, we conclude the chapter with section 4.5 in which we discuss
possible directions for future research.

4.1

Action-Angle Coordinates

In the treatment of weakly nonlinear resonators we often employ the van der Pol transformation in order to write the equations as a set of coupled first order differential
equations with small perturbations (see section 2.1). The van der Pol transformation
represents the solution to the linear unperturbed problem, and the effects of the perturbations are incorporated by allowing the constants of motion (e.g., amplitude and
phase) to vary in time. These variations are slow when the perturbations are small.
To handle a strongly nonlinear resonator we can also perform a change of coordinates
based upon the solution of the unperturbed system, which is nonlinear in this case.
For perturbations of strongly nonlinear conservative systems the appropriate transformation is to action-angle coordinates, as described in this section. This review of
action-angle coordinates follows from [50, 69]. The coordinate transformation leads
to equations for the action and angle of the resonator, which are analogous to the
amplitude and phase coordinates employed in the weakly nonlinear case. For the
complete solution it is necessary to know the solution of the unperturbed nonlinear
problem. In this section, however, we perform the transformation using general functions which contain the solution, and appeal only to properties that these functions
can be shown to posses. Thus, while we demonstrate the method, its application is
case specific. We provide an example with explicit calculations in section 4.4.
We begin by considering a nearly conservative nonlinear resonator described by
80

the equations of motion
q = p,
˙
(4.1)
p = −U (q) + F (q, p, t),
˙
where U (q) is the potential describing the nonlinear stiffness of the resonator and F is
a perturbation. As discussed above, we wish make a change of coordinates reflecting
the solution of the conservative subsystem described by the Hamiltonian

H=

p2
+ U (q).
2

(4.2)

We assume that the potential, U , is confining so that the system exhibits periodic
behavior. The conservative subsystem can also appropriately be called an oscillator
since it exhibits periodic motion. For a given energy, E = H, the system will oscillate
between two turning points, qa and qb , which are solutions of the equation E =
U (qa ) = U (qb ). In phase space, the trajectory of the system belongs to a family
of closed orbits parameterized by E. To achieve action-angle coordinates we want
to transform the coordinate and momentum of the system such that this family of
orbits becomes a family of circles, each of which is traversed at a constant, but
energy dependent, angular rate, Ω(E). This transformation must preserve the area
of the orbits so that the parameterized family remains ordered by the energy; such
a transformation is a canonical transformation. By virtue of this transformation, an
orbit will be transformed into a circle with a radius given by its area normalized by
2π. This radius is the action, I, defined by

I(E) =

1
2π

p dq.,

(4.3)

which is a monotonic function of E. We can therefore use I in place of E to param-

81

eterize the family of periodic orbits. To show that I is monotonic, we consider the
derivative,
∂I
1 qb ∂p
1 qb 1
1 tb
T
=
dq =
dq =
dt =
≥ 0,
∂E
π qa ∂E
π qa p
π ta
2π

(4.4)

where T is the period of the orbit. To obtain this result we have used ∂p E = ∂p H = p
and dq = pdt. Thus, for any orbit of finite non-zero period, I is monotonic and we
find that the resonator frequency is given by

ω(I) =

∂I −1
.
∂E

(4.5)

In order to ensure that the phase space transformation preserves area, we used a generating function for the transformation [50]. For the action-angle coordinate transformation, the generating function is

S0 (q, I) =

p(q, I) dq.

(4.6)

The function p(q, I) is obtained from equations (4.2) and (4.3) with E(I) = H. The
canonical coordinate transformation defined by a generating function of this type is
given by [50]
∂S0 (q, I)
,
∂q
∂S0 (q, I)
Θ =
,
∂I
p =

(4.7)
(4.8)

where q and p are the old coordinate and conjugate momentum and Θ and I are the
new coordinate and conjugate momentum. In principle, to perform this coordinate
change, it is necessary to evaluate and invert equations (4.3) and (4.8). However,
here we are looking at the general case, and so we do not need to do this, only note
82

that it is possible. In order to apply this transformation to equations (4.1), we make
use of the following identities,
−1
∂ 2 S0 (q, I)
∂p ∂E
∂q
∂Θ
∂p
ω
=
=
=
=
= ,
∂Θ I
∂q I
∂I∂q
∂I q ∂E q ∂I
p
∂p
1
∂q
1
= − U (q)
= − U (q),
∂Θ I
p
∂Θ I
ω
∂p
ω 1
∂q
=
− U (q)
.
∂I Θ
p p
∂I Θ

(4.9)
(4.10)
(4.11)

Now, we can rewrite equation (4.1) in terms of the new coordinates, Θ and I, by
using the functions p(Θ, I) and q(Θ, I), as follows
∂q ˙
Θ+
∂Θ
∂p ˙
p =
˙
Θ+
∂Θ
q =
˙

∂q ˙
p ˙
∂q ˙
I = Θ+
I = p,
∂I
ω
∂I
∂p ˙
1
ω U (q) ∂q
˙
I = − U (q)Θ +
−
∂I
ω
p
p ∂I

(4.12)
˙
I = −U (q) + F. (4.13)

˙
˙
Solving for I and Θ gives the resonator equation in action-angle coordinates
˙
I = (∂Θ q)F,

(4.14)

˙
Θ = ω(I) − (∂I q)F.
Note that when F = 0 this provides the desired form in which the phase space consists
of circles on which I remains constant and which are traversed at constant angular
rate ω(I). When F is small, we use this convenient solution as a basis for perturbation
calculations.

4.2

Microscopic Resonator Model

Before proceeding to treat equations (4.14), we pause to discuss what kind of terms
we might expect to find in the perturbation, F , for an oscillator. Generally, we
83

can split F up into two parts, one arising from the resonator’s interaction with the
environment and another arising from the feedback loop. In this section we consider
the first part. In section 4.3 we consider the complete closed loop dynamics.
In order to examine components of F that arise from the resonator’s interaction
with the environment, we must model the resonator and environment as a complete
system. We therefore review the theory of an oscillator coupled to a medium, following
Dykman and Krivoglaz [36]. The interested reader is directed to this reference and
its cited works for a more complete discussion.
Our microscopic resonator model is a nonlinear conservative oscillator coupled to
a bath modeled by a large number of harmonic oscillators. This model is appropriate
for a NEMS resonator where the nonlinear resonator of interest is a particular mode of
the complete device and substrate system. The thermal bath consists of the myriad of
other vibrational modes of this system. We assume the nonlinear oscillator is driven
by an external force and therefore it’s coordinate and momenta are not necessarily
small. On the other hand, the bath oscillators are assumed to have small amplitude
so that the we can drop cubic and higher order terms in the Hamiltonian. A change of
coordinates among the bath variables is always possible that will decouple the linear
bath subsystem, and so we also drop coupling between modes of the bath. Finally,
we assume that the coupling between the oscillator of interest and the bath is weak,
and so we keep only the linear terms in the bath coordinates. The Hamiltonian for
this combined system is given by

H=

p2
0
2

p2
k

+ U (q0 ) +
k

1 2 2
+ ωk qk + k qk HI (p0 , q0 )
2
2

− q0 f0 cos ωt,

(4.15)

where the singled-out nonlinear oscillator is described by coordinate q0 , momentum
p0 , and potential U (q0 ), and the bath oscillators are described by coordinates qk and
momenta pk . The function HI captures the dependence of the interaction energy on
84

the nonlinear oscillator’s coordinate and momentum. For now, we leave this function
unspecified. The equations of motion for this system can be written as
∂HI
q0 = p 0 +
˙
∂p0
k

∂HI
p0 = −U (q0 ) −
˙
∂q0
k
qk + ωk qk = − k HI (p0 , q0 ),
¨

k qk ,

k qk + f0 cos ωt,

(4.16)

(4.17)
(4.18)

where we have removed the bath momenta variables by the substitution pk = qk . We
˙
can remove the bath variables from equations (4.16) and (4.17) by formally solving
equation (4.18). The solution is
t

qk (t) = Ak cos(ωk t + φk ) −

k sin ω (t − τ ) H p (τ ), q (τ ) dτ,
0
I 0
k
ωk
0

(4.19)

where Ak and φk capture the initial conditions of the k th environmental oscillator.
These are unknown and so we model them as random variables distributed according
to the laws of equilibrium thermodynamics. Since the bath oscillators are uncoupled,
their initial conditions are independent random variables. The probability for the k th
bath oscillator to have amplitude Ak ∈ [A, A + dA] and phase φk ∈ [φ, φ + dφ] is
given by the Boltzmann distribution
2
2
ωk
ωk A2
P {Ak ∈ [A, A + dA], φk ∈ [φ, φ + dφ]} =
exp −
A dA dφ, (4.20)
2πkb T
2kb T
where Ak ≥ 0, φk ∈ [0, 2π], kb is Boltzmann’s constant, and T is the bath temperature. The phase, φk , is uniformly distributed, so equation (4.20) is independent of
φ. To apply equation (4.19) to equations (4.16) and (4.17) we must consider the sum

85

over k on k qk . We partition this term into two pieces. The first piece is the sum
over the homogeneous solutions, which gives a random process, ξ(t), parameterized
by the random variables Ak and φk . The second piece is the sum over the particular
solutions which gives a functional, L[p0 , q0 ], describing the retarded reaction of the
environment on the nonlinear oscillator. These effects are expressed as

ξ(t) ≡

k Ak cos(ωk t + φk ),

(4.21)

t 2
k sin ω (τ − t) H p (τ ), q (τ ) dτ.
0
I 0
k
ωk
0

(4.22)

k
L[p0 , q0 ] ≡
k

The environment is assumed to contain a large number of harmonic oscillators. Accordingly, we approximate the sum over k as an integral. Toward this end, we define
the spectral density of the bath, g(ω), as
2
k = ωg(ω)dω.
ωk
ω<ωk <ω+dω

(4.23)

As mentioned above, the sum over the homogeneous solutions of the bath oscillators,
ξ(t), represents a stochastic process. This process has zero mean, since the phases are
uniformly distributed over [0, 2π], and the autocorrelation function for this process is
given by
∞
ξ(t)ξ(t )

= kb T

g(ω) cos ω(t − t ) dω.

(4.24)

0

In the simple case where g is a constant, we have ξ(t)ξ(t ) = πgkb T δ(t − t ). The
power spectral density for ξ is then given by
∞
Sξ (ω) =

−∞

e−iωt ξ(t)ξ(0) dt = πkb T g(ω).

86

(4.25)

The sum over the particular solutions of the bath oscillators, equation (4.22), requires
more detailed treatment. First, we write the interaction term, HI , as
HI p0 (τ ), q0 (τ ) =

An (τ )eiΩn τ ,

(4.26)

n

where Ωn = nΩ are the harmonics of the expansion frequency Next we note that
when t − τ becomes sufficiently large the summation of sinusoids in equation (4.22)
becomes approximately zero. The time beyond which this occurs is the correlation
−1
time of the bath, tc . This correlation time is of the order ωc where ωc is the cutoff
frequency, i.e., the largest frequency of the bath oscillators. We now assume that the
coefficients An change very little over the bath correlation time, tc . Accordingly, we
expand An in a Taylor series about τ = t and keep only the first term,

HI p0 (τ ), q0 (τ )

˙
eiΩn τ (An (t) + An (t)(τ − t) + · · · ),

(4.27)

An (t)eiΩn τ .

=

(4.28)

n
≈
n

This assumption requires that the resonator relaxes slowly with respect to the bath
correlation time. After making this approximation in equation (4.22) we can evaluate
the integral over the time variable, τ . This gives
2
0
k A (t)eiΩn t
dt1 eiΩn t1 sin(ωk t1 )
n
ωk
−t
n,k
2
−i(Ωn −ωk )t 1 − e−i(Ωn +ωk )t
k A (t)eiΩn t 1 − e
=
−
.
n
2ωk
Ωn − ωk
Ωn + ωk
n,k

L≈

(4.29)

In order to further simply the expression for the functional L, we consider only times
much large than the correlation time of the bath, t

tc . This allows us to ignore

the initial transient of the bath and arrive at a “steady-state” resonator model. For

87

large t, we make the approximation (for t > 0)
1 − e−i(Ωn ±ω)t
1
≈ iπδ(Ωn ± ω) + v. p.
,
Ωn ± ω
Ωn ± ω

(4.30)

where v. p. stands for the principal value. Approximating the sum as an integral in
equation (4.29) and applying equation (4.30) gives
An (t)eiΩn t [2iΓn Ωn − Pn ] ,

L≈

(4.31)

n
where Γn and Pn are
1
πg(Ωn ),
4
∞ ω 2 g(ω)
Pn ≡ v. p.
dω.
0 Ω2 − ω 2
n
Γn ≡

(4.32)
(4.33)

This completes the development of the microscopic resonator model. We now rewrite
the resonator equations so that they can be viewed collectively. We also drop the
subscript ’0’ on the coordinate and momentum, since all other similar coordinates
have been removed. The resulting resonator equations are
∂HI
ξ(t) + L[p, q] ,
∂p
∂HI
p = −U (q) −
˙
ξ(t) + L[p, q] + f0 cos ωt,
∂q
Ω π/Ω
An =
H e−inΩt dt,
2π −π/Ω I

(4.35)

einΩt An 2inΩΓn − Pn .

(4.37)

q = p+
˙

L[p, q] ≈

(4.34)

(4.36)

n

1 2
As a simple example, suppose that U = 2 ω0 q 2 and HI = q, which is the simplest example that yields a nontrivial result. We suppose that q = ueiωt + c.c. and
88

p = iωueiωt + c.c. where c.c. stands for complex conjugate. For this example, the
expansion frequency for HI is Ω = ω and we have A1 = u, A−1 = u∗ and An = 0
for all other n. Applying this example to equations (4.34)-(4.37) yields a differential
equation u(t). Assuming that u changes slowly in time, we apply the method of
averaging. The resulting slow flow equation for u is

u = − Γ − i(ω − ω0 ) u −
˙
ˆ

if0 ie−iωt
+
ξ(t),
4ω
2ω

(4.38)

2
where ω0 is the renormalized resonator frequency ω0 = ω0 − P and we have dropped
ˆ
ˆ2
the subscripts on Γ and P since there is only one value for each, namely Γ1 and
P1 . This result is self consistent with our assumptions so long as the damping rate,
Γ, the detuning, ω − ω0 , the external excitation amplitude, f0 , and the noise ξ
ˆ
are sufficiently small. This is consistent with a high quality resonator forced near
its resonant frequency. The noise force ie−iωt ξ(t)/2ω can be approximated by a
complex white Gaussian noise with uncorrelated real and imaginary parts, each with
strength Γkb T /ω 2 , so smallness of Γ implies smallness of ξ. This is the same slow
flow equation one would obtain from the phenomenological resonator model
q + 2Γq + ω0 q = f0 cos ωt − ξ(t).
¨
˙ ˆ2

(4.39)

Thus we see that the coupling to the bath has the net effect of producing a frequency
shift, a damping force, and a noise force for the resonator.

1
1 2
As another example, suppose that U = 2 ω0 q 2 + 4 γq 4 and HI = q + αq 2 . Then
we employ the conventional van der pol transformation q = x cos ωt + y sin ωt, p =
−ωx sin ωt + ωy cos ωt. Applying this transformation to equations (4.34)-(4.37) gives
the equations for the quadrature coordinates x and y. These coordinates vary slowly in
2
time if the detuning ω− ω0 , the renormalized nonlinear coefficient, γ = γ − 3 α2 (2P0 +
ˆ
ˆ
89

P2 ), the damping coefficients, Γn , and the noise, ξ(t), are small. As in the previous
2
example, ω0 is the renormalized linear oscillator frequency, ω0 = ω0 − P1 . Assuming
ˆ
ˆ2
that these terms are small, we apply the method of averaging and obtain
3ˆ 2
γ
x = −Γ1 x − α2 Γ2 (x2 + y 2 )x − (ω − ω0 )y −
˙
ˆ
(x + y 2 )y
8ω
ξ(t)
+
sin ωt 1 + 2αx cos ωt + 2αy sin ωt ,
ω
f
3ˆ 2
γ
(x + y 2 )x + 0
y = −Γ1 y − α2 Γ2 (x2 + y 2 )y + (ω − ω0 )x +
˙
ˆ
8ω
4ω
ξ(t)
−
cos ωt 1 + 2αx cos ωt + 2αy sin ωt .
ω

(4.40)

(4.41)

These slow flow equations are consistent with those resulting from the phenomenological model
q + 2Γ1 q + 8α2 Γ2 q 2 q + ω0 q + γ q 3 + (1 + 2αq)ξ(t) = f0 cos ωt.
¨
˙
˙ ˆ2
ˆ

(4.42)

This model exhibits additive and multiplicative noise. It is known that the diffusion
rate of the phase of a self excited oscillator resulting from additive noise decreases
with increasing oscillator amplitude [72]. This motivates strong driving of the oscillator’s frequency selective element, eliciting nonlinear behavior. The above example
illustrates that how larger resonator amplitudes can bring multiplicative noise from
the coupling to the environment. The multiplicative noise can be thought of as a fluctuation in the frequency of the resonator. We show below, equation (4.84), that the
multiplicative noise produces phase diffusion in an oscillator that does not decrease
with oscillator amplitude while the oscillator is operated in the linear range.
In the analysis above we have neglected quadratic terms in the bath coordinates,
qk , in the interaction energy which couples the oscillator to the environment, see
equation (4.15). This was done to simplify the discussion in this dissertation. These
terms, however, can contribute to phase noise in self excited oscillators. Analysis
90

of a resonator model including these terms can be found in references [34, 36]. Our
analysis in the following section begins with a general phenomenological model, so the
simplified microscopic model employed in this section does not diminish the generality
of the subsequent results.

4.3

Oscillator Frequency Fluctuations

In this section, we consider general closed loop equations for a self-sustained nonlinear
oscillator. We neglect the dynamics of the feedback loop, modeling the feedback forces
on the resonator as instantaneous functions of the resonator state. As in the previous
section, the resonator is assumed to be lightly damped as a result of weak coupling
to a thermal bath. In the previous section we demonstrated that coupling to the
bath produces three types of forces acting on the resonator: conservative forces that
renormalize the resonator potential, damping forces which account for the loss of
energy into the bath, and stochastic forces which result from the summed motion of
the many degrees of freedom present in the bath. From this insight, we model our
self-excited oscillator phenomenologically as follows. The resonator is assumed to be
nearly conservative with a potential U (q). This contains the conservative forces due
to coupling to the bath. In addition, we add the friction forces L(p, q) which capture
the energy loss into the bath, the driving, or gain, forces G(p, q) which arise from
the feedback loop, and the fluctuating forces f (p, q, t) which arise from the thermal
motion of the bath and the noise processes produced by the feedback electronics.
Thus, the resonator equation of motion is

q = p,
˙

(4.43)

p = −∂q U + G(p, q) − L(p, q) + f (p, q, t).
˙

(4.44)

91

Note that we found that coupling to the bath can produce additional terms in the q
˙
equation when the coupling is a function of p. For this model, we assume that should
any such coupling be present, we can redefine p to be identically q and lump the
˙
resulting perturbations into f . We assume that L and G are small compared to the
characteristic value of the conservative force, |∂q U |, which implies small dissipation
and small gains required to maintain oscillation. The various effects lumped into f
are also assumed to be similarly small, and the stochastic processes found in f (p, q, t)
are taken to be stationary with zero mean. Stationarity is required for the oscillator
to have time-translational symmetry, and any nonzero mean may be moved from f
to the other terms in the equation.
Important for our analysis is the separation of time scales into the “fast” time
on the order of the vibration period, and slow times on the order of the vibration
decay rate. It is this slow scale on which small random perturbations accumulate.
As discussed in section 4.1, action-angle coordinates are the natural variables for the
treatment of this nearly conservative system, for which the action and angle will vary
on the slow time scale. Thus, we use a canonical coordinate transformation from (p, q)
to action-angle coordinates, (I, φ), expressed by functions p(φ, I), q(φ, I) that are 2π
periodic in φ. Without loss of generality, we set φ = 0, π at the turning points at
which p = 0. Accordingly, q is an even function in φ and p is odd in φ. To describe the
oscillator in action-angle coordinates, we use equation (4.14) with the perturbation,
F , from equations (4.43) and (4.44). This gives

˙
I = (∂φ q)(G − L + f ),

(4.45)

˙
φ = ω(I) − (∂I q)(G − L + f ).

(4.46)

The assertion that G, L, and f are weak perturbations on the resonator imply that all
terms in equations (4.45) and (4.46) are small compared to ω(I) during steady-state
92

operation. The oscillator action and frequency fluctuate about their mean values,
˙
I0 = I and Ω = φ respectively. It is assumed that these fluctuations are small
and that ω(I) remains close to Ω, i.e. any shift in the resonator frequency due to
the perturbative terms is also small. Equations (4.45) and (4.46) exhibit a separation
of time scales, since the angle, φ, varies rapidly while the action, I, and the time
deviation, x(t), defined by [3],

x(t) =

φ − Ωt
,
Ω

(4.47)

vary slowly. The time deviation x(t) is a quantity of primary interest for the applications of a resonator as a time standard, as it describes the deviation of time as
measured by the oscillator from the desired (true) time. A related quantity is the
time derivative of x(t), known as the fractional frequency fluctuation [3],

y(t) = x(t) =
˙

˙
φ−Ω
,
Ω

(4.48)

which is an instantaneous measure of the oscillator frequency deviation. Among the
common measures of oscillator stability [3] is the spectrum of fractional frequency
fluctuations, Sy (ν), which we presently use to quantify the frequency noise in our
oscillator model.
We begin by defining the phase deviation ϕ = φ − Ωt. Changing phase variables
from φ to ϕ in equations (4.45) and (4.46) gives

˙
I = (∂ϕ q)(G − L + f ),

(4.49)

ϕ = σ(I) − (∂I q)(G − L + f ),
˙

(4.50)

where σ(I) = ω(I) − Ω is the oscillator detuning, which is assumed to be small. With
the change to the phase deviation, the right hand sides of equations (4.49) and (4.50)
93

have become small, and are also periodic in time. Thus, the oscillator equations are
now in a form amenable to the method of averaging.
In appendix B we discuss how a near-identity coordinate transformation can be
used to push the time-periodic aspects of the system out to higher order. Subsequent
truncation of the series yields a simplified model that asymptotically describes the
system behavior. This method can be used to treat non-white noises, and therefore
can handle 1/f type noise (here f is frequency) so long as the perturbations remain
small. This requires that either a finite time interval be considered, or a low frequency
cut-off be imposed on 1/f type noises, which amounts to the same thing. We now
treat equations (4.49) and (4.50) with this method. First, we expand the stochastic
terms in a Fourier series in φ. We assume that the Fourier coefficients can be separated
into functions of action, I, and time, t,

(∂µ q)f =

A

(µ)
(I)eik(ϕ+Ωt) ζn (t),
n,k

µ = I, ϕ.

(4.51)

n,k
The near identity coordinate transformation operates on the use of a fictitious narrow band filter with impulse response h(t) and frequency response H(ω) centered
at ω = 0 with H(0) = 1. Moreover, we assume that the cutoff frequency for this
filtering function is less than Ω/2. By use of this filtering function, we construct the
transformed stochastic processes, ξϕ and ξI , by the integral
ξµ ≡

∞
(µ)
A eikϕ
h(t − t )ζn (t )eikΩt dt .
n,k
−∞
n,k

(4.52)

The processes ξµ have zero mean by our original assumption on f . The spectra of ξµ
can be calculated from the spectra of ζn (t). Let

ζn (t)ζm (t ) =

∞
1
Rn,m (ω)eiω(t−t ) dω.
2π −∞
94

(4.53)

Since we assumed the filtering function has a cutoff frequency less than Ω/2, the
transformed noise processes are wide-sense stationary. Letting Sµ,ν be the crossspectrum of ξµ and ξν , we find
∞
Sµ,ν (I, ω) =

ξµ (I, t)ξν (I, 0) e−iωt dt

−∞

(µ)
(ν)
A (I)A
(I)Rn,m (ω − kΩ)|H(ω)|2 .
n,k
m,−k

=

(4.54)

n,m,k
From equation (4.54) we can see how the nearly periodic motion of the oscillator
selects the noise spectra near the harmonics of its fundamental frequency. The noise
in these narrow bands is modulated down to baseband by the oscillator motion and
accumulates over long times, resulting in diffusion of the phase. A diagram illustrating
this phenomenon is shown in figure 4.2. In this figure a hypothetical noise spectrum
is plotted along with several narrow-band filtering functions. These functions are
given different heights illustrating that the magnitudes of the harmonics modulating
the noise typically decrease as the order of the harmonic increases. This effect is
(µ) (ν)
described in the functions A A
, which we expect to show a general trend of
n,k m,−k
decreasing with increasing k.
The formulation in appendix B illustrates how the noise terms (∂µ q)f also give
rise to a deterministic drift term. The noise-induced drift can likewise be written in
terms of the Fourier expansion coefficients, as follows

(µ)
∂A
(I) (µ) 
 (ϕ) m,−k
+ ikA A
An,k
n,k m,−k 
∂I

dµ =

n,m,k
∞
1 − |H(ω)|2
×
dω.
Rn,m (ω − kΩ)
2πiω
−∞

(4.55)

From equation (4.55) we can see that the noise-induced drift is a mean effect arising
from the noise at frequencies away from the oscillator harmonics. If the noise is taken
95

Schematic representation of the
oscillator transfer function

Hypothetical input
noise spectrum

0

1

2

3

4

5

6

7

ν/Ω
Figure 4.2: Illustration depicting how noise is filtered and down-converted. The noise
is selected near the oscillator harmonics, that is, in the shaded regions, weighted by
the Fourier coefficients in equation (4.51), and modulated down to baseband by the
oscillator motion.

to be white, then Rn,m is a constant and the integrand is odd. Thus, the integration
gives zero and this drift is zero.
The final term in the transformed equations comes from the deterministic part of
the oscillator equation. This is simply the average of the deterministic terms over one
period. The final form of equations (4.49) and (4.50) after the transformation is

˙
I ≈ (∂ϕ q)(G − L) + dϕ (I) + ξϕ (I, t),

(4.56)

ϕ ≈ σ(I) − (∂I q)(G − L) − dI (I) − ξI (I, t),
˙

(4.57)

where the overline means averaging over ϕ. We note that the right hand sides of
equations (4.56) and (4.57) are independent of ϕ. This is a consequence of the timetranslational symmetry of the oscillator. Without knowledge of a reference time, the
phase must have a zero eigenvalue and so it does not appear in equations (4.56) and
(4.57). We expect the oscillator to reach a steady state where the action fluctuates
about a mean value, denoted I0 . Assuming these fluctuations are small, we expand
96

equations (4.56) and (4.57) in I about I0 . Let
0 =

(∂ϕ q)(G − L) + dϕ (I) I=I ,
0

ω(I) − (∂I q)(G − L) − dI (I) I=I ,
0
∂
α =
(∂ϕ q)(G − L) + dϕ (I) I=I ,
∂I
0
∂
σ(I) − (∂I q)(G − L) − dI (I) I=I .
β =
∂I
0

Ω =

(4.58)
(4.59)
(4.60)
(4.61)

Keeping only leading order terms in the expansion, we approximate equations (4.56)
and (4.57) as

˙
I = −α(I − I0 ) + ξϕ (I0 , t),

(4.62)

ϕ = β(I − I0 ) − ξI (I0 , t).
˙

(4.63)

The fractional frequency fluctuation, defined in equation (4.48), can be expressed as
y = ϕ/Ω. To find its first two moments, we solve equation (4.62) and substitute the
˙
solution into equation (4.63). Naturally y has zero mean since we defined Ω to be
the mean frequency. We characterize the second moment of y by its power spectral
density. The spectrum of frequency fluctuations can be written as a function of the
spectra of the noises acting on the action and angle coordiantes, ξϕ and ξI , and is
given by

Sy (ω) =

β 2 Sϕ,ϕ (ω) 2αβ Sϕ,I (ω)
1
−
+
S (ω),
2 α2 + ω 2
2 α2 + ω 2 Ω2 I,I
Ω
Ω

(4.64)

where the spectra Sµ,ν are given in equation (4.54). From equation (4.64) we see
that the fluctuations that cause phase diffusion can be sorted into two groups. The
first group are noise sources that cause fluctuations in the resonator action, which
then couples to the phase, resulting in phase diffusion. The second group are noise

97

sources that directly affect the phase. The first group can be removed by setting
β = 0 such that the oscillator action and phase coordinates are decoupled. For an
oscillator employing a linear resonator as the frequency selective element, we have
σ (I) = 0. In addition, we note that ∂I q is an even function of φ. Therefore, the
action and phase coordinates will be uncoupled if G − L is an odd function of φ and
dI = 0. This is often the case, and even if it is not, the coupling will generally be
small since σ (I) = ω (I) is the dominant term in (4.61). This is one reason why it
is often favorable to use linear resonators in oscillator designs. However, nonlinear
resonators may also exhibit decoupling of the action and phase coordinates [38]. In
the next section we discuss an example and examine the phase noise characteristics
of a hypothetical oscillator constructed with such a resonator.
Before discussing the example, however, we will provide a summary of the necessary steps required to calculate the spectrum of frequency fluctuations. The necessary
tasks are as follows.
1. Obtain an ODE model for the oscillator.
2. Separate the model into i) a conservative subsystem with potential U , ii) a deterministic perturbation, G − L, and iii) a zero mean stochastic perturbation, f .
3. Calculate the action-angle transformation for the conservative subsystem to obtain
the orbit function q(Θ, I) and the frequency dispersion ω(I).
4. Calculate the Fourier coefficients of the noise terms (∂µ q)f for µ = Θ, I, implicitly
defined in equation (4.51).
5. Choose a filtering function, H(ω), e.g. a rectangular window of width min[ω(I)/2].
6. Calculate the noise induced drift using equation (4.55).
7. Average the deterministic perturbation over one cycle of the angle.
8. Calculate Ω, α, and β using equations (4.59-4.61).
9. Calculate the slow noise spectra using equation (4.54).
10. Calculate the spectrum of frequency fluctuations using equation (4.64).
98

4.4

The Biased Duffing Oscillator

In this section, we examine the phase noise in an oscillator employing a prototypical
nonlinear resonator modeled by Duffing’s equation with constant bias force. The exercise is an illustration of the results obtained in the previous section and demonstrates
how phase noise can be reduced by tuning the oscillator’s operating point. We begin
by considering the conservative subsystem that generates the periodic orbits along
which the perturbations are calculated. The potential of the conservative subsystem
is given by

U (q) =

1 4 1 2
q + q + λq,
4
2

(4.65)

where λ is the bias force. Equation (4.65) is a one-parameter family of resonator
potentials which brings another parameter to our optimization problem in addition
to those describing the feedback force. In order to do this, we must first calculate the
action-angle coordinate transformation. This requires the four roots of the equation
U (q) − E = 0, which are difficult to find in closed form for arbitrary λ. Thus, we let
a, b, c, and c∗ represent the four roots, to be found numerically, and proceed with
the analysis. Note that these roots are functions of the bias parameter, λ, and the
operating action, I. We let a and b be the real roots with b > a, q(0, I) = b, and
q(π, I) = a. The complex root are not of interest for the present analysis.
The conservative subsystem of our oscillator model is described by the differential
equation q + U (q) = 0. The solution to this equation is
¨

q(t) =
u =

(aB + bA) + (bA − aB) cn(u|m)
,
(A + B) + (A − B) cn(u|m)
AB
t,
2

(4.66)
(4.67)

where cn(u|m) is a Jacobi elliptic function with argument u and parameter m, and
99

2.0
λ=2
Frequency, ω

1.8

λ = 1.5

1.6

λ=1

1.4

λ = 0.5

1.2
1.0
0.0

λ=0
0.5

1.0

1.5

2.0

2.5

3.0

Action, I
Figure 4.3: Frequency dispersion for the biased Duffing resonator for several values
of the bias. Red dashed line is the locus of critical (zero dispersion) points.
A, B, and m are

A = |a − c|,

(4.68)

B = |b − c|,

(4.69)

m =

(a − b)2 − (A − B)2
.
4AB

(4.70)

The fundamental frequency of q(t) given in equation (4.66) is
√
π 2AB
ω(I) =
,
4K(m)

(4.71)

where K(m) is the complete elliptic integral of the first kind with parameter m [1, 18].
The frequency as a function of action is shown in figure 4.3 for various values of λ. The
key feature of the frequency dispersion curves is the extrema, called zero dispersion
points, which exist for sufficiently large λ. The red dashed line in figure 4.3 traces
the locus of these extrema. While this resonator exhibits points of zero dispersion,
it is just one example of a broad class of nonlinear resonators with this feature; for
other examples see [101]. The complete action-angle coordinate transformation can
be constructed from the solution (4.66) by assigning each point on the orbit an angle
100

momentum, p

momentum, p

8
6
4
2
0
-2
-4
-6
-8
-4 -3 -2 -1 0 1 2
coordinate, q
(a)

3

4

8
6
4
2
0
-2
-4
-6
-8
-4

-3

-2

-1
0
1
coordinate, q
(b)

2

3

Figure 4.4: Several closed orbits belonging to the potential (4.65) with (a) λ = 0 and
(b) λ = 4.

from 0 to 2π. We do this by relating the angle to time,

t=

2K(m)
π

2
φ.
AB

(4.72)

This gives the coordinate transformation
(aB + bA) + (bA − aB) cn(u|m)
,
(A + B) + (A − B) cn(u|m)
∂q(φ, I)
p(φ, I) = ω(I)
,
∂φ
2K(m)
u =
φ.
π
q(φ, I) =

(4.73)
(4.74)
(4.75)

Several example orbits are shown in figure 4.4, for λ = 0 and λ = 4. These figures
show that the one parameter family of potentials considered in this problem includes
orbits with a variety of shapes.

Before proceeding to discuss the other aspects of our oscillator model, we will
examine the function q(φ, I) in the limits of large and small action. For small action
101

U (q) can be approximated as a parabola. Thus

q(φ, I) ≈ q0 (λ) +

2I
cos φ,
ω0 (λ)

(4.76)

where ω0 (λ) is the frequency of the resonator in the limit I → 0. For small action
the frequency is independent of the action and so we have q ∝ I 1/2 and p ∝ I 1/2 .
For large action, the potential can be approximated by U ≈ q 4 /4, giving

q(φ, I) ≈

1/3
3π
I
cn
2K(0.5)

2K(0.5)
φ .
π

(4.77)

In this approximation we have q ∝ I 1/3 and ω ∝ I 1/3 , so p ∝ I 2/3 . From these
limiting cases for q(φ, I), we can develop scaling laws for noise sources that may
provide insight useful for locating an optimal operating point. For example, suppose
that f includes the term (q − q0 )a pb ζi (t). Such noise terms may arise from nonlinear
coupling, through the coordinate, to the environment or noise in the feedback loop.
By considering time symmetries of the response and the small and large action limits
of p and q, we can obtain some useful information about the noise processes ξϕ and
ξI that are ultimately responsible for diffusion of the oscillator phase. First, we recall
that q is evenin φ and p is odd in φ. It follows that
(I)
(I)
A
= (−1)b A ,
i,−k
i,k
(ϕ)
(ϕ)
A
= (−1)b+1 A .
i,−k
i,k

(4.78)
(4.79)

Next, by applying the scaling limits for p and q noted above we determine the following

102

scaling of the harmonic coefficients,

(I)
A
∝
i,k
(ϕ)
A
∝
i,k

I (a+b−1)/2

for small I

I (a+2b−2)/3 for large I
I (a+b+1)/2

for small I

I (a+2b+1)/3 for large I

,

(4.80)

.

(4.81)

Using the symmetry arguments, equations (4.78) and (4.79), the sum over k with
n = m in equation (4.55) gives zero. Therefore, a noise such as the one considered
here cannot directly induce a frequency shift, though such a shift may be possible
through correlations with other noises. Next we examine the scaling of the spectral
components, Sµ,ν , arising from the noise term considered. Applying the scaling laws
in equations (4.80) and (4.81) to equation (4.54) we find

Sϕ,ϕ ∝

Sϕ,I ∝

SI,I ∝

I a+b+1

for small I

I (2a+4b+2)/3 for large I
I a+b

for small I

I (2a+4b−1)/3 for large I
I a+b−1

for small I

I (2a+4b−4)/3 for large I

,

(4.82)

,

(4.83)

.

(4.84)

It is understood that additive noise in an oscillator employing a linear resonator produces phase noise that decreases as the inverse of the oscillator amplitude squared
[72]. This result is captured here as well. For small action the resonator is approximately linear and the action-angle coupling coefficient is approximately zero. The
spectrum SI,I dominates the spectrum of fractional frequency fluctuations. For an
additive noise a = b = 0 and we have SI,I ∝ I −1 . For large action multiplicative
noise will dominate additive noises. Moreover, Sϕ,ϕ will dominate SI,I unless the
103

action linear coefficient α is very large and/or the action-angle coupling coefficient β
is very small. Thus we again see, as we saw in equation (4.64), that is favorable to
operate the oscillator where the action and angle become decoupled, that is, where
β = 0. Operating in this parameter subspace, the remaining parameters should be
chosen to minimize SI,I . Note that a minimum w.r.t. action should exist if multiplicative noises of sufficiently high order, 2a + 4b − 4 > 0, exist. Otherwise, we should
raise the action as high as possible to reduce the phase noise.
To continue with our example, we now construct the remaining terms in our
oscillator model. We choose our loss and noise terms by appealing to the microscopic
resonator model discussed in section 4.2. Moreover, we will assume the “Ohmic
dissipation model” by taking g(ω) to be a constant value g for ω up to a cutoff
frequency ωc . This model is valid when the bath can respond quickly compare to the
motion of the resonator. The consequence of this assumption is that we take Pn = P
and Γn = Γ for all n. Next, we take
1
HI = (q − q0 ) + (q − q0 )2 ,
2

(4.85)

and assume those terms proportional to P have already been incorporated into the
potential given in equation (4.65). With this form of the interaction function, HI , the
loss of energy to the bath is given by the microscopic resonator formulation, equations
(4.16) and (4.17),
L = 2Γp 1 + 2(q − q0 ) + (q − q0 )2 .

(4.86)

The first term is the usual viscous friction, the averaged effect of which is linear in
the action, I,
Γ π ∂q
Γ
(∂φ q)L1 =
p dφ =
π −π ∂φ
π
104

p dq = 2ΓI.

(4.87)

The higher order terms, which cannot be averaged in closed form, give nonlinear
friction. This loss model does not shift the frequency since q − q0 and ∂I q are even in
φ and p is odd, rendering (∂I q)L = 0 by symmetry. The microscopic resonator model
also gives the multiplicative noise acting on the resonator due to the environment

fenv = ξenv (t) 1 + q − q0 .

(4.88)

For simplicity we assume that ξenv (t) is a white noise with strength Denv . The white
noise assumption is the simplest since the noise-induced drift terms dϕ and dI in the
averaged oscillator model, equations (4.56) and (4.57), are zero.

To compensate for the energy loss, the feedback loop must replenish the resonator’s
energy. This can be done in a variety of ways from simple positive feedback to
parametric forcing. The various methods may differ in the details, but they result
in a qualitatively similar average effect. We therefore choose to consider a saturated
feedback of the resonator velocity, with specific gain model of

G = g1 tanh g2 q(φ + θ, I) − q0 ,

(4.89)

where g1 is the feedback saturation level, g2 is a parameter used to capture the
rapid, but finite, rate of change of the saturation device in the feedback circuit, and
θ is a phase shift. This phase shift is necessary for the force G to do work on the
resonator. The maximum work is done when θ = ±π/2. This type of feedback could
be constructed by measuring the resonator motion, amplifying it, and then passing it
through a hard limiter. The physical implementation of the feedback loop must also
bring with it some noise. We therefore add to our noise model a term arising from
the feedback loop. A simple model is to assume that the feedback coefficients g1 and
105

θ fluctuate about their nominal values. This leads to the feedback noise model

ffb = ξfb,1 (t) tanh g2 q(φ + θ, I) − q0
g g
+ξfb,2 (t) 1 2 sech2 g2 q(φ + θ, I) − q0 p(φ + θ, I),
ω

(4.90)

where ξfb (t) represents the fluctuations of g1 and ξfb (t) represents the fluctuations
1
2
of θ. Both noises are assumed to be white and uncorrelated with each other and ξenv ,
and the noise strengths are taken to be Dfb,i , i = 1, 2.
This oscillator model possesses three tunable parameters: the feedback parameters, g1 , and θ, and the resonator (bias) parameter λ. All other parameters can be
thought of as given once a resonator element is selected. The feedback parameter g1
˙
is related to the action at the oscillator operating point by requiring I = 0. Thus,
we can replace this parameter with the operating action I0 . For simplicity, we set
θ = π/2 so that the feedback does the maximum work on the resonator. We are
thus left with two tunable parameters which we can vary over a reasonable range.
The action will be confined to values less than 10, where the nonlinearity becomes
extremely strong. The bias will be confined to values less than 2.5 which corresponds
with shifting the equilibrium position to q0 = −1.11, pushing the system again into
the strongly nonlinear region. For our numeric study we will let Γ = 5 × 10−4 which
corresponds with a quality factor of Q = 1000 when the resonator is unbiased and
undergoing small oscillations. This value of the quality factor is commonly obtainable
for N/MEMS devices [81]. For the slew rate parameter we set g2 = 500, and for the
noise strengths we set Dfb,1 = Dfb,2 = 10Denv .
Figure 4.5 shows the logarithm of the spectrum of fractional frequency fluctuations
evaluated at zero frequency, Sy (0), normalized by Denv . This quantity represents
the long term phase diffusion coefficient. Figure 4.6 shows the corresponding mean
oscillator frequencies. In both plots a dashed white line indicates the locus of zero
106

4.5

0

4
3.5

-1

3
-2

2.5
2

-3
-4

1.5
1
0

0.5

1
1.5
Bias, λ

Log10 Sy (0)/Denv

Log Action, Log10 (I)

5

2

Figure 4.5: The logarithm of the normalized long term phase diffusion coefficient,
log10 Sy (0)/Denv , as a function of operating action and bias force. Dashed white
line is the locus of zero dispersion (∂I ω(I) = 0) points. For interpretation of the
references to color in this and all other figures, the reader is referred to the electronic
version of this dissertation.

dispersion points. Figure 4.5 is remarkable in its complex topography. The key
feature is the narrow valley around the zero dispersion line for large action. Recall
that we noted, based upon the scaling results found above, that SI,I dominates in
the expression for Sy for small action. This is why the narrow valley around the zero
dispersion curve fades away for small action. Meanwhile, for large action, we observed
that Sϕ,ϕ should dominate. Accordingly, we see a reduction in the phase diffusion
coefficient near the zero dispersion line, where the phase is protected from fluctuations
in the action. In the parameter range considered the phase diffusion coefficient can
be reduced by at most 5dB by these tuning parameters. This is consistent with
observations made by Chris Burgner at the University of Santa Barbara (private
communication) on a real device. Improvement on the order of 10’s of decibels will
likely require attacking the noises at their source.
107

0

2.2
2

-1
1.8
-2

1.6

Frequency, Ω

Log Action, Log10 (I)

2.4

1.4
-3
1.2
-4

0

0.5

1
1.5
Bias, λ

2

1

Figure 4.6: Mean oscillator frequency as a function of operating action and bias
force. Dashed white line is the locus of zero dispersion (∂I ω(I) = 0) points. For
interpretation of the references to color in this and all other figures, the reader is
referred to the electronic version of this dissertation.

4.5

Outloook

Equation (4.64) is the primary result of this chapter. It can be employed as a tool to
guide the effort to engineer better N/MEMS oscillators. The equation can serve as a
cost function for oscillator optimization or to help target primary noise sources. In
this work, we demonstrated the utility of this expression with a hypothetical oscillator
model. Building a model matched to experimental observations of a real oscillator
may prove more informative. In addition, an interesting problem that could be investigated using this formulation is to determine the optimum potential for an oscillator.
Equation (4.64) evaluated at ω = 0 would provide the cost function for the optimization problem. It is likely that the answer will depend on the noise model and so an
appropriate model would need to be generated. Such an effort is beyond the scope of
this work, but could be the next step.

108

Chapter 5
Conclusions
In this work we examined several noise-induced phenomena that are to be found
in nano- and micro-electromechanical systems. Our aim has been to exploit these
phenomena, where possible, for system identification and sensing applications. The
exception being our discussion of nonlinear oscillators in chapter 4 where our goal was
to remove the effect of noise as much as possible. In general, however, our viewpoint
has been one of embracing the noise and nonlinearity that are an inevitable part of
miniaturization.
In chapter 2 we examined delayed subcritical pitchfork and saddle-node bifurcations in the presence of noise. This investigation was born out of an attempt to port,
to MEMS sensors, the bifurcation detection methods used in the 1970s to locate the
critical current in Josephson junctions. During the initial stages of our investigation
the issue of sweep rate become important, and it was not clear what it meant to
sweep fast or slow. In the present work we determined that “fast” means fast with
respect to the noise strength, and that this is indeed the situation in many MEMS
devices, leading to delayed bifurcations. The experimental demonstration of our analytic results met with limited success. It turned out that the sweeping experiments
took so long that drift in the device parameters became problematic. Despite efforts
109

to remove this effect by measuring the device frequency before each parameter sweep,
this effect may have influenced the data. In addition the sweep rate required for
the normal form to remain valid appears to have been so slow that we swept clean
though the local region. These issues arose in large part because these experiments
were performed before the analysis was complete. A more systematic experiment on
a MEMS device has yet to be done, and is left for future work.
In chapter 3 we examined the balanced dynamical bridge as a general use detector and a detector of non-Gaussian noise. This dynamic bridge sensing strategy
was inspired by the exponential dependence of the escape rate on device parameters
and the success of the Josephson bifurcation amplifier. As a general use detector,
we found that the bridge requires a long time to produce a measurement. The problem is that the system encodes the information in the occupation probabilities and
measuring these probabilities is time-consuming, requiring the record of many events.
The detection of non-Gaussian noise suffers from this same problem. However, using
the bridge in this manner can be slightly faster since the non-Gaussian noise can
increase the switching rate. Moreover, other methods of detecting the non-Gaussian
character of a stochastic process operate on similar principles, and so here the dynamic bridge becomes competitive. We discussed an example where we considered
shot noise measured by a one-dimensional bridge. Future work on this project may
include a comparison of the dynamical bridge sensitivity to other detection methods.
In addition, one may consider a dynamical bridge constructed from a parametric
resonator. The inherent symmetry in parametric resonance will automatically balance the occupation probabilities and create a device that is sensitive to symmetry
breaking perturbations such as direct forcing at the resonant frequency.
In chapter 4 we examined the spectrum of fractional frequency fluctuations in an
oscillator with a nonlinear frequency selective element. We ignored the dynamics of
the feedback loop elements, focusing on the resonator equation. We developed an
110

expression for the spectrum of fractional frequency fluctuations as a function of the
oscillator noise model and the nominal resonator orbit. We proposed that such an
expression could be useful to tune an oscillator to minimize the rate of phase diffusion.
We considered a prototype model based on the biased Duffing resonator with which to
demonstrate this oscillator tuning. We showed how the oscillator can be tuned using
the resonator zero-dispersion points to decouple the action and angle fluctuations.
The minimum phase diffusion constant w.r.t. the oscillator operating action was
found to be near the zero-dispersion points for sufficiently large action. The example
system showed reduction in the phase diffusion rate of up to 5 decibels by operating
at this point. Future work on this subject of oscillator design can be divided into two
separate projects. First, a detailed model of an existing oscillator can be constructed
in an effort to optimally tune that oscillator. Second, one could attempt to minimize
the phase diffusion constant by shaping the limit cycle. The tools we employed would
limit the class of limit cycles to those arising from a one degree-of-freedom potential
system. However, other perspectives would allow the consideration of more general
limit cycles [27, 80]. This problem naturally depends on the noise model employed,
and so these two suggested directions for future research are closely related.

111

Appendices

112

Appendix A
Stochastic Processes
In this appendix we discuss some of the tools for the analysis of dynamic systems
with random terms. This discussion provides some background for the material in
the body of this dissertation. The interested reader is directed to references [93, 61, 60]
for additional material.

A.1

The Kinetic Equation for a Stochastic Process

Consider a stochastic process, X(t), where X ∈ Rd . We’ll use X(t) to denote the
stochastic process itself and x to denote possible values that X(t) can take at various
times. Now, let ρ(x, t) be the probability distribution for X(t) at time t. Thus,
ρ(x, t)dx is the probability that X(t) = x. Also, let P (x, t + τ |x , t) be the probability
that X(t+τ ) = x given that X(t) = x ; this is known as the transition probability.
It is the probability that the process X(t) will make a transition from x to x over
the time [t, t + τ ]. A Markov process is one whose future value depends only on
the present, not the past. We can see that P (x, t + τ |x , t) takes only two arguments,
the value of X at the present time t and the value of X at the future time t + τ . So,
if X(t) is Markovian, then P (x, t + τ |x , t) contains the complete description of the
113

process. Using the transition probability, P , we can determine the evolution of the
probability distribution ρ. The probability for X(t + τ ) = x can be understood as
the the sum of the probabilities for the process, X, to transition to x over the time
[t, t + τ ] from all values x at time t. That is, ρ(x, t + τ ) is the probability to get to
x in the time interval [t, t + τ ] from anywhere. Thus, we can write

ρ(x, t + τ ) =

P (x, t + τ |x , t)ρ(x , t)dx .

(A.1)

We can make use of (A.1) to develop a partial differential equation for ρ(x, t).

We begin by expanding the integrand for small x − x . We will assume that x is a
d-dimensional real vector, x ∈ Rd . Let ∆ = x − x , then we can write the integrand
of equation (A.1) as

P (x, t + τ |x , t)ρ(x , t) = P (x + ∆, t + τ |x , t)ρ(x , t),
(−∆1 )n1 · · · (−∆d )nd ∂ n1 +···+nd
=
n
n
n1 ! · · · nd !
∂x1 1 · · · ∂x d
n1 ,··· ,nd
d
×P (x + ∆, t + τ |x, t)ρ(x, t),

(A.2)

(A.3)

where we have expanded x about x in a Taylor series and the sum for ni goes from
0 to ∞. Substitution of this into equation (A.1) gives

ρ(x, t + τ ) =

(A.4)
∂ n1 +···+nd
(−1)n1 +···+nd n
n ρ(x, t)
∂x1 1 · · · ∂x d
n1 ,··· ,nd
d
nd
n1
∆1 · · · ∆
d P (x + ∆, t + τ |x, t)d∆.
×
n1 ! · · · nd !
114

Now, the integral in equation (A.4) can be written as
d
i=0

n
∆i i
P (x + ∆, t + τ |x, t)d∆ =
ni !

d
i=0

1
[Xi (t + τ ) − xi ]ni
ni !

, (A.5)
Xi (t)=xi

which are the moments of the probability distribution of X(t + τ ) − X(t) given that
X(t) = x, weighted by 1/(n1 ! · · · nd !). The collection of these quantities describe the
probability distribution of the next step that the stochastic process, X(t), takes over
the time τ if it takes on the value x at time t. We now introduce the Kramers-Moyal
coefficients, D(n1 ···nd ) (x, t).
1
D(n1 ···nd ) (x, t) = lim
τ →0 τ

d
i=0

1
[X (t + τ ) − xi ]ni
ni ! i

.

(A.6)

Xi (t)=xi

These are the coefficients of the first order term in an expansion of equation (A.5) in
τ about τ = 0, so long as at least one ni is nonzero. When all ni are zero the integral
equation (A.5) is equal to 1 by the normalization of P . Thus, for small τ , we have
∂m
ρ(x, t)D(n1···nd ) (x, t)τ.
ρ(x, t + τ ) = ρ(x, t) +
(−1)m n
1 · · · ∂xnd
∂x1
n1 ,··· ,nd
d
m=n1 +···+nd =0
(A.7)

Next we subtract ρ(x, t) from both sides of equation (A.7) and divide by τ . In the
limit of τ → 0 we have the kinetic equation of the stochastic process X(t), given
by
∂ρ(x, t)
∂ n1 +···+nd (n1···n )
d (x, t)ρ(x, t).
=
(−1)n1 +···+nd n
n D
∂t
∂x1 1 · · · ∂x d
n1 ,··· ,nd
d
n1 +···+nd =0

(A.8)

We now further discuss the meaning of the Kramers-Moyal coefficients. There are d
115

coefficients for which only one of the integers ni = 1 and the rest are zero. Changing
(1)
to a more compact notation, we can write this sub set as Di ≡ D(0,0,...1,...,0)
(1)
where the 1 occurs in the ith place. In Di the superscript denotes the order of
the coefficient and the subscript denotes which ni = 1. This set of coefficients give
the mean step made by the stochastic process X(t) over a small time normalized by
the time. Thus these coefficients are called the drift vector or drift coefficient.
(2)
Similarly, we can consider the coefficients of order 2. We can label these as Dij ≡
D(0,0,...1,0,0,...0,1,...0) denoting the second order coefficient where ni = 1 and nj = 1
if i = j and ni = 2 if i = j. These coefficients describe the variance of the step made
by the stochastic process over a small time normalized by this time. Thus D(2) is
called the diffusion coefficient.
Now, equation (A.8) is pretty unwieldy when it represents an infinite series. Calculating the infinite number of coefficients D(n1 ···nd ) could be difficult, if not impossible. However, the Pawula Theorem says that the series D(n1 ···nd ) will terminate
at either the first or the second order (i.e. n1 + n2 + · · · nd = 1 or 2) or, if it does
not, it must have an infinite number of terms in it. Thus, one only has to check the
third order term to find out if the series terminates at second order or is infinite. This
holds under the condition that ρ(x, t) is to be always positive, which it must be, since
its a probability distribution.
As an example, consider the standard one-dimensional Brownian motion, a.k.a.
a Wiener process, defined by
D(1) = 0

D(2) = 1

D(n) = 0,

n ≥ 3.

(A.9)

In this case equation (A.8) and its solution are

ρ=ρ
˙

⇒

ρ(x, t) =

1
(x − x(0))2
exp −
.
4πt
4t
116

(A.10)

The case where the Kramers-Moyal coefficients are nonzero only at first and second order is very interesting and useful. Now we know one process that fits in that
category, Brownian motion, but what about a stochastic process that represents the
response of a dynamic system? As we will see in the next section, when a dynamic
system is forced by Gaussian white noise, it’s Kramers-Moyal expansion is also truncated at second order. So, all that is left to do is to calculate these coefficients for our
system of interest. To do this, all we have to do is solve our system equation over a
short time interval τ , and then take the limit according to (A.6). Solving our system
over a small time interval τ , however, is where the trouble comes in. It turns out that
the answer depends on how we interpret the noise. We will discuss this further in the
next two sections. We begin with a “physical” interpretation of the noise, assuming
that it is smooth on time scales much shorter than the time scale of the evolution of
our system.

A.2

The Langevin Equation

Suppose we have the system

˙
Xi = fi (X, t) + gij (X, t)ξj (t).

(A.11)

This is our Langevin equation. The time functions ξj (t) are assumed to be zero
mean random functions that are smooth on a sufficiently short time scale. This
smoothness comes from the ’physical’ interpretation, arguing that real quantities are
smooth. Smoothness implies that the standard rules of calculus apply. Now, we want
to integrate our equation over a short time τ with initial condition X(t) = x. So we
117

integrate (A.11) w.r.t. time
t+τ
Xi (t + τ ) − xi =

fi (X(t ), t ) + gij (X(t ), t )ξj (t ) dt ,

t

(A.12)

and apply an iterative approach, assuming f and g are sufficiently well behaved
functions. First we approximate f and g as

f( X(t ), t ) ≈ fi (x, t ),

(A.13)

gij (X(t ), t ) ≈ gij (x, t ).

(A.14)

Applying these approximations to equation (A.12) gives the zeroth order approximation,
t+τ

(0)
Xi (t + τ ) − xi =

t

fi (x, t ) + gij (x, t )ξj (t ) dt .

(A.15)

∂fi (x, t ) (0)
(X (t ) − xk ),
k
∂xk

(A.16)

Going to the next order, we take

fi (X(t ), t ) ≈ fi (x, t ) +

gij (X(t ), t ) ≈ gij (x, t ) +

∂gij (x, t )
∂xk

(X

(0)
(t ) − xk ),
k

(A.17)

to get the first order approximation
(1)
Xi (t + τ ) − xi =

t+τ
t

fi (x, t ) + gij (x, t )ξj (t )

t ∂f (x, t )
∂f (x, t )
i
fk (x, t ) + i
gkj (x, t )ξj (t )
∂xk
∂xk
t
∂gij (x, t )
∂gij (x, t )
+
fk (x, t )ξj (t ) +
gkl (x, t )ξj (t )ξl (t ) dt
∂xk
∂xk

(A.18)

+

118

dt .

This is a sufficient order for our purposes. Now, we take the expectation or ensemble
average, then we take f (x, t ) ≈ f (x, t) and g(x, t ) ≈ g(x, t), and find
t+τ
Xi (t + τ ) − xi

fi (x, t )dt

=
t
t+τ
+
t

(A.19)

t ∂gij (x, t )
gkl (x, t ) ξj (t )ξl (t ) dt dt + O(τ 2 ).
∂xk
t

Now, we assume that the noise processes ξi are wide sense stationary. This means
that ξi (t)ξj (t ) depends only on the time difference t − t . Moreover, we assume
that the processes have short memory compared to the rate of change of X. Thus,
ξi (t)ξj (t ) is a narrow function of t − t that goes to zero for t − t

τ . It is

reasonable, therefore, to make the approximation ξi (t)ξj (t ) = 2Dij δ(t − t ). The
δ-correlation implies that the power spectral density is constant, so the noise is white.
Unfortunately, things are not quite that simple. In equation (A.19) we integrate δ(t)
up to t = 0 from below and thus the result of the integration depends on how we
define the δ function. For example, we could assume

δ(t) = lim
h→0

θ(t + h) − θ(t − h)
,
2h

(A.20)

where θ(t) is the unit step function. In this case integration over half the δ-function
gives 1/2. Alternatively, we could assume

δ(t) = lim
h→0

θ(t) − θ(t − h)
.
h

(A.21)

In this case the integration gives zero. The choice we make for the definition of the
δ-function must reflect the noise properties. We will continue with our discussion

119

assuming the former case so the integration gives 1/2. Thus, we obtain

Xi (t + τ ) − xi

= fi (x, t)τ +

∂gij (x, t)
∂xk

gkl (x, t)Djl τ + O(τ 2 ).

(A.22)

Next, we consider the second moments of Xi (t + τ ) − xi .
[Xi (t + τ ) − xi ][Xj (t + τ ) − xj ] =
(A.23)
t+τ t+τ
gik (x, t )gkl (x, t ) ξk (t )ξl (t ) dt dt + · · ·
t
t
= 2τ gik (x, t)gjl (x, t)Dkl + O(τ 2 ).
The third and higher moments are all higher order in τ ,
[Xi (t + τ ) − xi ][Xj (t + τ ) − xj ][Xk (t + τ ) − xk ] = O(τ 2 ).

(A.24)

Thus we find that a dynamic system subject to white Gaussian noise is Markovian
with Kramers-Moyal coefficients
∂gij (x, t)
(1)
Di (x, t) = fi (x, t) +
gkl (x, t)Djl ,
∂xk
(2)
Dij (x, t) = gik (x, t)gjl (x, t)Dkl ,
D(n) (x, t) = 0, n ≥ 3.

(A.25)
(A.26)
(A.27)

Plugging these coefficients into equation (A.8) gives the Fokker-Planck equation.

A.3

Ito and Stratonovich Integrals

The stochastic process ξ(t) is not a proper function when it is truely δ-correlated
because it is discontinuous at every point. Rather than deal with this, mathematicians
prefer to deal with the Brownian motion or the Wiener process we saw earlier. And
120

so, rather than write a Langevin equation like (A.11), they write the stochastic
differential equation (SDE)

dXi = fi (X, t)dt + gij (X, t)djk dWk ,

(A.28)

where dWk is the increment of a Brownian motion. The noise strength and correlation
has been explicitly written in the matrix d. Note that, to agree with the above,
ddT = D. Since a Brownian motion is continuous, the increments dWk are well
defined. Now, if we integrate the SDE (A.28) with initial conditions X(t) = x we get
t+τ

t+τ
Xi (t + τ ) − xi =

fi (X(t ), t )dt +

t

gij (X(t ), t )djk dWk (t ). (A.29)

t

Applying the same iterative procedure as above, we find
(0)
Xi (t + τ ) − xi =

t+τ

t+τ
fi (x, t )dt +

t

t

gij (x, t )djk dWk (t ),

(A.30)

and
(1)
Xi (t + τ ) − xi =

t+τ
fi (x, t )dt + gij (x, t )djk dWk (t )

t
t

+
t
+
+
+

(A.31)

∂fi (x, t )
fk (x, t )dt dt
∂xk

∂fi (x, t )
gkj (x, t )djl dWl (t )dt
∂xk
∂gij (x, t )
∂xk
∂gij (x, t )
∂xk

fk (x, t )djl dt dWl (t )
gkl (x, t )djn dlm dWm (t )dWn (t )

121

.

As in the Langevin approach discussion, we take f (x, t ) ≈ f (x, t) and g(x, t ) ≈
g(x, t), and find
(1)
Xi (t + τ ) − xi = fi (x, t)τ + gij (x, t)djk
+

∂gij (x, t)
∂xk

τ
0

dWk (t )
τ

gkl (x, t)djn dlm

0

(A.32)

t
0

dWm (t )dWn (t ) + · · · .

We have pushed all terms that won’t contribute into the · · · and shifted the limits on
the integrals w.r.t. dW because the Brownian motion increment dW is stationary,
that is its statistics don’t depend on the the location of the time origin. The reason
why we needed to shift the limits of integration is because the definition of these
stochastic integrals is written with these limits. Now, these stochastic integrals can
be of either the Ito or Stratonovich type, depending on interpretation. These integrals
are defined as follows. Let 0 = t0 < t1 < t2 ... < tN −1 < tN = τ and ∆ = ti+1 − ti .
Then an integral of the Ito type of a function Φ(W (t), t) with respect to dW (t) is
defined by
N −1
Φ(W (ti ), ti ) W (ti+1 ) − W (ti ) , (A.33)
Φ(W (t), t)dW (t) = lim
∆→0
0
i=0
τ

II =

and an integral of the Stratonovich type is defined by
τ
IS =

Φ(W (t), t)dW (t)
N −1

=

(A.34)

0
lim
∆→0

Φ
i=0

W (ti+1 ) + W (ti ) ti+1 + ti
,
2
2

122

W (ti+1 ) − W (ti ) .

In both integrals, W (t) is a standard Wiener process or Brownian motion defined as

W (0) = 0,

W (t) = 0,

(A.35)

Var(W (t) − W (s)) = 2t − 2s.

Var(W (t)) = 2t,

(A.36)

τ
Now if we consider the term 0 dW (t ) in (A.32), we have Φ = 1 and so both Ito and
Stratonovich give the same answer
N −1
dW (t ) = lim
[W (ti+1 ) − W (ti )] = W (τ ).
∆→0
0
i=0
τ

(A.37)

So, the double intergal in (A.32) gives
t

τ

τ
dW (t )dW (t ) =

0

0

W (t )dW (t ).

(A.38)

0

What we need is the expectation of this integral. The Stratonovich definition gives
N −1
=
W (t )dW (t )
[W (ti+1 + W (ti ))][W (ti+1 ) − W (ti )]
0
S
i=0
N −1
N −1
1
2 (t
2 (t ) = lim
= lim
ti+1 − ti = τ.
W i+1 ) − W i
∆→0 2
∆→0
i=0
i=0
(A.39)
τ

1
lim
∆→0 2

The Ito definition gives
N −1

τ
W (t )dW (t )
0

=
I
=

lim
∆→0

W (ti )[W (ti+1 ) − W (ti )]

(A.40)

i=0

1
1
lim
W 2 (τ ) −
2
∆→0 2

123

N −1
i=0

[W (ti+1 ) − W (ti )]2 = 0.

Applying equations (A.37)-(A.40) to (A.32) we obtain our expression for X(t+τ )−x.

Xi (t + τ ) − xi = fi (x, t)τ +

∂gij (x, t)
∂xk

τ
gkl (x, t)djn dlm

0

Wm (t)dWn (t)

,
I,S
(A.41)

with
τ
0
τ
0

Wm (t)dWn (t)
Wm (t)dWn (t)

= 0,

(A.42)

= τ δnm .

(A.43)

I
S

Thus, we have
(1)
Di,Ito = fi (x, t),

∂gij (x, t)
(1)
Di,Strat = fi (x, t) +
gkl (x, t)Djl .
∂xk

(A.44)

It can also be shown that the diffusion coefficient, D2 , is the same for both. The
Stratonovich interpretation matches the Langevin approach of the previous section
with definition (A.20) and the Ito interpretation matches the Langevin approach
with definition (A.21). The choice of stochastic calculus to be used should be made
by considering the origin and character of the noise processes. In [61] there is a
chapter on this topic. Ito calculus appears better suited for numeric solutions of
stochastic differential equations [61] because the equations can be directly integrated.
Stratonovich calculus would require an implicit scheme. When performing simulations of Stratonovich stochastic systems it is usually necessary to write an equivalent
Ito SDE. To do this one would take the drift and diffusion coefficients from the
Stratonovich system and write an Ito SDE that gives rise to these same coefficients
using Ito calculus.

124

Appendix B
Stochastic Averaging
In this appendix we discuss the method of averaging applicable to dynamic systems
with random terms. This method should apply to systems with non-white and/or
non-Gaussian noise. The approach is meant to be as simple as possible, employing
only basic understanding of stochastic processes. No proof is given however, just an
explanation of the method and the thought processes behind it.
We are interested in equations of the form

xi =
˙

fi (x, t) + 1/2 gij (x, t)ξj (t),

(B.1)

where ξ represents some stochastic process, the functions f and g are periodic, or
at least have discrete spectrum, and

is a small parameter. Note that this scaling

is consistent with the microscopic resonator model in section (4.2). The noise ξ is
assumed to be continuous and differentiable so that conventional calculus rules apply.
In other words, equation (B.1) is to be interpreted in the Stratonovich sense. The
goal is to simplify this system by pushing the non-essential aspects of the system
to higher orders in

by means of a near-identity coordinate transformation. These

non-essential terms are then ignored by virtue of their smallness, recognizing that
125

the truncated equations will provide a good approximation over a long but generally
finite time scale. In addition to the above, we assume that the stochastic process ξ is
wide-sense stationary with zero mean, that is

ξi (t)

= 0,

(B.2)

= rij (τ ),

(B.3)

rij (τ )e−iωτ dτ.

(B.4)

ξi (t)ξj (t + τ )
and, we introduce the power spectrum of ξ,
∞
Rij (ω) =

−∞

For this averaging process only the second order statistics are needed, so we will not
specify ξ any further. As stated above, we assume the time variation of f and g is
made up of several discrete frequencies. Thus we write

fi (y, t) =
gij (y, t) =

n
n

(n)
Fi (y)eiΩn t ,

(B.5)

(n)
Gij (y)eiΩn t ,

(B.6)

(−n)
(n)
(−n)
(n)
with Fi
= (Fi )∗ , Gij
= (Gij )∗ , and Ω−n = −Ωn . Here we must make
a point of notation. We make use of the convention of summation over repeated
indices written as subscripts, and we often use i as an index. So, when i appears as
a subscript it is meant as an index. However, when it appears elsewhere it is meant
as the imaginary unit. We now proceed with the coordinate transformation.

126

B.1

A Near-Identity Coordinate Transformation

In order to simplify equation (B.1), we would like to change to a coordinate system
in which the system appears as simple as possible, given the smallness of . To do
this, we begin by defining the general near-identity coordinate transformation
xi = yi + 1/2 Ui (y, t) + Vi (y, t).

(B.7)

The functions U and V are left unspecified for the present. The problem of defining
these two functions is the subject of the next two sections. In this section, we simply
solve for the equation of motion of y with these general functions. We do this by
first taking the time derivative of (B.7) and equating it with (B.1) in which we have
replaced x with y via a Taylor series expansion
∂U (y, t)
∂Vi (y, t)
˙
˙
xi = yi + 1/2 i
˙
˙
yj +
˙
yj + 1/2 Ui (y, t) + Vi (y, t),
˙
∂yj
∂yj
∂gij (y, t)
Uk (y, t)ξj (t) + · · · .
= fi (y, t) + 1/2 gij (y, t)ξj (t) +
∂yk

(B.8)
(B.9)

From (B.8) and (B.9) we solve for y up to O( ).
˙

yi =
˙

1/2 g ξ − U +
˙
ij j
i

fi +

∂gij

∂Ui ˙
∂Ui
˙
U − Vi
Uk ξj −
gjk ξk +
∂yk
∂yj
∂yj j

.
(B.10)

Now, it remains to choose U and V such that equation (B.10) becomes as simple
as possible. The constraint on the choice is that the coordinate transformation,
equation (B.7), must remain well ordered. Thus, we must choose U and V such that
they remain bounded. It is by this condition that we identify the principle part of
equation (B.1) that cannot be reduced.
127

B.2

Stochastic Averaging Part I: O(

1/2

)

To begin, we consider the function U . The role of this function is to replace the
¯
stochastic process g(y, t)ξ(t) with something simpler, call it ξ(y, t). The process gξ
¯
is wide-sense cyclo-stationary. Thus, the most general form of ξ then is also a widesense cyclo-stationary process. We don’t go through the exercise here, but it can
be shown that the principle part of gξ, that is the part that must be contained in
¯
¯
ξ, need only be wide-sense stationary. This is done by taking ξ to be a general
wide-sense cyclo-stationary process and calculating U 2 . One then finds that the
¯
secular terms contain only the wide-sense stationary component of ξ and its (cyclostationary) correlation with ξ. We therefore need only seek a wide-sense stationary
¯
process ξ to replace gξ. We can construct such a process by defining our “slow noise”
¯
process, ξ, as a low-pass filtering of gξ. Thus, we take

¯
˙
Ui (y, t) = gij (y, t)ξj (t) − ξi (y, t),
∞
¯
h(t − t )gij (y, t )ξj (t ) dt ,
ξi (y, t) =
−∞

(B.11)
(B.12)

where h(t) is the time-domain filtering function, the filter’s ‘impulse response’. Note
¯
that equation (B.12) defines the simplified noise ξ as a functional of the original
¯
process, ξ. Thus, in principle, the statistics of ξ may be calculated for non-white
and/or non-Gaussian ξ. The definition of this coordinate transformation, however,
does not require a causal filtering function and so we are free to choose h(t) = 0 for
¯
t < 0. By a proper choice of the filter we can make ξ wide-sense stationary. For
example, we may choose an ideal low-pass filter with appropriate cut-off frequency.
Our filter may also be characterized by the frequency domain filtering function,
∞
H(ω) =

−∞
128

h(t)e−iωt dt.

(B.13)

For an ideal low-pass filter, this is simply a rectangular pulse, in ω space, centered
zero. We now wish to show that U remains bounded, at least with high probability.
2
To do this, we demonstrate that Ui is bounded for every i. For this to be so, the
tails of the probability distribution for U (y, t) should decay sufficiently fast and so the
probability for U to be very large is very small. In order to facilitate the calculation,
¯
it is convenient to first calculate the second moment of ξ and its correlator with ξ.
From equation (B.12) we find

¯
ξj (t)ξi (y, s)

=

¯
¯
ξi (y, t)ξj (y, s)

=

∞ (n)
1
G R (ω)H(Ωn − ω)eiω(t−s) eiΩn s dω,
2π −∞ ik kj
n
(B.14)

and
∞ (n) (m)
1
G G
Rk (ω)H(Ωn − ω)H(Ωm + ω)
ik j
n,m 2π −∞
×eiω(s−t) eiΩn t eiΩm s dω.

(B.15)

Equation (B.15) is of particular importance since it is the autocorrelation of the slow
¯
noise, ξ. From this equation we can see how the oscillations of g select the noise
around the frequencies Ωn . The oscillating function g modulates components of the
noise ξ down to baseband. These low frequency components of the modulated noise
are then selected by the low-pass filter. Moreover, we see that if H(Ωm − ω) and
¯
H(Ωn + ω) are orthogonal, except in the case m = −n, then the slow noise ξ is
wide-sense stationary. Returning to the ideal low-pass example, suppose H(ω) is
a rectangular window centered at ω = 0 and with width less than the minimum
distance between any pair of frequencies Ωn and Ωm , (m = n). In this case, only the
wide-sense stationary terms in equation (B.15) survive.
˙
It remains to show that our choice of U in equation (B.11) is valid. To do this we
129

will calculate U 2 and show that it remains bounded in time. Thus, from equation
(B.11), we take
t
Ui =

−∞

¯
gij (y, t )ξj (t ) − ξi (y, t ) dt ,

(B.16)

and calculating the variance gives
2
Ui

t
=

t

−∞ −∞

¯
gij (y, t )gik (y, s ) ξj (t )ξk (s ) − gij (y, t ) ξj (t )ξi (y, s )

¯
¯
¯
−gij (y, s ) ξj (s )ξi (y, t ) + ξi (y, t )ξi (y, s )

dt ds .

(B.17)

Applying equations (B.14) and (B.15) and moving into the frequency domain yields
2
Ui

=

∞
1
(n) (m)
dω Gij G
Rkj (ω) 1 − 2H(Ωm − ω)
ik
2π −∞
n,m

(B.18)

+H(Ωm − ω)H(Ωn + ω) ei(Ωn +Ωm )t
0

0
×

−∞ −∞

ei(Ωn +ω)t ei(Ωm −ω)s ds dt .

The time integrals over t and s can be understood as Fourier transforms of the unit
step,
0

1
eiat dt =
+ πδ(a).
ia
−∞

(B.19)

Now, the secular terms in equation (B.18) are those for which m = −n. In these terms
the poles and delta functions double up to make the integration over ω unbounded.
However, these terms are removed by the filtering on the condition that H(0) = 1.

130

This can be assumed without any restriction on the approach, in which case we find
2
Ui

ei(Ωn +Ωm )t (n) (m)
G G
Rkj (−Ωn ) 1 − H(Ωn + Ωm ) ,
i(Ωn + Ωm ) ij ik
n,m
m=−n
(B.20)
einΩt
(m) (n−m)
=
1 − H(nΩ)
Gij G
Rkj (−mΩ) . (B.21)
ik
inΩ
m
n=0
=

2
Thus, we find that Ui is given by a Fourier series and it is bounded so long as this
(n)
series converges. For example, if Gij diminish sufficiently fast with large n.

B.3

Stochastic Averaging Part II: O( )

With the choice of U (y, t) derived in the previous section, we now move onto the
selection of the deterministic term V (y, t). With our selection of U (y, t) equation
(B.10) now becomes

yi =
˙

1/2 ξ (y, t) +
¯
i

˙
fi (y, t) + ηi (y, t) + di (y, t) − Vi (y, t) + · · · , (B.22)

where η + d is the O( ) noise, η is the fluctuating part and d is the expectation, given
by

ηi (y, t) =
di (y, t) =

∂gij

∂Ui ¯
Uk ξj −
ξ − di (y, t),
∂yk
∂yj j
∂gij
∂Ui ¯
Uk ξj −
ξ .
∂yk
∂yj j

(B.23)
(B.24)

We have explicitly separated d from η because η will be dropped when we scale time
(in the next section), but d, the term of interest, will be retained.
In order to simplify equation (B.22) as much as possible, understanding that η
131

will be dropped, we take

¯
¯
˙
Vi (y, t) = fi − fi + di − di ,
∞
¯
h(t − t )fi (y, t ) dt ,
fi (y) =
−∞
∞
¯
h(t − t )di (y, t ) dt ,
di (y, t) =
−∞

(B.25)
(B.26)
(B.27)

where h is the filter function from the previous section. It is clear that V will remain
bounded since it is the integral of oscillating functions with their low frequency components removed by the filter. We are therefore satisfied that this choice of V is valid
¯
and proceed to calculate d. Evaluating the two parts of equation (B.24) gives

∂gij

U ξ
∂yk k j

(m)
∂Gij
1 − H(Ωn − ω) i(Ωn +Ωm )t
(n)
dω
e
,
=
G R j (ω)
k
∂yk
2πi(Ωn − ω)
−∞
n,m
(B.28)
∞

and

∂Ui ¯
ξ
∂yj j

(m)
∂Gij
(n)
dω
=
G R j (ω)H(Ωn − ω)
k
∂yk
n,m −∞
1 − H(Ωm + ω) i(Ωn +Ωm )t
×
e
.
2πi(Ωm + ω)
∞

(B.29)

Subtracting equation (B.29) from equation (B.28) and then filtering the difference
¯
gives the expression for d,
(m)
∂Gij
1 − H(Ωn − ω)
(n)
¯
dω
di (y, t) =
G R j (ω)H(Ωn + Ωm )
k
∂yk
2πi(Ωn − ω)
n,m −∞
1 − H(Ωm + ω) i(Ωn +Ωm )t
−H(Ωn − ω)
e
.
(B.30)
2πi(Ωm + ω)
∞

132

Equation (B.30) simplifies significantly if H(Ωn + Ωm ) = δm,−n . In this case we
have

¯
di (y, t) =

(−n)
∂Gij
1 − |H(ω)|2
(n)
dω
G R j (Ωn − ω)
.
k
∂yk
2πiω
−∞
∞

n

(B.31)

¯
If R j is even and real di (y, t) becomes

(−n)
∂G
1 − |H(ω)|2
(n)  ∞
 ij
¯
di (y, t) =
Im 
G 
dωR j (ω − Ωn )
.
k
∂yk
πω
−∞
n=1
∞



(B.32)
.

¯
From equation (B.31) or (B.32) it is easy to see that d is the contribution to the slow
dynamics produced by the noise away from the frequencies Ωn . For white noise, for
¯
which R j = const, the weights are all zero and so d is zero. In the more general
case, however, it may nonzero.

B.4

Time Scaling

With the coordinate transformation chosen, we scale time to put all important terms
at O(1). Thus, we introduce the slow time τ ,

τ = t.

(B.33)

Assuming we have removed all oscillating components by the coordinate transformation,the equation of motion for y now takes the form
dyi
¯
¯
¯
= fi (y) + di (y) + −1/2 ξi (y, −1 τ ) + ηi (y, −1 τ ) + O( 1/2 ).
dτ
133

(B.34)

¯
We replace the noise terms ξ and η with noises defined in the new time τ as follows.
Let
¯
ζi (y, τ ) ≡ −1/2 ξi (y, −1 τ ),

(B.35)

ηi (y, τ ) ≡ −1/2 ηi (y, −1 τ ).
ˆ

(B.36)

We introduce the factor of 1/2 to that ζ and η maintain the same spectral magnitude.
ˆ
The signal power, however, increases. This can be understood by imagining that the
noises appear to be ‘more white’ over the long time scale τ . To see this, consider the
spectrum, Sζ ,
∞
−1 τ
1
dω Sξ (y, ω )eiω
,
(B.37)
¯
2π −∞
∞
dτ ζ(y, t)ζ(y, t + τ ) e−iωτ ,
(B.38)
Sζ (y, ω) =
−∞
∞
∞
−1 ω −ω)τ
1
, (B.39)
dω dτ Sξ (y, ω )ei(
=
¯
2π −∞ −∞
∞
(B.40)
dω Sξ (y, ω )δ(ω − ω),
=
¯
−∞

ζ(y, t)ζ(y, t + τ )

=

Sζ (y, ω) = Sξ (y, ω).
¯

(B.41)

¯
Thus the spectrum of ζ is the stretched spectrum of ξ. It is similar for η. In this
way, we keep the spectral strength of the noise processes at a unit level, and let the
coefficients multiplying these processes capture their scale. Accordingly, the process
η becomes 1/2 η with the scaling of time, and is thus deemed small in comparison to
ˆ
the other terms. We therefore drop η .
ˆ

The resulting averaged equations are
dyi
¯
¯
= fi (y) + di (y) + ζi (y, τ ) + O( 1/2 ),
dτ
134

(B.42)

¯ ¯
where f +d are the deterministic parts of the approximate system, defined in equations
(B.26) and (B.31), and ζ is a non-Markovian stochastic process, for which the relevant
equations are (B.12), (B.15), (B.35), and (B.41).
Since our coordinate transformation remains well ordered, and the neglected terms
are small, a solution of (B.42) should remain close to our original system, (B.1), for
some time until the neglected terms in y build up. This time is expected to be
˙
O( −1/2 ). Moreover, again since the coordinate transformation is well ordered, a
stationary distribution of (B.42) is expected to be O( 1/2 ) close to a stationary
distribution of (B.1).

B.5

Approximation by a Markov Process

Now the averaged equation (B.42) contains a non-Markovian process, ζ, and so y is
also non-Markovian. However, in some cases we may be able to approximate y as
Markovian. If we suppose the cut-off frequency of the filter H(ω) used to generate
equation (B.42) is given by ωc , then ζ might be uncorrelated over times ∆ >> 2π /ωc .
For example, suppose that R(ω) is flat around the frequencies Ωn , so that ζ is a bandlimited white noise. In this case its autocorrelation is given by
ωc
sinc
ζi (τ )ζj (τ + ∆) ∝
π
Since

ωc ∆

.

(B.43)

is small, this is a sharp delta-like function which goes approximately to zero

for times of O( 1/2 ), over which the system state y has changed very little. We can
therefore approximate y as a Markov process. We can describe this Markov process by
its Kramers-Moyal coefficients. These coefficients are the moments of an infinitesimal
135

step of the process defined for a process z(t) as
1
.
D(n) = lim
(z(t + τ ) − z(t))n
n! τ
z(t)=z
τ →0

(B.44)

In this section we will calculate the first two coefficients. These are all that is required
for a Gaussian noise. For a non-Gaussian noise all coefficients are required for a
complete description. However, it may be that higher-order coefficients are higherorder in . This suggests that in the limit of small , the noise may become Gaussian
as well as white, which in turn suggests a type of of central limit result. This can
be viewed in a discrete sense, noting that for very small , each step in coordinate
results from the sum of many noise pulses, and the central limit theorem says such a
sum should have a Gaussian distribution.

Returning to our calculation, in our approximation of y as Markovian, we cannot
take the limit to zero as was done in equation (B.44), but rather we must take it
to a small time of O( 1/2 ). Thus, to approximate D(n) , we must first calculate
y(τ + 1/2 ∆) − y(τ ) given y(τ ) = y (0) . We do this by the method of successive
integration. We first integrate equation (B.42)

(0)
yi (τ + 1/2 ∆) − yi =
τ

τ + 1/2 ∆

¯
¯
fi (y(t )) + di (y(t )) + ζi (y(t ), t ) dt ,
(B.45)

and then expand the integrand in a Taylor series in y about y (0) . For example, for
¯
f,
¯
∂ f (y (0) )
(0)
¯
¯
(yj (t ) − yj ),
fi (y(t )) ≈ fi (y (0) ) + i
∂yj

(B.46)

and similarly for the other terms. We then insert expression (B.45) into equation
136

(B.46) and repeat the Taylor series expansion. This process yields

(0)
yi (τ + 1/2 ∆) − yi
=

1/2 ∆ f (y (0) ) + d (y (0) ) +
¯
¯
i
i

τ + 1/2 ∆
τ

τ

ζi (y (0) , t )dt

∂ζi (y (0) , t )
dt
ζj (y (0) , t ) + O( ).
∂yj
τ
t

dt

+

τ + 1/2 ∆

(B.47)

The first Kramers-Moyal coefficient, the drift coefficient is thus approximated by
(1)
Di
≈

1
y (τ + 1/2 ∆) − yi (τ )
1/2 ∆ i

(B.48)
y(τ )=y

(−n)
∂Gij
(n) ∞
¯
¯
= fi (y) + di (y) +
G
dωR j ( ω + Ωn )H( ω)
k
∂yk
−∞
n

(B.49)

1/2 ∆ω
1 + i 1/2 ∆ω − ei
.
×
2π 1/2 ∆ω 2
Since we asserted that R(ω) is approximately constant around the frequencies Ωn ,
we have
(−n)
∂Gij
(n)
G R j (Ωn )
k
∂yk
n
1/2 ∆ω
ωc /
1 + i 1/2 ∆ω − ei
.
×
dω
−ωc /
2π 1/2 ∆ω 2

(1)
¯
¯
Di
≈ fi (y) + di (y) +

(B.50)

Moreover, at the limits of integration, the integrand in equation (B.50) is O( 3/2 ).
We therefore extend the limits to ±∞ and the value of the integration is 1/2. Thus,
we find
(−n)
∂Gij
1
(1)
(n)
¯
¯
Di
≈ fi (y) + di (y) +
G R j (Ωn ).
k
2 n
∂yk

137

(B.51)

The second Kramers-Moyal coefficient is approximated by
1
(B.52)
(yi (τ + 1/2 ∆) − yi (τ ))(yj (τ + 1/2 ∆) − yj (τ ))
1/2 ∆
2
y(τ )=y
ωc / 1 − cos( 1/2 ∆ω)
1
(n) (−n)
≈
G G
Rk (Ωn )
.
(B.53)
2 n ik j
−ωc
π 1/2 ∆ω 2

(2)
Dij ≈

Again, we let the integration limits go to ±∞ and the integration evaluates to 1,
resulting in
1
(n) (−n)
(2)
G G
Rk (Ωn ).
Dij ≈
2 n ik j

(B.54)

¯
Note that when the noise in the original equation is white, then d = 0 and the drift
and diffusion coefficients are simply those belonging to the original equation averaged
over one period. Other cases are significantly more challenging.

138

Appendix C

Derivation of Equation (2.65)

In this appendix we derive equation (2.65), the stochastic normal form for the pitchfork bifurcation arising from parametric resonance. Our main source for this development is reference [62]. We start with the equation for a parametrically forced
resonator,
2
z + 2Γz + ωo (1 + δ1 + δ1 cos 2ωt)z + γ(1 + δ3 + δ3 cos 2ωt)z 3 = ξ(t),
¨
˙
ξ(t)ξ(t ) = 2Dδ(t − t ).

ξ(t) = 0,

(C.1)
(C.2)

In principle this equation could include all terms of the form (δn + δn cos 2ωt)xn for
n = 0, 1, 2, 3, but only the n = 1, 3 terms survive the averaging, and so only these are
included here. We first transform this equation into a form amenable for stochastic
averaging. We begin by rescaling time as

τ = Γt,
ξ(τ )ξ(τ )

= 2DΓδ(τ − τ ),
139

(C.3)
(C.4)

which gives the rescaled equation
2
Γ2 z + 2Γ2 z + ωo (1 + δ1 + δ1 cos (2ωτ /Γ))z
+γ(1 + δ3 + δ3 cos (2ωτ /Γ))z 3 = ξ(τ ),

(C.5)

where (·) = d(·)/dτ . Next, the system is transformed into scaled Cartesian coordinates (a, b) that rotate with ω/Γ, defined via

z =
z

=

8ωΓ
(a cos(ω/Γτ ) + b sin(ω/Γτ )),
3|γ|

(C.6)

8ωΓ
(−aω/Γ sin(ω/Γτ ) + bω/Γ cos(ω/Γτ )).
3|γ|

(C.7)

Applying this substitution and solving for a and b gives




 a 
 =

b





3|γ|  − sin(ω/Γτ ) 


8ω 3 Γ3
cos(ω/Γτ )

2
− 2Γ2 z + ω 2 z − ωo (1 + δ1

+δ1 cos (2ωτ /Γ))z − γ(1 + δ3 + δ3 cos (2ωτ /Γ))z 3 + ξ(τ ) . (C.8)
In appendix B we showed that the slow time drift and diffusion coefficients for a
system like (C.8) can be approximated by their time-average. The averaged drift and

140

diffusion coefficients are
¯ (1) = −a −
D1

2 δ ω2
2
ω 2 − ωo − δ1 ωo
+ 1 o
2ωΓ
4ωΓ

¯ (1) = −b +
D2

1
+ sgn(γ)δ3 b a2 + b2 ,
3
2
2 δ ω2
ω 2 − ωo − δ1 ωo
− 1 o a − sgn(γ)a(a2 + b2 )
2ωΓ
4ωΓ

b + sgn(γ)b(a2 + b2 )

− sgn(γ)δ3 a
¯ (2)
D11
¯ (2)
D12
¯ (2)
D21
¯ (2)
D22

=

3D|γ|
,
16ω 3 Γ2

(C.9)

(C.10)

5 2
a + b2 ,
3
(C.11)

= 0,

(C.12)

= 0,

(C.13)

=

3D|γ|
.
16ω 3 Γ2

(C.14)

To simplify the expressions for the drift coefficient, we define the nondimensional
parameters
2
ω 2 − ωo (1 + δ1 )
,
2ωΓ
2
δ1 ωo
β =
.
4ωΓ

α =

(C.15)
(C.16)

The parameter α is the renomalized detuning, accounting for the shift in the resonant
frequency from δ1 . β is the effective parametric pump strength and |β| = 1 is the
minimum parametric pump required to induce parametric resonance. For the next
several steps it is convenient to work with the Langevin description of the system.
We therefore write a Langevin system that gives the drift and diffusion coefficients
above. In addition, we now use the overdot to denote derivatives w.r.t. τ . The

141

averaged Langevin equations are
1
a = −a − (α + β)b + sgn(γ)b(a2 + b2 ) + sgn(γ)δ3 b a2 + b2 + ξ1 (τ ), (C.17)
˙
3
5 2
˙
a + b2 + ξ2 (τ ),(C.18)
b = −b + (α − β)a − sgn(γ)a(a2 + b2 ) − sgn(γ)δ3 a
3
3D|γ|
ξi = 0,
ξi (τ )ξj (τ ) =
δij δ(τ − τ ),
(C.19)
8ω 3 Γ2
which describe the evolution of (a, b) in convenient form. Next, we reduce the dimension of this system by restricting it to the slow invariant manifold that exists near
the bifurcation point. We begin by applying a linear change of coordinates to the
eigendirections. The linear deterministic part of this system is

a = −a − (α + β)b,
˙

(C.20)

˙
b = (α − β)a − b,

(C.21)

β 2 − α2 and eigenvectors

which has eigenvalues λ± = −1 ±
1
v+ = (−
β
1
v− = (
β

β(β + α)
,
2
β(β + α)
,
2

β−α
),
2β

(C.22)

β−α
).
2β

(C.23)

This motivates the linear coordinate transformation,

u =
v = −

β
2β(β + α)

a+

β
2β(β + α)

a+

β
b,
2(β − α)
β
b,
2(β − α)

142

a=

1
β

b=

β(β + α)
(u − v),(C.24)
2
β−α
(u + v),
2β

(C.25)

that will render the linear part of the equation in normal form. The full equations in
these coordinates are given by

u = (−1 −
˙

sgn(γ)

β 2 − α2 )u −

(3αβ + δ3 (α2 + β 2 + 3αβ))u3 (C.26)

β 2 − α2
−(6αβδ3 + 6α2 (1 + δ3 ) + 3β 2 (1 + δ3 ))u2 v
3β

+(3α2 δ3 + 3β 2 δ3 + 9αβ(1 + δ3 ))uv 2
−β(2αδ3 + 3β(1 + δ3 ))v 3 +
v = (−1 +
˙

β
2β(β + α)

sgn(γ)

β 2 − α2 )v −

β 2 − α2
−(3α2 δ3 + 3β 2 δ3 + 9αβ(1 + δ3 ))u2 v

ξ1 (τ ) +

β
ξ (τ ),
2(β − α) 2

β(2αδ3 + 3β(1 + δ3 ))u3

(C.27)

3β

+(6αβδ3 + 6α2 (1 + δ3 ) + 3β 2 (1 + δ3 ))uv 2
−(3αβ + δ3 (α2 + β 2 + 3αβ))v 3 −

β
2β(β + α)

ξ1 (τ ) +

β
ξ (τ ).
2(β − α) 2

As is typical, the linear terms are uncoupled, but at the expensive of more complicated
nonlinear terms. We will need the deterministic center manifold u = h(v), determined
from the equation u|h = ∂h v|h with ξ1 and ξ2 set to zero. This gives, to leading
˙
˙
∂v
order,

u=

sgn(γ)(3β + δ3 (2α + 3β)) 3
v .
2 − α2 ) − 6 β 2 − α2
12(β

143

(C.28)

Next we write the Fokker-Planck equation for the system in (u, v) coordinates. For
ease of notation, let

u = F1 (u, v) + g1 ζ1 (τ ) + g2 ζ2 (τ ),
˙
ˆ
ˆ

(C.29)

v = F2 (u, v) − g1 ζ1 (τ ) + g2 ζ2 (τ ),
˙
ˆ
ˆ

(C.30)

g1 = β
ˆ
g2 =
ˆ

3D|γ|
1
,
3 Γ2 β(β + α)
32ω
3D|γ|
β
,
3 Γ2 β − α
32ω

(C.31)

(C.32)

ζ1 (τ )

=

ζ2 (τ ) = ζ1 (τ )ζ2 (τ ) = 0,

(C.33)

ζ1 (τ )ζ1 (τ )

=

ζ2 (τ )ζ2 (τ ) = 2δ(τ − τ ),

(C.34)

and the Fi ’s are the deterministic part of the equations. Moreover, by defining
g1 = g1 + g2 ,
ˆ2 ˆ2
g2 = 2(ˆ2 − g1 ).
g 2 ˆ2

(C.35)
(C.36)

the Fokker-Planck equation for the probability distribution of u and v, denoted as P,
is given by
∂P
∂
∂
∂ 2P
∂ 2P
∂ 2P
= − (F1 P) −
(F2 P) + g1
+ g2
+ g1
.
∂τ
∂u
∂v
∂u∂v
∂u2
∂v 2

(C.37)

To restrict the system to be near the suspended center manifold, we assume that P
is confined to a narrow strip around u = h(v), specifically, P is taken to be a narrow
tube distribution with mean along the deterministic center manifold and Gaussian

144

cross section, as follows,

P = p(u, τ |v)P (v, τ ),
p(u, τ |v) =

Λ(v)
exp − Λ(v)[u − h(v)]2 .
π

(C.38)
(C.39)

To formulate the FP equation for P we will need the partial derivatives of p(u, τ |v)
up to second order. These are
∂p
∂u
∂p
∂v
∂ 2p
∂u2
∂ 2p
∂u∂v
∂ 2p
∂v 2

= −2Λ(u − h)p,
=

(C.40)

Λ
p + (u − h)p(2Λh − Λ (u − h)),
2Λ

= −2Λp + 4Λ2 (u − h)2 p,

(C.41)
(C.42)

= 2Λh p + (u − h)p 2ΛΛ (u − h)2 − 3Λ − 4Λ2 h (u − h) , (C.43)
Λ2
Λ2
Λ
p + (u − h)p 6Λ h − (u − h)
p+
2Λ
Λ
4Λ2
+(u − h)[2Λh − (u − h)Λ ]2 + 2Λh − (u − h)Λ .
(C.44)

= −2Λh 2 p −

Substitution of the form of P given in equation (C.38) into the Fokker-Planck equation
gives

p

∂P
∂τ

∂
∂p
∂
∂p
(F1 P ) − F1 P
− p (F2 P ) − F2 P
+
(C.45)
∂u
∂u
∂v
∂v
∂ 2P
∂ 2P
∂p ∂P
∂ 2p
∂p ∂P
∂ 2p
p
g1 + p
g1
+2
+P
+2
+P
∂u ∂u
∂v ∂v
∂u2
∂u2
∂v 2
∂v 2
∂ 2P
∂p ∂P
∂p ∂P
∂ 2p
+ p
+
+
+P
g .
∂u∂v ∂u ∂v
∂v ∂u
∂u∂v 2

= −p

145

To get an equation for P (v, τ ), we integrate this equation over u, treating p(u, τ |v)
like a δ(u − h(v)). This yields
∂P
∂τ

∂F1
∂
Λ
∂ 2P
∂P
Λ
−
(F2 P ) −
F2 P − 2Λg1 P + g1
+ g1
∂u
∂v
2Λ
Λ ∂v
∂v 2
Λ2
Λ
−2Λh 2 g1 P −
g P + 2Λh g2 P
.
(C.46)
g1 P +
2Λ 1
4Λ2
u=h
−P

=

Next, we assume Λ = Λ0 + Λ1 v + Λ2 v 2 , and use scaling arguments to greatly simplify
this equation. For notational ease, we also express F1 and F2 as
F1 = a1 u + b1 u3 + c1 u2 v + d1 uv 2 + e1 v 3 ,

(C.47)

F2 = a2 v + b2 u3 + c2 u2 v + d2 uv 2 + e2 v 3 .

(C.48)

The scaling employed is motivated by proximity to the bifurcation point and weak
noise, and is given by

v → v,

(C.49)

a2 → 2 a2 ,

(C.50)

gi → 4 gi ,

(C.51)

Λi → −4 Λi .

(C.52)

Keeping terms up to 2 gives
∂P
∂τ

=

2

−

∂
∂ 2P
[(a2 v + e2 v 3 )P ] + g1
∂v
∂v 2

+P

− 2Λ0 g1 − a1 − 2 Λ1 g1 v

+ 2 (−d1 v 2 + 6Λ0 a3 g2 v 2 − 2Λ2 g1 v 2 ) .

(C.53)

To ensure that the P equation represents a conservation of probability along the
center manifold, it is necessary that all terms proportional to P vanish. Accordingly,
146

we take
a
Λ0 = − 1 ,
2g1
Λ1 = 0,

(C.54)
(C.55)

d
6Λ0 a3 g2
− 1.
2g1
2g1

Λ2 =

(C.56)

The resulting Fokker-Planck equation is given by
∂P
∂
=−
∂τ
∂v

a2 v + e 2 v 3 P + g 1

∂ 2P
,
∂v 2

(C.57)

which corresponds with the Langevin equation
v = a2 v + e 2 v 3 +
˙

√

g1 ξ(τ ).

(C.58)

We scale this according to
√
e2
q=
v,
2

(C.59)

resulting in
q = a2 q + 4q 3 +
˙

e2 g1
ξ(τ ).
4

(C.60)

Substitution of the values of a2 , e2 , and g1 into the above gives the final normal
form, expressed as

q = (−1 +
˙

β 2 − α2 )q + 4q 3 +

ξ = 0,

Dγ αβ 2 (3+2δ3 )+δ3 β(β+α)2
ξ(τ ),
64ω 3 Γ2
(β 2 −α2 )3/2
(C.61)

ξ(τ )ξ(τ ) = 2δ(τ − τ ).

147

(C.62)

Bibliography

148

BIBLIOGRAPHY

[1] Milton Abramowitz and Irene A. Stegun. Handbook of Mathematical Functions:
with Formulas, Graphs, and Mathematical Tables. Dover, 1965.
[2] J.S. Aldridge and A.N. Cleland. Noise-enabled precise measurements of a duffing nanomechanical resonator. Phys Rev Lett, 94:156403, 2005.
[3] D. Allan, H. Hellwig, P. Kartaschoff, J. Vanier, J. Vig, G.M.R. Winkler, and
N.F. Yannoni. Standard terminology for fundamental frequency and time
metrology. In Frequency Control Symposium, 1988., Proceedings of the 42nd
Annual, pages 419–425, June 1988.
[4] Vinay Ambegaokar and B. I. Halperin. Voltage due to thermal noise in the dc
josephson effect. Phys Rev Lett, 22:1364–1366, 1969.
[5] A. D. Armour. Current noise of a single-electron transistor coupled to a nanomechanical resonator. Phys. Rev. B, 70:165315, 2004.
[6] A. D. Armour, M. P. Blencowe, and Y. Zhang. Classical dynamics of a
nanomechanical resonator coupled to a single-electron transistor. Phys Rev
B, 69:125313, 2004.
[7] D. A. Bagrets and Yu. V. Nazarov. Full counting statistics of charge transfer
in coulomb blockade systems. Phys. Rev. B, 67:085316, 2003.
[8] Nils Berglund and Barbara Gentz. Pathwise description of dynamic pitchfork
bifurcations with additive noise. Probability Theory and Related Fields, 122:341–
388, 2002.
[9] Nils Berglund and Barbara Gentz. Noise-Induced Phenomena in Slow-Fast
Dynamical Systems: A Sample-Paths Approach. Springer, 2006.
[10] Lora Billings, M. I. Dykman, Marie McCrary, and Ira B. Schwartz. Switching
barrier scaling near bifurcation points for non-gaussian noise. Phys. Rev. Lett.,
104:140601, September 2010.
149

[11] Lora Billings, Mark I. Dykman, and Ira B. Schwartz. Thermally activated
switching in the presence of non-gaussian noise. Phys. Rev. E., 78:051122,
2008.
[12] Lora Billings, Ira B. Schwartz, Marie McCrary, A. N. Korotkov, and M. I.
Dykman. Switching exponent scaling near bifurcation points for non-gaussian
noise. Phys. Rev. Lett., 104:140601, 2010.
[13] G. Binnig, Ch. Gerber, E. Stoll, T. R. Albrecht, and C. F. Quate. Atomic resolution with atomic force microscope. Europhysics Letters, 3:1281–1286, 1987.
[14] G. Broggi, A. Colombo, L. A. Lugiato, and Paul Mandel. Influence of white
noise on delayed bifurcations. Physical Review A, 33:3635–3637, 1986.
[15] Eyal Buks and Bernard Yurke. Mass detection with a nonlinear nanomechanical
resonator. Physical Review E, 74:046619, 2006.
[16] C. B. Burgner, K. L. Turner, N. J. Miller, and S. W. Shaw. Parameter sweep
strategies for sensing using bifurcations in mems. In Proceedings of the 13th
Hilton Head Solid State Sensors and Actuators Conference, June 6-10, SC,
USA, 2010.
[17] Christopher Burgner, William Snyders, and Kimberly Turner. Control of mems
on the edge of instability. In IEEE Transducers, Beijing, China, June, 2011.
[18] Paul F. Byrd and Morris D. Friedman. Handbook of Elliptic Integrals for Engineers and Physicists. Springer-Verlag / Berlin, 1954.
[19] J. Carr. Applications of Centre Manifold Theory. Springer; 1 edition, 1981.
[20] H. B. Chan, M. I. Dykman, and C. Stambaugh. Paths of fluctuation induced
switching. Phys Rev Lett, 100:130602, 2008.
[21] H. B. Chan, M. I. Dykman, and C. Stambaugh. Switching-path distribution in
multidimensional systems. Physical Review E, 78:051109, 2008.
[22] Shui-Nee Chow and Hao-Min Zhou. An analysis of phase noise and fokkerplanck equations. J. Differential Equations, 234:391–411, 2007.
[23] A. N. Cleland and M. R. Geller. Superconducting qubit storage and entanglement with nanomechanical resonators. Phys Rev Lett, 93:070501, 2004.
[24] A. N. Cleland and M. R. Geller. Mechancial quantum resonators. Electronic Properties of Novel Nanostructures: XIX International Winterschool/Euroconference on Electronic Properties of Novel Materials. AIP Conference Proceedings, 786:396–400, 2005.
[25] A. N. Cleland and M. L. Roukes. Noise processes in a nanomechanical resonators. Journal of Applied Physics, 92:2758–2769, 2002.
150

[26] Andrew N. Cleland. Foundations of Nanomechanics. Springer-Verlag, 2003.
[27] Alper Demir, Amit Mehrotra, and Jaijeet Roychowdhury. Phase noise in oscillatiors: A unifying theory and numerical methods for characterization. IEEE
Transactions on circuits and systems - I: Fundamental Theory and Applications,
47:655–674, 2000.
[28] Michel H. Devoret, Daniel Esteve, John M. Martinis, Andrew Cleland, and
John Clarke. Resonant activation of a brownian particle out of a potential well:
Microwave-enhanced escape from the zero-voltage state of a josephson junction.
Phys. Rev. B., 36(1):58–73, July 1987.
[29] O. K. Dudko, A. E. Filippov, J. Klafter, and M. Urbakh. Beyond the conventional description of dynamic force spectroscopy of adhesion bonds. PNAS
September 30, 100:11378–11381, 2003.
[30] M. Dykman. Quantum measurements with dynamically bistable systems. In
Visarath In, Patrick Longhini, and Antonio Palacios, editors, Applications of
Nonlinear Dynamics, volume 35 of Understanding Complex Systems, pages 367–
378. Springer Berlin / Heidelberg, 2009.
[31] M. I. Dykman. Poisson-noise induced escape from a metastable state. Phys.
Rev. E, 81:051124, 2010. (in press).
[32] M. I. Dykman, B. Golding, L. I. McCann, V. N. Smelyanskiy, D. G. Luchinsky,
R. Mannella, and P. V. E. McClintock. Activated escape of periodically driven
systems. Chaos, 11:587–594, 2001.
[33] M. I. Dykman, M. Khasin, J. Portman, and S.W. Shaw. Spectrum of an oscillator with jumping frequency and the interference of partial susceptibilities.
Phys. Rev. Lett, 105:230601, 2010.
[34] M. I. Dykman and M. A. Krivoglaz. Classical theory of nonlinear oscillators
interacting with a medium. Phys. Stat. Sol. (b), 48:497–512, 1971.
[35] M. I. Dykman and M. A. Krivoglaz. Theory of fluctuational transistions between
stable states of a nonlinear oscillator. Sov. Phys. JETP, 50:30–37, 1979.
[36] M. I. Dykman and M. A. Krivoglaz. Theory of nonlinear oscillator interacting
with a medium. Sov. Sci. Rev. A Phys., 5:265–442, 1984.
[37] M. I. Dykman, C. M. Maloney, V. N. Smelyanskiy, and M. Silverstein. Fluctuational phase-flip transitions in parametrically driven oscillators. Phys Rev E,
57:5202–5212, 1998.
[38] M. I. Dykman, R. Mannella, P. V. E. McClintock, S. M. Soskin, and N. G.
Stocks. Noise-induced narrowing of peaks in the power spectra of underdamped
nonlinear oscillators. Phys Rev A, 42:7041–7049, 1990.
151

[39] Albert Einstein. Investigations on the theory of the Brownian Movement. Dover
Publications, 1956.
[40] K. L. Ekinci, X. M. H. Huang, and M. L. Roukes. Ultrasensitive nanoelectromechanical mass detection. App. Phys. Lett., 84:4469 – 4471, 2004.
[41] K. L. Ekinci and M. L. Roukes. Nanoelectromechanical systems. Review of
Scientific Instruments, 76:061101, 2005.
[42] K. L. Ekinci, Y. T. Yang, and M. L. Roukes. Ultimate limits to inertial mass
sensing based upon nanoelectromechanical systems. Journal of Applied Physics,
95:2682–2689, 2004.
[43] Kilho Eom, Harold S. Park, Dae Sung Yoon, and Taeyun Kwon. Nanomechanical resonators and their applications in biological/chemical detection: Nanomechanics principles. Physics Reports, 503:115–163, 2011.
[44] Mykhaylo Evstigneev. Statistics of forced thermally activated escape events out
of a metastable state: Most probable escape force and escape-force moments.
Phys. Rev. E, 78:011118, 2008.
[45] X.L. Feng, C.J. White, A. Hajimiri, and M.L. Roukes. A self-sustaining
ultrahigh-frequency nanoelectromechanical oscillator. Nature Nanotechnology,
3:342–346, 2008.
[46] Richard P. Feynman and Albert R. Hibbs. Quantum Mechanics and Path Integrals: Emended Edition (Dover Books on Physics). Dover Publications, 2010.
[47] T. A. Fulton and L. N. Dunkleberger. Lifetime of the zero-voltage state in
josephson tunnel junctions. Phys Rev B, 9:4760, 1974.
[48] Anupam Garg. Escape-field distribution for escape from a metastable potential
well subject to a steadily increasing bias field. Phys. Rev. B, 51(21):15592–
15595, June 1995.
[49] Cevat G¨k¸ek. Tracking the resonance frequency of a series rlc circuit using a
o c
phase locked loop. Control Applications, 2003. CCA 2003. Proceedings of 2003
IEEE Conference on, 1:609 – 613, 2003.
[50] Herbert Goldstein. Classical Mechanics, Second Edition. Addison-Wesley, 1980.
[51] A. C. R. Grayson, R. S. Shawgo, A. M. Johnson, N. T. Flynn, Li Yawen, M. J.
Cima, and R. Langer. A biomems review: Mems technology for physiologically
integrated devices. Proceedings of the IEEE, 92:6–21, 2004.
[52] J. Guckenheimer and P. Holmes. Nonlinear Oscillations, Dynamical Systems,
and Bifurcations of Vector Fields. Springer-Verlag New York, 1983.
152

[53] S. Gustavsson, R. Leturcq, B. Simovi˘, R. Schleser, T. Ihn, P. Studerus, K. Enc
sslin, D. C. Driscoll, and A. C. Gossard. Counting statistics of single electron
transport in a quantum dot. Phys Rev Lett, 96:076605, 2006.
[54] Peter Hanggi. Escape from a metastable state. Journal of Statistical Physics,
42:105–148, 1986.
[55] Tero T. Heikkil¨, Pauli Virtanen, G¨ran Johansson, and Frank K. Wilhelm.
a
o
Measuring non-gaussian fluctuations through incoherent cooper-pair current.
Phys. Rev. Lett., 93:247005, 2004.
[56] B. Huard, H. Pothier, Norman O. Birge, D. Esteve, X. Waintal, and J. Ankerhold. Josephson junctions as detectors for non-gaussian noise. Ann. Phys.
(Leipzig), 16(10-11):736–750, 2007.
[57] G. Hummer and A. Szabo. Kinetics from nonequilibrium single-molecule pulling
experiments. Biophysical Journal, 85:5–15, 2003.
[58] L. D. Jackel, W. W. Webb, J. E. Lukens, and S. S. Pei. Measurement of the
probability distribution of thermally excited fluxoid quantum transitions in a
superconducting ring closed by a josephson junction. Phys Rev B, 9:115, 1974.
[59] Kalvis M. Jansons and G. D. Lythe. Stochastic calculus: Application to
dynamic bifurcations and threshold crossings. Journal of Statistical Physics,
90:227–251, 1998.
[60] N.G. Van Kampen. Stochastic Processes in Physics and Chemistry. Elsevier
Academic Press, 2007.
[61] Peter E. Kloeden and Eckhard Platen. Numerical Solution of Stochastic Differential Equations. Springer, 3rd edition, 1999.
[62] E. Knobloch and K. A. Wiesenfeld. Bifurcations in fluctuating systems: The
center-manifold approach. J. of Stat. Phys., 33(3):611–637, 1983.
[63] R. H. Koch, G. Grinstein, G. A. Keefe, Yu Lu, P. L. Trouilloud, and W. J. Gallagher. Thermally assisted magnetization reversal in submicron-sized magnetic
thin films. Phys Rev Lett, 84:5419, 2000.
[64] H. A. Kramers. Brownian motion in a field of force and the diffusion model of
chemical reactions. Physica VII, 7(4):284–304, 1940.
[65] I. N. Krivorotov, N. C. Emley, A. G. F. Garcia, J. C. Sankey, S. I. Kiselev, D. C.
Ralph, and R. A. Buhrman. Temperature dependence of spin-transfer-induced
switching of nanomagnets. Phys. Rev. Lett., 93:166603, 2004.
153

[66] Vijay Kumar, J. William Boley, Yushi Yang, Hendrik Ekowaluyo, Jacob K.
Miller, George T.-C. Chiu, and Jeffery F. Rhoads. Bifurcation-based mass
sensing using piezoelectrically actuated microcantilevers. App. Phys. Lett., 98:1,
2011.
[67] Juhani Kurkij¨rvi. Intrinsic fluctuations in a superconducting ring closed with
a
a josephson junction. Phys. Rev. B, 6(3):832–835, August 1972.
[68] M. D. LaHaye, J. Suh, P. M. Echternach, K. C. Schwab, and M. L. Roukes.
Nanomechanical measurements of a superconducting qubit. Nature, 459:960–
964, 2009.
[69] L. D. Landau and E. M. Lifshitz. Mechancics, Third Edition: Course of Theoretical Physics Volume 1. Elsevier Academic Press, 2008.
[70] Q. Le Masne, H. Pothier, Norman O. Birge, C. Urbina, and D. Esteve. Asymmetric noise probed with a josephson junction. Phys. Rev. Lett., 102:067002,
January 2009.
[71] Sangwoo Lee, Sangjun Park, Jongpal Kim, Sangchul Lee, and Dong
il (Dan) Cho. Surface/bulk micromachined single-crystalline-silicon microgyroscope. Microelectromechanical Systems, Journal of, 9:557–567, 2000.
[72] Thomas H. Lee and Ali Hajimiri. Oscillator phase noise: A tutorial. IEEE
Journal of Solid-State Circuits, 35:326–336, 2000.
[73] R. Lifshitz and M. C. Cross. Nonlinear dynamics of nanomechanical and micromechanical resonators. In Reviews of Nonlinear Dynamics and Complexity.
Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim, Germany., 2009.
[74] Ron Lifshitz and M. L. Roukes. Thermoelastic damping in micro- and nanomechanical systems. Phys Rev B, 61:5600–5609, 2000.
[75] Chang Liu. Foundations of MEMS. Pearson Education, 2006.
[76] A. Lupascu, S. Saito, T. Picot, P. C. De Groot, C. J. P. M. Harmans, and J. E.
Mooij. Quantum non-demolition measurement of a superconducting two-level
system. nature physics, 3:119–123, 2007.
[77] I. Mahboob, C. Froitier, and H. Yamaguchi. A symmetry-breaking electromechanical detector. App. Phys. Lett., 96:213103, 2010.
[78] Francois Mallet, Florian R. Ong, Agustin Palacios-Laloy, Franois Nguyen,
Patrice Bertet, Denis Vion, and Daniel Esteve. Sinle-shot qubit readout in
circuit quantum electrodynamics. Nature Physics, 5:791–795, 2009.
[79] C. Meunier and A. D. Verga. Noise and bifurcations. Journal of Statistical
Physics, 50:345–375, 1988.
154

[80] Jeff Moehlis. On the precision of noisy oscillators. Department of Mechanical
Engineering, University of California, Santa Barbara, CA, Preprint, 2012.
[81] P. Mohanty, D. A. Harrington, K. L. Ekinci, Y. T. Yang, M. J. Murphy, and
M. L. Roukes. Intrinsic dissipation in high-frequency micromechanical resonators. Phys Rev B, 66:085416, 2002.
[82] F. Mohd-Yasin, D. J. Nagel, and C. E. Korman. Noise in mems. Measurement
Science and Technology, 21:012001, 2010.
[83] Martin H. M¨ser, Michael Urbakh, and Mark O. Robbins. Statistical mechanics
u
of static and low-velocity kinetic friction. In I. Prigogine and Stuart A. Rice,
editors, Advances in Chemical Physics, volume 126, pages 187–272. John Wiley
& Sons Inc., 2003.
[84] E. B. Myers, F. J. Albert, J. C. Sankey, E. Bonet, R. A. Buhrman, and D. C.
Ralph. Thermally activated magnetic reversal induced by a spin-polarized current. Phys Rev Lett, 89:196801, 2002.
[85] N. Sri Namachchivaya. Stochastic bifurcation. Applied Mathematics and Computation, 38:101–159, 1990.
[86] N. Sri Namachchivaya and Y. K. Lin. Method of stochastic normal forms.
International Journal of Non-Linear Mechanics, 26:931–943, 1991.
[87] Anatoly Neishtadt. On stability loss delay for dynamical bifurcations. Discrete
and Continuous Dynamical Systems Series S, 2:897–909, 2009.
[88] Tom´ˇ Novotn´. Josephson junctions as threshold detectors of full counting
as
y
statistics: open issues. Journal of Statistical Mechanics: Theory and Experiment, 2009:P01050, 2009.
[89] L. A. Oropeza-Ramos, C. B. Burgner, and K. L. Turner. Robust micro-rate
sensor actuated by parametric resonance. Sensors and Actuators A: Physical,
152(1):80–87, 2009.
[90] H.W.C. Postma, I. Kozinsky, A. Husain, and M.L. Roukes. Dynamic range
of nanotube- and nanowire-based electromechanical systems. Applied Physics
Letters, 86:223105, 2005.
[91] Michael V. Requa and Kimberly L. Turner. Precise frequency estimation in
a microelectromechanical parametric resonator. App. Phys. Lett., 90:173508,
2007.
[92] Jeffery F. Rhoads, Steven W. Shaw, and Kimberly L. Turner. Nonlinear dynamics and its applications in micro- and nanoresonators. Journal of Dynamic
Systems, Measurement, and Control, 132:034001, 2010.
155

[93] H. Risken. The Fokker-Planck Equation: Methods of Solution and Applications.
Springer, 2nd edition, 1996.
[94] M. L. Roukes. Nanoelectromechanical systems face the future. Physics World,
14:25–31, 2001.
[95] Y. Sang, M. Dub´, and M. Grant. Thermal effects on atomic friction. Phys Rev
e
Lett, 87:174301, 2001.
[96] I. Serban, M. I. Dykman, and F. K. Wilhelm. Relaxation of a qubit measured
by a driven duffing oscillator. Phys. Rev. A, 81:022305, 2010.
[97] I. Siddiqi, R. Vijay, M. Metcalfe, E. Boaknin, L. Frunzio, R. J. Schoelkopf, and
M. H. Devoret. Dispersive measurements of superconducting qubit coherence
with a fast latching readout. Phys Rev B, 73:054510, 2006.
[98] I. Siddiqi, R. Vijay, F. Pierre, C. M. Wilson, M. Metcalfe, C. Rigetti, L. Frunzio,
and M. H. Devoret. Rf-driven josephson bifurcation amplifier for quantum
measurement. Phys Rev Lett, 93:207002, 2004.
[99] J. Sieber, A. Gonzalez-Buelga, S. A. Neild, D. J. Wagg, and B. Krauskopf.
Experimental continuation of periodic orbits through a fold. Phys. Rev. Lett.,
100:244101, 2008.
[100] V. N. Smelyanskiy, M. I. Dykman, H. Rabitz, and B. E. Vugmeister. Fluctuations, escape, and nucleation in driven systems: Logarithmic susceptibility.
Phys Rev Lett, 79:3113–3116, 1997.
[101] S. M. Soskin, R. Mannella, and P. V. E. McClintock. Zero-dispersion phenomena
in oscillatory systems. Physics Reports, 373:247–408, 2003.
[102] Michael J. Stephen. Noise in a driven josephson oscillator. Physical Review,
186:393–397, 1969.
[103] N. G. Stocks, R. Mannella, and P. V. E. McClintock. Influence of random
fluctuations on delayed bifurcations: The case of additive white noise. Phys
Rev A, 40:5361–5369, 1989.
[104] N. G. Stocks, R. Mannella, and P. V. E. McClintock. Influence of random
fluctuations on delayed bifurcations. ii. the cases of white and colored additive
and multiplicative noise. Phys Rev A, 42:3356–3362, 1990.
[105] Eugene V. Sukhorukov and Andrew N. Jordan. Stochastic dynamics of a josephson junction threshold detector. Phys Rev Lett, 98:136803, 2007.
[106] J. Z. Sun, J. C. Slonczewski, P. L. Trouilloud, D. Abraham, Ian Bacchus,
W. J. Gallagher, J. Hummel, Yu Lu, and G. Wright. Thermal activationinduced sweep-rate dependence of magnetic switching astroid. App. Phys. Lett.,
78(25):4004–4006, 2001.
156

[107] Masuo Suzuki. Anomalous fluctuation and relaxation in unstable systems. Journal of Statistical Physics, 16:477–504, 1977.
[108] Masuo Suzuki. Scaling theory of non-equilibrium systems near the instability
point. iii. Progress of Theoretical Physics, 57:380–392, 1977.
[109] J. Tobiska and Yu. V. Nazarov. Josephson junctions as threshold detectors for
full counting statistics. Phys Rev Lett, 93:106801, 2004.
[110] M. C. Torrent and M. San Miguel. Stochastic-dynamics characterization of
delayed laser threshold instability with swept control parameter. Phys Rev A,
38:245–251, 1988.
[111] M. C. Torrent, F. Sagu´s, and M. San Miguel. Dynamics of sweeping through an
e
instability: Passage-time statistics for colored noise. Phys Rev A, 40:6662–6672,
1989.
[112] A. V. Tsukanov. Nanoelectromechanical systems and quantum information
processing. Russian Microelectronics, 40:254–267, 2011.
[113] Michael Urbakh, Joseph Klafter, Delphine Gourdon, and Jacob Israelachvili.
The nonlinear nature of friction. Nature, 430:525–528, 2004.
[114] D. F. Urban and Hermann Grabert. Feedback and rate asymmetry of the
josephson junction noise detector. Phys. Rev. B, 79:113102, 2009.
[115] R. H. Victora. Predicted time dependence of the switching field for magnetic
materials. Phys Rev Lett, 63:457–460, 1989.
[116] R. Vijay, M. H. Devoret, and I. Siddiqi. Invited review article: The josephson
bifurcation amplifier. Review of Scientific Instruments, 80:111101, 2009.
[117] D. Vion, M. Gotz, P. Joyez, D. Esteve, and M. H. Devoret. Thermal activation
above a dissipation barrier: Switching of a small josephson junction. Phys Rev
Lett, 77(16):3435–3438, 1996.
[118] Philip S. Waggoner and Harold G. Craighead. Micro- and nanomechanical
sensors for environmental, chemical, and biological detection. Lab on a Chip,
7:1238–1255, 2007.
[119] W. Wernsdorfer, M. Murugesu, A. J. Tasiopoulos, and G. Christou. Fieldsweep-rate dependence of the coercive field of single-molecule magnets: A classical approach with applications to the quantum regime. Phys. Rev. B, 72:212406,
2005.
[120] W. Wernsdorfer, E. Bonet Orozco, K. Hasselbach, A. Benoit, B. Barbara, N. Demoncy, A. Loiseau, H. Pascard, and D. Mailly. Experimental evidence of the
nel-brown model of magnetization reversal. Phys Rev Lett, 78:1791–1794, 1997.
157

[121] Y. T. Yang, C. Callegari, X. L. Feng, K. L. Ekinci, and M. L. Roukes.
Zeptogram-scale nanomechanical mass sensing. Nano Letters, 6:583586, 2006.
[122] H. Zeghlache and Paul Mandel. Influence of noise on delayed bifurcations. Phys
Rev A, 40:286–294, 1989.

158