MEASUREMENTS USING NOISE-DRIVEN NONLINEAR MICRO-RESONATORS
By
Pavel M. Polunin

A THESIS
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
Mechanical Engineering - Master of Science
2013

ABSTRACT
MEASUREMENTS USING NOISE-DRIVEN NONLINEAR
MICRO-RESONATORS
By
Pavel M. Polunin
In this work we discuss new measurement applications of nonlinear micro-resonators that
are subject to stochastic forcing in combination with primary periodic excitation. First,
we describe how the noise-driven switching of a nonlinear oscillator subjected to harmonic
excitation with bistable response can be used as a new type of sensor for measuring parameter
changes in the system; we refer to this as the balanced dynamic bridge. For small noise we
develop a predictive theory that describes the dynamical behavior of these systems on the
oscillator time scale, as well as on the characteristic time scale of the switching. A general
theory of activated escape in the presence of a Gaussian noise allows one to compute the
switching rates between the two stable states, and the manner in which these rates are related
to system parameters. We discuss the high sensitivity of the bridge with respect to changes
of system parameters and derive expressions describing the precision of the method and the
time required to perform experimental measurements to a given precision. In addition, we
discuss application of the dynamical bridge as a sensor of non-Gaussian noise in the presence
of additional weak non-Gaussian stochastic forcing. We show how non-Gaussian noise can
be detected and estimated by measuring the long-time occupation probability ratio of the
bistable states of the system. We conclude our discussion with two examples: they describe
non-Gaussian noise estimation in a one-dimensional system and in the nonlinear Mathieu
resonator.

ACKNOWLEDGMENTS

I would like to express all my gratitude and respect to my research advisor Professor Steven
W. Shaw. Apart from his professional guidance and supervision, Professor Shaw has been
my excellent teacher and instructor in all issues that I ever had during the time I worked
on this thesis. I am really glad that fortune has given me the chance to work with this
outstanding person. Additionally, I would especially like to thank Professor Mark I. Dykman
for very fruitful discussions and valuable comments on my work process. Thanks also to
committee members who oversaw this work: Professor Brian Feeny and Professor Ranjan
Mukherjee. I owe my inspiration to my former lab-mate Nick Miller for a great number
of useful discussions and his example of successful junior researcher; I am sure it will play
important role in my future life as a PhD student. Furthermore, I thank my fellow labmates for providing a friendly atmosphere in the lab that has made me feel welcome and
comfortable. These include Scott Strachan, Brendan Vidmar, Venkat Ramakrishnan, Smruti
Panigrahi, Abhisek Jain and Xing Xing. I also thank to my friends, Andrey Maslennikov and
Oleksii Karpenko, for their support and making my work more enjoyable. Finally, I would
like to thank the funding agencies that made this work possible. The National Science
Foundation under grant CMMI-0900666 and the Defense Advanced Research Agency under
the DARPA-MTO-DEFYS project.

iii

TABLE OF CONTENTS

List of Figures
Chapter 1

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

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

v

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

Chapter 2 The Balanced Dynamical Bridge . . . . . .
2.1 Theory of activated escape with Gaussian excitation . .
2.2 The Duffing Resonator as a Balanced Dynamic Bridge
2.3 Sensitivity of the Dynamic Bridge . . . . . . . . . . . .
2.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

7
8
13
23
31

Chapter 3 Non-Gaussian Noise Detection . . . .
3.1 Switching in the Presence of Non-Gaussian Noise
3.2 Detection Technique . . . . . . . . . . . . . . . .
3.3 An Example: Shot Noise Measurement in 1D . . .
3.4 Noise in Parametrically Excited Resonators . . . .
3.5 Summary . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

33
34
39
42
47
60

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

62

Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

65

Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

69

Chapter 4

Conclusions

iv

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

LIST OF FIGURES

Figure 1.1

Figure 1.2

Figure 2.1

Figure 2.2

Figure 2.3

Figure 2.4

Shift in natural frequency of a linear resonator, as a tool for estimating parameter changes. (For interpretation of the references to color
in this and all other figures, the reader is referred to the electronic
version of this thesis.) . . . . . . . . . . . . . . . . . . . . . . . . . .

3

Nonlinear resonator frequency response for different amplitudes of
excitation. Dashed lines represent responses that are below the noise
level while responses in solid lines exceed the noise level, and exhibit
nonlinearity. This demonstrates a case where a linear dynamic range
does not exist. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

Solution to Eq. (2.27). The solid line represents the function g with
Ω = 4.7 while the dashed line represents β = 8.5. The fixed points
correspond to the intersection points, which provide the steady state
values of |u|2 . Generally the outer two solutions are dynamically
stable while the central solution is unstable. . . . . . . . . . . . . . .

17

The phase portrait of a Duffing resonator depicted in the u phase
plane Eqs. (2.26a) and (2.26b). The solid lines represent the unstable
manifolds of the saddle, which originate at the saddle and terminate
at the stable fixed points, while the dashed line designates the stable
manifolds of the saddle, representing separatrices between the two
basins of attraction. The parameter values are Ω = 4.7 and β = 8.5.

18

The (Ω, β) parameter space of the Duffing resonator, depicting the
bistable region. The solid branches enclose the region of bistability
and represent the bifurcation curves given in Eq. (2.30). . . . . . . .

19

The heteroclinic solution to Eqs. (2.31a) to (2.31d) projected onto
(ur , ui ) plane. The solid lines represent noisy dynamics and the
dashed lines are the noise-free saddle manifolds. The parameter values are Ω = 4.7 and β = 8.5. . . . . . . . . . . . . . . . . . . . . . .

20

v

Figure 2.5

Figure 2.6

The heteroclinic solution to Eqs. (2.31a) to (2.31d) projected onto
(pr , pi ) plane. The parameter values are Ω = 4.7 and β = 8.5. . . . .

21

The (Ω, β) parameter space for the Duffing resonator. The solid lines
represent the boundaries of the bistable region while the dashed line
depicts the balance curve. . . . . . . . . . . . . . . . . . . . . . . . .

22

Figure 2.7

The occupation probabilities as a functions of the action change, ∆S. 27

Figure 2.8

The lower bound on the number of measurements required to obtain
an accurate estimate of the parameter change using the balanced
dynamical bridge. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

28

The lower bound on the required total measurement time as a function of the shift in the occupation probability. The dimension of the
vertical axis, Tm rbp , represents the number of the switching events
that would happen if the system stays at the balance point for the
time Tm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

30

Figure 2.9

Figure 3.1

Schematic one-dimensional representation of the potential for a bistable
nonlinear system, Eq. (3.1). . . . . . . . . . . . . . . . . . . . . . . . 35

Figure 3.2

The change in the probability ratio due to exposure the balanced
dynamic bridge to a non-Gaussian perturbation. . . . . . . . . . . .

40

The double-well potential Eq. (3.24); q1 and q2 denote positions of the
stable fixed points, while qS refers to the saddle-point. Poisson pulses
with positive area g change the switching rates, such that r1 > r2 ,
therefore shifting the probability from q1 towards q2 . . . . . . . . . .

43

Figure 3.3

Figure 3.4

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

Figure 3.5

Dependence of the amplitude of steady-state oscillations on the normalized detuning parameter Ω. Parameter α is taken to be equal
√
17. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

50

Phase portrait of the nonlinear Mathieu resonator. . . . . . . . . . .

52

Figure 3.6

vi

Figure 3.7

Figure 3.8

Figure 3.9

Figure A.1

Shaded area represents the bistable region of the nonlinear Mathieu
resonator. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

53

Solutions to Eq. (3.43) are represented by the intersection points of
the solid and the dashed lines. The left-hand side of Eq. (3.44) as a
function of |uo |2 is drawn by the solid line, while β 2 is drawn by the
√
dashed line. Parameters values: α = 17, Ω = 0, ψ = 30◦ , β = 1. . .

54

The numerical solution to the four-dimensional system of equations
corresponding to the nonlinear Mathieu resonator depicted as a projection onto: (a) the (ur , ui ) plane, and (b) the (pr , pi ) plane. The
noise-free dynamics is exactly the same as shown in Fig. 3.6(b), and
is shown in black dotted lines. Red and blue lines are the projections
of the noisy switching trajectories from different stable states to the
saddle. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

58

Schematic representation of a quantum-mechanical system that has
more than one possible trajectory between points A and B. Different
paths have different probabilities of realization, and these paths can
significantly deviate one from another. . . . . . . . . . . . . . . . . .

66

vii

Chapter 1
Introduction
It is well known that noise can have a significant effect on system response in many applications, and this is especially the case when the characteristic size of the system of interest
is very small. Here we use the term noise to designate various types of stochastic processes
that can occur in systems. These include small fluctuations of charge on capacitors in electric
circuits, tiny mass or stiffness fluctuations in mechanical systems, thermal vibrations, and
so on. In most systems these noise processes have a negative effect on system response, and
consequently, engineers search for ways to eliminate them or reduce their effects. However,
this task becomes more challenging for systems at micro- and nano- length scales. In fact,
for many systems, these random processes are an essential part of system dynamics and must
be taken into consideration in system analysis and design.
In this work we discuss vibrating micro/nano-electro-mechanical systems (M/NEMS)
that operate in a nonlinear regime and are subject to stochastic effects in addition to resonant excitation. These systems are relevant to a number of existing and potential applications, including the counting statistics of charge transport systems, and employment as
1

elements in small-scale frequency generators. Of particular interest are systems that exhibit
nonlinear responses, for which we give a detailed analysis of their dynamics and show that
noise appearing in these systems can be constructively employed in particular applications,
specifically: (i) for the measurement of parameter changes in the system, which can be linked
to the system environment for purposes of sensing, and (ii) for the detection of specific types
of non-Gaussian random processes. In this chapter we give a brief review of M/NEMS, discuss how noise and random processes are related to M/NEMS, and conclude the chapter by
providing more specific motivation for the present work.
Having more than 50 years of development and enhancement [13], M/NEMS have found
their place in many areas of modern human life. These include applications in microelectronics, medicine, aerospace and automotive engineering, and so on. One classification of
M/NEMS divides most existing systems into two major groups: sensors like inertial, pressure
and flow sensors, to name a few, and devices used for signal processing, for example filters,
oscillators, switches, and relays. In this study we focus our attention on micro-mechanical
resonators, which have applications in both groups.
It is well known that linear systems are much easier to understand and analyze than
their nonlinear counterparts, cf.[16]. So, there is no surprise that researchers initially worked
with simple linear models of resonators for different types of measurement and processing
purposes. Essentially, the concept of M/NEMS measurements using linear resonance is the
following: suppose we have a micro-resonator that is coupled to a environmental parameter of
interest in such a way that changes of this parameter will affect the resonant frequency of the
system, see Fig. 1.1. Then, by calibrating and measuring these shifts in the natural frequency
of the resonator, we can determine the associated change of the parameter of interest. For
2

Response

example, many mass-sensing devices and humidity sensors work on this principle [12, 26].

Ω0

Ω
Figure 1.1: Shift in natural frequency of a linear resonator, as a tool for estimating parameter
changes. (For interpretation of the references to color in this and all other figures, the reader
is referred to the electronic version of this thesis.)

However, new fabrication techniques are making it possible to fabricate M/NEMS of
smaller and smaller dimensions and, as mentioned previously, noise effects can come into
play, and are often an essential component of the system response. Noise leads to corruption
of the response signal, and, as a consequence, reduces the bandwidth of the system. In
some cases, the ambient noise signal of a device can even overcome the signal containing
measurement information, and the device loses its usefulness. The way around this obstacle
is to increase the amplitude of system excitation in order to reach a desired signal-to-noise
ratio. However, increasing the system input signal inevitably leads to a transition to a
nonlinear response regime, where linear model approximation fail, see Fig. 1.2. In this
situation one must consider the effects of nonlinearities and noise in the response, that is, a
linear dynamic range does not exist.
Nonlinear behavior of M/NEMS resonators involves many interesting aspects, including
so-called bistability, that is, the property of the system to have multiple stable operating
3

Response

Noise level

Detuning
Figure 1.2: Nonlinear resonator frequency response for different amplitudes of excitation.
Dashed lines represent responses that are below the noise level while responses in solid lines
exceed the noise level, and exhibit nonlinearity. This demonstrates a case where a linear
dynamic range does not exist.

regimes [21]. It turns out that noise in these cases can cause switching of the micro-resonator
response between these stable operating regimes, cf. [2, 3, 8, 9, 20]. The rates of these
transitions depend exponentially on the height of the potential barrier between stable states
and on the noise strength, according to the Arrhenius law [1, 25],

r ∝ exp −

R
,
D

(1.1)

where R is the height of potential barrier and D is intensity of the noise. In the weak
noise limit, i.e., when D is small, even small changes in the activation energy can produce
significant changes in switching rates, which are reflected in system behavior. This sensitivity
is a primary motivation for the present investigation. Specifically, we are interested in the
possibility of using switching in noisy bistable systems for purposes of sensing. The significant
part of the background and motivation to this work lies in Ph.D. thesis of N.J.Miller, [20].
In this work we consider two topics that are closely connected and represent novel alternatives to contemporary measurement techniques. First is the application of noise-driven
4

bistable systems for parametric sensing, which we refer to as the balanced dynamical bridge.
In Chapter 2 we consider the bridge application for the measurement of changes in system
parameters due to environmental changes. The theory of activated escape forms the basis of
the analysis of the bridge, and this is described first. We then consider a Duffing resonator
as a model for our subsequent discussion and analyze its switching dynamics in detail. We
derive and analyze a model describing transitions of the system between its stable states and
investigate it numerically. Finally, we turn to the switching dynamics on a long time scale
and derive expressions that describe the measurement precision of the bridge, as well as the
time needed for measurements.

In Chapter 3 we examine the second application of M/NEMS, specifically, use of the
balanced dynamical bridge as a detector of non-Gaussian noise. The presence of an additional
non-Gaussian modulation to a bistable system switching under Gaussian noise requires a
substantial change in the theory of the activated escape, and we address this modification
at the beginning of Chapter 3. We find that this non-Gaussian forcing is projected onto
the so-called susceptibility of the system on switching trajectories in the phase space. This
projection affects the switching rates of the resonator in each direction differently, thus
providing a shift of the population ratio. Using these results we derive expressions that
allow us to determine the statistical properties of the non-Gaussian noise from changes in
transition rates. Two examples of non-Gaussian detection are presented. The first example
considers a simple one-dimensional system where analytical predictions can be completed
in closed form, and the second example deals with a nonlinear Mathieu resonator, which
is substantially more complex and requires more numerical calculations. The results from
these examples allow us to draw some general conclusions about the effectiveness of using
5

switching in noisy bistable M/NEMS for measurement and detection purposes.

6

Chapter 2

The Balanced Dynamical Bridge

In this chapter we discuss a new measurement paradigm which we refer to as balanced
dynamic bridge. This sensing application can be realized in bistable systems that are driven
by a white Gaussian noise such that the system randomly switches between the two stable
states. Measurement is achieved by relating the ratio of the occupation probabilities of the
two states to system properties that are affected by the system environment. It is shown
that this measurement approach is most sensitive when one operates near a point where the
occupation probabilities are equal. It is also shown that this method has high resolution,
but longer measurement time, when compared to more standard measurement techniques.
In Section 2.1 we describe a general model and theory for switching dynamics of bistable
systems in the presence of white Gaussian noise. In Section 2.2 we consider the specific
example of a resonantly driven Duffing oscillator. In Section 2.3 we describe a general theory
for measurement sensitivity and measurement time for the balanced bridge and show results
derived from the Duffing example. The Chapter is closed with a Summary in Section 2.4.
7

2.1

Theory of activated escape with Gaussian excitation

In this section we review a general theory of noise-activated escape in the presence of a
Gaussian noise. Originated by Kramer in his famous work [17], this theory has been extended
significantly since that time and has found its place in the many applications.
We start with a general system those dynamics can be expressed in the state variable
form as follows,
ˆ
˙
u = K(u) + f (t),

(2.1)

where u ∈ Rn is the state vector of the system of interest, K is the deterministic nonlinear
ˆ
vector field, and f ∈ Rn is a stationary Gaussian noise vector with characteristic strength D.
It is assumed that form of K gives rise for the system to experience at least two stable steadyˆ
states, and that all components of the noise vector, f , are delta-correlated and independent
from each other. For further convenience we transform the system of interest to a discrete
form,
ui+1 − ui
ˆ
= K(ui ) + f i ,
∆t

(2.2)

using the following time-slicing

∆t = ti+1 − ti ,

i = 1, 2, . . . , N.

(2.3)

ˆ
Since f is the vector of the independent white Gaussian noise processes, its auto-correlation
function is simply
D
ˆ ˆT
δ I,
f if j =
∆t ij
8

(2.4)

where δij denotes the Kronecker delta function and I is the identity matrix.
Using these introductory concepts we can formulate the main results of this section,
namely we can compute the probability for the system to make a transition from a region
near a stable fixed point in the configuration space to another region, near another fixed
point, typically a saddle point, during a given time interval. Since the noise is considered to
be small, such transitions rare, and in order to be realized they require quite specific time
evolutions of the stochastic driving. In the case when system dynamics is represented in the
ˆ
ˆ
form of Eq. (2.2), this time evolution of f is described in terms of the {f i }. Because of the
ˆ
ˆ
δ-form of the autocorrelation function for f , we have N independent noise vectors f i . Each
of these vectors, in turn, has Gaussian probability distribution given by,

ˆ
p(f i ) =

where σi =

1

exp −

2
2πσi

ˆ
|f i |2
,
2
2σi

(2.5)

ˆ
D/∆t is constant for all elements f i .

During the time evolution of the system, N noise bursts occur in a series and we are
interested in particular sequence of these bursts which will force the system to move between
points of interest in the configurations space. Following Jacobs [15], the joint probability
distribution of these independent processes is the product of individual probability densities
for each noise vector
N

ˆ
ˆ ˆ
ˆ
p({f }) = p(f 1 , f 2 , . . . , f N ) =

ˆ
p(f i ).

(2.6)

i=1

Applying Eq. (2.5) to the expression for the joint probability distribution, we find

1
ˆ
p({f }) = CN exp −
2D
9

N

ˆ
|f i |2 ∆t ,
i=1

(2.7)

where CN is a coefficient that absorbs all exponential prefactors.

Now it is possible to calculate the probability that the system will come to point u2 at
time t2 being near point u1 at time t1 . This can be done by using the standard relationship
between the cumulative probability function and its density distribution, yielding [15]

P (1 → 2) = P (u2 , t2 ; u1 , t1 ) =

...

ˆ
ˆ
ˆ ˆ
ˆ
p({f })δ(u[u1 , {f }) − u2 ]df 1 df 2 . . . df N , (2.8)

where we have introduced δ-function term reflecting the fact that the sequence {f } must
bring the system from point u1 to point u2 . By changing the variables of integration from
ˆ
df i to dui , and using Eqs. (2.2) and (2.7), we obtain the following expression for the total
probability, P (1 → 2),

P (1 → 2) = CN

u2
u1

1
det(J) exp −
2D

N
i=1

2
ui+1 − ui
− K(ui ) ∆t δ(u − u2 )Du, (2.9)
∆t

ˆ
where J is the Jacobian of the [f → u] transformation and Du designates the path integral
between points u1 and u2 . Appendix A describes more details of this development. Since
the Jacobian J describes the transformation between different system parameters, it remains
fixed for a given system and, as a result, it is independent of the path the system follows
between given points in the phase space. So, we can also include its determinant in the
integral prefactor, and by taking the limit as ∆t → 0, we come to the final form for the
transition probability expressed in terms of the system dynamics

P (1 → 2) = C

u2
u1

exp −

t2
1
˙
|u − K(u)|2 dt Du.
2D t
1

10

(2.10)

It is useful to note here that there are an infinite number of trajectories that the system can
follow as it undergoes the transition between two given points in the state space. However,
the probability of the realization of a particular trajectory depends strongly on the path, and
so there exists a probability distribution for the possible trajectories. Then, the usual way
to calculate the total probability for for the system of interest to switch from one state to
another is to integrate this probability density distribution over the whole range of possible
trajectories. That is exactly the role of the path integral in Eq. (2.10); it simply integrates
the probability distribution over all possible time evolutions of the system between two given
points.

The probability distribution discussed above possesses one remarkable property which
significantly facilitates the calculation of the total transition probability. According to previous work of Dykman and collaborators, cf.[9, 2, 3], this probability distribution is sharply
peaked about so-called most probable path. Using this property, we can approximate the
total transition probability by the expression in the integrand of Eq. (2.10) evaluated along
the most probable path, as follows,

P (1 → 2)

C exp −

t2
1
˙
|u − K(u)|2 dt
.
2D t
MPP
1

(2.11)

Hereafter we assume that we only work with the most probable path. According to the
path integral formulation, the transition probability P (1 → 2) depends exponentially on the
action, calculated along the most probable path, namely,

P (1 → 2) ∝ exp
11

S
− 1→2
D

(2.12)

Then, by comparison Eqs. (2.11) and (2.12) one can see that the action S1→2 can be expressed
in the following form
1 t2
˙
|u − K(u)|2 dt.
2 t
1

S1→2 =

(2.13)

On the other hand, it is well known from the classical mechanics [11] that mechanical action
is related with the Lagrangian of the system as

S=

t2

˙
L(q, q, t)dt,

(2.14)

t1

So, the Lagrangian for the general system, Eq. (2.1), reads

1
˙
˙
L(u, u) = |u − K(u)|2
2

(2.15)

along with the generalized momenta defined as

p=

∂L
˙
= u − K(u),
˙
∂u

(2.16)

which is exactly the vector of white Gaussian noise processes defined in Eq. (2.1). Now, we
only have to specify the equations of motion for the generalized momenta in order to get a
complete description of the system dynamics. It can be done in the following way,

p=−
˙

∂H
∂L
∂K
=
=−
p,
∂u
∂u
∂u

(2.17)

where H is the associated Hamiltonian of the system. Eqs. (2.1) and (2.17) are the full set of
equations of motion that completely describe the system dynamics in the extended R2n state

12

space, which contains the noise-free dynamics and the realizations of noise. The solution to
these equations describing the heteroclinic trajectory between two stable states represents
the most probable path for the system in this space.

2.2

The Duffing Resonator as a Balanced Dynamic Bridge

In this section we discuss the harmonically excited Duffing resonator as an example of a
nonlinear system that exhibits bistability under certain conditions, and its operation as a
dynamic bridge when subjected to additive noise. We start from the model of the Duffing
resonator, subject to both periodic and stochastic excitation,

2
q + 2Γq + ω0 q + γq 3 = h cos ωt +
¨
˙

√
Df (t),

(2.18)

where q is the resonator coordinate, ω0 is its natural frequency, Γ is its damping coefficient,
γ is the nonlinearity parameter, h is the amplitude of the harmonic forcing, ω is the driving
√
frequency, and Df (t) is white Gaussian noise of strength D. It is assumed that the
Gaussian noise has zero mean, f (t) = 0, and is delta-correlated, f (t)f (t ) = 2Dδ(t − t ).
We also assume that |ω − ω0 |
lightly damped, Γ

ω, ω0 , that is, the system is driven near resonance and is

ω0 . In order to simplify the analysis, it is convenient to switch from

coordinate q to a slowly-varying complex amplitude u in a rotating frame, achieved by the
well known van der Pol transformation,

q(t) =

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

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

(2.19a)
(2.19b)

Using this transformation, we can obtain the following expressions for the first and second
derivatives of the primary variable q, namely

q(t) =
˙

2ωΓ
iω(ueiωt − u∗ e−iωt ),
3|γ|

q (t) = −ω 2
¨

(2.20a)

2ωΓ iωt
(ue + u∗ e−iωt ) + 2iω
3|γ|

2ωΓ iωt
ue .
˙
3|γ|

(2.20b)

Substitution of q and its derivatives as functions of u, u∗ , u, and u∗ into Eq. (2.18), elimination
˙
˙
of fast oscillating terms of the deterministic part of the equation, and rescaling of time and
parameters yields the following complex first order differential equation for the slow-varying
amplitude u,
u = −(1 + iΩ)u + i sgn(γ)|u|2 u − i
˙

ˆ
β − if (τ ),

(2.21)

where (˙) is understood to be the derivative with respect to the rescaled time τ , defined below.
Additionally, in order to make the equation fully non-dimensional, we have introduced the
following quantities,

t = Γτ,

(2.22a)

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

Ω=

(2.22b)
(2.22c)
(2.22d)

Since the spectrum of the white noise f (t) has constant magnitude over infinite frequency
range, and because we are interested in switching dynamics of the system, we have not
averaged the noise term, and thus we retain the complex exponential factor in Eq. (2.22d).
14

Instead of trying to extract its mean value and neglect fast oscillating components, we use
ˆ
ˆ
Euler’s formula for the complex exponent to separate real, fr (τ ), and imaginary, fi (τ ), parts
of the complex noise function, resulting in

ˆ
f (τ ) =

3|γ| ˆ
ˆ
(fr (τ ) − ifi (τ )),
8ω 3 Γ3

(2.23)

ˆ
ˆ
ˆ
ˆ
where fr (τ ) = f (τ ) cos(ωt) and fi (τ ) = f (τ ) sin(ωt). These projections of the original white
Gaussian noise onto two orthogonal quadratures makes them approximately independent
white noise processes with following properties [9],

ˆ
ˆ
fr (τ ) = fi (τ ) = 0,

(2.24a)

ˆ
ˆ
ˆ
ˆ
ˆ
fr (τ )fr (τ ) = fi (τ )fi (τ ) = Dδ(τ − τ ),

(2.24b)

ˆ
ˆ
fr (τ )fi (τ ) = 0,

(2.24c)

ˆ
where D is the nondimensional noise intensity of each quadrature, defined as,

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

(2.25)

Equation (2.21) describes the time evolution of the complex slowly-varying amplitude u.
In the absence of the stochastic excitation, it is possible to find stationary solutions of the
resonator by setting u = 0, cf.[16]. In order to determine these stationary amplitudes, we
˙
15

first separate real, ur , and imaginary, ui , parts of u and get the following dynamical system,

ur = −ur + Ωui − sgn(γ)(u2 + u2 )ui ,
˙
r
i
ui = −ui − Ωur + sgn(γ)(u2 + u2 )ur −
˙
r
i

(2.26a)
β.

(2.26b)

Now, by setting ur = ui = 0 and taking a convenient linear combination of the resulting
˙
˙
equations we obtain,

g(|u|2 , Ω) = |u|6 − 2 sgn(γ)Ω|u|4 + (1 + Ω2 )|u|2 = β,

(2.27)

where |u| is the magnitude of u and function g is introduced for convenience. Equation (2.27)
is a cubic algebraic equation in terms of |u|2 which, in general, provides three stationary
solutions to the system of Eqs. (2.26a) and (2.26b). Since we now work with |u|2 , we are
only interested in the real roots of Eq. (2.27). The number of real roots of any cubic equation
depends on the values of the coefficients of the polynomial, which in this case depend on the
values of Ω and β. Fig. 2.1 shows a sample of the cubic function g versus |u|2 and a sample
level line of β.

There are two qualitatively different types of system dynamics. Depending on the values
of Ω and β, the system of interest generically has either one or three fixed points in the u
phase plane. In the latter case, two fixed points are stable while the third one is an unstable
saddle point[16]. A system that exhibits this type of dynamics is called bistable. The typical
phase portrait of a bistable Duffing resonator is depicted in Fig. 2.2, where the bistability
of the system and the attendant domains of attraction are clearly shown. The basins of
attraction of the two stable fixed points are delineated by separatrices, that is, the stable
16

20

g u2 , Β

15

10

5

0
0

1

2

3

4

5

6

7

u2
Figure 2.1: Solution to Eq. (2.27). The solid line represents the function g with Ω = 4.7
while the dashed line represents β = 8.5. The fixed points correspond to the intersection
points, which provide the steady state values of |u|2 . Generally the outer two solutions are
dynamically stable while the central solution is unstable.

manifolds of the saddle point.
Since the concept of the balanced dynamical bridge is based on the switching dynamics
between the system’s stable states, it is reasonable to determine the region in the (Ω, β)
parameter space where the Duffing resonator has two stable solutions. In order to find this
region of bistability, we can look at Fig. 2.1 and note that for a fixed value of Ω, which
determines the form of the cubic function g, there are two values of β where a stable and
an unstable solution merge. These are saddle-node bifurcations [14], and it is clear that the
derivative of the cubic function g with respect to |u|2 is equal to zero at these bifurcation
points, i.e.,
3|u|4 − 4 sgn(γ)Ω|u|2 + (1 + Ω2 ) = 0.

(2.28)

Equation (2.28) is quadratic in terms of |u|2 , and we can solve it obtaining the following
17

1
0.5
0
ui

-0.5
-1
-1.5
-2
-2.5
-2

-1

0
ur

1

2

Figure 2.2: The phase portrait of a Duffing resonator depicted in the u phase plane
Eqs. (2.26a) and (2.26b). The solid lines represent the unstable manifolds of the saddle,
which originate at the saddle and terminate at the stable fixed points, while the dashed
line designates the stable manifolds of the saddle, representing separatrices between the two
basins of attraction. The parameter values are Ω = 4.7 and β = 8.5.
values of |u|2 at the bifurcation points as a function of Ω,

1
|u|2 = {2 sgn(γ)Ω ±
B
3

Ω2 − 3}.

(2.29)

In order to obtain the bifurcation curves in (Ω, β) parameter space, one substitutes the values
of |u|2 into Eq. (2.27), which yields the saddle-node parameter condition for β,
B

βB =

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

Ω2 > 3.

(2.30)

The bifurcation curves obtained from Eq. (2.30) are depicted in a convenient form in Fig. 2.3.
18

0.4

Β

3

0.3
0.2
0.1
0.0
0.0

0.1

0.2

0.3

0.4

2

Figure 2.3: The (Ω, β) parameter space of the Duffing resonator, depicting the bistable
region. The solid branches enclose the region of bistability and represent the bifurcation
curves given in Eq. (2.30).
By defining the bistable region in the (Ω, β) parameter space, we have completely described the deterministic part of the system dynamics. In order to understand the switching
dynamics of the resonator, we use the results of Section 2.1, namely, Eqs. (2.1) and (2.17)
capturing the noise dynamics. Thus, the equations of motion governing the optimal switching
trajectories are four-dimensional and have the following form,

ur = −ur + Ωui − sgn(γ)(u2 + u2 )ui + pr ,
˙
r
i

(2.31a)

ui = −ui − Ωur + sgn(γ)(u2 + u2 )ur −
˙
r
i

(2.31b)

β + pi ,

pr = [1 + 2 sgn(γ)ur ui ]pr + [(u2 + 3u2 ) sgn(γ) − Ω]pi ,
˙
r
i

(2.31c)

pi = [Ω − sgn(γ)(u2 + 3u2 )]pr + [1 − 2 sgn(γ)ur ui ]pi .
˙
r
i

(2.31d)

Solutions to this system of equations represent the most probable path the system follows
in the extended state space, where the u’s are the system states and the p’s represent noise
19

1
0.5
0

ui

-0.5
-1
-1.5
-2
-2.5
-3
-2

-1.5

-1

-0.5

0

0.5

1

1.5

2

ur
Figure 2.4: The heteroclinic solution to Eqs. (2.31a) to (2.31d) projected onto (ur , ui ) plane.
The solid lines represent noisy dynamics and the dashed lines are the noise-free saddle
manifolds. The parameter values are Ω = 4.7 and β = 8.5.

quadratures. It is important to note that the noise-free dynamics are invariant, that is, if one
starts with zero initial conditions for the p’s, they remain zero for all time, and the system
simply follows the deterministic equations of motion. Thus, fixed points of the deterministic
system are also fixed points of the full system. However, the stability of these points can
be different; in fact, fixed points that are stable in the deterministic system become saddle
points in the full system, and heteroclinic trajectories that connect deterministic fixed points
to one another, via a noisy trajectory, are an essential feature of the switching dynamics.
The two optimal switching trajectories that go between the two deterministically stable
foci, one each way, each consist of two parts. The first part is induced by the noise and is
20

1.5

1

pi

0.5

0

-0.5

-1

-1.5
-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

p

r

Figure 2.5: The heteroclinic solution to Eqs. (2.31a) to (2.31d) projected onto (pr , pi ) plane.
The parameter values are Ω = 4.7 and β = 8.5.
represented by the solution of Eqs. (2.31a) to (2.31d) that goes from a deterministically stable
focus to the deterministic saddle; such a trajectory starts and ends with the p’s both zero,
but must extend into the noise states. The second half of a switching trajectory is essentially
noise-free and goes along the corresponding unstable manifold of the saddle, asymptotically
approaching the region of other deterministically stable focus. Of course, once the solution
gets to the saddle point, there is a 50/50 chance of it going towards either deterministically
stable focus. Thus, the two complete switching trajectories are each composed of the joining
of two heteroclinic solutions, the first of which involves all four states, while the second is
essentially deterministic. The challenge to determining these trajectories is that of finding
saddle-to-saddle heteroclinic solutions of a four-dimensional system, since the deterministic
part is simple to find. By using a shooting method, these heteroclinic solutions to the four
21

dimensional equations can be obtained, a sample of which is depicted in Figs. 2.4 and 2.5,
depicted by two dimensional projections.
The results shown in Figs. 2.4 and 2.5 were obtained for a special case when the action
calculated along the two switching trajectories has the same value in both directions, which
requires a constraint on parameters Ω and β. In this case we call the system a balanced
dynamic bridge, which is especially useful for the measurement purposes, as described below.
Additional numerical calculations show that there exists a so-called balance curve, that is,
a locus of points in the (Ω, β) parameter space, where the system is in the balanced bridge
configuration. This curve is depicted as the dashed line in Fig. 2.6.

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

0.25

0.30

2

Figure 2.6: The (Ω, β) parameter space for the Duffing resonator. The solid lines represent
the boundaries of the bistable region while the dashed line depicts the balance curve.

The results presented above set the stage for using the Duffing resonator with additive
Gaussian noise as a balanced bridge. In the next section we discuss some general features of
the bridge as a measurement tool.
22

2.3

Sensitivity of the Dynamic Bridge

In this section we discuss the quantitative aspects of the balanced dynamical bridge as a
measurement tool. The idea is to use the occupation probability ratio of the two stable
states, and its dependence on system and input parameters, as the source of information.
Since this ratio depends on the switching rate, we can use the results derived above to develop
a theoretical basis for such a sensor. Specifically, we investigate the sensitivity properties of
the bridge and the time required to perform measurements to a given level of precision.
Section 2.2 describes the dynamics of the resonator on the time scale of the order of one
switching event. In order to analyze the performance of this measurement techniques, we
need to consider a longer time scale since the probability ratio is a statistical quantity of
the system accumulated over many switching events. Considering the transitions that the
system makes between the two stable states on a long time scale, it is convenient to use a
master equation governing the dynamics of a two-level system, given by [23],

˙
P1 = −r1 P1 + r2 P2 ,

(2.32a)

˙
P2 = r1 P1 − r2 P2 ,

(2.32b)

where Pi is the probability to find the system in ith state, also known as the occupation probability, and ri is the switching rate out of the ith state, for i = 1, 2. Since there are just two
states that system can occupy, the occupation probabilities must satisfy the normalization
condition, namely,
P1 + P 2 = 1 .

(2.33)

Hereafter, we focus our attention only on P1 , assuming that P2 and its dynamics can be found
23

using this condition and the master equation. The time-dependent solution to Eqs. (2.32a)
and (2.32b) can be obtained by using standard techniques [7, 16] and reads,

P1 (t) =

r2
1 − e−(r1 +r2 )t
r1 + r2

+ P1 (0)e−(r1 +r2 )t .

(2.34)

From this expression it is clear that the characteristic time constant of the bridge is quite
large, namely tr ∼ (r1 + r2 )−1 .
The concept of the dynamic bridge as a sensor is to track changes in occupation probabilities that vary as a parameter of interest changes. In order to measure the occupation
probability, P1 , we have to have some quantity x in the system that takes on the value
X1 when the resonator is in the first state and X2 when the resonator is in the second
state. Suppose we take N measurements of the system states. The time between successive
measurements has to be long enough for the system to relax; particularly, it has to be on
the order of tr , in order to statistically capture the switching dynamics. Taking this into
account, we can conclude one requires a time on the order of N tr , with N sufficiently large
(as described below) in order to perform the measurement with sufficient accuracy.
Having x as the test quantity, we can write the following expression for the estimator of
P1 [19],
ˆ
P1 =

1
N (X1 − X2 )

N

(xi − X2 ).

(2.35)

i=1

ˆ
The mean and the variance of the estimator P1 are then

ˆ
¯
P1 = P 1 ,

(2.36a)

¯
¯2
ˆ1 − P1 )2 = P1 − P1 ,
¯
(P
N

(2.36b)

24

where
¯
P1 =

r2
r1 + r2

(2.37)

is the steady-state value of P1 . Equation (2.36b) shows that the uncertainty of the probability
measurement is proportional to N −1/2 , which proves that the accuracy of the measurement
technique is higher when more measurement experiments are performed.
Now, suppose that some parameter λ of the resonator changes by an amount ∆λ. Ex¯
¯
pansion of P1 in a Taylor series yields that the system sensitivity is described by dP1 /dλ.
In the following calculation we assume that the prefactor of the switching rate depends on λ
weakly, and we neglect this dependence in the weak noise limit. As a result, the corresponding change in the occupation probability is caused by the change of the activation energy in
both directions and is given by,
¯
r1 r2
∂(S1 − S2 ) S1 − S2 ∂D
dP1
=
−
,
2
dλ
∂λ
D
∂λ
D(r1 + r2 )

(2.38)

where Si is the activation barrier that system has to overcome in order to escape from ith
basin of attraction, i = 1, 2. This expression provides the link between the general analysis
of the switching dynamics and the measurement sensitivity.
Equation (2.38) shows that the operating point where the system exhibits the highest
sensitivity with respect to the parameter change is located on the balance curve, as shown
in Fig. 2.6 for the Duffing equation. In this case r1

r2 and S1

S2 , and the sensitivity

term simplifies to,
¯
d P1
dλ

1 ∂(S1 − S2 )
,
D
∂λ

(2.39)

which shows that the dynamic bridge can be extremely sensitive as the noise intensity D
25

decreases.

After a sufficiently long time after a system parameter changes, the switching process
relaxes and new occupation probabilities are established in the system. The corresponding
shift in the occupation probability P1 is proportional to the parameter change, namely,

¯
∆P1 =

¯
dP1
∆λ.
dλ

(2.40)

In order to detect this change with the bridge, it is necessary that the uncertainty in the
measurement process must be of the order of the probability shift that we want to measure,
implying that,
ˆ
¯
¯
(P1 − P1 )2 ∼ ∆P1 .

(2.41)

This condition allows one to determine the number of measurements, N , required to detect
the shift in the parameter value, assuming that the dynamic bridge is operating on the
¯ =
balance curve, so that P1 ∼ 1/2. When the parameter λ changes, it causes a corresponding
change in the action, ∆S, and using the fact that the switching rate depends exponentially
on the path action
ri ∝ exp −

Si
2D

(2.42)

and Eq. (2.37) one finds,
¯
¯
P1 + ∆ P1 =

∆S
e 2D
∆S
− ∆S
e 2D + e 2D

.

(2.43)

¯
In order to obtain an expression for ∆P1 , we subtract the value of 1/2 of equilibrium prob26

¯
ability P1 from the RHS of Eq. (2.43) and simplify the result to

¯
∆P1 =

1
tanh
2

∆S
.
2D

(2.44)

Expressing ∆S from Eq. (2.43), we have

1

P2

0

P1

SD

Figure 2.7: The occupation probabilities as a functions of the action change, ∆S.

∆S = D ln

¯
1 + 2∆P1
¯ .
1 − 2∆P1

(2.45)

ˆ
However, the estimated value of the action change, ∆S, is slightly different from the ∆S
shown in Eq. (2.45). This deviation stems from the uncertainty in the occupation probability
that arises from the finite number of measurements; see Eq. (2.36b). Up to the first order
expansion, the estimation of the change in action reads,

ˆ
∆S = D ln

¯
1 + 2∆P1
¯
1 − 2∆P1

−

4D
ˆ
¯
(∆P1 − ∆P1 ).
¯2 − 1
4∆P1

(2.46)

With these results in hand, we can formulate a condition for a target level of precision of the
ˆ
measurement process. The measurement is deemed sufficiently accurate if ∆S
27

∆S, which

is only possible if the second term of the expansion in Eq. (2.46) is much smaller than the
first one, or
ˆ
¯
(∆P1 − ∆P1 )2

1
¯2
− ∆P1 ln
4

¯
1 + 2∆P1
¯ .
1 − 2∆P1

(2.47)

Additionally, following from Eq. (2.36b), the uncertainty in the occupation probability is
proportional to (4N )−1/2 at the balance operating point, which gives one the opportunity
to obtain the expression for the required number of measurements, given by,

N

2

1
¯2
− ∆P1 ln
4

¯
1 + 2∆P1
¯
1 − 2∆P1

The RHS of Eq. (2.48), with an equality in place of

−2

.

(2.48)

to indicate the boundary, is depicted

100
80

N

60
40
20
0

0.4

0.2

0.0
∆P1

0.2

0.4

Figure 2.8: The lower bound on the number of measurements required to obtain an accurate
estimate of the parameter change using the balanced dynamical bridge.

on Fig. 2.8. These curves, which are symmetric about zero, exhibit some interesting feature
that demonstrate the special utility of the bridge, as discussed further below.
28

The attendant total measurement time, Tm , is proportional to the product of the number
of measurements N and the characteristic relaxation time of the resonator, tr , which can be
expressed in the following way,

tr ∼ (r1 + r2

)−1

=

−1
rbp

−1
¯ 2 −1
∆S
− ∆S
−1 1 + 4∆P1
2D + e 2D
e
= rbp
,
¯2
1 − 4∆P1

(2.49)

where rbp is the switching rate at the balance operating point. Then the expression for the
lower bound on the required total measurement time reads as,

1
Tm ∼
rbp

1 3 ¯2
¯4
− ∆P1 − ∆P1 ln2
4 4

¯
1 + 2∆P1
¯
1 − 2∆P1

−1

,

(2.50)

and this result is depicted in Fig. 2.9, which has the same general features as the required
number of measurements N , which are discussed next.

To sum up, we have described and proved that in general the measurement time is
the longest time scale. In particular, it needs to be much longer than the period of the
resonator oscillations, and much longer than its mean switching time. However, as Fig. 2.9
clearly shows, when the parameter λ changes only very slightly, it causes an infinitesimally
small change in the occupation probabilities, due to the linear dependence (see Eq. (2.40))
and, consequently, the measurement process will take a very long time to detect this tiny
deviation, because the variance of the estimated probability has to be less or equal to the
¯
infinitesimal change in P1 . As Eq. (2.36b) dictates, this variance is inversely proportional to
√
N , and it only can converge to zero when N → ∞. This results in the unboundedness of
¯
¯
N and Tm as δ P1 → 0. The situation is quite similar when ∆P1 approaches the value of 1/2,
where clearly the occupation probabilities are nearly constant and therefore insensitive to
29

80

Tmrbp

60

40

20

0
0.4

0.2

0.0
∆P1

0.2

0.4

Figure 2.9: The lower bound on the required total measurement time as a function of the
shift in the occupation probability. The dimension of the vertical axis, Tm rbp , represents the
number of the switching events that would happen if the system stays at the balance point
for the time Tm .

parameter changes. Equation (2.36b) says that the occupation probability variance becomes
¯
very small as P1 → 1. In this case the probability goes to its boundary value, and we
call this process the probability saturation. So, the probability distribution for the estimated
occupation probability will collapse to a sharp peak in the close vicinity of this boundary, and
¯
¯
in the limit P1 → 1, it approaches a delta-function, δ(P1 − 1). Additionally, the information
about the values of the occupation probabilities is based on the relative number of particular
outcomes of the system, X1 and X2 , and one would need a great number of experiments
in order to detect the difference between the actual value of the occupation probability and
¯
its saturation value. In fact, in the limit P1 → 1, the number of measurements becomes
unbounded. These extremes are consistent with the trends shown in Figs. 2.8 and 2.9.
¯
However, between these limits there is a range of ∆P1 , and thus a range of changes in λ, for
30

which the measurement time can be relatively short. This is the remarkable feature of the
dynamic bridge as a measurement tool.

2.4

Summary

The analysis presented in this chapter was originated by the idea to employ the noisy switching dynamics of bistable systems for purposes of sensing. The applicability of the dynamic
bridge requires only that one couple the quantity of interest with the system and be able
to measure shifts in the occupation probabilities. After presenting an introduction into the
general theory of activated escape, we analyzed the Duffing resonator subject to the nearresonant periodic driving as an example of a bistable system, with additive white Gaussian
noise to induce switching. In Section 2.2 we discussed the switching dynamics of the resonator on the time scale related to the relaxation time of the system. In the next section we
quantified how bistable systems can be used as a measurement tool for detecting changes in
system parameters, examining sensitivity and measurement time.
We determined the sensitivity limits of the proposed paradigm and showed that the
dynamic bridge can be extremely sensitive to shifts in system parameter values over a certain
range. This allows one to detect and estimate small changes in the parameter of interest,
although the method can require relatively long measurement times. However, we determined
that, in an intermediate range of parameter shifts and the attendant occupation probabilities,
the total measurement time is comparatively short. Additionally, the proposed technique
allows one to track changes in large variety of system parameters, but only in those that
change the natural frequency of the system. These results are quite promising since they
provide confidence that this measurement paradigm can be competitive with other methods,
31

for example, the measurement of the shift in the natural frequency.

32

Chapter 3
Non-Gaussian Noise Detection
In Chapter 2 we discussed the balanced dynamic bridge as a measurement tool that can
be employed in bistable systems with additive stochastic forcing. However, there is another
possible application of this class of stochastic nonlinear systems, specifically, these systems
can be used to detect the presence of non-Gaussian noise acting on the system, and estimate
its statistical properties. In many applications this information, about an unknown stochastic
processes, is difficult to access. The possibility to extract this information from the system
response is important and relevant, for example, in the problems of the electron transport
through a nanowire when one considers the current flow as a series of discrete random pulses
having their own statistical properties.
In this chapter we develop a theory that describes how information about a non-Gaussian
noise process can be obtained using the balanced bridge. The main idea is that non-Gaussian
noise will change the occupation probabilities in a manner that can be quantified, whereas an
additional additive Gaussian noise will simply shift the rates, not the balance. In Section 3.1
we give an introduction to the theory of the activated escape in the presence of additive
33

non-Gaussian stochastic excitation, and in Section 3.2 we describe the general detection
technique. Section 3.3 provides an example of the estimation of the statistics of a nonGaussian noise process using a simple first-order system for which the calculations can be
done explicitly. We then discuss more complex example of a nonlinear Mathieu resonator, and
apply the detection procedure to this system. Section 3.5 closes the chapter by highlighting
important aspects of the analysis.

3.1

Switching in the Presence of Non-Gaussian Noise

In Section 2.1 we discussed in detail the theory of activated escape in the presence of additive
white Gaussian noise. This noise affects the system by causing random switching between
stable states. In the present chapter, we extend this theory in order to capture another
interesting aspect of the system dynamics that arises under the influence of an additional
non-Gaussian modulation [10]. This extra noise source modifies the usual switching dynamics
of the system of interest in some way, and its understanding is essential for the discussion of
the proposed technique for estimating the properties of the non-Gaussian noise.
We start the discussion with a general nonlinear system of the form,

˙
q = K(q) +

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

(3.1)

where q ∈ Rn is the state variable vector, K ∈ Rn is the nonlinear vector field, f ∈ Rn
is the vector of white Gaussian noise processes with effective strength D, and G ∈ Rn×m
is the coefficient matrix of the non-Gaussian noise vector ξ ∈ Rm . It is assumed that ξ
is independent of f . Hereafter, we also assume that the Gaussian noise f is weak, and
34

the strength of the non-Gaussian modulation ξ is such that we consider the non-Gaussian
excitation as a perturbation to the main switching process. Suppose that in the absence of
all stochastic terms the deterministic system has two stable fixed points, q A and q B , and one
unstable fixed point, q S , whose separatrix is the boundary between the basins of attraction
for q A and q B , as depicted in Fig. 3.1.

Uq

qA

qS

qB

q

Figure 3.1: Schematic one-dimensional representation of the potential for a bistable nonlinear
system, Eq. (3.1).

In the weak noise limit, the system motion is nearly concentrated about one of the
stable fixed points, say q A , for some time. However, due to the stochastic excitation the
response can jump from the one stable state to another. For small noise these events are rare.
Suppose also that we know the switching path that system follows during these transitions
when ξ(t) = 0, and we also assume that the realization of the non-Gaussian stochastic
vector, ξ(t), is also known. Using Eq. (3.1), one can express the white Gaussian noise in the
following way,
˙
f (t) = D−1/2 q − K(q) − G(q)ξ(t) .

(3.2)

The corresponding probability for such a noise realization can be expressed through its
35

probability density functional [1, 10],

P [f (t)] ∝ exp −

1
2

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

(3.3)

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

(3.4)

where F is the inverse pair correlator of f (t) and we assume that f possesses time reversal
ˆ
ˆ
ˆ
symmetry so that F(t) = F(−t) = F T (t). It follows from Eq. (3.2) that the integrand of
Eq. (3.3) is inversely proportional to the characteristic noise strength, D, and in the weak
noise limit, the probability for the transition to occur is sharply peaked about the most
probable switching path, see Section 2.2. This path minimizes the integral expression in
Eq. (3.3), where f is given by Eq. (3.2).
The switching rate depends exponentially on the action calculated along the specific path
between any two points of the configuration space. However, due to the sharpness of the
probability ditribution about the most probable path, one can approximate the expression
for the mean switching rate assuming that system always follows the most probable switching
path. So, for ξ = 0, the mean switching rate reads,

r0 ∝ exp [−S0 /D] ,

(3.5)

where the action S0 should be calculated using Eqs. (3.2) and (3.3) in the following way,

S0 = min

1
2

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

.

(3.6)

˙
˙
In the above equation (q − K)t is meant to imply q(t) − K(q(t)), and the minimum is taken
36

over all paths that the system can follow from q A to q B .
Since f is a white Gaussian noise vector, its autocorrelation matrix reads,

f (t)f T (s) = Iδ(t − s),

(3.7)

where I denotes the identity matrix. Further, the inverse of the correlation function simpliˆ
fies, giving F(t) = Iδ(t), and therefore the action S0 becomes,

S0 = min

1
2

˙
|q − K|2 dt .

(3.8)

Equation (3.8) has exactly the same form as Eq. (2.13), so the canonical momentum p reads,

˙
p ≡ q − K(q).

(3.9)

At this point it is convenient to introduce the so-called auxiliary Hamiltonian, given by,

1
H = |p|2 + pT K(q) ,
2

for which the conjugate momentum p is the optimal noise realization multiplied by

(3.10)
√

D.

As mentioned in Section 2.2, the switching path that minimizes Eq. (3.8) is the composed
of two parts. The first part is the heteroclinic trajectory going from q A to q S , which is
necessarily noise-driven, so that p = 0 along this path. Once the system reaches the noisefree saddle point q S , there is a 50-50 chance of the system going to either of the noise-free
fixed points, which it does by tracking either branch of the noise-free unstable manifold of q S .
Therefore, the second part of a switching trajectory is the noise-free heteroclinic trajectory
37

going from q S to q B , on which p = 0. This second part of the switching trajectory does not
contribute to the value of the switching probability, and so we ignore it for the remainder of
the discussion.

Now, consider the case with the non-Gaussian noise, that is, with ξ = 0. Assuming that
f and ξ are independent noise processes and ξ is small compared to f , the mean switching
rate can be expressed, cf.[1], as

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

(3.11)

where the subscript ξ indicates that the exponential term should be averaged over all possible
realizations of ξ. Using Eq. (3.2) we can expand S[ξ] in a Taylor series, which up to the first
order gives the following result,

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

S[ξ] = S0 +

(3.12)

where χi is so-called the susceptibility of the system to a perturbing non-Gaussian process
as it follows out of ith state, and S0 is the action calculated along the unperturbed path.
From Eqs. (3.2) and (3.3), the susceptibility has the form,

χT (t) = −
i

˙
q i (t ) − K(q i (t ))

T

ˆ
F(t − t)G(q i (t))dt .

(3.13)

Following Billings et al., [1], the susceptibility χi can be considered as a direction along which
the additional non-Gaussian modulation is projected, for a given unperturbed response q i (t).
This projection changes the generalized work done on the system as it follows a trajectory,
38

which for present purposes is the optimal switching trajectory.
As a result of this additional noise, the activation energy of the system changes, but because of the weakness of the white Gaussian noise, the total correction of the mean switching
rate from the non-Gaussian noise will have different order of magnitude. The new switching
rate in the presence of this perturbing noise can be expressed as,

r ∝ Ar0 ,

(3.14)

where the prefactor A is the term that describes the correction of the switching rate r0 when
ξ = 0. Using Eq. (3.12), one can obtain the following form for this correction prefactor,

A=

exp −

1
D

χT ξ i dt

= Φξ [iχi /D],

(3.15)

ξ

where Φξ is the characteristic functional of ξ, cf.[10].
In the case when the non-Gaussian modulation is much weaker than the primary white
Gaussian noise, the whole procedure of the calculation the switching rate reduces to the
calculation of the correction prefactor A, and this is the essential aspect of the non-Gaussian
noise detection technique.

3.2

Detection Technique

In this section we construct a theoretical approach that allows one to detect the presence of
a small non-Gaussian driving term and determine its statistical properties. Suppose we have
a system described by Eq. (3.1), and it is assumed that we have the ability to turn ξ on and
off. By turning on/off the non-Gaussian noise we are able to remove its mean by examining
39

Uq

Uq

q

q

(a) ξ(t) = 0.

(b) ξ(t) = 0.

Figure 3.2: The change in the probability ratio due to exposure the balanced dynamic bridge
to a non-Gaussian perturbation.
the mean fluctuation of the coordinate q about the fixed points q A or q B and applying a
required bias that compensates for the mean of ξ, essentially recalibrating the system.
In general, the addition of the perturbing noise ξ changes the transition rates, see Section 3.1, in both switching directions, and it shifts the relative occupancy ratio in the system.
The ratio of occupation probabilities for the system to be in state 1 and 2, respectively, is
given by,
(0)

P1
r
A r
= 2 = 2 2 ,
P2
r1
A1 r(0)
1
(0)

where ri

(3.16)

is the switching rate out of the ith state when ξ = 0. If the bridge is operated
(0)

(0)

at a balance point before exposing it to the non-Gaussian perturbation, we have r1 = r2 ,
and the resulting probability ratio is equal to the ratio of the correction prefactors, given in
Eq. (3.15). It is convenient to introduce the random variable,

zi =

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

(3.17)

where χi (t) is given by Eq. (3.13), and zi is the increment in the action due to the per(i)
turbation. Let κn be the nth cumulant of zi . From Eq. (3.15), Ai is simply the moment

40

generating function for zi , and ln Ai is the cumulant generating function [15]. So, one has
∞

ln Ai =
n=1

(i)

κn
.
n!(−D)n

(3.18)

Note, if ξ is Gaussian, then zi is also Gaussian, and κn = 0 for all n > 2. Substitution of
Eq. (3.18) into Eq. (3.16) yields,

P1
= exp
P2

∞
n=1

(2)

(1)

κn − κn
n!(−D)n

.

(3.19)

At this point it is assumed that ξ has already zero mean because we eliminated the effects of
(1)

(2)

its mean by applying the required bias. If the second cumulants κ2 and κ2 can be balanced,
a condition considered in more detail below, they do not contribute to the exponential term in
Eq. (3.19), and in this case the population ratio depends only on the higher order cumulants,
that is, on the non-Gaussian nature of ξ. The second cumulant of zi can be expressed as
follows,
(i)

κ2 =

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

Using Eqs. (3.13) and (3.20) it can be shown that κ2

(2)

and κ2

(3.20)

are equal if

ˆ
dsds F (s − t)G(q i (s)) ξ(s)ξ T (s )
(3.21)
×GT (q

i (s

ˆ
ˆ
))F (s − t ) = cF (t − t ),

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

41

processes, and G ∝ I. If the mean has been removed, and if the second cumulants cancel,
then the population ratio of the bridge reflects 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.3. Note that while the second cumulants drop out from Eq. (3.19), they remain in
Eq. (3.18), and thus they change the switching rates, but not the occupancy ratio. 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.

3.3

An Example: Shot Noise Measurement in 1D

In this section we consider an example of the non-Gaussian noise detection using a onedimensional system for which the calculations can be carried out explicitly. These results
are relevant to the Duffing resonator described in Chapter 2 near the cusp of the system
bistability, that is, near Ω−2 = 1/3 (see Fig. 2.6), where the system dynamics effectively
collapses onto a one-dimensional slow manifold, [14]. On the balance curve, the system behaves according to the classic overdamped symmetric double well quartic potential. Suppose
that such a system is subjected to a shot noise with unknown parameters, and the goal is
to estimate them by measuring the switching statistics. Following [10], the shot noise is
characterized by the pulse area, g, and the pulse rate, ν. Thus, the equation describing the
dynamics of the system is

q = −U (q) +
˙

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

(3.22)

where q is the scalar state variable, f (t) is a white Gaussian noise with f (t)f (t ) = δ(t−t ),
and ξ(t) is a Poisson process. In the context of the Duffing resonator, we assume that the
shot noise acting on the resonator is modulated by the resonator driving signal so that it
becomes a simple Poisson process on the slow manifold. Since the mean value of the Poisson
noise
ξ = νg,

(3.23)

we introduce the constant driving term −νg in Eq. (3.22) in order to make the bridge sensitive
only to higher-order statistics, as discussed in Section 3.2. Here we assume that, in spite the
fact that ν and g are unknown quantities, their product can be measured and compensated
for, which is possible comparing the response near one of the fixed points with ξ(t) turned
off and on and measuring a shift in the mean of q. The system is assumed to move in a
symmetric double-well potential, as shown in Fig. 3.3, described by

1
1
U (q) = − q 2 + q 4 .
2
4

g

q1

(3.24)

0

qs

q2

Figure 3.3: The double-well potential Eq. (3.24); q1 and q2 denote positions of the stable
fixed points, while qS refers to the saddle-point. Poisson pulses with positive area g change
the switching rates, such that r1 > r2 , therefore shifting the probability from q1 towards q2 .
Following Feynman, [10], the cumulant generating function for the zero-mean Poisson
43

process, ξ − νg, is given by

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

ln Ai = ν

(3.25)

Now, Eq. (3.18) still holds, and equating like powers of D gives,

(i)

κ1 = 0,

(3.26a)

(i)

κn = νg n

dtχn (t),
i

n > 1.

(3.26b)

The susceptibility χi can be understood as a direction for the Poisson noise projection
computed along the heteroclinic solution of the auxiliary Hamiltonian system. For this
system, the auxiliary Hamiltonian reads,

H=

p2
− pU (q),
2

p = q + U (q),
˙

(3.27)

and the heteroclinic trajectories are those that 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.28)

that is, the time-reversed deterministic trajectories. For the potential given in Eq. (3.24),
these solutions are,

(h)

q± (t) = ± 1 + e2t
44

−1/2

(3.29)

and the attendant susceptibilities are given by,

(h)

χ1,2 = −2q± (t) = ±2e2t 1 + e2t
˙

−3/2

.

(3.30)

Figure 3.4: Change in the probability ratio due to addition of the zero mean shot noise,
Eq. (3.59). Data are results from Monte-Carlo simulations with D = 0.04 and ν = 0.5. The
inset shows the same quantities over a larger range, demonstrating the limitations of the
approximate theory.
Using equation (3.28), the integrals in Eq. (3.26b) can be evaluated explicitly. Applying
the result to Eq. (3.19) yields,

P1
= exp −ν
P2

n=1

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

(3.31)

where Γ(n) is the Euler Gamma function. Figure 3.4 shows ln(P2 /P1 )/ν as a function of g/D
along with results from Monte-Carlo simulations. The inset in the figure shows the same
quantities as the primary axes over a larger range. The Monte-Carlo simulations were carried
45

out with parameter values D = 0.04 and ν = 0.5. The weak noise limit assumed by the
analysis requires D/(2∆U )

1, where ∆U is the barrier height. We also assumed that the

non-Gaussian noise is weak compared to the Gaussian noise. This is true if νg 2
parameters used give D/(2∆U ) = 0.08 and require g/D

D. The

7. Figure 3.4 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 ξ and yields 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 the following explicit expression for g 2 as a function of the measured
population ratio, the mean of the shot noise, and the strength of the Gaussian noise,

g2

715 2
=−
D +
32

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

(3.32)

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. For positive g, ξ is also positive. Moreover, since the
system is overdamped, pulses pushing in the positive direction can only decrease P1 while
increasing P2 . Thus, ln P1 /P2 < 0 and ξ > 0 for g > 0 and the solution gives a positive
result. On the other hand, when g < 0, we have ln P1 /P2 > 0 and ξ < 0, and we still
obtain a positive result. The sign of g is found from the sign of ln(P1 /P2 ).
Equation (3.32) can be considered the final result for this example, from which g can be
determined. Use of this equation requires measurement of the quantities D, ξ , and P1 /P2 .
The Gaussian noise strength D and the mean of the non-Gaussian noise ξ can be measured
from the quasi-steady-state distribution of the drift vector when the system stays in the
neighborhood around a stable fixed point. The Gaussian noise strength can be determined
46

from the width of the distribution about the fixed point, and the mean is determined from
the force required to rebalance the distribution, as described above. To reiterate, with ξ
turned off, the mean of this quasi-steady-state distribution is measured. With ξ turned on,
this mean will 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 latter quantity is required in
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 is sampled over a sufficiently long period, and
the probabilities are estimated by use of Eq. (2.35). From these three measurements, the full
counting statistics can found, in theory. Of course, there are limitations to the measurement
resolution, etc., that limit the amount of information that can be attained.

3.4

Noise in Parametrically Excited Resonators

In this section we consider an example of non-Gaussian noise detection in a more complicated
system by considering a nonlinear Mathieu resonator as a model. This model is relevant to
parametrically excited MEMS, which have received a lot of attention, especially in terms
of sensor development [24, 5, 4, 2]. Here we consider the situation in which the noisefree system exhibits a pair of stable steady-state responses that are identical except for a
phase shift relative to the excitation. These arises from a time translation symmetry in the
response [6], which is broken if the system is subjected to external resonant excitation[22, 24].
In the present case we consider random switching between these states induced by additive
Gaussian noise, which has been examined experimentally in MEMS [2, 3]. In order to detect
the non-Gaussian noise it is modulated by resonant external (additive) excitation, and thus
this additional noise is tied with symmetry-breaking of these responses. The equation of
47

motion for the resonator with these properties is given by,

2
q + 2Γq + (ω0 + h cos ωt)q + γq 3 =
¨
˙

√
Df (t) + ξ(t) cos

ω
t+ψ ,
2

(3.33)

where q represents the system configuration, Γ is a coefficient of linear dissipation, ω0 is the
resonator natural frequency, h is the amplitude of the parametric forcing with frequency ω,
γ is a nonlinearity parameter, f is a white Gaussian noise, and ξ is a Poisson process that is
modulated by the harmonic of frequency ω/2.

For convenience we separate the mean and alternating parts of ξ in the following way

ξ(t) = ξ + ξa (t),

(3.34)

where ξa (t) is the alternating part of the Poisson process. Suppose that the system driving
is resonant, i.e. ω = ω0 + σ, with σ

ω0 , which leads one to the following convenient form

for the detuning,
2
ω0 =

ω2
+ δ,
4

where δ = −ω0 σ.

(3.35)

Similarly to what was done in Section 2.2 for the Duffing equation, we introduce the Van
der Pol transformation,

q=

ωΓ
iωt
u exp
+ c.c.,
3γ
2

(3.36a)

q=
˙

ωΓ ω
iωt
ui exp
+ c.c.,
3γ 2
2

(3.36b)

which, when applied to the equation of motion, results in the following form of the equation
48

governing the complex amplitude u,

u = −u + iΩu + iα¯ + i|u|2 u − i
˙
u
−i

3γD
iωτ
f (τ/Γ) exp −
3 Γ3
2Γ
ω

3γ
ωτ
iωτ
ξ (t) cos(
+ ψ) exp −
+ β(sin ψ − i cos ψ),
3 Γ3 a
2Γ
2Γ
ω

(3.37)

where we have switched to the non-dimensional time τ = Γt. Also we introduce the set of
non-dimensional parameters Ω, α, and β, defined as,

δ
,
ωΓ
h
α=
,
2ωΓ
3γ ξ
,
β=
ω 3 Γ3 2
Ω=

(3.38a)
(3.38b)
(3.38c)

where Ω captures the detuning, α describes the strength of the parametric excitation, and
β is the bias induced by the mean value of the non-Gaussian noise.

We first briefly review the dynamics of the nonlinear Mathieu resonator subject to deterministic forcing only; for details we refer the reader to [21]. First, we define u = ur + iui ,
where ur and ui are real and imaginary parts of the complex slowly-varying amplitude u.
Then, the noise-free dynamics of the system obeys the following equations,

ur = −ur − Ωui + αui − ui (u2 + u2 ),
˙
r
i

(3.39a)

ui = −ui + Ωur + αur + ur (u2 + u2 ).
˙
r
i

(3.39b)

Fixed points of Eqs. (3.39a) and (3.39b), which represent near-harmonic periodic responses
of the original system at frequency ω/2, can be found by setting ur = ui = 0, which has
˙
˙
49

three possible solutions,

A1,2 =
φi1 =

α2 − 1,

−Ω ±
1 −1
sin
2

1
,
α

A3 = 0,

φi2 =

1 −1
sin
2

(3.40a)
1
α

+ π,

i = 1, 2.

(3.40b)

Thus, depending on the value of Ω one can have either one, three or five possible solutions
to Eqs. (3.39a) and (3.39b).

These steady-state amplitudes are depicted in Fig. 3.5 as a

u2
5

3

1
10

8

6

4

2

2

4

Figure 3.5: Dependence of the amplitude of steady-state oscillations on the normalized
√
detuning parameter Ω. Parameter α is taken to be equal 17.

function of the frequency detuning parameter Ω. The case of a hardening system is shown,
whereas the case of a softening system is identical with the sign of Ω reversed. Of course,
the zero solution always exists, but can be unstable near the resonance. When an instability
zone exists, the nontrivial responses branch out from the zero solution. For each non-zero
amplitude, in turn, there exist two corresponding phases, Eq. (3.40b), and these fixed points
are degenerate in the sense that they have the same amplitude but differ in phase by π.
On Fig. 3.5 these responses coincide, so we see only two nontrivial branches, while a phase
50

portrait, see Fig. 3.6(a), shows that there are indeed two distinct fixed points for each
nontrivial branch. From this picture it is also clear that these two points are physically
indistinguishable from one another, and switching involves relatively rapid changes in the
phase of response by π. From this symmetry, the Mathieu resonator with additive Gaussian
noise will be automatically balanced, and from Section 3.2 it is seen that any zero-mean
stochastic perturbation does not shift the operating point of the dynamic bridge away from
the balance curve; this type of noise simply affects the switching rate.

Figure 3.5 suggests that there exist three regions in (α, Ω) parameter space with distinct
√
types of response. The first region, Ω > α2 − 1, where only the trivial response exists and
√
is stable. The second region, |Ω| < α2 − 1, has three fixed points with zero stable and a
√
symmetric pair of stable solutions. The third region, Ω < − α2 − 1, has five fixed points,
one at zero and two symmetric pairs of nontrivial responses, one pair stable (at larger
amplitude) and one pair unstable.

Since the dynamic bridge operates using a bistable

system, the naturally balanced response in the second parameter region described above is
√
of interest here, and its boundaries in (α, Ω) parameter space are |Ω| < α2 − 1. This
region is depicted in Fig. 3.7. Note that typically the region depends on the level of damping,
that is, on Γ, but the scaling used here eliminates that parameter, reducing the number of
√
parameters for the deterministic system to two. The third region, Ω < − α2 − 1, can also
be treated as a bistable one if we consider the transitions between stable zero and non-trivial
solutions. However, in this case the system is not symmetric, requiring prior balancing,
which is essentially the case of the Duffing oscillator.

Now, we include in Eqs. (3.39a) and (3.39b) the term representing the mean value of the
51

3
2

ui

1
0
-1
-2
-3
-1

-0.5

0
ur

0.5

1

(a) Symmetrical phase portrait corresponding to Eqs. (3.39a)
and (3.39b).

3
2

ui

1
0
-1
-2
-3
-1

-0.5

0
ur

0.5

1

(b) Breaking the symmetry in the phase portrait due to a constant
(non-zero) mean value of the Poisson perturbation. β = 1 and
ψ = 30◦ .
Figure 3.6: Phase portrait of the nonlinear Mathieu resonator.

52

10
5
0

1

3

5

7

9

Α

5
10
Figure 3.7: Shaded area represents the bistable region of the nonlinear Mathieu resonator.

non-Gaussian process, β. The corresponding equations read,

ur = −ur − Ωui + αui − ui (u2 + u2 ) + β sin ψ,
˙
r
i

(3.41a)

ui = −ui + Ωur + αur + ur (u2 + u2 ) − β cos ψ.
˙
r
i

(3.41b)

Then, applying the condition for the fixed points, uro = uio = 0, we find
˙
˙

uro = β

sin ψ − {α − (Ω + |uo |2 )} cos ψ
,
1 − α2 + (Ω + |uo |2 )2

(3.42a)

uio = β

{α + (Ω + |uo |2 )} sin ψ − cos ψ
,
1 − α2 + (Ω + |uo |2 )2

(3.42b)

where |uo |2 = u2 + u2 . One can combine Eqs. (3.42a) and (3.42b) to obtain an expression
ro
io
of the form,
β 2 = g(|uo |2 , α, Ω, ψ),
53

(3.43)

where

g(|uo |2 , α, Ω, ψ) =

|uo |2 [1 − α2 + (Ω + |uo |2 )2 ]2
.
1 + α2 + (Ω + |uo |2 )2 − 2α sin 2ψ + 2α(Ω + |uo |2 ) cos 2ψ

(3.44)

The resonator phase which is defined as φ = tan−1 [uio /uro ], in turn, has the following form,

φ = tan−1

{α + (Ω + |uo |2 )} sin ψ − cos ψ
.
sin ψ − {α − (Ω + |uo |2 )} cos ψ

(3.45)

Solutions to Eq. (3.43) are graphically depicted in Fig. 3.8. From this plot we can conclude

g uo 2, Α, , Ψ , Β2

15

10

5

0
0

1

2

3

4

5

uo 2
Figure 3.8: Solutions to Eq. (3.43) are represented by the intersection points of the solid
and the dashed lines. The left-hand side of Eq. (3.44) as a function of |uo |2√ drawn by the
is
2 is drawn by the dashed line. Parameters values: α = 17, Ω = 0, ψ =
solid line, while β
30◦ , β = 1.

that when β = 0, i.e., when there is no noise is the system, we recover the noise-free dynamics
perfectly. In this case the amplitude of the degenerate fixed points is represented by the single
point where the solid line touches (in fact, is tangent to) the horizontal axis. If we introduce
54

the non-Gaussian noise with non-zero mean value, the position of all fixed points are shifted
in both amplitude and phase, see Eqs. (3.43) and (3.45). This bias term, proportional to β,
breaks the symmetry, or splits the degeneracy, in the nonlinear Mathieu resonator, as shown
in Fig. 3.6(b). Now it is clear why we preserve the non-zero mean of the non-Gaussian
excitation, and do not compensate for it as was the case in Section 3.3, since it causes the
system operate in a desirable way. The only constraint that we put on the bias is that is
should be sufficiently small to maintain the resonator in the bistable regime.

Additionally, we see that this symmetry breaking changes the heteroclinic trajectories
corresponding to the primary switching due to the Gaussian excitation. Strictly speaking,
as we increase the mean value of the non-Gaussian excitation, the system behaves more like
the Duffing resonator. This fact can be used as a method for the estimation of ξ , which
we will need subsequently.

We next examine what can be done in terms of the detection and the estimation of
the statistics of the non-Gaussian input. To begin, we write the equations describing the
dynamics of the system in the following form,

u = K(u) + f (t) + G(u)ξ(t),
˙

(3.46)

where


2
 −ux − Ωuy + αuy − uy (ux

K(u) = 

−uy + Ωux + αux + ux (u2
x

f (t) = −

ωτ
Γ



+ u2 ) + β sin ψ 
y
+ u2 ) − β cos ψ
y

(3.47)



3γD  sin

 f (τ ),
3 Γ3 
ω
ωτ
cos Γ
55

,

(3.48)


G(u) = −



3γ 1 0

,
ω 3 Γ3
0 1


ξ(t) = ξa (t) cos

(3.49)



ωτ
2Γ 

ωτ
 sin
+ψ 
.
2Γ
ωτ
cos 2Γ

(3.50)

According to the previously stated assumption, the non-Gaussian stochastic process is much
weaker than the primary Gaussian input. Only in this case can we apply the perturbation
theory to the system and treat ξ as a perturbation. Now, it is noted that Equation (3.46)
has exactly the form of Eq. (3.1) that was introduced in Section 3.1, which implies that
we can apply the line of reasoning described in Sections 3.1 and 3.2. Then, the first-order
correction to the activation energy reads,

Si =

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

(3.51)

where the system susceptibility is given by,

χT (t) = −
i

1
2

ˆ
˙
dt (ui − K i (u))T F (t − t )G.

(3.52)

Since f (t) is the vector of white Gaussian noise processes, the inverse of its correlation matrix
is simple to obtain and is given by,
4ω 3 Γ3
ˆ
F (t − t ) =
Iδ(t − t ),
3γ

56

(3.53)

which significantly simplifies the form of the system susceptibility, as follows,

χT (t) =
i

4ω 3 Γ3
˙
(ui − K i (u))T .
3γ

(3.54)

Further, by using results of Section 3.1 and, in particular, Eq. (3.15), it is possible to determine the attendant correction prefactor to the switching rate, Ai , which is found to be,

exp −

A(i) =

(i)

χT (t)
i
ξ(t)
D

= exp −

χef f (t)

ξ

D

ξa (t) ,

(3.55)

where we have introduces the so-called effective susceptibility, which reads,

(i)

χef f (τ ) =

4ω 3 Γ3
ω(τ − τ0 )
(i)
{(ur − Kr (u(i) )) sin
˙
3γ
2Γ

(3.56)

ω(τ − τ0 )
ω(τ − τ0 )
(i)
+ (ui − Ki (u(i) )) cos
˙
} cos(
+ ψ) τ .
0
2Γ
2Γ
Each switching event between the coexisting stable states in unperturbed system is the same
in its dynamics, but occurs at random moments of time. We know that the mean escape
rate depends exponentially on the action along the most probable path between two states.
As a result of such switching, the effect of the Poisson noise on the transitions depends on
the relative phase between the Poisson process and the unperturbed heteroclinics. In order
to take into account all possible cases of such interactions, we have to average the effective
susceptibility over the whole range of these relative phases, which amounts to averaging over
ξ in Eq. (3.55). Integration of Eq. (3.56) is straight-forward and gives

(i)

χef f (τ ) =

ω 3 Γ3
(i)
(i)
[−(ur − Kr (u(i) )) sin ψ + (ui − Ki (u(i) )) cos ψ].
˙
˙
3γ

57

(3.57)

3
2

ui

1
0
-1
-2
-3
-1

-0.5

0
ur

0.5

1

(a)
3
2

ui

1
0
-1
-2
-3

-1.5

-1

-0.5

0
ur

0.5

1

1.5

(b)
Figure 3.9: The numerical solution to the four-dimensional system of equations corresponding
to the nonlinear Mathieu resonator depicted as a projection onto: (a) the (ur , ui ) plane, and
(b) the (pr , pi ) plane. The noise-free dynamics is exactly the same as shown in Fig. 3.6(b),
and is shown in black dotted lines. Red and blue lines are the projections of the noisy
switching trajectories from different stable states to the saddle.

58

Introducing this effective susceptibility allows one to consider the nonlinear Mathieu resonator from the point of view of one-dimensional system discussed in Section 3.3. In particular, one just needs to replace χ by χef f in the corresponding expression for the characteristic
functional for ξa
(i)

∞

ln A(i) = −ν

1−

χef f (t)g
D

−∞

(i)

− exp

χef f (t)g
D

dt,

Then we can get expressions for the cumulants as it is done in Section 3.1

(i)

κ1 = 0,
(i)
κn

=

(3.58a)

νg n

∞
−∞

(i)

[χef f (t)]n dt,

n > 1.

(3.58b)

The expression for the logarithm of the probability ratio which is directly measured reads

P
ln 1 =
P2

∞
n=2

(2)

(1)

κn − κn
=ν
n!(−D)n

(2)
T
n
0 ([χef f (t)]
n
g

∞

(1)

− [χef f (t)]n )dt

n!(−D)n

n=2

.

(3.59)

The last expression in Eq. (3.59) is an infinite polynomial in powers of g, where the mean
pulse rate, ν, can be expressed in terms of the known mean of the Poisson noise, ξ , as
follows,
ν=

ξ(t)
.
g

(3.60)

Finally, we can write an expression containing only g as an unknown,

P
ln 1 = ξ(t)
P2

∞

(2)
T
n
0 ([χef f (t)]
n−1
g

(1)

− [χef f (t)]n )dt

n!(−D)n

n=2

59

,

(3.61)

which, after making a reasonable truncation, leads one to a polynomial of finite size that
can be solved for g, and then for ν, using Eq. (3.60).

3.5

Summary

In the present chapter we developed the theoretical background that provides a means of
detecting and characterizing random excitation of a non-Gaussian type acting on a balanced
dynamical bridge. From this perspective, Eq. (3.19) is the main result of the analysis. The
theory developed is applicable to some engineering problems, for example, the estimation of
the statistics of the current fluctuations in electro-mechanical devices whose models exhibit
the bistability in some range of system parameters. In the analysis we do not restrict the
nature of the non-Gaussian noise in any manner; the only assumptions made in Section 3.2
relate to the intensity of the additional driving and its independence from the Gaussian noise
excitation that causes the underlying switching process of the bridge system. Specifically, the
strength of the non-Gaussian noise is assumed to be much smaller than the corresponding
intensity of the Gaussian forcing, so that it can be treated essentially as a perturbation.
In conclusion, we considered two examples of the estimation of unknown parameters of
a Poisson process in systems with different complexity, namely, a simple one-dimensional
model and a parametrically excited resonator model. According to the results obtained, we
can report that when system is operated near the cusp, its dynamics reduces in order and the
whole analysis can be done analytically in closed form. This observation is import because it
may help significantly with the design of systems used for the estimation of some stochastic
processes.
Looking at contemporary methods of non-Gaussian noise detection, it should be noticed
60

that the dynamical bridge offers an attractive alternative for this purpose in some cases.
For instance, one of the conventional methods that works with Josephson junctions involves
conducting two distinct experiments and measuring the difference in the statistics of currents
going through the junction in both directions [18]. The bridge, in contrast, benefits in the
measurement time since it requires a only single experiment involving estimation of the shift
in the switching rates.

61

Chapter 4
Conclusions
In this work we discussed new measurement applications of nonlinear micro/nano-electromechanical systems (M/NEMS) that are subject to a stochastic forcing in combination with
some periodic excitation. Since noise processes of different kinds are an important part of
many M/NEMS, due to their small dimensions, our primary goal in this study was to exploit
these noise sources in order to explore new types of M/NEMS applications.
In Chapter 2 we examined the concept of the balanced dynamic bridge as a new measurement paradigm in noise-driven micro-resonators. This analysis is an attempt to employ
bistable M/NEMS subject to a white Gaussian excitation as a sensor that would allow one
to measure changes in some parameters of interest with high accuracy. The general results
of this analysis are very promising and potentially important. In particular, we found that
a bistable micro-resonator is extremely sensitive in the limit of the weak noise driving. The
main idea is summarized as follows: First, we tune the dynamic bridge in such a way that
the occupation probabilities for the system to be in each stable state are equal to each other.
In order to measure changes of the particular parameter of interest, we must then couple it
62

somehow to the bistable resonator in a way that changes in this parameter affects a system
variable which is easy to access through the measurement process. Then, we have to perform a sufficient number of measurements on this variable, and the number of measurements
should be estimated a priori, based on a reasonable range of expected changes of the parameter. Finally, we need to measure the shift in the occupation probabilities, from which we
estimate the associated parameter change by means of Eq. (2.40).
As analysis shows, this new measurement approach can take quite long time to give
sufficient precision. However, for a range of changes of the parameter, it requires a comparatively small number of measurements. Because of this, we suppose that the balanced
dynamic bridge can be effectively used as a measurement tool in M/NEMS where stochastic
processes play an essential role.
In Chapter 3 we discussed the possible application of noise-driven M/NEMS, namely, the
detection of a non-Gaussian noise acting on the system of interest, in addition to the white
Gaussian excitation that causes the switching. We exploit again the balanced dynamic
bridge as a starting configuration of the resonator driven by the noise. General analysis
reveals a remarkable feature of the additional non-Gaussian modulation as it acts on the
system. This extra forcing affects the switching rates in both directions, but generally with
a bias, providing a shift in the ratio of the occupation probabilities. We found that the
non-Gaussian random process, as projected onto the system susceptibility, results in higher
order statistics that change the switching rates of the bridge differently for the opposite
switching directions. In fact, while the first and the second moments/cumulants simply
renormalize the effective properties of the Gaussian driving term, moments/cumulants of
order three and higher distinguish the non-Gaussian noise. As a result, by looking at the
63

changes in the occupation probability ratio, we can determine if there is a non-Gaussian
noise source in the system of interest. We constructed a theoretical approach which allows
one to extract information about the non-Gaussian process if one can guess its type, e.g.,
Poisson or telegraph. The analysis was supported by two examples that showed the potential
usefulness and effectiveness of the proposed technique.

64

APPENDIX

65

A

Path Integral Formulation

In this Appendix we review Feynman’s path integral formulation [10], since it is an important
mathematical tool that is used the calculation of the transition probability. It should be noted
that the path integral generally appears in systems that are quantum-mechanical. However,
this approach is also very useful in the study of stochastic dynamic systems.
We start with the assumption that the system can follow different paths in the configuration space as it undergoes transitions from a point A at time tA to a point B at time tB ,
as depicted in Fig. A.1. For each of these possible trajectories there exists a probability of

Figure A.1: Schematic representation of a quantum-mechanical system that has more than
one possible trajectory between points A and B. Different paths have different probabilities
of realization, and these paths can significantly deviate one from another.

66

its realization that can be expressed in the form,

P [A → B](k) ∼ exp

tB

i

tA

˙
L(xk (t), xk (t), t)dt ,

where L is the Lagrangian of the system of interest [11],

(A.1)

is Planck’s constant, and k an

index corresponding to the possible trajectory.

The total probability that the system will make a transition between the two points is
computed by performing so-called time-slicing, meaning that we divide the time interval
[tA , tB ] into n equal discrete steps,
t − tA
,
n= B

= ti+1 − ti .

(A.2)

As a result, the coordinate vector x also becomes discrete, and for each time ti there is a
corresponding xi , and the system path can now be represented by the set {xi }. Consequently,
in order to account for all possible trajectories that go through the points in the set {xi },
we perform the following integration,

P [A → B] =

...

exp

i

tB
tA

˙
L(xi (t), xi (t), t)dt dx1 dx2 . . . dxn ,

(A.3)

where P [A → B] is the total probability of the system making a transition between the
given points, and integration is done over the entire range of possible values at points xk
of the configuration space. Now, we take the limit
67

→ 0 so that the time-slicing becomes

infinitesimally short, yielding Feynman’s formulation of the path integral,
B

P [A → B] =

exp
A

where Dx =

tB

i

tA

∞
i=1 dxi .

68

˙
L(x(t), x(t), t)dt Dx,

(A.4)

BIBLIOGRAPHY

69

BIBLIOGRAPHY
[1] L. Billings, M. I. Dykman, and I. B. Schwartz. Thermally Activated Switching in the
Presence of non-Gaussian Noise. Phys. Rev. E., 78:051122, 2008.
[2] H. B. Chan, M. I. Dykman, and C. Stambaugh. Paths of Fluctuation Induced Switching.
Phys Rev Lett, 100:130602, 2008.
[3] H. B. Chan, M. I. Dykman, and C. Stambaugh. Switching-Path Distribution in Multidimensional Systems. Physical Review E, 78:051109, 2008.
[4] H. B. Chan and C. Stambaugh, Phys. Rev. Lett. 99, 060601, 2007.
[5] B.E.DeMartini, J.F.Rhoads, K.L.Turner, J.Moehlis, S.W.Shaw, Linear and Nonlinear
Tuning of Parametrically Excited MEM Oscillator, Journal of Microelectromechanical
Systems 16, 310-318.
[6] B.E.DeMartini, J.Moehlis, K.L.Turner, J.F.Rhoads, S.W.Shaw, and W.Zhang, Modeling
of Parametrically Excited Microelectromechanical Oscillator Dynamics with Application
to Filtering, Paper PID134530, IEEE Sensors 2005, 4th IEEE Conference on Sensors,
Irvine, CA, 2005.
[7] R. C. Dorf, R. H. Bishop. Modern Control Systems. Prentice Hall, 12th edition, 2010
[8] M. I. Dykman and M. A. Krivoglaz. Classical theory of nonlinear oscillators interacting
with a medium. Phys. Stat. Sol. (b), 48:497-512, 1971.
[9] M. I. Dykman and M. A. Krivoglaz. Theory of Nonlinear Oscillator Interacting with a
Medium. Sov. Sci. Rev. A Phys., 5:265-442, 1984.
[10] Richard P. Feynman and Albert R. Hibbs. Quantum Mechanics and Path Inte- grals:
Emended Edition (Dover Books on Physics). Dover Publications, 2010.
70

[11] H. Goldstein, C. Poole, and J. Safko. Classical Mechanics, Third Edition. AddisonWesley, 2002.
[12] M. Gurjar and N. Jalili. Toward Ultrasmall Mass Detection Using Adaptive Self-Sensing
Piezoelectrically Driven Microcantilevers. IEEE/ASME Transactions on mechatronics ,
vol. 12, no. 6, December 2007.
[13] H. van Heeren, P. Salomon. MEMS: Recent Developments, Future Directions. Loughborough University, 2007.
[14] J. Guckenheimer and P. Holmes. Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields. Springer-Verlag New York, 1983.
[15] K. Jacobs. Stochastic Processes for Physicists: Understanding Noisy Systems. Cambridge University Press, 2010
[16] Hassan K. Khalil. Nonlinear Systems. Prentice Hall, 3rd edition, 2002.
[17] H. A. Kramers. Brownian motion in a field of force and the diffusion model of chemical
reactions. Physica VII, 7(4):284-304, 1940.
[18] Le Masne, Q., Pothier, H., Birge, N. O., Urbina, C., and Esteve, D., 2009. Asymmetric
Noise Probed with a Josephson Junction. Phys. Rev. Lett., 102, January, p. 067002.
[19] A. Leon-Garcia. Probability and Random Processes for Electrical Engineering. AddisonWesley, 1993.
[20] N.J.Miller, Noise in Nonlinear Micro-Resonators, Ph.D. Thesis, Michigan State University, 2012.
[21] A. H. Nayfeh, D. T. Mook. Nonlinear Vibrations. John Wiley & Sons, 1995.
[22] J. F. Rhoads and S. W. Shaw. The Impact of Nonlinearity on Degenerate Parametric
Amplifiers. Applied Physics Letters, 96(23): 234101, 2010.
[23] H. Risken. The Fokker-Planck Equation:
Springer, 2nd edition, 1996.

Methods of Solution and Applications.

[24] D.Ryvkine and M.I.Dykman. Resonant symmetry lifting in a parametrically modulated
oscillator. Phys. Rev. E 74, 061118 (2006)
71

[25] 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.
[26] W. Zhang and K. L. Turner. A Mass Sensor Based on Parametric Resonance. Proceedings Solis State Sensor, Actuator and Microsystems Workshop, Hilton Head Island,
South Carolina, June 6-10, 2004, pp.49-52.

72