EXPLOITING INTERNAL RESONANCE IN MEMS FOR SIGNAL PROCESSING
APPLICATIONS
By
Brian Scott Strachan

A DISSERTATION
Submitted
to Michigan State University
in partial fulfillment of the requirements
for the degree of
Mechanical Engineering – Doctor of Philosophy
Electrical Engineering – Doctor of Philosophy
2017

ABSTRACT
EXPLOITING INTERNAL RESONANCE IN MEMS FOR SIGNAL PROCESSING
APPLICATIONS
By
Brian Scott Strachan
This research focuses on the development and analysis of predictive models for frequency
converters and frequency generators that are based on micro-electro-mechanical-system
(MEMS) technology. In contrast to applications in which nonlinearity is sought to be
avoided, frequency conversion and frequency generation necessarily involve nonlinear
processes, and while many existing technologies are available for realizing these operations, MEMS technology offers a potentially advantageous combination of size, power
requirements, and noise characteristics.
This dissertation describes a series of investigations related to MEMS frequency conversion and generation, including: (i) an analytical investigation of a class of passive
multi-stage frequency dividers, (ii) the design and realization of this behavior in a MEMS
device, (iii) the development of a model for nonlinear modal interactions in closed loop
MEMS and (iv) the development of a computational method for optimizing their nonlinear resonant response through shape optimization. Items (ii) and (iii) were carried out in
close collaboration with experimental groups at the University of California at Santa Barbara and Argonne National Labs, respectively. Item (iv) was carried out in collaboration
with the topology optimization group at the Technical University of Denmark.
The subharmonic frequency divider is based on a class of mechanical structures with
nonlinearly coupled high Q vibration modes with sequential 2:1 internal resonances, for
which sequential parametric resonances are used to transfer energy from a high frequency
mode down to lower frequency modes. We analyze the normal form for this subharmonic
resonance cascade and predict the system response based on system and driving signal

parameters. We then show how to design and experimentally implement this subharmonic cascade in MEMS, and we demonstrate frequency division by a factor of eight.
The frequency generator model is based on a closed loop oscillator in which the resonator element has vibration modes with 1:3 frequency ratio and nonlinear intermodal
coupling. Experimental observations have shown that the oscillator phase noise performance is significantly improved when operating in a coupled mode regime, in which a
flexural mode is nonlinearly coupled to a torsional mode. The device is characterized by
comparing its measured open loop response against a model based on 1:3 internal resonance, demonstrating good agreement. The closed loop version of the model is analyzed
with a focus on how noise sources are filtered through the system into phase noise. This
model predicts the significant drop in phase noise observed when operating with internal
resonance. This predictive model provides a basis for future designs that take full advantage of this nonlinear behavior, which has potential for commercialization in the growing
area of MEMS oscillators.
Lastly, we describe the development of a computational tool that allows one to tailor
the nonlinear resonant response of mechanical structures using a combination of normal
forms and structural optimization tools. This approach is used to improve a device’s
nonlinear modal coupling by nearly an order of magnitude. Such tools will be important
for the continuing development of MEMS that utilize nonlinear resonant behavior.
In summary, it is shown that internal resonance, in addition to offering interesting
dynamic behavior, can be used to improve the performance of signal processing devices.
This work also demonstrates that devices that use internal resonance can be analyzed
with generic dynamic models, thereby providing a basis for understanding fundamental
device characteristics and future design development.

Copyright by
BRIAN SCOTT STRACHAN
2017

For Austin, again.

v

TABLE OF CONTENTS
LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii
LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
CHAPTER 1 INTRODUCTION . . . . . . . . . . . . . . . . .
1.1 Frequency Dividers . . . . . . . . . . . . . . . . . . . .
1.2 Frequency Generators . . . . . . . . . . . . . . . . . . .
1.3 Internally Resonant Systems . . . . . . . . . . . . . . .
1.4 Summary of Topics Considered and Results Obtained

ix

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

1
1
3
6
7

CHAPTER 2 SUBHARMONIC CASCADE FREQUENCY DIVIDER
2.1 The Model and The Averaged Equations . . . . . . . . . . . .
2.2 The Semi-Infinite and Long Chains . . . . . . . . . . . . . . .
2.3 Activating the Cascade . . . . . . . . . . . . . . . . . . . . . .
2.4 Period Doubling Cascade Accumulation . . . . . . . . . . . .
2.5 Numerical Simulations . . . . . . . . . . . . . . . . . . . . . .
2.6 Other Phenomena . . . . . . . . . . . . . . . . . . . . . . . . .
2.6.1 Multistability . . . . . . . . . . . . . . . . . . . . . . .
2.6.2 Small Forcing and Small Detuning . . . . . . . . . . .
2.7 Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

9
11
15
17
19
21
25
25
26
27

EXPERIMEN. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .

.
.
.
.

.
.
.
.

29
29
33
38

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

40
42
42
46
46
49
50
52
55
56
59

CHAPTER 3
3.1
3.2
3.3

A MEMS FREQUENCY DIVIDER:
TAL RESULTS . . . . . . . . . . . .
Designs . . . . . . . . . . . . . . . . . . . .
Experimental Results . . . . . . . . . . . .
Outlook . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.

.
.
.
.
.

DESIGN AND
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .

.
.
.
.
.

CHAPTER 4 COUPLED MODE FREQUENCY GENERATOR . . . . . . .
4.1 The Open Loop Analysis . . . . . . . . . . . . . . . . . . . . . . . . .
4.1.1 The Open Loop Model . . . . . . . . . . . . . . . . . . . . . .
4.1.2 The Open Loop Frequency Sweeps . . . . . . . . . . . . . . .
4.1.3 The Open Loop Bifurcation Diagram . . . . . . . . . . . . . .
4.2 The Closed Loop Analysis . . . . . . . . . . . . . . . . . . . . . . . .
4.2.1 Overview of Oscillators with Nonlinear Resonator Elements
4.2.2 The Closed Loop Oscillator Model . . . . . . . . . . . . . . .
4.2.3 The Deterministic Operating Frequency . . . . . . . . . . . .
4.2.4 The Oscillator Frequency Fluctuations . . . . . . . . . . . . .
4.3 Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

CHAPTER 5 DESIGN OPTIMIZATION FOR NONLINEARITY . . . . . . . . . . . 64
5.1 Characterization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
5.2 Optimization and sensitivity analysis . . . . . . . . . . . . . . . . . . . . . . . 67
vi

5.3
5.4

Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.3.1 Single-mode resonator . . . . . . . . . . . . . . . . .
5.3.2 Coupled-Mode Resonator with Internal Resonance
Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

CHAPTER 6

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

71
71
77
82

CONCLUSION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

APPENDIX . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95

vii

LIST OF TABLES
Table 3.1

Device parameters: beam lengths, designed and measured resonance
frequencies, and measured quality factors (in vacuum). Direct excitation of each mode is achieved through the actuation beam shown in
Figure 3.3, with signal Vdc = 90 V and Vac = 3−10 V. . . . . . . . . . . . . 33

Table 4.1

Model coefficients as measured from the open loop experimental results.

viii

45

LIST OF FIGURES
Figure 1.1

Schematic diagram of a frequency generator showing the essential
components. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 2.1

A mechanical implementation of the subharmonic cascade, consisting of rigid bars with elastic hinges at their bases, and coupled by
springs with linear stiffnesses. The energy is down-converted from
the high frequency beam, parametrically driven by a source at approximately twice its natural frequency, down the chain to a frequency
of Ω/2 N at the terminal beam. For small stiffnesses of the coupling
springs the bar displacements ui are roughly equal to the system
modal coordinates qi , that is, the system modes are localized in the
individual beams. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

Figure 2.2

Analytical predictions for the response regimes in (α, F ) space; the
destabilization boundaries αi and Fi for i = 1, 2 and ∞ are shown. The
fully active regime is shaded grey and the partially activated regime
+ corresponds to the stable
is shaded light grey. The curve labeled req
nontrivial response for the infinite chain. Parameter values for all
simulations are as follows, unless specified otherwise: ζ = 0.03, γ =
0.075, δ = 0.064, and β = 0.008. . . . . . . . . . . . . . . . . . . . . . . . . 20

Figure 2.3

Simulation of a six resonator chain, showing the sequential activation of the six elements when started with small initial conditions.
Resonators 1-6 are shown from top to bottom. The thick lines indicate results of the amplitudes obtained by simulating the averaged
equations, and the underlying fast oscillations are from simulations
of the full equations of motion. The settling time of resonator j is
proportional to Q/2 j−1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

Figure 2.4

Steady state amplitudes of a six resonator cascade for α = 0(> α∞ )
for various forcing amplitudes. The fully activated response is achieved
− , which is depicted as the grey shaded region. The
when F > F∞ = req
top two examples (circles and triangles) converge to the infinite lat+ , which is denoted by the row of asterisks, with the
tice amplitude req
amplitude of the final resonator given by a slightly lower value, r N
. The partially activated solutions occur for F1 < F < F∞ denoted
by the light grey region; the diamonds show a sample response. The
trivial response occurs for F < F1 and is shown by the unshaded
region; the squares show a sample response. . . . . . . . . . . . . . . . . 24

ix

4

Figure 2.5

Frequency response at F = 0.6 for six resonators, computed from the
averaged equations. The resonator amplitudes are the thick curves
whose darkness indicates place along the chain, that is, the first resonator is light grey and the 6th resonator is black. The vertical lines
are the predicted existence region boundaries, and the equal amplitude solution for the infinite chain is shown as a black dashed line.
The light dashed line is the activation amplitude for r j given in terms
of r j−1 , given in Eqn (2.11); where it crosses the r j−1 response line r j
emerges from the trivial response. . . . . . . . . . . . . . . . . . . . . . . 25

Figure 2.6

Activation boundaries in (α, F ) space computed from the averaged
equations. The black lines correspond to the activation boundaries
determined from simulations of the original equations of motion for
a six resonator cascade. The red dotted lines correspond to the analytical approximations for the boundaries of the first two resonators and
the infinite lattice. In general, there is good agreement except near the
bottom of the Arnold tongue where modes interact in a complicated
manner. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

Figure 2.7

Activation boundaries with small forcing in the limit of zero damping. While the case with zero backcoupling, β = 0, has the infinite
chain solution as a subset for every other region, the case with backcoupling, β = 0, suggests an aphysical scenario in which the entire
chain is activated before second resonator activates. At very low forcing (F < 0.1), this incorrectly suggests that the infinite chain will activate before the first resonator activates, indicating a breakdown of
the applicability of the infinite chain results. . . . . . . . . . . . . . . . . 27

Figure 3.1

Schematic model of a mechanical implementation of a multi-stage
frequency divider with four localized modes and coupling spring elements [53]. Driving element x5 parametrically excites element x4 ,
which parametrically excites element x3 , and so on down to element
x1 , with frequency division by 2 at each stage. Resonant response
amplitudes are dictated by system nonlinearities. . . . . . . . . . . . . . 30

Figure 3.2

(a) A three stage “candy cane” frequency divider cascade and its linear mode shapes. (b) A four stage “T-bar” frequency divider cascade
and its linear mode shapes. For each device, note that the beam ends
are fixed and that the modes of interest have resonant frequency ratios of 2:1; mode shapes are computed using COMSOL Multiphysics R .

Figure 3.3

31

SEM micrograph of the device, which spans 400 ¯m x 400 ¯m, with
uniform depth of 10 ¯m and feature size of 1.85 ¯m. . . . . . . . . . . . . . 34

x

Figure 3.4

The first three vibrational modes from a COMSOL Multiphysics R finite
element model. Corresponding modal frequencies are given in Table 3.1. 34

Figure 3.5

Normalized modal parametric instability boundaries for modes 1, 2,
and 3 with frequencies normalized by modal number n. Divide-byeight operation is achieved in the shaded overlap region. . . . . . . . . . 36

Figure 3.6

Measured amplitudes extracted from the spectrum at 21 , 14 , and 18 of
the drive frequency at a drive amplitude of 30 Vac . Displacement
values, taken from the envelope of the large amplitude solution, are
measured using the LDV, with the laser directed near the point of
maximum displacement of each beam. Due to modal coupling, the
mode 3 (first to activate) amplitude drops when modes 1 and 2 are active. 37

Figure 3.7

Measured response for each mode (beam) for 824.6 kHz drive signal at 34 Vac , in the fully cascaded operating mode. Displacement
values, taken from the envelope of the large amplitude solution, are
measured using the LDV. . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

Figure 3.8

SEM micrographs of two alternative frequency divider designs: the
Tbar (left) and ellipse-coupler (right) designs. The orange block represents the drive electrode. The first three modes are concentrated in
the red, green, and aqua colored beams, respectively. . . . . . . . . . . . 39

Figure 4.1

The clamped-clamped flexural-torsional modal coupling device. The
fundamental flexural mode has a resonant frequency near 65kHz and
the fundamental torsional mode has a resonant frequency near 198kHz.
This device is driven and sensed through an interdigitated comb drive
mechanism [2]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

Figure 4.2

Experimental (blue) and model (red) results of the open loop frequency sweep response for various drive amplitudes (20mV, 60mV,
200mV, 600mV), noting that the feedthrough capacitance modifies
the experimental amplitude results for larger forcing levels. For the
model results the dark red curve segments indicate a stable response
and the light red segments indicate an unstable response. . . . . . . . . . 47

Figure 4.3

Predicted saddle node (green curves, SNB) and Hopf bifurcations
(red curves, HB) as a function of drive frequency and drive amplitude, indicating bifurcations from both stable (dark) and unstable
(light) response branches. The experimentally measured SNB and HB
are denoted by (x) and (o), respectively. The bottom plot is a blowup
of the upper plot in the region of interest. . . . . . . . . . . . . . . . . . . 60

Figure 4.4

The feedback loop circuit producing a self-sustained oscillation. . . . . . 61

xi

Figure 4.5

The mean (noise-free) operating frequency as a function of feedback
amplitude. (Left) The main graph shows predictions from the noisefree model. Left inset shows corresponding experimental results where
red indicates sweeping upward in the amplitude and black indicates
sweeping downward. The right inset shows a blow up of the frequency as a function of amplitude in the internally resonance regime
for different levels of mode two cubic (Duffing) nonlinearity; the green
point indicates zero slope, that is, a point of zero dispersion. Model
parameters for this figure and subsequent figures are ω2 /ω1 ≈ 3.067, Q1 =
104 , Q2 = 104 , γ1 = 0.01, β = 0.00015. . . . . . . . . . . . . . . . . . . . 62

Figure 4.6

This figure shows two measures of the loop phase noise predicted
from the model as the forcing amplitude is varied. The three traces
from gray to black represent different levels of the mode two cubic
nonlinearity, as in Figure 4.5. The upper plot is a measure of the ratio of the contribution of the total (modes one and two) amplitude
noise to those of the total phase noise, and the lower plot is measure
of the ratio of the contributions of mode two to those of mode one.
We see that in all three cases there is a significant drop in the contributions of the amplitude noise and the overall mode one noise when
in the internal resonance operating regime. The inset shows experimentally measure standard deviation of the frequency fluctuations as
the amplitude is swept up (red) and down (black), showing the drop
in frequency fluctuations when operating in the internal resonance
regime. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

Figure 5.1

Initial (uniform width) design and the first two linear vibration modes
obtained using COMSOL Multiphysics modal analysis. Left: linear
vibration mode 1, right: linear vibration mode 2, with ω2 ≈ 2ω1 .
Color indicates vibration amplitude. . . . . . . . . . . . . . . . . . . . . . 68

Figure 5.2

Optimized design for maximizing the absolute value of the essential
modal coupling coefficient and the two coupled linear vibration modes
obtained using COMSOL modal analysis. Left: linear vibration mode
1, right: linear vibration mode 2 and ω2 ≈ 2ω1 . The color in online
version indicates the vibration amplitude. . . . . . . . . . . . . . . . . . . 68

Figure 5.3

Flowchart of the structural optimization routine for nonlinear dynamic response. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

Figure 5.4

Initial (uniform width) design of a clamped-clamped beam and its
linear vibration mode. Color indicates vibration amplitude. . . . . . . . 72

xii

Figure 5.5

Optimized design for maximizing the cubic nonlinearity of a clampedclamped beam and its linear vibration mode. Color indicates vibration amplitude. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73

Figure 5.6

Optimized design for minimizing the cubic nonlinearity of a clampedclamped beam and its linear vibration mode. Color indicates vibration amplitude. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73

Figure 5.7

Evolution of the objective function and shapes encountered during
the maximization process. The vertical axis is the absolute value of
the objective function divided by its initial value and the horizontal
axis is the iteration number. . . . . . . . . . . . . . . . . . . . . . . . . . . 74

Figure 5.8

Evolution of the objective function and shapes encountered during
the minimization process. The vertical axis is the absolute value of
the objective function divided by its initial value and the horizontal
axis is the iteration number. . . . . . . . . . . . . . . . . . . . . . . . . . . 75

Figure 5.9

Evolution of the objective function and shapes encountered during
the optimization process. The vertical axis is the absolute value of
the objective function divided by its initial value and the horizontal
axis is the iteration number. . . . . . . . . . . . . . . . . . . . . . . . . . . 80

Figure 5.10 Left: Typical time simulation with added envelope curves showing
beats, for the final design with initial first mode amplitude of 7.5 ×
10−11 (using eigenvectors normalized by the mass matrix). Right:
Data points and fitted exponential curves for the beat period, normalized by the period of the first linear mode, versus the initial amplitude of the first mode, for designs corresponding to iteration numbers 0 (initial design, top curve), 30 (middle curve), and 250 (final
design, bottom curve). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
Figure 5.11 (a) Evolution of the eigenfrequencies of linear vibration modes 1 and
2 encountered during the optimization process. (b) Evolution of the
pumping ratio and the mode localization ratio encountered during
the optimization process. . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
Figure .1

(a) The mechanical device along with the two mechanical modes and
the eigenfrequencies as computed in COMSOL. (b) The layout with
the red region being the mechanical device, the dark green indicating the drive and sense electrodes for mode 1, and the light green
electrodes allow for modal tuning and coupling bias. . . . . . . . . . . . 94

xiii

CHAPTER 1
INTRODUCTION
In the realm of signal processing, nonlinear behavior is generally viewed as unwanted
since it produces spurious harmonics and undesirable system characteristics. As such,
much work has been done to operate in the linear regime or in a regime that does not
elicit nonlinear phenomena. There are, however, useful devices that operate in a fundamentally nonlinear manner, and the successful design of these devices requires the designer to understand and utilize the nonlinearity in an advantageous way. These include
devices such as photonic lasers, superconducting Josephson junctions, electronic diodes,
and mechanical nonlinear vibration absorbers [70, 30, 57, 52]. In the field of signal processing, there are two such components that are found in almost every electronic device:
frequency dividers and clocks.
The aims of the present work are to develop and examine these type of devices based
on internally resonant modal interactions. Below, we provide an overview of frequency
dividers, clocks, internally resonant dynamical systems, with MEMS as the technology
for the device implementation.

1.1

Frequency Dividers

Frequency converters are essential elements in modern portable electronic devices. Frequency down-conversion and up-conversion, respectively known as frequency division
and multiplication, are important to many processes including power conversion, clock
distribution, frequency synthesis, radio frequency instrumentation, and optical link frontends. The performance of converters is measured by factors that include size, cost, frequency coverage, frequency resolution, power consumption, spurious harmonic content
(including phase and amplitude noise performance), and stability [3]. Frequency con-

1

version can be passive or active, and can be achieved in electrical, optical, or mechanical
domains, or combinations thereof, and the strategy chosen for frequency conversion ultimately depends on performance specifications, power requirements, size, and cost.
In the context of portable devices, electronic frequency converters are the most prevalent due to their cost, tunability, integrability, and wide operating range. There are several types of electronic dividers. Digital dividers are useful for low frequencies due to
their wide band of operation but are limited at higher frequencies due to energy cost
and that they are limited by the analog-to-digital converter clock signal. Resonator based
electronic dividers, such as microstrip or cavity based, are narrowband devices useful
for higher frequencies (> 10Ghz), but not for lower frequencies due to size constraints.
The successful high frequency analog topologies are injection locked, regenerative frequency dividers, and parametric dividers [4, 58, 9] Injection locking, on the one hand,
requires two self-sustained oscillation circuits with commensurate resonant frequencies.
A weak coupling between the two resonators will, just like Hugyen’s clocks, synchronize
the clocks if the coupling to noise ratio is large enough [73, 50]. Regenerative dividers,
on the other hand, place an over damped or under damped resonator in a feedback loop
with an amplifier in the feedback and a mixer that mixes the input and feedback signal.
The limitation of these devices is the quality factor and amplifier speed, although these
devices have been shown to work in the low Terahertz range [64]. Parametric dividers
are generally passive devices that are extremely low power consuming devices and also
have the best phase noise characteristics since they do not have amplifiers. The parametric divider topology has also been utilized with amplifiers in order to increase their
operating regime [23, 24].
For both second order regenerative dividers and parametric dividers, a nontrivial
steady state operation occurs when a resonator with natural frequency ωo has a reactive element (stiffness or capacitance) modulated periodically with a frequency Ω ≈ 2ωo ,
and the system responds with frequency Ω/2 ≈ ωo . The archetypal equation govern2

ing this behavior, when the pump is the only input, is the nonlinear Mathieu equation,
˙
x¨ + ωo x/Q
+ ωo2 (1 + δ cos(Ωt)) x + γx3 = 0, in which Q represents dissipation and γx3
captures nonlinear effects. The threshold for instability, in terms of the modulation amplitude δ, depends on Q and the nearness of the excitation frequency to twice the natural
frequency, as measured by a detuning parameter α = (Ω − 2ωo )/(2ωo ) [23, 40]. The
effective Q factor can also be increased based on the amplifier (if there is one) [51].
Although electronic parametric frequency dividers are useful due to their power savings, their size is a prohibitive dilemma for small devices. MEMS resonators have a
much smaller footprint and a much higher Q than electronic resonators at the same frequency, and parametric resonance has already been demonstrated in a variety of MEMS
devices [63, 47, 11]. To date, with the exception of a patent [34], the author is unaware
of any effort to utilize a MEMS 2:1 parametric resonance device for frequency division,
although there has been a general effort to replace many timing application electronic
devices with MEMS devices [43]. Frequency dividers are generally used in timing distribution circuits. Since these dividers are slaved to the clock, we turn our attention to the
subject of the frequency generator.

1.2

Frequency Generators

Frequency generators are devices that produce a periodic signal that is used for time
keeping. A clock is just another name for a frequency generator. The dynamics of these
systems can be represented by mathematical models possessing a limit cycle of a given
period. The precision of these systems is affected by the inevitable fluctuations of the period of oscillations. Since these devices are not slaved to a reference timing system, they
have time translational symmetry, so that any source of noise will cause drift of the phase,
resulting in phase noise. The goal of clock designers is to identify and remove any systematic drift which usually entails understanding the source of the noise and eliminating

3

Resonator
Limiter

Amplifier

Phase Shifter
Figure 1.1 Schematic diagram of a frequency generator showing the essential
components.
the clock sensitivity to the noise. Different types of clocks have different phase noise characteristics. For example, an atomic clock, which is often a quartz resonator that is coupled
to the resonant frequency of atomic transitions, has a very good long term time stability
but the short term time stability is not as good as a purely quartz based clock. A majority of popular single mode clocks can be modeled by the schematic shown in Figure 1.1,
where a resonator element is in a closed loop with an amplifier that replaces the energy
dissipated by other loop elements, a phase shifter, and an amplitude limiter. The amplifier and phase shifter destabilize the system, resulting in oscillation whose amplitude is
fixed by the limiter. The resulting nominal frequency depends on all the loop elements
parameters and can be tuned by varying them. A more generic model commonly used
lumps all elements in the feedback loop as an impedance Z = ZR + jZ I . In this case, the
amplifier and amplitude limiter affect ZR and the phase shifter affects Z I [66]. This linearization impedance approach works well when the resonator element can be regarded
as linear. When a mechanical resonator is used, its time-dependent displacement is typically represented by a motional series RLC circuit parallel to a feedthrough capacitor C,
that is, the effective capacitance that occurs between voltages across the input and output leads of the transducer. The response is then a second order impedance, where the
applied voltage is converted to a current [25]. Usually, the system is linearized about an
operating point and analyzed from there. Generally, clocks are made to be as linear as
possible and are made to utilize a single second order filter as the frequency selective el-

4

ement. Due to the linearity constraint, they are engineered to be physically large enough
in order to have a large enough signal to noise ratio when transducing the mechanical
signal to an electronic signal. Furthermore, other “parasitic” modes of the oscillator are
engineered to be far away from the primary mode and all of its harmonics, to keep it
as isolated as possible. As a result, the resonator models considered in these circuits are
primarily single mode and linear.
Integration of a resonator into a circuit topology to create a self-sustaining feedback
loop requires specific types of circuits [35]. Due to the fact that increasing the number
of elements in the circuit will produce additional noise, it is desirable to minimize the
number of elements in the device. Thus, a popular topology used with a mechanical
(quartz) resonator is the Pierce oscillator circuit, since it uses only one active element,
necessary for counteracting the energy loss in the circuit. There are additional sources of
noise in these systems, such as thermo-mechanical noise, internal temperature gradients,
and ambient environment fluctuations, and Nyquist-Johnson noise [5, 17]. In this work,
we do not consider the details of these noise sources, but rather how to optimize the
system based on assumed noise profiles. That is, we take a design engineering point of
view.
An important goal in this field is to make extremely small, high frequency resonators
with good phase noise characteristics [1]. Since the desired size is small, since frequency
is inversely related to the size of the device, the devices may need to be driven into their
nonlinear regime in order to overcome measurement noise, and generally to achieve good
signal to noise characteristics. Generally, a distinctive feature of operation of an oscillator
at an amplitude where the resonator element behaves nonlinearly is that amplitude noise
will contribute to the overall phase noise, due to the amplitude-frequency dependence of
nonlinear resonators; this does not happen in the linear range of the resonator [13]. As
a result, one of the objectives for engineering a generator based on a nonlinear resonator
element is to make it as near to a linear system as possible, in terms of phase noise, which
5

is done by eliminating this amplitude to phase (AF) noise coupling. This amplitude noise
contribution to the phase noise can be eliminated for certain resonators by operating at
a so-called “sweet spot,” that is, where ∂I/∂ω ≈ 0 (a point of zero dispersion), thus locally mimicking with a linear oscillator [36, 18]. Another approach is to adjust the phase
shifter so that the amplifier noise is orthogonal to the vector that produces phase diffusion [71, 26]. All of these MEMS clock models are based on a single mode response of the
resonator element, but recent experimental work has shown that the phase noise performance of a given device can be significantly improved when the resonator responds with
two nonlinearly coupled modes. Predictive modeling of these systems requires analysis
of internally resonant systems, to which we now turn.

1.3

Internally Resonant Systems

Internal resonance occurs in systems when nonlinear coupling terms become strong given
an alignment of the modal resonant frequencies. Internal resonances in mechanical systems have been the subject of many investigations over the past fifty years, and some
internally resonant MEMS have been proposed and built in recent years. The nonlinear
coupling in these systems is necessarily passive and can thus be derived from a potential.
For two modes, the following potential displays a variety of different nonlinear coupling
terms
1
1
V = ω12 x12 + ω22 x22 + β 1 x1 x22 + β 2 x13 x2 + β 3 x12 x22 + . . .
2
2

(1.1)

The β 1 coupling term provides the so-called autoparametric resonance phenomena and is
important in the 1:2 case, and the β 2 term provides the essential modal coupling in the 1:3
case. The 1:1 case is considerably more difficult since all of the terms are relevant. The β 3
term is the so-called “diffusive coupling” term, which should be acknowledged as possibly being important in all of the interactions. Fortunately, standard perturbation methods
provide a systematic means of determining which coefficients are essential for different
6

resonance conditions [39], and the theory of normal forms gives this a mathematically
rigorous foundation [20]. Certain mechanical structures, such as cables and flexural and
torsional beams and plates, exhibit internal resonance and have been explored with nonlinear coupling between their various modes [27, 32, 48, 38, 41, 7]. Detailed investigations
of the response of systems composed of many nearly-identical resonators with parametric excitation, that is, with 1:1:· · · :1 internal resonances have been considered with both
nearest-neighbor [31] and global [10, 72, 37] coupling. Several micro mechanical devices
have been designed for a 2:1 internal resonance [68, 65] and 3:1 internal resonance [2].

1.4

Summary of Topics Considered and Results Obtained

In this work we discuss two main topics. The first topic is an investigation into a novel
class of frequency dividers, namely, a cascade of nonlinearly coupled micro-mechanical
resonators. This cascade is analyzed in Chapters 2 and 3. Chapter 2 considers an idealized
model for the cascade, based on the normal form for sequential 1:2 internally resonant
elements. Perturbation methods are used to analyze the system and determine viable
operating points in the driving signal parameter space and how these depend on system
parameters. In Chapter 3, we describe how one realizes such a cascade in mechanical
frames using coupled beam elements, suitable for MEMS implementation. We provide
simple designs rules for these systems and show different design configurations. Experimental results, carried out in collaboration with our UCSB partners, and demonstrating
a frequency division of eight, are described for one of these designs.
The second topic is an investigation of the class of frequency generators that utilize
internally resonant modal interactions to improve phase noise performance. This work
is based on close collaboration with an experimental group at Argonne National Labs. In
chapter 4, we develop a generic model for this class of problems, including small noise
sources, and develop a method for predicting how these noises are filtered through the

7

nonlinear system response. We demonstrate that certain operating points, in particular
those where the modal coupling becomes strong, can result in significant phase noise
reduction. We develop a detailed model based on an Argonne device with 1:3 internal
resonance, and consider its open and closed loop response. The model is shown to capture the experimentally measured responses of the device, and it is used to predict the
limits of phase noise reduction using internal resonance.
In Chapter 5 we describe computational tools that allow one to determine the normal
form coefficients for a given nonlinear resonance from a nonlinear finite element model.
These tools are used with an optimization algorithm to develop structures with tailored
nonlinear response, specifically, either maximized or minimized normal form coefficients.
This work forms the basis on which future generations of devices will be developed.
The dissertation closes in Chapter 6 with a discussion of the results, their significance,
and directions for future work in this area.

8

CHAPTER 2
SUBHARMONIC CASCADE FREQUENCY DIVIDER
In this chapter, we propose a novel frequency conversion method based on nonlinear
dynamics that exploits the robustness of parametric resonance. The primary approach
is to use a subharmonic resonance cascade, which involves a chain of internally resonant
subsystems of a form that allows energy exchange between the subsystems. Although detailed investigations of the response of systems composed of nearly-identical resonators
with parametric excitation have been considered with both nearest-neighbor [31] and
global [10] coupling, the resonance structure considered here has only been analyzed
for a small number of modes (up to three) [42]. As discussed in the introduction, this
investigation is motivated by providing an alternative to electronic frequency dividers
that are used in many devices [16, 22]. Most of the work presented here is taken from the
following published papers [54]. In the following chapter we consider realizations of this
cascade in MEMS hardware.
In order to physically intuit the dynamics of interest, consider the system shown in
Figure 2.1, consisting of N hinged beams (resonators) coupled by springs with linear stiffness. The first beam is parametrically excited and, for a range of driving frequencies and
amplitudes, the first resonator will experience parametric resonance; its response will
parametrically excite the next resonator through the coupling spring and so on down the
cascade. The inter-beam springs provide nonlinear coupling of the desired form via finite deformation kinematics. Note that the coupling is two way, since the j + 1 resonator,
when it is active, necessarily back-acts on resonator j, thereby providing two-way coupling, and this coupling is resonant in both directions.
In this chapter, we first present the generic equations of motion that model such a
subharmonic resonant cascade system, from which we derive averaged equations that
approximate the slowly varying amplitudes and phases of the system elements. These
9

u1

u0

u3

u2

Figure 2.1 A mechanical implementation of the subharmonic cascade, consisting of rigid
bars with elastic hinges at their bases, and coupled by springs with linear stiffnesses.
The energy is down-converted from the high frequency beam, parametrically driven by
a source at approximately twice its natural frequency, down the chain to a frequency of
Ω/2 N at the terminal beam. For small stiffnesses of the coupling springs the bar
displacements ui are roughly equal to the system modal coordinates qi , that is, the
system modes are localized in the individual beams.
equations are essentially the normal form for this type of resonance. For the present
study we consider special conditions on the parameters that provide a simple form that
is amenable for analysis and offers some explicit results. We start the analysis by considering a semi-infinite chain, considering a special solution in which all resonators are
active, helping set the stage for further analysis. This approach also offers approximate
closed form expressions for the steady state response of most elements of the semi-infinite
chain, which are shown to be very useful approximations for many elements in chains of
finite length. We then examine the conditions on the forcing amplitude and frequency for
nontrivial steady state operation of the first few resonators, specifically computing how
elements of the chain become activated as the input parameters are varied. Using the
results from these analyses, we delineate the regions in excitation parameter space over

10

which the system is expected to produce a complete subharmonic cascade. This parameter space is then given more resolution using numerical simulations, particularly in the
small forcing amplitude case, and other interesting phenomena are discussed. The analytical approximations are confirmed using simulations of a sample system. The chapter
closes with some thoughts about applications and extensions of these results.

2.1

The Model and The Averaged Equations

For the generic cascade model we start with a Hamiltonian (per unit mass) for a system
of N resonators, described in terms of vibrational modal amplitudes q j and conjugate
momenta p j , each with cubic nonlinear stiffness and quadratic nearest-neighbor coupling,
specifically given by,
N

H=

1 2 1 2 2 1 4
p j + ω j q j + νj q j + κ j q j q2j+1 ,
2
2
4
j =0

∑

(2.1)

where ω j and νj are, respectively, the linear natural frequency and mass-normalized cubic
stiffness coefficient for the jth resonator, and the κ j ’s are the mass-normalized nonlinear
coupling coefficients. The resonator nonlinearity will limit the amplitudes of vibration
when undergoing subharmonic resonance, and the coupling nonlinearity is required for
energy transfer between resonators. Note that many nonlinear coupling terms may be
present, but only one is present in the normal form, that is, only one is needed to capture
the energy transfer and the others can be removed by coordinate transformations [20].
In order to account for dissipation, we introduce resonator viscous damping factors, ζ j
directly into the modal differential equations, via generalized forces. We also assume
that nonlinear effects and damping are relatively small, consistent with many MEMS resonators [48]. The boundary conditions on the resonator chain consist of: (i) a source with
q0 = F sin (Ω0 t + φ0 ), which parametrically excites the first resonator with frequency

11

Ω0 = 2ω1 (1 + α), where α describes small detuning from resonance; 1 and (ii) a grounded
terminal condition with q N +1 ≡ 0.
The cascade occurs when each resonator j is parametrically excited by resonator j −
1 via the coupling, which requires that the resonators be tuned such that ω j−1 /ω j ≈

2. Conditions for complete activation of all N resonators must account for the effects
of damping, detuning, and quadratic backaction from resonator j + 1. The mechanical
system depicted in Figure 2.1 offers the desired features, and is modeled by the given
Hamiltonian.
To simplify the model we rescale time using the natural frequency of the first resonator, τ = ω1 t, and the displacement according to x j = (

νj /ω j )q j , thus normalizing

the cubic nonlinearities and the resonators’ amplitudes. Linear frequencies of adjacent
resonators are set to be close to 2:1 ratios, as required for the cascade to occur. For the
present model we make some additional assumptions that simplify the analysis, by making the non-dimensional resonator equations have identical form. It is important to note
that these assumptions are taken for convenience in the analysis, but are not required for
the desired cascade to occur. These additional assumptions are: cubic nonlinearity and
quadratic coupling coefficients for adjacent resonators are taken to be close to 4:1 ratios,
that is,
νj
νj − 1

=

κj
κ j −1

=4

(2.2)

and the modal damping values are nearly the same, that is, ζ j ≈ ζ ∀ j. 2 These assumptions render the nondimensional parameters for each resonator to be conveniently
1 Note

that one can use direct excitation of the first resonator, but the present form is
taken for convenience.
2 Note that the 4:1 ratio for resonator cubic nonlinearity holds for fixed-fixed beams
with identical uniform cross-sections where the lengths are varied to account for the 2:1
change of modal frequency.

12

identical, resulting in equations of motion,
x¨ j + 22− j ζ x˙ j + 22(1− j) x j = f j

(2.3)

8
f j = 22(1− j) −σj x j − γx3j − 4δx j−1 x j − 4βx2j+1
3
for j = 0, · · · , N + 1
where the non-dimensional parameter definitions are: nonlinearity, γ = 3/8; forward
coupling, δ = 2κ j /(

νj ω j ); backward coupling, β, where for passive coupling β = δ/2;

and detuning α = (Ω j − ω j )/ω j , where Ω j = Ω0 /2 j . Note that σj << 1 represents deviations from the 2:1 natural frequency ratio conditions. We further simplify the model for
our present investigation by taking σj = 0 ∀ j = 1, that is, we make the resonators perfectly tuned relative to one another. This assumption is also not required for the cascade
to occur, but it simplifies the analysis which, for present purposes, focuses on the general
features of the cascade. Of course, σj = 0 will need to be incorporated in system design
studies to account for variabilities arising from fabrication tolerances and other effects.
Note that this approach does allow for tuning of the input amplitude F, via x0 , and input
frequency Ω0 , via α.
Due to the nature of the excitation and the internal resonances, xi = 0 ∀i satisfies the
equations of motion, but its stability depends on the input and system parameters. The
response of interest will be one that is dynamically stable and has all resonators operating
with non-zero amplitudes. To investigate all possible responses and their stability we turn
to a perturbation analysis.
To approximate the system response, we use a transformation to rotating coordinates, the van der Pol transformation, by defining amplitude and phase coordinates via
x j (t) = r j (t) cos(Ω j t + φj (t)) and x˙j (t) = −r j (t)Ω j sin(Ω j t + φj (t)). Applying a first order perturbation expansion about the harmonic solution, averaging over the period of
the lowest frequency resonator, and ignoring higher order terms, we obtain equations approximating the slowly varying amplitude and phase of each resonator [40]. Under the
13

stated assumptions and transformations, the averaged equations for the slowly varying
amplitudes and phases are given by,
r˙ j = 21− j −ζr j − δr j−1 r j sin Θ j + βr2j+1 sin Θ j+1
r j φ˙ j = 21− j [−αr j + δr j−1 r j cos Θ j

(2.4)

+ βr2j+1 cos Θ j+1 + γr3j ]

where Θ j = φj−1 − 2φj . The array end conditions are r0 = F, φ0 = 0, and r N +1 =

0. Under the given assumptions about the parameters, the averaged equations become
identical for each resonator, except for a factor of 21− j for resonator j. This implies that the
time scale of resonator j is slower than that of j − 1 by a factor of two, resulting in slower
transient behavior as one moves down the chain. This is expected, since the resonator Q

values are the same and the frequencies of adjacent resonators scale by a factor of 1/2.
To find the steady state solutions, we set r˙ j = φ˙ j = 0, and solve for the resulting
amplitudes and phases. The stability of the steady-state responses can then be determined
from the averaged equations from the Jacobian in the usual manner. However, these
equations do not conveniently uncouple, and no closed form solution is possible, except
for some special responses (for example, the trivial response). However, in the case with
no back coupling, β = 0, one can sequentially solve the equations in closed form. While
this case is not of primary interest, it is considered in Section 2.4 below, since it allows for
a calculation of interest.
We begin by considering a special case that helps set the stage for a more full investigation, namely, a special type of fully activated response for an infinite length chain.
This is followed by an analysis that considers sequential activation of the chain as input
parameters are varied.

14

2.2

The Semi-Infinite and Long Chains

Instead of solving the full system of 2N equations, which requires numerics, we consider
a special response wherein each of the resonators in a semi-infinite chain has the same
amplitude. It is shown below that this solution provides a useful approximation for the
response of many of the central elements in finite length chains. By setting r j = req for
all j, and assuming a semi-infinite lattice of resonators, we can solve for the steady state
amplitudes from the averaged equations. Regarding the phases, it can be shown that
sin Θ j = sin Θ j+1 and cos Θ j = cos Θ j+1 for an equal amplitude solution, yielding 3
possible steady-state amplitudes: the trivial solution and
± =
req

A±

α2
A2 − 2 −
γ

( β + δ)ζ
γ( β − δ)

2

,

(2.5)

where A = [(δ + β)2 + 2αγ]/(2γ2 ). It is notable that these solutions are independent
± appear in a saddle-node bifurcation
of F since we assumed F = r0 = req . Solutions req
+ and an unstable response for r − .
as parameter vary, resulting in a stable response for req
eq

The trivial response is always stable. While these equal amplitude responses have explicit
solutions, due to the symmetry of the solutions there are 2 N possible phase solutions for a
lattice of N resonators. The phases for the steady state solutions must satisfy Θ j − Θ j+1 =
2πn for n ∈ Z, recalling that Θ j = φj−1 − 2φj .

It will be seen that for finite length chains, the steady-state amplitude for the bulk of

the resonators, not including the first few or the last (j = N) resonator, is approximately
+ . Whether the bulk of the chain is activated or quiescent depends
equal to either zero or req
− , as follows: If F > r − , F < r − , or F = r − , then the
on the magnitude of F relative to req
eq
eq
eq
+ , r ≈ r0 , or r ≈ r − , respectively, where the r − case is
bulk of the resonators r j ≈ req
eq
eq
eq
j
j

unstable.
Returning to the amplitude solutions, we can analytically derive the conditions for
± exist and are real. In particular, we are interested
which the two nontrivial solutions req

15

± exist. By considering their form it is found that r ±
in the region of (α, F ) for which req
eq

exist when α ≥ α∞ , where the bifurcation value of α is given by,
α∞ =

γζ 2
( β + δ )2
−
.
4γ
( β − γ )2

(2.6)

− , the entire cascade becomes active.
We can then apply the assertion that when F ≥ req

In later sections we will utilize the notation that F∞ is the critical forcing that makes all
− . This
the resonators destabilize from the 0 solution, where we have shown that F∞ = Feq

provides an approximate activation condition for the excitation amplitude that depends
on the other system parameters.
For long chains, the behavior of the terminal resonator N is completely determined
by the response of the N − 1 resonator, since it is parametrically driven by the N − 1 resonator from one side and the other side is grounded, r N +1 = 0. As such, r N is described
by the nonlinear Mathieu equation with excitation provided by r N −1 from the infinite
chain solutions, which has the trivial solution as well as
r±
N =

δr N −1 2 ζ 2
− 2.
γ
γ

α
±
γ

(2.7)

These solutions are identical to the solutions of the first resonator when r2 = 0 (given
in Equation 2.10 below), except that the effective forcing is r N −1 in this case, in place
+ as in Equation (2.7) and
of F. In order to approximate r N , we approximate r N −1 ≈ req

choose the stable solution in the desired parameter space, r +
N . The result is found to be

quite accurate for chains even with only a few elements, as demonstrated in the examples
simulated below.
Since low-power usage is of interest in applications, it is relevant to consider how a
semi-infinite chain can respond with equal amplitudes when each resonator is dissipative,
here modeled with viscous friction. For an entirely activated semi-infinite chain, we can
show that that the total power dissipated is finite, as follows. The calculation of the total

16

loss through dissipation is given by,
N

Pdiss =

1 T 2− j
∑ T 0 (2 ζ x˙ j )x˙ j dt .
j =1

(2.8)

+ , the
Assuming a semi-infinite chain where each element responds with amplitude req

result is found to be,
Pdiss =

2 2 + 2
ζΩ r
,
7 0 eq

(2.9)

− < F < r + , the actual
since the infinite sum is convergent. We would expect that for req
eq

power used would be slightly less than Pdiss since the first few resonators have ampli+ , and yet the chain is fully activated. If one’s goal is to limit the power
tudes less than req

consumption of the device, then one can either decrease the damping or decrease steady
+ . From Equation 2.5, the latter can be achieved by either decreasing
state amplitudes, req

the coupling or by increasing the cubic nonlinearity of the resonators.

2.3

Activating the Cascade

Here we derive conditions related to activation of the cascade as input parameters are varied. We begin by considering the system with parameter conditions such that the fully
trivial response is the only stable response, regardless of initial conditions. As the excitation amplitude or frequency is slowly varied, the first step in activation of the cascade occurs when the driven resonator (j = 1) experiences a subharmonic instability, specifically,
a period-doubling bifurcation of the system, which corresponds to a pitchfork-bifurcation
in the averaged equations. As parameters are varied further, the second resonator activates through another similar bifurcation. This occurs at a point at which the first resonator, which is driving the second, reaches a critical amplitude, and this scenario is then
repeated along the entire chain. For these activations, we append the additional subscript
to denote the “k-bifurcation”, that is, when rk becomes nonzero. As such, r j,k is the amplitude of the jth resonator at the k-bifurcation and αk , for example, is the value of the
17

detuning parameter at the k-bifurcation. In the desired fully activated response, the zero
solution for each resonator is unstable so that all N resonators are active. In fact, the
only way the entire chain can be activated is by sequential activation of the resonators,
since there is no mechanism for stimulating isolated internal elements in the chain. In
this section, we obtain expressions for the critical values of the excitation parameters that
correspond to this sequence of instabilities, and explicit expressions for the first two instabilities are derived.
We begin by considering a partially activated chain, for which ri = 0 for i ≤ j and
r = 0 for > j. Since resonator j + 1 is inactive, the response of resonator j is governed by
the simple nonlinear Mathieu equation for which the excitation is provided by resonator
j − 1. Its nontrivial steady-state response amplitude, where r j+1 = 0, can be found from
the averaged equations and is given by,

rj =

δr j−1 2

α
±
γ

γ

−

ζ2
.
γ2

(2.10)

From this expression we find a condition for the j-bifurcation by based on the amplitude
r j−1,j ,
r j−1,j =

ζ 2
α 2
+
.
δ
δ

(2.11)

Note that while useful, this expression does not give an explicit condition for the jbifurcation in terms of the excitation parameters. However, using this result we can obtain
such conditions for the first two bifurcations, as follows.
For the activation condition of the first resonator, we consider Equation (2.10) with
j = 1, recognizing that r0 = F is the excitation amplitude. One can solve r1 in terms of F,
which is the usual solution of the Mathieu equation. The point where r1 becomes nonzero
yields a condition that can be solved for a relationship between the critical values of the
excitation parameters (amplitude and frequency) and the system parameters. This initial
activation is the 1-bifurcation, and is simply that for a simple Mathieu equation. For the
18

next activation condition, the 2-bifurcation, we use the fact that r1 is known from Equation
(2.10) in closed form in terms of the input and system parameters, so long as r2 = 0. We
again use Equation (2.10), this time with j = 2, r2 = 0, and the known expression for
r1 to determine the 2-bifurcation condition in closed form. These calculations result in
the following expressions for the parameter conditions for the j-bifurcations for j = 1, 2,
which are expressed in terms of both the force amplitude and frequency,
α1 =
F1 =
α2 =
F2 =

(δF )2 − ζ 2

(2.12)

α 2
ζ 2
+
δ
δ

1
δ2
−
2γ 2γ

δ4 − 4γ2 ζ 2 + 4γ2 ζ 2

2
ζ 2
α γ
+
− ( F1 )2 .
δ
δ
δ

(2.13)

( α1 )2
γ2

(2.14)
(2.15)

These conditions describe boundaries in the (α, F ) input parameter space, shown as the
indicated solid lines in Figure 2.2, which delineate where each of the first two resonators
become active. Note that bifurcation conditions for j > 2 cannot be obtained explicitly
in this manner since no closed form expression is available for the r j ’s when two or more
resonators are active.

2.4

Period Doubling Cascade Accumulation

In Section 2.3 we found conditions for the critical input values for destabilizing the jth
resonator, which corresponds to a period-doubling bifurcation. Explicit expressions were
obtained for α j and Fj for j = 1, 2. Due to the fact that the infinite cascade is fully activated
at finite values of the input parameters, there must be an accumulation of period-doubling
bifurcations. Period-doubling cascades and the attendant accumulation ratios for the amplitudes and parameter values at bifurcation points are well understood for uni-modal
maps [6]. The universal nature of these sequences was described by Feigenbaum, who
19

1.4
1.2

(α1 ,F1 )

+
req

(α2 ,F2 )

Forcing

1.0
0.8

(α∞ , F∞ )

0.6
0.4
0.2
0.0
-0.10

-0.05

0.00
Detuning

0.05

0.10

Figure 2.2 Analytical predictions for the response regimes in (α, F ) space; the
destabilization boundaries αi and Fi for i = 1, 2 and ∞ are shown. The fully active
regime is shaded grey and the partially activated regime is shaded light grey. The curve
+ corresponds to the stable nontrivial response for the infinite chain. Parameter
labeled req
values for all simulations are as follows, unless specified otherwise: ζ = 0.03, γ = 0.075,
δ = 0.064, and β = 0.008.
computed values for the accumulation ratios [19]. It is important to note that the present
system differs from those satisfying Feigenbaum’s “universal” conditions, since here degrees of freedom are added at each instability, and the resulting response is not chaotic.
In all these cases, however, accumulation ratios of the following form can be determined,
a p = lim

j→∞

p j − p j −1
p j +1 − p j

(2.16)

where p represents some parameter value or response amplitude of the system.
For the present system we take p j as the square of the forcing parameter, F2 , and
formulate the accumulation ratio problem for the general cascade. By taking zero backcoupling, β = 0, we can obtain an explicit solution for the accumulation constant for F2
in terms of the system parameters. While not completely general, this special solution
demonstrates that the accumulation occurs and that the accumulation ratio is not universal, but rather depends on the system parameters. For β = 0 amplitude r j−1 depends on
20

r j , but not on r j+1 , and the steady state averaged equations give
r2j−1 =

ζ 2
α γ 2 2
+
− rj .
δ
δ
δ

(2.17)

The expressions for F1 and F2 are, even in this special case, given in Equations (2.13) and
Equation (2.15), respectively. The following pattern emerges emerges for all Fj ,
ζ 2
+
δ

( Fj )2 =

2 2
α γ
−
Fj−1
.
δ
δ

(2.18)

We can then find a F2 by writing Fj2 and Fj2+1 in terms of Fj2−1 . Since we are interested in
− , which is the known critical forcing
the accumulation ratio as j → ∞, we set Fj−1 = req

amplitude for destabilizing the entire cascade, thus yielding the accumulation ratio for F2
in the infinite limit,
1

a F2 =
1−

1+4

αγ
δ2

−

γ2 ζ 2
δ4

.

(2.19)

This validity of this result is confirmed by simulations of the system with β = 0.

2.5

Numerical Simulations

In this section, we perform a set of numerical simulations in order to validate and extend
the analysis. A set of N = 6 resonators is considered, which is sufficiently large for
the infinite lattice approximation to be useful. We take the resonator nonlinearity to be
hardening, that is, γ > 0. For a softening nonlinearity, the results are simply flipped in
terms of the sign of α.
First, the averaged equations and the original equations of motion are numerically
integrated and the results compared. As seen in Figure 2.3, the amplitudes from the
averaged equations, ri (t), match the transient and steady state response envelopes of the
original system of equations quite well in terms of decay times, transient beat frequencies,
and amplitudes, albeit with time shifts for the transient beating response of the resonator
21

1
0
−1
1
0
−1
1
0
−1
1
0
−1
1
0
−1

1
0
−1
0

5000

10000
Time (t)

15000

Figure 2.3 Simulation of a six resonator chain, showing the sequential activation of the
six elements when started with small initial conditions. Resonators 1-6 are shown from
top to bottom. The thick lines indicate results of the amplitudes obtained by simulating
the averaged equations, and the underlying fast oscillations are from simulations of the
full equations of motion. The settling time of resonator j is proportional to Q/2 j−1 .
furthest down the chain. These time shifts arise due to variations in the startup transients,
which are very sensitive to the initial conditions selected. The steady-state conditions are
very accurately captured.
This numerical simulation, as well as others carried out but not shown here, confirm
the use of the averaged equations as a basis for analytical predictions of the system response, and confirms our intuitive predictions about the system dynamics. It is important to note that the periodicity of each resonator is on the slowest timescale of the system,
namely, on the period of the last resonator in the chain. This is not visible at the resolution
shown, since the low frequency amplitude modulations, which arise from back coupling,
are quite subtle, and the steady-state response each resonator is dominated by Ω j .
Figure 2.4 shows the steady state amplitudes of all resonators, ri for i = 1 − 6, for
4 different values of the forcing strength F, for the case of zero detuning, α = 0. Re-

22

+ , the finite lattice is well approximated by the infinite lattice; this
call that when F = req

corresponds to the solution depicted by the row of asterisks in the figure. The critical am− = r − is depicted by the row of asterisks that differentiates the fully activated
plitude F∞
eq

(shaded) and partially activated (unshaded) regions. The resonant solutions described
by the circle and triangle symbols correspond to two cases with F ≥ F∞ , for which the
entire chain is activated; we refer to this as a fully activated response. The two solutions
represented by diamond and square symbols correspond to cases with F < F∞ . The solution with diamonds shows a partially activated chain, which occurs when the forcing is
strong enough to activate only the first two resonators, since F2 < F < F3 . For F slightly
lower, shown as squares in the Figure, not even the first resonator is activated, resulting
in a fully trivial response. The simulation results shown also confirm the assertion that
+ for all
for F ≥ F∞ , the steady state amplitude profile is well approximated by ri = req

but the first few and the final resonators. As discussed above, the amplitude of the final
resonator, i = N, indeed shows a slight decrease in amplitude from the infinite chain
amplitude, and matches its predicted value, given in Equation (2.7), very well.
The system frequency response, determined by sweeping the detuning parameter α
for fixed levels of forcing, is also of interest. As seen in Figure 2.5, as the frequency is
increased, the resonators sequentially activate, and the amplitude profile of many ele+
ments of the lattice, specifically, all except for the first few, are well approximated by req

when the chain is fully activated. Also, these responses are quite insensitive to the forcing
− ; recall that this follows since r + is independent of F. Note
strength, so long as F ≥ req
eq
√
+ has the form r + ≈ α − a + b where a and b are constants, which dictates how
that req
eq

the response amplitude depends on α, an observation confirmed by Figure 2.5. Also, it
is interesting to note that the critical values α j at which individual resonators activate accumulate onto α∞ as j becomes large. In fact, this accumulation must be geometric and,
although in Section 2.4 we only discussed the accumulation ratio for F2 for the β = 0 case,
due to its closed form availability, an accumulation ratio must exist for α as well, and for
23

Steady State Amplitude

1.5
rN
+
req
1

Full Activation Regime
−
F∞ , req
F2
F1

Partial Activation Regime
Trivial Regime

0
F

1

2

3
4
Resonator #

5

6

Figure 2.4 Steady state amplitudes of a six resonator cascade for α = 0(> α∞ ) for various
− , which
forcing amplitudes. The fully activated response is achieved when F > F∞ = req
is depicted as the grey shaded region. The top two examples (circles and triangles)
+ , which is denoted by the row of asterisks,
converge to the infinite lattice amplitude req
with the amplitude of the final resonator given by a slightly lower value, r N . The
partially activated solutions occur for F1 < F < F∞ denoted by the light grey region; the
diamonds show a sample response. The trivial response occurs for F < F1 and is shown
by the unshaded region; the squares show a sample response.
the more interesting case with backcoupling.
Lastly, Figure 2.6 shows numerically calculated (in solid black) activation boundaries
for each of the six resonators, along with the analytical predictions (in dashed red) for the
first two resonators and for the infinite cascade . Note that this figure focuses on the lower
central portion of the region shown in Figure 2.2, since this is where the most interesting
features occur. As expected, the analytical solutions correctly predict the first and second
resonator activation boundaries. The figure also shows that for all but a small range
of forcing amplitudes the infinite lattice activation boundary is a good predictor for the
finite lattice. However, due to secondary bifurcations, discussed below, the infinite lattice
approximation becomes invalid near the bottom of the shaded region. Other features of
this Figure as also described below.
24

2
α1
Amplitude

1.5

α2 α∞

−α1

+
req

Multistable

1
−
req

0.5
0
-0.06

r1
-0.04

r2 r3..6
-0.02

0
0.02
Detuning

0.04

0.06

Figure 2.5 Frequency response at F = 0.6 for six resonators, computed from the averaged
equations. The resonator amplitudes are the thick curves whose darkness indicates place
along the chain, that is, the first resonator is light grey and the 6th resonator is black. The
vertical lines are the predicted existence region boundaries, and the equal amplitude
solution for the infinite chain is shown as a black dashed line. The light dashed line is
the activation amplitude for r j given in terms of r j−1 , given in Eqn (2.11); where it
crosses the r j−1 response line r j emerges from the trivial response.

2.6

Other Phenomena

In this section, we consider two other types of phenomena that occur in this system:
multistability, and the behaviors that arise near the bottom of the shaded region in Figure
2.6. We do not pursue an investigation of these issues in this thesis, since they are not
relevant to the applications we have in mind, although they are no doubt interesting.

2.6.1

Multistability

So far, we have focused on the activation solution from the near zero initial conditions
and have defined this region to be localized in the (α, F ) parameter space as a generalized
type of Arnold tongue. However, as one would expect from the single resonator case,
there does exist a response regime for α > −α1 (for the hardening case), where multiple

+ solution (which approximates
stable solutions occur, namely, the zero solution, the req

the fully activated chain), and partially activated solutions. Simulated frequency sweeps

25

0.8

Forcing

0.7

(α1 ,F1 )

(α∞ ,F∞ )

(α2 ,F2 )
(α3...6 ,F3...6 )

0.6

Multistable

0.5
0.4
0.3
0.06

0.04

0.02

0
0.02
Detuning

0.04

0.06

Figure 2.6 Activation boundaries in (α, F ) space computed from the averaged equations.
The black lines correspond to the activation boundaries determined from simulations of
the original equations of motion for a six resonator cascade. The red dotted lines
correspond to the analytical approximations for the boundaries of the first two
resonators and the infinite lattice. In general, there is good agreement except near the
bottom of the Arnold tongue where modes interact in a complicated manner.
have confirmed the existence of the fully active solution and the zero solution in the multistable regime.

2.6.2

Small Forcing and Small Detuning

As indicated in Figure 2.6, there is a region where the infinite lattice solution breaks down
as a predictor for the activation of the entire chain. The main dynamics observed in this
region involve quasi-periodic energy transfer between the first two resonators, indicating
the presence of a secondary bifurcation not captured by the present analysis. In these responses, simulations show that the first resonator approaches an amplitude large enough
to excite the second resonator, but when it transfers its energy to the second resonator, the
first resonator amplitude decays, and does so to an amplitude where it no longer excites
the second resonator, at which point it begins to gain amplitude, starting the cycle over
again. This energy exchange phenomena does not occur for zero back-action (β = 0).
In addition to the energy exchange phenomena, Figure 2.7 shows that in the limit of
26

0.4
(α∞ ,F∞ ),β = 0
(α∞ ,F∞ ),β = 0

Forcing

0.3

0.2

(α2 ,F2 )
(α1 ,F1 )

0.1
0.0
-0.020

-0.015

-0.010
-0.005
Detuning

0.000

0.005

Figure 2.7 Activation boundaries with small forcing in the limit of zero damping. While
the case with zero backcoupling, β = 0, has the infinite chain solution as a subset for
every other region, the case with backcoupling, β = 0, suggests an aphysical scenario in
which the entire chain is activated before second resonator activates. At very low forcing
(F < 0.1), this incorrectly suggests that the infinite chain will activate before the first
resonator activates, indicating a breakdown of the applicability of the infinite chain
results.
small damping, ζ = 0, the (α∞ , F∞ ) curve, corresponding to full activation, intersects with
both the (α1 , F1 ) and (α2 , F2 ) curves near the bottom of the Arnold tongue. This presents
an impossible situation, since the first and second resonators must destabilize before the
rest of the cascade activates. This observation indicates a breakdown of the applicability
of the infinite chain results. This may be an important region in the parameter space
to consider for a given physical implementation of a cascade, since the effective forcing
strength depends on many design variables and is often relatively small.

2.7

Outlook

The results of this investigation demonstrate that cascades of subharmonic resonances
will occur in systems designed with the appropriate forms of internal tuning and cou27

pling. The infinite chain approximation was shown to be very useful for predicting the
fully activated response for chains of finite lengths, even as small as the six element example considered. It was also found that the analytical perturbation results are very effective
in predicting the transient and steady-state response of these systems, thereby providing
a useful tool for investigating system parameter dependence, which will be a significant
aid for device design. Also, while it was found that the response of this system can be
quite complicated, especially near the bottom of the Arnold tongue, there exist large open
sets of system and input parameter values over which the system behaves quite simply,
specifically, in the desired mode of a robust frequency dividing cascade.
Ongoing work on this class of systems includes the following topics: robustness to
parameter variations such as the effects of inter-resonator mistuning; design realization
of MEMS devices with the desired properties; the possibility of using this arrangement
for frequency up-conversion; and issues related to practical performance measures of
subharmonic cascades, including phase noise characteristics and power requirements.

28

CHAPTER 3
A MEMS FREQUENCY DIVIDER: DESIGN AND EXPERIMENTAL RESULTS
As discussed in the Chapter 2, the frequency divider cascade combines the benefits of
cascading, internal resonance, and mechanical coupling in a single micro-device. The operation is based on nonlinear modal coupling and exploits the robustness of parametric
resonance. This narrow-band approach uses a subharmonic resonance cascade in a chain
of internally resonant subsystems with specific coupling that allows energy exchange between successive divide-by-two stages. In this chapter we develop realistic devices based
on the generic model of a mechanical implementation of such a system consisting of four
stages, with the frequencies of consecutive elements (vibration modes) having 2:1 ratios,
as shown in Figure 3.1.
This work was done in close collaboration with Kamala Qalandar and Kim Turner
of the University of California at Santa Barbara. The author was primarily responsible
for the MEMS preliminary design and simulation, while the UCSB group fabricated and
set up the testing of the device. The experimental testing of the device was done jointly.
Much of the second half of this chapter is verbatim from the publication [46].

3.1

Designs

The design objectives for the subharmonic cascade are: a single mechanical structure, localized modes each with a modal nonlinearity and high Q, the desired modal coupling,
simple modal frequency tuning, and capacitive based parametric driving of the highest
mode in the cascade. The single mechanical structure is for simplicity and to guarantee the passive coupling between the modes. Due to the large number of modes in the
subharmonic cascade, in order for a conceptual design (as opposed to, for example, an
optimization based design), we need localized linear modes. Thin flexural-based beams

29

x4
x5
x2

x3

x1

Figure 3.1 Schematic model of a mechanical implementation of a multi-stage frequency
divider with four localized modes and coupling spring elements [53]. Driving element
x5 parametrically excites element x4 , which parametrically excites element x3 , and so on
down to element x1 , with frequency division by 2 at each stage. Resonant response
amplitudes are dictated by system nonlinearities.
that each can represent a mode will provide the mode localization. The modal nonlinearity is necessary for amplitude saturation. The high Q is necessary for several reasons: a
large response is usually necessary for a MEMS based readout in order to overcome the
measurement noise floor and the effect of noise is weaker in a large amplitude solution,
and, since we do not want to drive the system hard, we need to be able to destabilize
the trivial solution for small forcing levels. Modal coupling nonlinearities come about for
many reasons, and they often will have a different form depending on the geometry of

30

Mode 1
Mode 1

Mode 2

Mode 2
Mode 3

Mode 3

Mode 4

(b)

(a)

Figure 3.2 (a) A three stage “candy cane” frequency divider cascade and its linear mode
shapes. (b) A four stage “T-bar” frequency divider cascade and its linear mode shapes.
For each device, note that the beam ends are fixed and that the modes of interest have
resonant frequency ratios of 2:1; mode shapes are computed using COMSOL
Multiphysics R .
the device; in this case, due to the desired 2:1 parametric coupling, we need to ensure
that the dominant modal coupling nonlinearity is of the quadratic form. In addition, due
to fabrication tolerances, the proper alignment of the resonant frequencies will likely not
occur in the nominal equilibrium state of the device so the resonators need to be tunable.
The tuning must be small so that the integrity of the desired geometry which produces
the desired nonlinearity is not affected. Finally, the mode localization allows for driving
the first mode and being nearly tangent to all other modes. In this section we review
two designs and discuss, qualitatively, their operation. We then look at the eigensystem,
namely the linear modes, as calculated in COMSOL Multiphysics R software. We then
simulate the excitation force in COMSOL Multiphysics R and discuss the results.
In practice, the schematic in Figure 3.1 is not feasible in a MEMS device for high frequency operation. Instead, more realistic implementations of the device are shown in
Figure 3.2. Since the modes are localized on a particular beam, we will discuss the dynamics in terms of beam modes.
As seen in Figure 3.2 (a), the "candy cane" cascade is composed of the several beams

31

whose lengths increase by factors near

√

2 since the resonant frequencies need to be in

a 2:1 ratio. It is clear that the eigenmodes are relatively localized, although there are
small deflections in the connecting beams, which are semi-circular springs. The coupling
springs are relatively weak so that the modes are spatially localized in the beams, but the
coupling is also sufficiently strong to provide parametric pumping from one stage to the
next. The eigenmodes have a mode shape similar to a beam with one fixed end and one
end with a spring. The jth mode is excited by the j − 1 beam through a modulation of the
jth beam’s length. This modulation of a length is the source of the parametric excitation.
Since the coupling is passive, there is also back coupling present, which implies that a
lengthening of the jth happens twice per j mode oscillation and each lengthening has the
same effect on the j − 1 beam. As shown in 3.2 (b), the "Tbar" cascade is also composed
√
of several beams whose length increase by factors near 2. In contrast to the candy cane
model, the Tbar divider is composed of modes that are nearly symmetric. Except for the
first beam, the end that is connected to the previous mode’ main beam is free to rotate,
with a grounded spring attached. Due to the ability to rotate the jth mode contains some
small anti-symmetric deflections in the j − 1 beam. Even so, the individual modes can be
approximated as fixed-pinned beams, since the jth mode is still concentrated on the jth
beam. For each of these designs, each main modal beam will undergo midline stretching
for large deflections, which results in a hardening cubic nonlinearity for each mode.
At this point, the devices satisfy the three of the desired criteria: a single mechanical
structure, localized modes with a modal nonlinearity, and the desired modal coupling.
The other criteria are satisfied via the surrounding electronic components, fabrication,
and packaging, which are standard features of MEMS devices. In this design, we use
electrostatic fields in both devices for driving, tuning, and readout. The devices are driven
by exciting the first mode of the structure. In the candy cane cascade, the first beam can
be driven directly at the frequency near ω1 , with a fluctuating potential across a small gap
along the length of the beam. The first beam can also be driven parametrically, by quasi32

statically exciting the 0th beam near the frequency 2ω1 . The Tbar divider can only be
excited directly, in this case, although the design can be slightly altered to also be driven
parametrically.

3.2

Experimental Results

We realize the "candy cane" design, as shown in the SEM micrograph in Figure 3.3, using
a single-crystal silicon structure fabricated with deep reactive ion etching (DRIE) and
released with HF vapor. The microbeam elements are uniformly 10 ¯m in depth, 1.85 ¯m
wide, and the beam lengths Ln are determined to satisfy the required 2:1 frequency ratios.
Table 3.1 shows measured device parameters and characteristics. Note that only modes
1-3, which are localized in the response of beams 1-3, of the device shown in Figure 3.3
are considered since, in this particular device they exhibit the necessary modal frequency
ratio alignments. As such, beam 5 is directly actuated and it directly drives beam 4 which
acts as the parametric driver to beam 3, the highest frequency mode in this cascade. Finite
element modeling, shown in Figure 3.4, is used to refine the basic device design and
provide the desired modal frequencies and ratios.

Mode 1
Mode 2
Mode 3

Beam
Length
[µm]
376
263
183

Design
Frequency
[kHz]
105
210
420

Measured
Frequency
[kHz]
103
206
413

Quality
Factor
Q
858
1035
1315

Table 3.1 Device parameters: beam lengths, designed and measured resonance
frequencies, and measured quality factors (in vacuum). Direct excitation of each mode is
achieved through the actuation beam shown in Figure 3.3, with signal Vdc = 90 V and
Vac = 3−10 V.
Actuation for both characterization and operation is provided through the beam shown
at the top left in Figure 3.3, which is driven using capacitive plate actuation with an applied voltage Vin = Vdc + Vac sin(ωd t). The electrostatic force f e is approximated by an
33

Drive
Electrode
Actuation
Beam

Mode 2 Response Beam

Mode 1 Response Beam

Mode 3 Response Beam

Figure 3.3 SEM micrograph of the device, which spans 400 ¯m x 400 ¯m, with uniform
depth of 10 ¯m and feature size of 1.85 ¯m.

Mode 3

Mode 2

Mode 1

Figure 3.4 The first three vibrational modes from a COMSOL Multiphysics R finite
element model. Corresponding modal frequencies are given in Table 3.1.

34

effective parallel plate capacitance and assuming that Vac << Vdc , the force is approximately linear in Vac . The velocity and phase of the in-plane motion of a point on each
beam are detected using a Polytec MSA-400 Laser Doppler Vibrometer (LDV). To direct
the laser beam perpendicular to the in-plane motion, 45 ◦ angled mirrors are etched into
anchored regions adjacent to the beams using FIB milling[62]. Modal parameters are
characterized by sweeping ωd for different levels of Vac and measuring the response of
individual beams, which requires a separate run for each mode due to the use of a single
point laser. The semi-circular coupling beams allow direct excitation of each mode individually through the actuation beam. The modal frequencies and quality factors listed in
Table 3.1 are measured at a pressure of 350 mTorr with drive signal of Vdc = 90 V and
Vac = 3−10 V.
In operation, a measurable response of mode n is observed when the amplitude of
the drive signal is sufficiently large and the drive frequency is sufficiently close to 2ωn .
This induces a parametric resonant response in the corresponding beam. The regions in
which these responses occur are mapped in the (Vac , ωd ) parameter space by sweeping
these input parameters for the three modes of interest (Figure 3.5). The instability regions
observed in Figure 3.5 are not the usual single mode parametric resonance regions due to
the fact that the modal responses include the dynamics of other modes. Specifically, mode
n is parametrically driven by mode n+1, which, due to the complicated geometry, is in
turn driven by both parametric and direct excitation from the drive signal. The various
mode n instability wedges are dictated by these subtle effects. For example, modes 3 and
1 are parametrically driven by modes 4 and 2, respectively. Since both beams 4 and 2
are perpendicular to the drive beam, the drive signal will have a dominant parametric
force and a much smaller direct forcing term on beams 4 and 2. In contrast, mode 2 is
parametrically driven by mode 3, which has a parallel orientation to the drive beam and,
as a result, the drive signal has a much larger direct force on beam 3 than either beams 4
or 2. As such, the effective parametric forcing on modes 3 and 1 is considerably smaller
35

than that on mode 2 which is why the second instability region is much larger than the
other two instability regions.

40

Drive Amplitude VAC [V]

Mode 3
30

20
Mode 1
Mode 2

10

0
0.996

0.998

1

1.002

1.004

Normalized Drive Frequency
Figure 3.5 Normalized modal parametric instability boundaries for modes 1, 2, and 3
with frequencies normalized by modal number n. Divide-by-eight operation is achieved
in the shaded overlap region.

To determine the input conditions for which the entire cascade is activated, the required drive frequencies of the individual modal activation regions are normalized by
the modal frequency and modal number n. The entire cascade activates when the drive
amplitude Vac and frequency are in the region delineated by the intersection of the three
individual modal activation regions, shown shaded in Figure 3.5, resulting in frequency
division by 2, 4 and 8 in sequential beam elements.
36

Frequency of Concurrent Large Amplitude
Responses of Modes 1, 2, and 3 [kHz]
fin /2

Measured Displacement Amplitude [µm]

0.4

412.25

413.25

414.25
Inside
Mode 3
Wedge

0.2
0

824
fin /4

825

826

206.25

827

828

206.75

829
207.25

0.6
0.3
Inside Mode 2 Wedge
0
824

825

fin /8
103.05
0.4
0.2
0
824

826

827

828

829

827

828

829

103.25
Inside
Mode 1
Wedge
825

826

Drive Frequency of Vac [kHz]

fin

Figure 3.6 Measured amplitudes extracted from the spectrum at 12 , 14 , and 81 of the drive
frequency at a drive amplitude of 30 Vac . Displacement values, taken from the envelope
of the large amplitude solution, are measured using the LDV, with the laser directed near
the point of maximum displacement of each beam. Due to modal coupling, the mode 3
(first to activate) amplitude drops when modes 1 and 2 are active.

37

Figure 3.6 shows the measured vibration amplitudes from beams 1-3 as the drive frequency is swept up across the range of interest. As the drive frequency crosses the left
boundary of the instability region for mode n, a jump to a finite amplitude response of
that mode (beam) is observed; conversely, the amplitude drops to the stable zero solution
as the drive frequency crosses the right boundary. There is no observed bi-stability in this
sweep. The response amplitudes are determined by modal frequency detuning mechanisms, such as mechanical and electrostatic nonlinearities captured by parameters γn in
the model, as well as by the modal interactions described by the κn terms. Note that the
observed drop in the mode 3 amplitude between 825 kHz and 826 kHz results from the
activation of mode 1 in that range. Figure 3.7 shows the scaled FFT of the responses of the
three beams in full cascade operation when the system is driven through the actuation
beam with a signal of 34 Vac at 824.6 kHz, which is well inside the full cascade activation
region. In this response, vibration amplitudes reach up to 70% of the gap size.

3.3

Outlook

The results show demonstration of a multi-stage micro-mechanical frequency divider
which takes an input signal and provides division by 2, 4 and 8 in different output elements. Ongoing work includes testing of packaged devices, experimental characterization of system nonlinearities, development of devices optimized for nonlinear response,
and consideration of phase noise and power usage in comparison with more traditional
devices. The phase noise of all passive, open loop, divider topologies will track the input
frequency at very low frequencies, with a drop of 6dB per divide-by-two stage. A variety
of designs that exhibit this type of frequency division in the UHF range are also in development, which are expected to have reduced power requirements when compared with
digital-based dividers operating in that range. A SEM micrograph of a couple of these
devices are shown in Figure 3.8.

38

0.8

Displacement [µm]

Mode 2
0.6
Mode 1
0.4
Mode 3
0.2

Drive
Signal

0
0

200

400
600
Frequency [kHz]

800

Figure 3.7 Measured response for each mode (beam) for 824.6 kHz drive signal at 34 Vac ,
in the fully cascaded operating mode. Displacement values, taken from the envelope of
the large amplitude solution, are measured using the LDV.

Figure 3.8 SEM micrographs of two alternative frequency divider designs: the Tbar (left)
and ellipse-coupler (right) designs. The orange block represents the drive electrode. The
first three modes are concentrated in the red, green, and aqua colored beams,
respectively.

39

CHAPTER 4
COUPLED MODE FREQUENCY GENERATOR
In this chapter, our ultimate goal is to understand and describe how various noise sources
contribute to phase noise in closed loop oscillators when the electromechanical frequency
selective element operates with internal resonance, specifically, with two nonlinearly coupled vibration modes. The closed loop produces a self-sustained oscillation, which can
be used as a frequency generator, that is, a clock. The performance of a clock is related to
measure of the frequency fluctuations during operation. Note that this is also related to
the phase fluctuations or the phase noise, but we will use the frequency fluctuations as the
main matric for clock performance in this chapter. The preferred operation of such oscillators has the resonator element maintaining vibration in its linear regime, since nonlinear
operation converts amplitude noise into frequency noise [18]. However, aggressive size
and power requirements require smaller resonator elements, which necessitates consideration of nonlinear resonator effects. There have been two recent experimental demonstrations of reducing frequency fluctuations in oscillators based on nonlinear resonator
behavior, one using a so-called parametric feedback oscillator in which the output has its
frequency doubled and pumped back into the loop [67], and a device in which nonlinear
activation of a secondary mode through a 1:3 internal resonance results in significant drop
in frequency fluctuations [2]. In both cases, the main effect appears to be that the oscillator
improves, in terms of frequency fluctuations, when the resonator element is coupled to a
signal that is off from the main oscillator frequency. Although much recent work has been
directed towards understanding how single mode nonlinear resonator frequency generators work [36, 28], coupled mode oscillators and the frequency fluctuation reduction
observed in them have not been adequately theoretically explained and understood.
The contribution of the present work focuses on modeling and predictive analysis,
done in close collaboration with David Czaplewski and Daniel Lopez from the Nanofab40

rication and Devices group at Argonne National Labs (ANL) who first observed internal resonance in a closed loop system [2], using the frequency selective resonator device
shown in Figure 4.1. Relevant experimental results from ANL are interspersed throughout the chapter to demonstrate the validity and utility of our modeling and analytical
efforts.
The analytical work is based on a normal form two-mode model that captures the experimental results observed in [2]. We begin with an analysis of the open loop response of
a system with a 1:3 frequency ratio for the two modes, corresponding to the fundamental
flexural and torsional modes in the ANL device. Open loop frequency sweep techniques
and bifurcation diagrams are used to develop a quantitative model with parameter values
that provide a good match with the experimental open response. Since the ANL device
follows closely the dynamics of the 1:3 internal resonance model, we have confidence in
the subsequent investigation of the closed loop system.
The closed loop model with the coupled mode resonator element provides a means of
predicting the spectrum and the statistics of frequency fluctuations. The analysis provides
a good qualitative match with experimental results, but a quantitative match is not possible without better models of the noise sources in the experimental device. The model also
takes into account the important effect of feedthrough capacitance. A key to understanding the system response is a nonlinear version of anti-crossing, where the crossing occurs
at the point of perfect tuning for the internal resonance. We conclude the chapter with a
discussion of how one might improve the frequency fluctuation performance in this class
of devices.

41

4.1
4.1.1

The Open Loop Analysis
The Open Loop Model

In this section, we provide a derivation and analysis of a two-mode model for the open
loop response of the MEMS device under consideration when the primary mode is subjected to an external harmonic drive. The mechanical device, shown in Figure 4.1, is a
clamped-clamped beam structure composed of three thin beams, and comb-finger paddles on each side that are used for transduction, namely, one for actuation and one for
sensing. The parallel beams are designed to reduce the internal strain which is associated
with thermo-elastic damping [33], which can cause additional noise terms. This device
acts as a relatively simple clamped-clamped beam with a central lumped mass (the comb
structures). This structure has many modes, but we are interested in the fundamental
flexural mode at 65kHz, the primary mode of interest, and the fundamental torsional
mode at slightly more than three times that frequency, at 198kHz. The stiffnesses of these
modes are dominated by an elastic potential, with a small contribution from the electrostatic fields from the transduction mechanisms. Thus, we model the two vibration modes
and their interaction using a potential, which will generally contains all possible cubic
and quartic modal coupling terms. However, only some of these terms will be resonant,
depending on the frequency ratio of the modal coupling, are they are determined by the
normal form [20]. Here, we observe a 1:3 frequency ratio, which suggests the following
for the interaction potential

Vint = αx12 x22 + βx13 x2

(4.1)

If one assumes symmetry in the electrostatic fields and the mechanical structure, the
modes are pure flexure and torsion, and the only mechanism for a 1 : 3 modal interaction
is through the αx12 x22 term, which provides only a very narrow 2:3 parametric resonance

42

Figure 4.1 The clamped-clamped flexural-torsional modal coupling device. The
fundamental flexural mode has a resonant frequency near 65kHz and the fundamental
torsional mode has a resonant frequency near 198kHz. This device is driven and sensed
through an interdigitated comb drive mechanism [2].
tongue [8]. A more viable model that matches the experiments includes a bias that creates asymmetry in the equations of motion, which can be represented by a simple linear
coupling of flexure and torsion. This force can originate physically from many sources
including the vertical asymmetry in the electrostatic transduction forces, from pretension
in the flexural beam producing a vertical force on the structure, from vertical asymmetries due to the crystal orientation, or from other physical asymmetries. The resulting
potential, with terms for the isolated modes and their coupling, has the form
2

V=

1

1

∑ 2 ωi2 q2i + 4 γi q4i

i =1

+ β q q31 q2

(4.2)

where the q1 are modal coordinates. These coordinates are dominated by flexural and
torsional motions of the device, but are combinations of the two, due to the asymmetry.
The open loop internally resonant response is governed by dynamic variables Si , i =
1, 2, which represent modal amplitudes scaled by a some complex transduction gain, that
is, Si = gi qi , where Si is the output signal and has units volts (V). From the Hamiltonian
whose potential is described by Equation (4.2), the model for the 1:3 internal resonance is

43

then given by
S¨1 +2Γ1 S˙ 1 + ω12 S1 + α1 S22 S1 + γ1 S13 + 3β 1 S12 S2 = F cos(Ωt)

(4.3)

S¨2 +2Γ2 S˙ 2 + ω22 S2 + α2 S12 S2 + γ2 S23 + β 2 S13 = 0

(4.4)

which includes coefficients for the individual modal damping, Γi , frequencies ωi , and
Duffing nonlinearities γi , along with coefficients for nonlinear dispersive coupling α1 and
α2 and modal coupling β i , and external drive of amplitude F and frequency Ω. We refer
to the forward and back coupling coefficients as β 2 and β 1 , respectively. While important
in some applications [44], dispersive coupling is not required in the present analysis, so
we take α1 = α2 = 0 in the subsequent development. The complex amplitudes associated
with the Si are more relevant to experimentally measured quantities, reflecting the amplitude and phase, or the quadratures, of the system response. As such, we utilize averaging
to obtain a more useful dynamic system.
The complex amplitude equations are found by expressing S1 = s1 exp(iΩt) + c.c.,
S2 = s2 exp(i3Ωt) + c.c.,S˙ 1 = iΩs1 exp(iΩt) + c.c., and S˙ 2 = i3Ωs2 exp(i3Ωt) + c.c.. Differentiating these and utilizing the constraints on (S1 , S˙1 ) and (S2 , S˙2 ), we obtain conditions S¨1 = −Ω2 S1 + 2iΩs˙ 1 exp(iΩt) and S¨2 = −9Ω2 S2 + i6Ωs˙ 2 exp(i3Ωt). Substitution
of these into Equations (4.3) and (4.4) and integrating over the timescale of the slowest
period of oscillation, T = 2π/Ω, we obtain the complex amplitude equations of interest:



α
i
3
iF
∗(2)
s˙ 1 = −Γ1 − i δω1 − 1 |s2 |2  s1 +
3β 1 s1 s2 + i
γ1 |s1 |2 s1 −
(4.5)
Ω
2Ω
2Ω
4Ω
s˙ 2 =

−Γ2 − i δω2 −

α2
| s1 |2
6Ω

β i
3
s2 + 2 s31 + i
γ2 |s2 |2 s2
6Ω
6Ω

(4.6)

where

(Ω2 − ω12 )
2Ω
(9Ω2 − ω22 )
δω2 =
6Ω
δω1 =

44

(4.7)
(4.8)

represent the frequency detuning of the drive relative to the fundamental mode and also
from the internal resonance condition. These model equations provide the basis for the
analysis of the open loop system. Periodic steady states in which the primary mode responds at frequency Ω and the secondary mode responds are 2Ω are given by the constant
solutions of Equations (4.5) and (4.6), denoted as (s¯1 , s¯2 ). The steady state amplitudes
and phases are determined in the usual manner from the complex amplitudes. Furthermore, the local stability of a steady state (s¯1 , s¯2 ) is determined in the usual manner using
linearization of Equations (4.5) and (4.6) about (s¯1 , s¯2 ). As parameters are varied, bifurcations of (s¯1 , s¯2 ) occur, and these correspond to qualitative changes in the steady-state
response. Details about the method of averaging for harmonically driven systems can be
found in many standard text on nonlinear vibrations, for example [20].
From the experimentally measured open loop response of the system, we utilize frequency sweep data and associated bifurcation diagrams in order to derive coefficients in
equations (4.5) and (4.6) that provide a match between the model and the measured response. The results are found in Table 4.1. Here we review the model predictions and
how they compare with the experimental results.
coeff
ω1
ω2
Q1
Q2
γ1
γ2
β1
β2
F1 (20mV )

values
64580 · 2π
3 · 66533 · 2π
7.2 · 105
4 · 105

units
s −1
s −1
-

2.45 · 1015

V −2 s −2

2 · 1016

V −2 s −2

0

3 · 1010

3.198 · 103

V −2 s −2
V −2 s −2
Vs−2

Table 4.1 Model coefficients as measured from the open loop experimental results.

45

4.1.2

The Open Loop Frequency Sweeps

In Figure 4.2, four upward frequency sweeps are shown, with different drive amplitudes,
with experimental data in blue and the model predictions in red. The experimental data
was taken by fixing the drive amplitude and sweeping the drive frequency, and recording
the mode 1 complex amplitude, s1 , as a function of these drive parameters. Mode 1 is
deeply into the non-linear regime before any modal interaction occurs. As a result, the
lowest drive amplitude (Vin = 20mV), provides a sweep where mode 1 is governed by
a single mode Duffing model. For modest drive amplitudes, Vin = 60mV and Vin =
200mV, mode 1 is affected by mode 2, saturating the upper frequency limit of the Duffing
(upper) branch response. For these two drive levels, the model predicts that the upper
branch solution continues after a small gap near the interaction frequency, although such
a gap cannot be found by one-way frequency sweeps. For the frequency range plotted, the
model predicts that the continued branch is unstable and stable for the drive amplitudes
Vin = 60mV and Vin = 200mV, respectively. Note that experimental observations can
be made only for stable responses, and while unstable responses exist and important to
the overall frequency response, they cannot be directly observed. Interestingly, for Vin =
600mV, the solution is able to jump the frequency gap and continue on the stable upper
branch solution. The reason for this is related to the structure of the solution stability, and
is better understood by looking at bifurcation diagram, shown subsequently.

4.1.3

The Open Loop Bifurcation Diagram

In addition to the frequency response approach, we also consider the structure of the bifurcations as both the drive amplitude and frequency are varied. The resulting bifurcation
diagrams clearly show the effects of modal coupling and form an important part of the
present analysis. In the single mode case, the system experiences only stable saddle node
bifurcations (SNB), that is SNB in which a stable and an unstable response merge and

46

�  �� ��

�  �� ��

������

������

������

���� � ���������

���� � ���������

������

������

������

������

������

������

������

������

������

������

�� ���

�� ���
�� ���
������� ���������

�� ���

�� ���

������

�  ��� ��

�� ���

�� ���

�� ���
������� ���������

�� ���

�� ���

�� ���

�� ���

�� ���

�� ���

�  ��� ��
������

������
������

������

������

������

���� � ���������

���� � ���������

������

������
������

������

������

������
������

������
������

�� ���

�� ���

�� ���
�� ���
������� ���������

�� ���

�� ���

������

�� ���

�� ���

�� ���
������� ���������

Figure 4.2 Experimental (blue) and model (red) results of the open loop frequency sweep
response for various drive amplitudes (20mV, 60mV, 200mV, 600mV), noting that the
feedthrough capacitance modifies the experimental amplitude results for larger forcing
levels. For the model results the dark red curve segments indicate a stable response and
the light red segments indicate an unstable response.
annihilate one another. However, modal interactions caused by internal resonance gives
rise to both stable and unstable SNBs and Hopf bifurcations (HB), where unstable bifurcations are associated with changes that involve only unstable response branches. We
first review these bifurcations and then examine the bifurcation diagram for the present
model.
We refer to the frequency sweeps, in Figure 4.2, to help explain the bifurcations. The
formal criteria for identifying these bifurcations are reviewed in the Appendix, which are
related to changes in the signs of the real parts of the eigenvalues of the model linearized
about a steady state operating point. The appearance and disappearance of solutions is
the hallmark of SNBs. As a result, when traversing the upper response branch (for any of
47

the frequency sweeps where Vin < 600mV), when the response jumps off the branch and
decays to a lower branch response, this is due to a SNB that annihilates the upper branch.
The SNB can also create solutions, such as the solutions after the “gap”, as seen when
Vin > 20mV. The HBs, on the other hand, appear in this system when solutions continue
to exist but switch stability as a complex conjugate pair of eigenvalues cross the imaginary
axis. There is a HB on the upper branch after the SNB bifurcation since that branch turns
from unstable (Vin = 60mV) to stable (Vin = 200mV). The hallmark of a (supercritical)
HB is that when the solution branch loses stability, a stable oscillatory response is created.
This is a preliminary explanation for how the solution jumps the “gap" for Vin = 600mV,
since after the initial SNB, as the frequency is increased, the response is first attracted to
the oscillatory response, then returns to the upper branch solution after the HB restabilizes the upper branch. (Note that here an “oscillatory” response corresponds to a beating
type response in the time domain, since it is the slowly varying amplitude and phase that
oscillate.)
In addition to considering the four frequency sweep responses, we can combine the
bifurcation data from many sweeps at different drive amplitudes into one plot, resulting
in the bifurcation diagram of interest. Utilizing the parameters in Table 4.1, we see that
the model predicts a bifurcation structure, as shown in Figure 4.3, consisting of stable and
unstable SNB (solid and translucent green curves), and stable and unstable HB (solid and
translucent red curves). Bifurcations from a stable response branch are shown as solid
curves while those from unstable response branches, which are not experimentally observable, are shown as translucent curves. For these bifurcation diagrams, the horizontal
axis is the frequency of operation and the vertical axis is the forcing voltage, so that, for a
given open loop frequency sweep, a horizontal line is traced out. The lower green branch
is the single mode standard Duffing SNB where the stable upper branch solution disappears by merging with an unstable response branch. The upper green branch is also a
standard Duffing SNB, where the lower stable branch disappears by merging with the un48

stable branch, as one sweeps in the other direction. The x and o points indicated in Figure
4.3 are the measured points where the experimental device undergoes a SNB sweeping
up or a HB sweeping up (jumping the gap) or sweeping down. Of course, only those bifurcations involving a stable response branch are observed. From this diagram is it clear
that the modal interactions begin to occur above the AC drive voltage Vin = 35mV, and
the SNB occurs near ω2 , rather than at the valued predicted by the single mode Duffing
model. This saturation in the SNB frequency is caused by the internal resonance and is
captured by the nearly vertical segment of the solid green SNB line, which provides a
close match with the experimental results. The SNB trend near ω2 is sensitive to the nonlinearity of the mode 2. The model here has a linear mode 2, for ease of calculation, but
the data shows a slight softening trend.
As one approaches Vin = 500mV, as seen in the frequency sweep case for Vin =
600mV, the resonator is able to jump the gap and maintain on the upper amplitude response branch. This is because of the HB structure that appears at these higher voltages,
as predicted.
With this successful model for the open loop system in hand, and confidence in the
type of modal coupling it experiences, we turn our attention to a closed loop oscillator
with this device as the frequency selective element. The goal is to understand how internally resonant coupled mode operation improves frequency fluctuation performance in
these systems, and estimate the limits one can achieve using this mode of operation.

4.2

The Closed Loop Analysis

In this section, we investigate the dynamics of a closed loop oscillator that employs an
internally resonant coupled mode system for its frequency selective element, with a focus
on the system frequency stability characteristics. We develop a model and analytical results that explain the significant frequency fluctuation improvements reported in [2] and

49

[67] and allow one to predict how this effect depends on system parameters and noise
characteristics. The device described in [2] is used to benchmark the analytical treatment,
and the results are used to explore the limits of this approach for frequency stabilization.
The experimental results are incorporated into the discussion to provide qualitative verification of the analytical predictions. Quantitative frequency fluctuation predictions are
not available, due to the lack of quantitive measures of the various noise sources in the
Argonne device, and the simplicity of the closed loop model, which ignores feedthrough
capacitance and some nuances of the phase shifter.
The oscillator model is based on a simple oscillator loop consisting of the standard
components, namely, an amplifier, and phase shifter, an amplitude limiter, and a resonator element, with the novel feature being that the resonator operates with two modes
of vibration that are internally resonant and nonlinearly coupled, as discussed above. In
the closed loop system the amplitude is fixed by the limiter, and the frequency of operation is set by this amplitude, the phase shift, and the resonator parameters. The dynamics
of the deterministic system are used to compute the operating conditions, and the effects
of small sources of noise in the resonator modes and the oscillator loop are determined by
linearizing about these operating points. The predictions show that indeed one can tune
the response to achieve dramatic reduction in the frequency fluctuations of the oscillator,
which is ultimately limited only by the intrinsic noise of mode 2.

4.2.1

Overview of Oscillators with Nonlinear Resonator Elements

The open loop analysis given above provides important insight into the behavior of the
closed loop system. The main difference between open loop and closed loop operation
is the locking of the phase operating point; in open loop this phase takes on a value determined by the system dynamics, while in closed loop it is manipulated by the phase
shifter. With the single mode linear open loop frequency response in mind, if the feedback phase is π/2, then the oscillator will operate at the peak of the response; this is also
50

true for the single mode nonlinear (Duffing) response. For a linear resonator, fluctuations
in the drive amplitude will cause fluctuations in the amplitude of the response, but since
the peak occurs at the same frequency for all amplitudes, these fluctuations will not affect the frequency. However, for a nonlinear resonator, both the response amplitude and
frequency depend on the drive amplitude, so that amplitude fluctuations are converted
into frequency fluctuations. This is the intuition behind the so-called amplitude to frequency (AF) fluctuations which often dominate the performance of oscillators that use a
nonlinear resonator element, and why linear resonator elements are preferred for oscillator loops. Note that this AF conversion occurs for any feedback phase, not just one that
sets the operating point at the peak.
We continue to use the open loop bifurcation diagram, Figure 4.3, to phenomenologically explain what happens in the closed loop. We are interested in what happens to
the peak operating frequency, Ω peak as we vary F, as would be the case for a feedback
phase of π/2. While obtain the peak is easy in the linear regim, experimentally obtaining the peak in the nonlinear regime is difficult; it is far easier to obtain the bifurcations.
Fortunately, the peak operating point for a nonlinear resonator is very close to the SNB
that occurs on the upper branch of the Duffing response curve. We will subsequently use
ΩSN , denoting the driving frequency where the SNB occurs, in place of Ω peak when in
the nonlinear regime. When the system is operated in the single mode nonlinear regime,
such as with Vin < 35mV in Figure 4.3, the AF conversion will be nearly identical to
that of the single mode resonator. The AF conversion is measured local to the operating
point by dΩSN /dF. As such, considering the Duffing SNB branch in Figure 4.3, flatter
response curves correspond to poor clock performance, while steeper response curves
relate to better performance, due to their AF characteristics. For example, the AF conversion degrades as one operates further into the single mode non-linear regime. However,
as Vin continues to increase, the resonator begins to operate in the coupled mode regime,
and the bifurcations begin to have a very steep slope, suggesting that dΩSN /dF ≈ 0 for
51

Vin = 600mV. This is analagous to the linear case when dΩ peak /dF = 0. This is known as
a “sweet spot” operating condition, where AF conversion is significantly reduced [36, 18].
In addition, there is a considerable range of operating points, from 50mV to 1V, where
dΩSN /dF is very small. Simulations of the ensuing analysis show this to be the case
for the present system. The unique feature of the present model is that the sweet spot
is generated by modal coupling, in contrast to the single mode approaches previously
considered [36, 18].
Using sweet spot operation for cancellation of the AF conversion is only one effect
of the nonlinear modal interaction. Another consequence of the interaction is that the
frequency stability of the (undriven) mode 2 dominates the clock performance when operating in this regime. This is due to the mode 2 frequency ultimately controlling the
frequency at which the modal interaction takes place; the interaction frequency is not determined at all by the fluctuations of the mode 1 frequency. The consequences of this
are considerable: one can evade electronic and feedthrough capacitance noise completely,
something that cannot be done in a single mode device.
In the next section, we provide an analysis of this closed loop system, showing how the
coupled mode response improves frequency stability. This model and predictive analysis
provide an essential step towards designing oscillators that take advantage of nonlinear
modal coupling from internal resonance.

4.2.2

The Closed Loop Oscillator Model

By assuming that there is no feedthrough capacitance, and that the amplifier, phase shifter,
and amplitude limiter are ideal, that is, their characteristics are static and independent of
frequency and amplitude (an idealization made for simplicity in the analysis), small adjustments to the open loop model equations can represent the closed loop model. From
Equations (.1) and (.2), we first set Ω → ω1 and φ1 → φ f , where φ f is the feedback phase,
dictated by the phase shifter. Here, the equations assume that |Ω − ω1 |/ω1 is small. In
52

contrast to the open loop equations, this closed loop model will cause a small deviations
from the actual solution since ω1 is not the actual operating frequency, whereas Ω is the
operating frequency for the open loop. As long as ΩCL , or the closed loop operating frequency, is near ω1 , this error will be small. Also, the forcing amplitude F is now the set
by the amplitude limiter, rather than by the amplitude of an external AC drive. In both
cases, F is the amplitude of the drive signal as seen from the perspective of the resonator.
Now, it is important to note that we are ignoring two aspects of the actual system in our
model: phase shifter dependence on operating frequency and feedthrough capacitance.
As such, we only expect qualitative agreement with experimental results.
Due to the qualitative nature of this model, we revert to looking at the modal equations
with the unscaled potential for simplicity. This is in contrast to the open loop characterization, which accounts for transduction gains for the modal amplitudes. Under these
assumptions, the model equations of motion for closed loop operation have the form

q¨1 + 2Γ1 q˙ 1 + ω12 q1 + γ1 q31 + 3βq21 q2 = F (q1 ) + ξ a1 (t) + ξ m1 (t)q1
q¨2 + 2Γ2 q˙ 2 + ω22 q2 + γ2 q32 + βq31 = ξ a2 (t) + ξ m2 (t)q2

(4.9)
(4.10)

with qi as the modal coordinates and noise sources, ξ ai and ξ mi , representing the additive and mulitplicative noises for the ith mode, respectively. As in the open loop model,
we assume that the feedback drive, F (q1 ), acts only on mode 1. While the drive clearly
acts on mode 2 at some level, it can be ignored since it is non-resonant, specifically at 1/3
of the resonant frequency of that mode. The drive is also dependent on the phase of q1
through the phase shifter, and the limiter sets the feedback amplitude at F. The intrinsic and feedback noises are represented, to first order, as the additive and multiplicative
noises, and stem from a variety of sources, including amplifier noise, feedthrough capacitance noise, intrinsic thermo-mechanical fluctuations, extrinsic thermal drift, and measurement noise. It is assumed that the overall noise additive and multiplicative noises are
53

small relative to the feedback force. These noise sources ultimately act to produce oscillator frequency fluctuations, the focus of the present investigation. There are additional
assumptions that are usually made on the noise sources for the sake of applying analytical tools, such assuming white noise with a flat power spectrum. While this is necessary
for traditional methods of averaging [55, 49], recent work has suggested that the results
are more broadly applicable to noises with a generic spectrum [12, 36].
Since we are interested in the frequency fluctuations over a long time period, it is
convenient to convert these equations to complex amplitude equations using the rotating
wave approximation, q1 = a1 exp(iω1 t) + c.c. and q˙ 1 = ω1 ( a1 exp(iω1 t) − c.c), and ignoring fast fluctuating terms (equivalent to the averaging procedure carried out for the
open loop system):
iFa1 iφs
3
i
e + Ξ1 (t, a1 )
3β( a1∗ )2 a2 + i
γ1 | a1 |2 a1 −
2ω1
2ω1
4ω1 | a1 |
(9ω12 − ω22 )
βi 3
1
a˙ 2 = − Γ2 a2 − i
a2 +
a1 + i
γ | a |2 a + Ξ2 (t, a2 )
6ω1
6ω1
2ω1 2 2 2

a˙ 1 = − Γ1 a1 +

(4.11)
(4.12)

which, as mentioned above, are similar to the open loop Equations (4.5) and (4.6), with
the already noted modifications. Note that the noise terms are converted in the averaging process, resulting in the terms Ξ1 and Ξ2 , which are considered in more detail below.
Also, the averaging here is only a first order averaging, ignoring the stochastic drift terms.
Since the noises are small, its effects are small, and including them in the analysis detracts
from the simplicity in describing how the internal resonance affects the noise characteristics of the system.
The complex modal amplitudes are represented as ai = ri exp(iφi )/2. To solve for the
steady state operating point, we transform the system to a four dimensional system with
r1 , r2 , φ1 , and φ2 as the state variables. As it typical for internal resonances, the dynamic
system actually depends on the modal amplitudes, ri , and only the phase difference, φ2 −
3φ2 , described here in the three-dimensional state vector R = ( r1 , r2 , φ2 − 3φ2 ). With

54

this notation the model equations can be expressed as
R˙ = f ( R) + g( R)Ξ,

(4.13)

where the noises acting on the modal amplitudes and phases are represented as Ξ =

[Ξr1 Ξr2 Ξφ1 Ξφ2 ] T , which consist of components near the resonant frequencies and DC
parts of their spectra, respectively. Note here that f ( R) is a three-to-three mapping and
g( R) is a 4 × 3 matrix.
The deterministic steady state operating point about which the noise acts is found by
¯ (Fluctuations in the the steady
setting the noise terms to zero and solving f ( R¯ ) = 0 for R.
state operating point caused by noise contribute only at higher orders due to the small
noise assumption, and are therefore neglected.) The local stability of these steady states
is determined by considering the local coordinate δR = R − R¯ and linearization. The

linearized model that includes the behavior of noise near R¯ is given by:

˙ = A( R¯ )δR + B( R¯ )Ξ(t)
δR

(4.14)

Here, the matrices A and B are the Jacobians for the deterministic and stochastic forces,
¯ The determinrespectively, that is, they are the linearizations of f ( R) and g( R) about R.
istic local stability of R¯ is dictated in the usual manner by the real parts of the eigenvalues
of the Jacobian matrix A. We next consider the deterministic operation of the oscillator,
which sets the stage for the analysis for the noise response, which is considered at the end
of the chapter.

4.2.3

The Deterministic Operating Frequency

The deterministic closed loop operating frequency, ΩCL , which is the output of interest, is obtained by substitution of R¯ into the φ˙ 1 equation and evaluating the expression
ΩCL = ω1 + φ˙ 1 . Plots of the noise-free steady state operating frequency as a function
of feedback amplitude F, from the model and for the experimental device, are shown in
55

Figure (4.5), showing the qualitative success of the model. Figure (4.5) indicates that in
single mode operation the frequency changes as a function of amplitude, as determined
by the cubic nonlinearity of mode 1. When ΩCL ≈ ω2 /3, however, the internal resonance stabilizes the oscillator frequency as a function of feedback amplitude, resulting in
an additional response branch that coexists with the single mode response over a range
of amplitudes beyond the ΩCL ≈ ω2 /3 point. On this branch the resonator frequency is
stabilized to near 1/3 that of the mode 2 frequency and therefore depends on the mode
2 cubic nonlinearity. This is indicated by the three different curves in the inset, which
represent different levels of modal softening in mode 2. This type of mode 2 dependency,
as mentioned in the previous section, arises from an anti-crossing behavior, known to occur in systems with linear coupling. Here, that notion is generalized to nonlinear modal
coupling. As shown in the figure, when the mode 2 cubic nonlinearity has an opposite
sign to that of mode 1, and with sufficient strength, the “sweet spot” appears, where
dΩCL /dF = 0, which is analogous to the dΩSN /dF = 0 condition in the bifurcation diagram of the open loop system. Note that there is a range of amplitudes/frequencies with
two stable operating points, representing single mode and coupled mode responses, and
the behavior will exhibit hysteresis as the amplitude F is varied up and down. The coupled mode regime is only reached by increasing F from a small value, while sweeping
down from a large value will experience only the single mode branch, except for a small
region near the anti-crossing point. We now turn our attention to the noise response.

4.2.4

The Oscillator Frequency Fluctuations

The eigenvalues of A not only determine the local stability of the operating point, they
also dictate how the system is perturbed by small levels of noise. The real parts of the
eigenvalues provide a relaxation time and the imaginary parts of the eigenvalues set the
frequency at which the slow complex amplitude will spiral around the steady state operating point. It is not surprising, then, that the noise analysis relies heavily on A. Specif56

ically, the behavior of closed loop frequency fluctuations due to the noises are found by
looking at the spectrum of the model linearized for small noise near the deterministic op˙ which is the variation of
erating point. We define the output for the linear system by δφ,
the oscillator frequency from its deterministic value, ΩCL . Using Equation (4.14) for the
˙ = CδR + DΞ where C and D map
underlying dynamics, this output is governed by δφ
the local deterministic (R) and noise (Ξ) dynamics onto the local frequency behavior. This
is solved in the frequency domain in conjunction with Equation (4.14), yielding the transfer function between the noise spectra and the spectrum of the frequency. The resulting
spectrum is given by

Sy ( ω ) = M ∗ T ( ω ) SΞ ( ω ) M ( ω )

(4.15)

where M = C (iω − A)−1 B + D. The coefficients of M = [mr1 mr2 mφ1 mφ2 ] represent
the weighting of how each of the noise sources, given by the elements of Ξ, map onto the
output frequency.
To illustrate the predictive capabilities of this model, we focus on the how the M coefficients behave as a function of the feedback forcing amplitude, F, since the actual strengths
of the spectrum SΞ , are highly system dependent, and not known for the Argonne device.
While the model is capable of predicting the spectrum, to show the frequency stabilization due to internal resonance we consider only the long term effects of frequency stability by evaluating M at ω = 0. In Figure 4.6, we consider two revealing quantities that
demonstrate what happens when the system enters into the internal resonance regime:
the contribution ratio of amplitude noise to frequency fluctuations, and the contribution
ratio of mode 2 noise to mode 1 noise.
The amplitude noise to phase noise ratio is obtained by taking the collective noise
weights from the modal amplitudes, |mr1 | and |mr2 |, which is the AF contribution, and
comparing them to the collective noise weights from the modal phases, |mφ1 | and |mφ2 |.
In the single mode regime, as expected, mode 1 is the dominant mode and the contribu57

tion of the mode 1 AF noise dominates. In fact, the AF contribution is several orders of
magnitude larger than the phase noise contribution, and increases as a function of amplitude due to the increasing effect of the nonlinearity. This trend is expected based on
our observations from the open loop, as discussed in the previous section. In the coupled
mode regime, however, the AF noise contribution drops by several orders of magnitude
and becomes comparable to the contribution from the phases. Experimentally, the improvement in frequency fluctuations qualitatively follows this trend as shown in the inset
in Figure 4.6. The gray curves represent the various amounts of the mode 2 nonlinearity,
as in Figure 4.5, which has little effect on the elimination of the amplitude to frequency
conversion. Although one might expect the AF contribution to become zero at the “sweet
spot”, as indicated by the green point, there is a small AF contribution from mode 2.
The contribution of the two modal noises to the overall frequency stability of the device is measured by taking the collective noise weights of mode 2, |mr2 | and |mφ2 |, and
comparing it to the collective noise weights of mode 1, |mr1 | and |mφ1 |. Not surprising,
outside of the internal resonance regime, the mode 1 noise dominates the contribution
to the frequency stability due to the large mode 1 AF contribution. However, as the system enters into internal resonance operation, we see that the contributions from the two
modes become comparable. The most interesting case is that when one approaches the
“sweet spot”, the stability is dominated by mode 2. This does not occur unless there is a
“sweet spot”.
The consequences of this shift in the noise filtering of the system are only made clear
when you consider the sources of the noise. Since mode 1 is the driven mode, the measurement, feedthrough, and feedback noises are all lumped into Ξr1 and Ξφ1 . If mode
2 is perfectly insulated from the electronics, then none of these noises will contribute to
the mode 2 noises. As such, perhaps the most important contribution provided by this
class of oscillator is that it has is the ability to evade all of the noises associated with the
feedback loop. Consequently, it is possible to obtain an oscillator whose stability is gov58

erned by thermal (both internal and external) noises acting on mode 2. If the device was
ovenized or temperature compensated, then one can operate at the fundamental intrinsic
thermal limit, or thermo-mechanical noise limit, of the oscillator [5]. This assumes, however, that the mode 2 noises from other sources are smaller than the thermo-mechanical
noise, which would require careful engineering. In any case, this approach offers a new
avenue for oscillator development.

4.3

Outlook

This chapter presented modeling and analysis of a frequency generator whose internal
resonant characteristics stabilize its frequency fluctuations by two to three orders of magnitude. At the outset of this chapter, we characterized an open loop model that produces excellent quantitative agreement with the open loop experimental results obtained
at ANL. We then expanded the resonator model to embody the dynamics of a simplified
closed loop circuit. This closed loop circuit model shows excellent qualitative agreement
with the closed loop experimental results in terms of both deterministic response and
noise characteristics. Further work in this area includes alternative designs, such as that
described in the Appendix, that allow for more thorough characterization of both modes,
and more detailed examinations of the noise sources and feedback impedance in order
to provide a better quantitative match. These will be crucial for development of devices
optimized to take advantage of couple mode response.
In terms of the more general impact on oscillator development, this work demonstrates that the decrease in frequency fluctuations due to internal resonance comes from
two sources: (i) a decrease in the AF noise conversion and (ii) an evasion of the feedback
noise. It may be possible to design devices utilizing internal resonance that outperform
traditional oscillators when the noise is dominated by either of these two sources.

59

������������
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

��

��� (�)

�

���

��-�

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
● ● ●
●�
�
��
●
���●
● �
�
�
●
●
�
���
● �
●
●
��� ●
� ●

����
�� ●
●
���●
●●
�●
�●●
�●●
●
●
�●●●●●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
�●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

��-�
�� ���

�� ���

�� ���
�� ���
������� ���������

�� ���

�� ���

������������

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
� ●●●
●
●
●
●
●
●
●
���
●
●
●
●
●
�●●●●●●●●●●●��
●
●
●
●
●
� ����●●�●
●
●
●
�
●
●
●
●
●
●
●
●
�
●
●
●
●
●
●
●
●
● ��
●
●
����●●●
●
●●�
●
●
●
●
●
●
●
�
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
�
●
●
●
●
●
●
●
●
●
�●●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● �●●●●●●●●●�
●
●
●
●
●
●
●
●
●
●
�
●
●
●
●
�����●● ●
●
�
●
●
●
●
●
●
●
●
��●●
�
●
●
●
●
●
●● ●
●
�
●
���●
●
●
●
●
●
●
●
��● ●
�
●
●
● ●
●
●
●
●
● �
��● ● ●
●
●
●
●
●
●●
●
��●
●
● ●
●
●
●
�
●
●
●●
●
●
● ●
�●
●
●
●
●
●
●
●
�● ●
●
●
�
●
●
●
●
�●
●
●
●
●
● ●
●
● ●
�●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
�
��
●
●
●
●
●
●
●
●
●
●
● ●
●
●
● �
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
�●
�
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●●● ●
●
● �
●
●
●
●
●
●
●
●
●
●
●
●
�
�
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
�●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

�

��� (�)

����

����
����

��-�

�� ���

�� ���
�� ���
������� ���������

�� ���

�� ���

Figure 4.3 Predicted saddle node (green curves, SNB) and Hopf bifurcations (red curves,
HB) as a function of drive frequency and drive amplitude, indicating bifurcations from
both stable (dark) and unstable (light) response branches. The experimentally measured
SNB and HB are denoted by (x) and (o), respectively. The bottom plot is a blowup of the
upper plot in the region of interest.
60

Limiter(F)	
  

Amplifier	
  
Phase	
  Shi/er	
  

Figure 4.4 The feedback loop circuit producing a self-sustained oscillation.

61

2

=0
2

=

2

=

1

2
1

( )
Figure 4.5 The mean (noise-free) operating frequency as a function of feedback
amplitude. (Left) The main graph shows predictions from the noise-free model. Left
inset shows corresponding experimental results where red indicates sweeping upward
in the amplitude and black indicates sweeping downward. The right inset shows a blow
up of the frequency as a function of amplitude in the internally resonance regime for
different levels of mode two cubic (Duffing) nonlinearity; the green point indicates zero
slope, that is, a point of zero dispersion. Model parameters for this figure and
subsequent figures are ω2 /ω1 ≈ 3.067, Q1 = 104 , Q2 = 104 , γ1 = 0.01, β = 0.00015.

62

|m+r1 | + |mr2 |
|m+ 1 | + |m 2 |

|mr2 | +
|m 2 |
+
|mr1 | +
+ |m 1 |

Amplitude Noise Regime

Phase Noise Regime
Mode 2 Regime
-

Mode 1 Regime

Figure 4.6 This figure shows two measures of the loop phase noise predicted from the
model as the forcing amplitude is varied. The three traces from gray to black represent
different levels of the mode two cubic nonlinearity, as in Figure 4.5. The upper plot is a
measure of the ratio of the contribution of the total (modes one and two) amplitude
noise to those of the total phase noise, and the lower plot is measure of the ratio of the
contributions of mode two to those of mode one. We see that in all three cases there is a
significant drop in the contributions of the amplitude noise and the overall mode one
noise when in the internal resonance operating regime. The inset shows experimentally
measure standard deviation of the frequency fluctuations as the amplitude is swept up
(red) and down (black), showing the drop in frequency fluctuations when operating in
the internal resonance regime.

63

CHAPTER 5
DESIGN OPTIMIZATION FOR NONLINEARITY
Results from the previous chapters describe applications where an internal resonance is
utilized in a MEMS application. All the devices were designed in an ad hoc manner, with
the required frequency ratios and modal coupling, but without any means of quantifying
system parameters at the design stage. The purpose of this chapter is to describe a method
for using finite element models to determine the linear and nonlinear modal parameters
that arise from the kinetic energy and the elastic potential, and to use this characterization
to optimize features of the nonlinear system response. This work will be important for
next generation devices, and has applicability to devices that operate using one and two
mode nonlinear response. This work was done in close collaboration with Suguang Dou
and Jakob Jensen of the Technical University of Denmark.
The optimization routine and single mode example were done primarily by our DTU
collaborators, and are included here for completeness. The author was primarily responsible for the general formulation of the normal forms from the nonlinear finite element
models and the coupled mode example. Much of this chapter is verbatim from the forthcoming publication [14].

5.1

Characterization

In order to demonstrate the methodology we consider systems for which the dominant
nonlinear effects arise from non-inertial conservative effects and we neglect the nonlinear inertial terms that might be important in, for example, cantilevered structures. This
includes planar frame structures for which mid-plane stretching is the primary nonlinearity. For convenience in implementation, the derivation is formulated in matrix-vector
form instead of tensor form. The full set of coefficients for individual mode nonlinear

64

stiffness and for nonlinear modal coupling are first derived for the Hamiltonian of the
system, and then the attendant essential coefficients in the reduced-order model are obtained in a straightforward manner.
The characterization theory is derived for general finite element models with quadratic
and cubic nonlinearities arising from nonlinear strain-displacement relations. With the finite element discretization of the continuous structure, the Hamiltonian for the system,
i.e. the sum of kinetic energy T and potential energy U, can be expressed as
Ne
1
1
H = T + U = u˙ TMu˙ + ∑
2
2
e =1

Ve

Tσ

dV

(5.1)

where M is the mass matrix of the finite element model, u is the global vector of nodal
displacements,

and σ are element-wise vectors of strain and stress components, respec-

tively, and Ve indicates that the volume integration is performed within the element e.
Assuming a linear strain-stress relation, the potential energy U can be written as
Ne
1
1
T
TC
σ dV = ∑
U= ∑
2 Ve
2 Ve
e =1
e =1
Ne

dV

(5.2)

where C is a symmetric constitutive matrix. We may now divide the strain into a linear
and nonlinear part as

= B0 ue +

1
B (ue ) ue
2 1

(5.3)

where B0 is the linear strain-displacement matrix and B1 (ue ) is a function of the elementwise vector of nodal displacements ue . Substituting Eq. (5.3) into Eq. (5.2) and using the
symmetry of C, the potential energy U can be divided into three parts, representing its
expansion at leading orders, as
U (2) =
U (3) =
U (4) =

Ne

1
(ue ) T B0T C B0 ue dV,
2 Ve
e =1

∑
Ne

1
(ue ) T B0T C (B1 (ue )) ue dV
2 Ve
e =1

∑
Ne

1
(ue ) T (B1 (ue )) T C (B1 (ue )) ue dV
8 Ve
e =1

∑

65

(5.4)

The first step in setting up the modal equations is to solve the corresponding linear eigenproblem

( ω 2 M − K) Φ = 0

(5.5)

where K is the linear stiffness matrix of the finite element model, to yield a set of modal
frequencies and mode shapes (ω p , Φ p ); these are normalized w.r.t. M such that Φ Tp MΦ p =
1 and Φ Tp KΦ p = ω 2p . We then express the displacements as a superposition of Nm linear
modes as
ue =

Nm

∑

p =1

q p Φep

(5.6)

where q p are the modal coordinates and Φep are the element-wise mode shape vectors
extracted from the global vector. Substituting Eq. (5.6) into Eq. (5.4), we obtain the strain
energy expressed in terms of mode shapes and modal coordinates as
U (2) =

Nm Nm

∑∑

i =1 j =1

U (3) =

(2)

αij qi q j ,

Nm Nm Nm

∑∑∑

i =1 j =1 k =1

U (4) =

(3)

αijk qi q j qk

Nm Nm Nm Nm

(4)

∑ ∑ ∑ ∑ αijkl

i =1 j =1 k =1 l =1

(5.7)

qi q j q k q l

(2)

where the linear modal coupling coefficients αij and the non-linear modal coupling co(3)

(4)

efficients, αijk and αijkl , are explicitly expressed as
Ne

(2)

1
(Φie ) T B0T C B0 Φej dV
2
V
e
e =1

(3)

Ne

1
(Φie ) T B0T C B1 (Φej ) Φek dV
2
Ve
e =1

(4)

Ne

αij =
αijk =
αijkl =

∑

∑

1
(Φie ) T
8
Ve
e =1

∑

B1 (Φej )

T

(5.8)

C B1 (Φek ) Φel dV

Since we use the linear eigenmodes normalized with respect to the mass matrix, we have
(2)

(2)

αij = 0 for i = j and αii

= ωi2 /2. With these modal coupling coefficients, we can now
66

write the Hamiltonian for the system in Eq. (5.1) in modal coordinates (qi , pi ) (where
pi = q˙ i ) as

H=

Nm

∑

i =1

Nm Nm Nm
Nm Nm Nm Nm
1 2 1 2 2
(3)
(4)
pi + ωi qi +∑ ∑ ∑ αijk qi q j qk +∑ ∑ ∑ ∑ αijkl qi q j qk ql
2
2
i =1 j =1 k =1
i =1 j =1 k =1 l =1

(5.9)

where the relation 12 p2i = 21 ( pi Φi ) T M( pi Φi ) has been used. With the Hamiltonian of the
system, it is straightforward to derive a set of ordinary differential equations in modal
coordinates as
q¨ p + 2ξ p ω p q˙ p + ω 2p q p +

Nm Nm

∑∑

i =1 j =1

p

gij qi q j +

Nm Nm Nm

p

∑ ∑ ∑ hijk qi q j qk =

i =1 j =1 k =1

f p (t)

(5.10)

where p = 1, . . . , Nm , q p is the modal coordinate corresponding to Φ p , and we have
introduced the modal force f p obtained from a projection of the load vector f in the full
finite element model onto the pth mode, i.e., f p = Φ Tp f, as well as modal damping ratios
(damping as a ratio of critical damping) expressed as ξ p . The modal coupling coefficients
p

p

gij and hijk are
p

(3)

(3)

(3)

gij = α pij + αipj + αijp
p
hijk

(4)
(4)
(4)
(4)
= α pijk + αipjk + αijpk + αijkp

(5.11)

It is noted that the differential equation representation of Eq. (5.10) with the modal coupling coefficients in Eq. (5.11) is equivalent to the formulation obtained using the principle of virtual work in previous studies [61].
In this work we apply the framework to structures modeled by beam elements [15].
However, it should be noted that the characterization and optimization procedure in this
paper works independent of the choice of elements.

5.2

Optimization and sensitivity analysis

In this section we present the general formulation for optimizing an objective function
that depends on non-linear coefficients of interest for a given model. For example, we can
67

Figure 5.1 Initial (uniform width) design and the first two linear vibration modes
obtained using COMSOL Multiphysics modal analysis. Left: linear vibration mode 1,
right: linear vibration mode 2, with ω2 ≈ 2ω1 . Color indicates vibration amplitude.

Figure 5.2 Optimized design for maximizing the absolute value of the essential modal
coupling coefficient and the two coupled linear vibration modes obtained using
COMSOL modal analysis. Left: linear vibration mode 1, right: linear vibration mode 2
and ω2 ≈ 2ω1 . The color in online version indicates the vibration amplitude.
maximize/minimize the hardening behaviour in a clamped-clamped beam by maximizing/minimizing the coefficient of the cubic non-linearity.
For generality, we consider an objective function c that may be an explicit function
of the nonlinear coefficients, as well as the eigenvectors and associated eigenvalues, and

68
1

formulate our optimization problem as the following minimization problem:
min
ρe

subject to (s.t.):

(3)

(4)

c(ω p , Φ p , αijk (Φ p ), αijkl (Φ p ))
he = hmin + ρe (hmax − hmin )
0 ≤ ρe ≤ 1

(element-wise beam thickness) (5.12)
(normalized design variable)

where the subscripts i, j, k, l, p = 1, . . . , Nm . The optimization problem is subjected to
a set of constraints associated with the beam shape parametrization: he is the thickness
of a beam element which is bounded by [hmin , hmax ] in the optimization and ρe is the
normalized design variable bounded between 0 and 1. In practice, the lower bound hmin
is dictated by fabrication tolerances, and the upper bound hmax is used to keep the beam
relatively slender. For coupled -mode resonators, such as the one treated in the second
example of this paper, we impose an additional constraint such that the ratio of two associated eigenvalues stays in a sufficiently small neighbourhood of n1 ω p1 = n2 ω p2 , where,
for internal resonance conditions, n1 and n2 are selected integers and p1 and p2 are the
orders of the two modes of interest.
For efficient structural optimization we will use a gradient-based approach. The sensitivity of the objective function can be calculated by using direct differentiation [45, 29],
but a more efficient approach for many design variables is the adjoint method [59] where
we merely need to solve Nm groups of adjoint equations. To derive the adjoint equation,
the objective function is first rewritten with adjoint variables λ p and η p as
(3)

(4)

c = c(ω p , Φ p , αijk (Φ p ), αijkl (Φ p )) +

Nm

∑

p =1

λ Tp (ω 2p MΦ p − KΦ p ) + η p (Φ Tp MΦ p − 1)
(5.13)

where it is noted that the terms in the two sets of parentheses that appear in the appended
term both vanish identically. Differentiation of the objective function with respect to de-

69

sign variable ρe is then expressed as
(3)
 ∂αijk

(3)
Nm Nm Nm
Nm ∂α
∂c
ijk dΦ p 
+
+

T dρe 
(3)
∂ρ
e
p=1 ∂Φ p
i =1 j=1 k=1 ∂αijk



Nm
∂c
∂c dΦ p
∂c dω p
dc
=
+ ∑
+
∑∑∑
∑
dρe
∂ρe p=1 ∂Φ Tp dρe
∂ω p dρe


(4)
(4)
Nm Nm Nm Nm
∂c ∂αijkl Nm ∂αijkl dΦ p Nm
T dΦ p + Φ T ∂M Φ
+∑
+∑ ∑ ∑ ∑
 + ∑ η p 2Φ p M

p
p
T
(4)
∂ρe
dρe
dρe
∂ρe
p=1 ∂Φ p
p =1
i =1 j=1 k =1 l =1 ∂αijkl

+

Nm

∑

p =1

λ Tp

ω 2p M − K

dΦ p
dω p
∂K
∂M
+ 2ω p MΦ p
+ ω 2p
−
dρe
dρe
∂ρe
∂ρe



Φp
(5.14)

dω p

dΦ p

Since the direct computation of dρ and dρ is computationally expensive, the adjoint
e
e
dω p
dΦ p
variables are selected such that the coefficients of dρ and dρ vanish. This leads to the
e

e

adjoint equations in terms of λ p and η p as
dc
+ ω 2p M − K λ p + 2MΦ p η p = 0
dΦ p
∂c
+ 2ω p Φ Tp M λ p
=0
∂ω p

(5.15)
(5.16)

where we use the symmetry of M and K, and
(3)

(4)

Nm Nm Nm
Nm Nm Nm Nm
∂c ∂αijk
∂c ∂αijkl
∂c
dc
=
+∑ ∑ ∑
+∑ ∑ ∑ ∑
(4) ∂Φ p
dΦ p
∂Φ p i=1 j=1 k=1 ∂α(3) ∂Φ p
i =1 j=1 k =1 l =1 ∂αijkl
ijk

(5.17)

In this case, the sensitivity of the objective function writes
(3)

(4)

Nm Nm Nm
Nm Nm Nm Nm
dc
∂c
∂c ∂αijk
∂c ∂αijkl
=
+∑ ∑ ∑
+∑ ∑ ∑ ∑
(4) ∂ρe
dρe
∂ρe i=1 j=1 k=1 ∂α(3) ∂ρe
i =1 j=1 k=1 l =1 ∂αijkl
ijk

+

Nm

∑

p =1

λ Tp

ω 2p

It is noted that the two quantities
ing element-wise quantities

∂αijk

(e)
∂Φ p

∂M
∂K
−
∂ρe
∂ρe

∂M
Φ p + η p Φ Tp
Φp
∂ρe

(3)
(4)
∂α
ijkl
ijk
∂Φ p and ∂Φ p
∂β ijkl

∂α

and

(e)
∂Φ p

(5.18)

are assembled from the correspond-

. For convenience of computational imple-

70

mentation, the adjoint equations are expressed in matrix form as


 

dc
2
2MΦ p  λ p 
 ωpM − K
 dΦ p 


  = −
∂c
2ω p Φ Tp M
0
ηp
∂ω

(5.19)

p

Based on the objective function and calculated sensitivities, the update of design variables is performed using the mathematical optimization tool called the method of moving asymptotes (MMA) [56], which solves a series of convex approximating subproblems.
The algorithm has proven efficient for large-scale structural optimization. A new system
analysis is then performed with the updated design variables. These steps are repeated
until the design variables no longer change within some prescribed small tolerance. To
summarize the procedure, a flow chart of the proposed optimization is displayed in Fig.
5.3. We now demonstrate the approach with examples.

5.3

Examples

Two examples are presented. The first is to maximize/minimize the cubic nonlinearity of the fundamental mode in a nonlinear resonator consisting of a clamped-clamped
beam, a common element in MEMS applications. The second example shows how one
can maximize the essential modal coupling nonlinearity in a T-bar frame with a 2:1 internal resonance, a structure that has been proposed as a MEMS frequency divider.

5.3.1

Single-mode resonator

A crucial feature of a nonlinear resonator is the hardening and softening behaviour associated with a given vibration mode. For a lightly damped single-mode resonator, we
will focus on the hardening and softening behaviour of its free responses, i.e., without
damping and external loads. This problem has been investigated from a structural dynamics/finite element perspective; see, for example, [60]. For its applicability in general
structures, including symmetric structures like clamped-clamped beams and asymmetric
71

Initial design

Solve one eigenvalue problem Eq. (2.5)
Calculate modal coupling coefficients using Eqs. (2.8) and (2.11)
Evaluate the objective function in Eq. (2.13)
Solve one adjoint equation Eq. (2.20)

update model

Calculate all sensitivities using Eq. (2.19)
Update of design variables using MMA
No

Stop?
Yes
Final design

Figure 5.3 Flowchart of the structural optimization routine for nonlinear dynamic
response.

Figure 5.4 Initial (uniform width) design of a clamped-clamped beam and its linear
vibration mode. Color indicates vibration amplitude.
structures like arches or shells, the model development includes both quadratic and cubic
terms. For single-mode resonator, the reduced order model based on a single linear mode
is written as
q¨ p + ω 2p q p + g pp q2p + h ppp q3p = 0
p

72

p

(5.20)

Figure 5.5 Optimized design for maximizing the cubic nonlinearity of a
clamped-clamped beam and its linear vibration mode. Color indicates vibration
amplitude.

Figure 5.6 Optimized design for minimizing the cubic nonlinearity of a
clamped-clamped beam and its linear vibration mode. Color indicates vibration
amplitude.
p

(3)

(4)

p

with g pp = 3α ppp and h ppp = 4α pppp . The frequency-amplitude relation is derived as
ωNL = ω p (1 + Γ a2 )

(5.21)

where a is the amplitude of the corresponding linear mode, and the effective coefficient Γ
is
p

p

2
5 ( g pp )
3 h ppp
−
Γ=
8 ω 2p
12 ω 4p

(5.22)

A more accurate model for the hardening/softening behaviour of these structures can
be created with the dynamics projected onto a single nonlinear normal mode and the
corresponding equation is then written as
p
p
p
R¨ p + ω 2p R p + ( A ppp + h ppp ) R3p + B ppp R p R˙ 2p = 0
p

(5.23)

p

with the new coefficients A ppp and B ppp given as explicit functions of nonlinear modal
l and hl , [60]. Based on Eq. (5.23), the frequency-amplitude relacoupling coefficients gij
ijk

tion is derived as
ωNL = ω p (1 + Γ∗ A2 )
73

(5.24)

15

final design

iteration=30

Objective function

10

iteration=20

5
iteration=10

inital design

0

0

10

20

30

40

50
60
Iteration number

70

80

90

100

Figure 5.7 Evolution of the objective function and shapes encountered during the
maximization process. The vertical axis is the absolute value of the objective function
divided by its initial value and the horizontal axis is the iteration number.
where A is the amplitude of the corresponding nonlinear normal mode in curved, normal
coordinates, and the effective coefficient Γ∗ is
p

Γ∗ =

p

p

3 ( A ppp + h ppp ) + ω 2p B ppp

(5.25)

8 ω 2p

as approximated by averaging methods. Note that in this case Γ∗ is not only linked to the
p

p

pth linear mode explicitly through ω p and h ppp and implicitly through g pp which is in
p

p

A ppp and B ppp , but also linked to other linear modes, whose contributions are taken into
p

p

account through A ppp and B ppp .
74

inital design
1

0.9
iteration=2

Objective function

0.8

0.7
iteration=5
0.6

0.5
iteration=10
0.4

0.3

final design

0.2

0

10

20

30

40
50
Iteration number

60

70

80

Figure 5.8 Evolution of the objective function and shapes encountered during the
minimization process. The vertical axis is the absolute value of the objective function
divided by its initial value and the horizontal axis is the iteration number.
We find that for the flexural mode of a clamped-clamped beam with geometric nonlinp

earity from mid-plane stretching, the dominating coefficient is h ppp and we can therefore
simplify our optimization problem considerably by approximating Γ∗ as:
p
(4)
3 h ppp
3 α pppp
∗
=
Γ ≈Γ=
8 ω 2p
2 ω 2p
p

(5.26)

since g pp ≡ 0, due to the symmetry. Based on the coefficient Γ∗ in Eq. (5.26) and omitting

75

the constant factor, the objective function is now selected as
(4)

min

c=±

ρe

α pppp

(5.27)

ω 2p

where the plus/minus sign corresponds to minimizing/maximizing the hardening behaviour, respectively. This simplification is possible since we consider the fundamental
mode of a clamped-clamped beam, which has a strictly hardening nonlinearity, that is
(4)

α pppp > 0. Substituting the objective function c in Eq. (5.27) into Eq. (5.18), its sensitivity
with respect to design variables ρe is obtained as
(4)

1 ∂α pppp
dc
=± 2
+ λ Tp
dρe
ω p ∂ρe

ω 2p

∂M
∂K
−
∂ρe
∂ρe

Φ p + η p ΦTp

∂M
Φp
∂ρe

(5.28)

with adjoint variables λ p and η p solved from adjoint equation in Eq. (5.19) with
(4)

dc
1 ∂α pppp
=± 2
,
dΦ p
ω p ∂Φ p

(4)

α pppp
∂c
= ∓2 3
∂ω p
ωp

(5.29)

where
(4)

∂α pppp
∂Φep

=

1
(Φe ) T
2 Ve p

B1 (Φep )

T

C B1 (Φep )

dV

(5.30)

The beam has a fixed length L of 300 µm and a fixed out-plane width of 20 µm. The initial design has a uniform in-plane thickness of 4 µm, and is discretized with 400 beam
elements as described in [15]. During shape optimization the in-plane thickness h is varied to tailor the cubic nonlinearity in the reduced order model. We set hmin =2 µm,
and hmax =6 µm. The material properties are assumed for Si, that is, mass density
ρ = 2329 kg/m3 , and Young’s modulus E = 170 GPa. The vibration modes of the initial
design and two optimized designs are shown in Fig. 5.4, Fig. 5.5 and Fig. 5.6, respectively. Evolution of the objective function and shapes obtained during the evolution are
shown in Fig. 5.7 and Fig. 5.8. In these optimizations, the objective function is increased
by a factor of 13 and reduced by a factor of 4, respectively. The optimized designs are
in accordance with the results in [15], obtained using the incremental harmonic balance
76

method, where we found the nonlinear strain energy due to mid-plane stretching reaches
its local maximum around x = 14 L and x = 43 L, which is precisely where the optimized
structures are altered most significantly relative to their general thickness. Furthermore,
the eigenfrequency of the first flexural mode decreases during optimization of maximizing the cubic nonlinearity, and increases during optimization of minimizing the cubic
nonlinearity.This follows from the fact that the structure is made generally thinner when
maximizing the nonlinearity, so that the critical sections can be made relatively thick,
and the opposite trend occurs when minimizing the nonlinearity. It should be emphasized that, no guarantee can be made that the obtained designs are global optima and
the found solutions will in general depend on the chosen initial design. However, repeated optimization runs reveal that the formulation is quite robust and we believe only
a marginal improvement in performance could be possibly found with another design.

5.3.2

Coupled-Mode Resonator with Internal Resonance

We present an example of two-to-one internal resonance for which the normal form has
a single important inter-modal coupling term. This example is motivated by a MEMS
frequency divider that makes use of internal resonance, for which a general theory is presented in [53]. The structure considered is also similar to other proposed MEMS devices
developed for filtering [68, 69]. In these structures two vibrational modes can exchange
energy during free vibration. For convenience, we will refer to modes p1 and p2 as modes
1 and 2, respectively, in our discussion. The divider device consists of two localized
modes with ω p2 ≈ 2ω p1 , which provides a division by two in frequency, as measured
in mode 1, when mode 2 is driven near its resonance. The term “localized mode” indicates that the dominant vibration associated with the mode occurs in a localized part of
the structure, even though the entire structure is generally involved in the modal vibration. For instance, for the T-bar structure shown in Fig. 5.1, the vibration of mode 1 is
localized in the vertical beam and the vibration of mode 2 is localized in the horizontal
77

beam. In operation, a harmonic load with frequency close to ω2 is applied to drive mode
2 into resonance. When the amplitude of mode 2 is sufficiently large, it will induce the
vibration of mode 1 due to the parametric pumping, wherein the transverse vibration of
the horizontal beam provides an axial force in the vertical beam, which in turn induces
transverse motion of the vertical beam when the frequency of the horizontal beam is approximately twice that of the vertical beam. The response of the structure is governed by
a model in which these two modes are coupled through resonant terms. In application, n
such elements are linked to achieve division by 2n ; division by eight has been experimentally achieved [46]. The reduced order dynamic model with its conservative dynamics
projected onto two nonlinear normal modes in curved, normal coordinates is written as
1 + g1 ) R R + ( A1 + h1 ) R3 + B1 R R
˙2
R¨ 1 + ω12 R1 + ( g12
21 2 1
111
111 1
111 1 1
1 R
˙ 2 R1 + B1 R2 R˙ 2 R˙ 1 = 0
+ ( A1212 + A1122 + h1122 + h1212 + h1221 ) R22 + B122
2
212
2 R2 + ( A2 + h2 ) R3 + B2 R R
˙2
R¨ 2 + ω22 R2 + g11
222 2 2
222
222 2
1

+

( A2112 +

2 R
˙2
A2211 + h2112 + h2121 + h2211 ) R21 + B211
1
p

R2 +

2 R R
˙
B112
1 1

R˙ 2 = 0

(5.31)

(5.32)

p

where the explicit expressions of Aijk and Bijk are given in [60]. The key term of interest is
1 , g1 and g2 , since they
that associated with the essential modal coupling coefficients g12
21
11

are the terms in the normal form that describes the resonant nonlinear coupling terms
that promote energy exchange between the modes. These modal coupling coefficients
can also be observed from the governing equation with its dynamics projected onto two
linear modes, which is written as
1 + g1 )q q + g1 q2 + h1 q3 + others = f (t )
q¨1 + ω12 q1 + ( g12
1
21 1 2
11 1
111 1
2 q2 + g2 q2 + h2 q3 + others
q¨2 + ω22 q2 + g11
222 2
22 2
1

(5.33)

= f 2 (t)

1 , g1 and g2 , between which there is an
where the most important coefficients are g12
21
11
1 + g1 = 2g2 . The reason for this exact relation is that the coefficients
exact relation g12
21
11

of q1 q2 and q21 in Eq. (5.33) arise from differentiation of the term βq21 q2 in the Hamilto(3)

(3)

(3)

nian H with respect to q1 and q2 , where β = α112 + α121 + α211 . In fact, this example is
78

particularly attractive for demonstrating the present approach since energy exchange between the modes can be described (to leading order) by this single nonlinear term. Also,
as observed in the equation for q1 in the reduced order model, β is essentially the amplitude of the parametric excitation provided to mode 1 from mode 2, given by 2β 1 q2 q1 , and
is therefore related to the stability region of the linear parametric resonance of q1 . The
q21 term in the equation for q2 captures the back-coupling of the driven mode onto the
driving mode, as required for passive coupling, which is beneficial in frequency dividers
[46]. The combined effect of this coupling is the possibility of energy transfer between the
modes when there is a two-to-one internal resonance.
The optimization problem for this resonance is therefore formulated as
min
ρe

s.t.
where

(3)

(3)

(3)

c = − α112 + α121 + α211

(5.34)

|ω1 /ω2 − 1/2| ≤

= 0.001 and additional constraints are imposed as in Eq. (5.12). The sensitivity

is computed with Eq. (5.18) written as


(3)
(3)
(3)
2
∂α
∂α
∂α
dc
= S ·  112 + 121 + 211  + ∑ λ Tp
dρe
∂ρe
∂ρe
∂ρe
p =1

ω 2p

∂K
∂M
−
∂ρe
∂ρe

Φ p + η p Φ Tp

∂M
Φp
∂ρe
(5.35)

(3)

(3)

(3)

where S = −sign α112 + α121 + α211 . It is noted that there are two groups of adjoint
variables, λ1 and η1 , λ2 and η2 , corresponding to the two vibration modes. These two
groups of adjoint variables are solved using adjoint equation in Eq. (5.19) with p = 1, 2.
For a specific example the length of the horizontal beam is taken to be 300 µm and the
length of the vertical beam is taken to be 195.5 µm, so that ω2 ≈ 2ω1 .The lengths of the
two beams are fixed during the optimization. The initial in-plane thickness is uniformly
4 µm along both beams, and the in-plane thickness is bounded between 2 µm and 6 µm
during the optimization. The material properties are the same as in example 1. An opti-

79

mized design and its two important vibration modes are displayed in Fig. 5.2. Evolution
of the objective function and optimized designs over iterations are displayed in Fig. 5.9.

final design

10
iteration=50
9

8

Objective function

7

6

iteration=30
intial design

5

4

3

2

1
0

50

100
150
Iteration number

200

Figure 5.9 Evolution of the objective function and shapes encountered during the
optimization process. The vertical axis is the absolute value of the objective function
divided by its initial value and the horizontal axis is the iteration number.
In order to demonstrate the effects of tuning and optimizing the coupling nonlinearity
in this example, we carry out simulations of the structure for different shapes encountered during the optimization iteration process. With this internal resonance and zero
damping, free vibrations of the system will exhibit energy exchange between the modes
in a beating type response, whose features depend crucially on the magnitude of the coupling coefficient and the initial conditions, as well as the other system parameters. For
80

each structural shape and initial amplitude of the starting mode, the exchange of energy
has a particular beat period; this period will be shorter for higher starting amplitudes and
for larger values of the coupling coefficient, since both these effects enhance the nonlinear
modal coupling. The left panel of Figure 5.10 shows a typical response obtained from the
finite element model of the final (optimized) shape, obtained by initiating the response
using the second linear mode at a moderate amplitude. The right panel of Figure 5.10
shows the normalized beat period for several values of the initial energy for three different designs, that is, for three different values of the coupling coefficient. The predicted
trends are evident as both the coupling coefficient and initial amplitude are varied. Note
that it would be also worthwhile to compare results from the finite element model with
simulations and analysis of the reduced two mode model, but this requires a more detailed study of the effects of all nonlinear coefficients that vary during the optimization
process.
−11

x 10

100

mode 1
mode 2
Normalized Beat Period

6
Modal Amplitude

4
2
0
−2
−4

80
60
40
20

−6
0

1

2
Time (s)

3

0

4
−5

x 10

0.5

1
1.5
Mode 1 Initial Amplitude

2
−10

x 10

Figure 5.10 Left: Typical time simulation with added envelope curves showing beats, for
the final design with initial first mode amplitude of 7.5 × 10−11 (using eigenvectors
normalized by the mass matrix). Right: Data points and fitted exponential curves for the
beat period, normalized by the period of the first linear mode, versus the initial
amplitude of the first mode, for designs corresponding to iteration numbers 0 (initial
design, top curve), 30 (middle curve), and 250 (final design, bottom curve).

81

Evolution of the eigenfrequencies of linear vibration modes 1 and 2 encountered over
iterations of the optimization process is displayed in Fig. 5.11(a). Other measures of interest for this system are (i) the degree of spatial energy localization in the vertical beam of
the first vibration mode (note that from symmetry the second vibration mode has perfect
localization) and (ii) the effectiveness of the horizontal beam in parametrically pumping
the vertical beam in the second mode. The localization of the first mode is measured
by the maximum transverse amplitude of the horizontal beam divided by the maximum
transverse amplitude of the vertical beam. Likewise, the pumping effectiveness of the
horizontal beam in the second mode is measured by the ratio of the transverse vibration
at the midspan of the horizontal beam to the maximum transverse vibration of the same
beam, which occurs near the quarter spans. The results in Fig. 5.11(b) show that the localization ratio decreases during optimization, which indicates improved localization, and
the pumping ratio increases, which indicates enhanced coupling of the two modes. The
intuition of the optimized design is that the axial deformation along the vertical beam
increases with a larger end mass and a thinner cross section, both of them occur in the
optimized design.

5.4

Outlook

This work provides a framework for shape optimization of a mechanical structure based
on nonlinear coefficients of a reduced order model based on normal forms [20]. The first
step is to extract the nonlinear coefficients of a reduced order model from a finite-element
model of the mechanical structure. These nonlinear parameters are then available to use
as objective functions for optimization routines. In the examples, we considered the shape
optimization for minimizing and maximizing the Duffing nonlinearity associated with
the fundamental mode of a clamped-clamped beam, and the in-plane flexural vibration
of a frame with 1:2 internal resonance, in which the quadratic modal coupling coefficient

82

2.5

0.3

2

0.25

1.5

0.2
Ratio

Frequency (MHz)

(a)

1

(b)

pumping ratio
localization ratio

0.15
0.1

0.5
0

mode 1
mode 2
0

50

100
150
Iteration number

0.05
0

200

0

50

100
150
Iteration number

200

Figure 5.11 (a) Evolution of the eigenfrequencies of linear vibration modes 1 and 2
encountered during the optimization process. (b) Evolution of the pumping ratio and
the mode localization ratio encountered during the optimization process.
is maximized. This approach creates designs that alter the relevant nonlinear coefficients
by about an order of magnitude.
While this work focuses on mechanical nonlinearities that arise from a potential, future work is geared to adding nonlinearities due to inertial effects and electrostatic nonlinearities, as well as other relevant physics. In addition, although this work focused
on maximizing the nonlinearity, future work will focus more on minimization of nonlinearity, thereby increasing the linear range of devices. This approach will also be useful
when designing dynamic MEMS devices such as frequency dividers, oscillators, and disk
resonator gyroscopes.

83
1

CHAPTER 6
CONCLUSION
As electronic and MEMS devices are pushed towards miniaturization, they are often
driven in the nonlinear o‘perating regime in order to obtain a good signal to noise ratio.
Most previous studies of nonlinear behavior in MEMS resonators have focused on single
mode operation. In this work, we demonstrate the use of nonlinear modal coupling in
MEMS frequency generators and converters.
In Chapter 2, we examined the extent to which a conceptual model of a passive Nstage frequency divider could lead to a stable solution. This conceptual model consists
of a chain of parametrically coupled resonator modes whose consecutive modes are arranged in a 2:1 frequency ratio. Since the energy is allowed to flow forwards and backwards in the cascade, one must consider the fully coupled response of the entire cascade,
as opposed to merely considering sequentially driven individual stages. We found that
the system could be tuned so that the cascade would settle to a steady state operation
wherein all stages are active.
In Chapter 3, we developed MEMS designs that exhibit the frequency division characteristics that were described by the models in Chapter 2. The designs consist of perpendicular beams connected by springs, ultimately providing a network of localized modes
with quadratic nonlinear modal interactions. The designs were fabricated and tested by
our collaborators at University of California at Santa Barbara. The main results were
a device that achieved successful division by 8. The devices did not exactly match the
theory in Chapter 1 due to the parameter ratios not being as specified for the equal amplitude solution, and the modal coupling coefficients could not be quantified for the device.
There is further work to be done in the design and characterization of these devices. Also,
additional design topologies were fabricated, but have yet to be tested.
In Chapter 4, we investigated the role of internal resonance in stabilizing the frequency
84

fluctuations of a frequency generator. This work was motivated by the experimental observations from Argonne National Labs, as described in [2]. The internal resonance is
associated with a modal frequency ratio of 1 : 3, where the fundamental mode is the
driven mode. The resonator was characterized through an open loop investigation, and
we found excellent agreement between the model and the experiments, in terms of coupled mode frequency response bifurcation diagrams. The closed loop model, which is
more complex than the open loop, provides good a qualitative match with the experiments, and clearly indicates the model is sufficient to capture the phenomena for the frequency stabilization. To achieve quantitative agreement, one must have specific knowledge about the noise sources, which is lacking for the experimental device. The main
result of this work is that the model shows that a frequency generator can move beyond
the limits imposed by the measurement and feedback noise in the main loop, due to the
nonlinear coupling to the second mode, which is nearly insulated from the actuation and
detection and thus is not subjected to the noise felt by the driven mode. If the device is
operated where there is a strong modal coupling, the precision of the clock is dependent
on the noise associated with the secondary mode, which is limited by its own thermomechanical noise. Future work is necessary to reach this metric, which may be feasible,
particularly by investigating systems with a 1:2 modal frequency ratio.
In Chapter 5, we developed and applied a tool that allows for shape optimization
for nonlinear response, where the cost function is associated with nonlinear phenomena,
specifically coefficients in a normal form for the resonance of interest. The first step was
to extract the required normal form coefficient (or combination thereof) from a nonlinear
computational model. This coefficient is maximized or minimized using a gradient-based
optimization algorithm, by changing the geometry of the beam elements, specifically,
the local beam thickness. Two examples wer considered: the Duffing nonlinearity in a
clamped-clamped beam and the modal coupling coefficient for a 1:2 internal resonance of
a T-frame structure, of the type developed as a frequency converter. The results show that
85

one can increase/decrease the coefficient of interest by about an order of magnitude. The
results from the modal coupling example were verified by examining modal energy exchange during transient motion using time-simulations of finite element models. These
tools open up the possibility for shape optimization based on dynamic nonlinearities.
Future work in this area includes extending the method to physics beyond mechanics,
such as electrostatics and thermal quantities, and implementation in a variety of MEMS
resonator applications.
In summary, the body of work described in this dissertation indicates that, in addition
to providing rich dynamic behavior that is of academic interest, nonlinear modal coupling
in MEMS holds promise for improving the performance of a variety of devices relevant
to signal processing technologies.

86

APPENDIX

87

APPENDIX

In this Appendix, we provide supplementary material to support the results presented
in Chapter 4. We begin with some solutions for a simple two resonator interaction with
nonlinear modal coupling and review how to identify various bifurcations that appear in
the open loop system as the drive strength and frequency are varied. Finally, we provide
an alternative resonator design for future studies of internal resonances, one that has
certain advantages over the Argonne device.

A.1

The 1 : 3 Modal Amplitude Equations

In this section, we derive the averaged equations, look at a specific case where the second
mode is linear and has faster dynamics than the first mode, and obtain the closed form
solution for the single mode Duffing resonator.

A.1.1

The Averaged Equations

The most basic equations of motion that capture the phenomena of interest, expressed in
terms of the flexural and torsional degrees of freedom coordinates, Si , are:

S¨1 +2Γ1 S˙ 1 + ω12 S1 + α1 S22 S1 + γ1 S13 + 3β 1 S12 S2 = F cos(Ωt)

(.1)

S¨2 +2Γ2 S˙ 2 + ω22 S2 + α2 S12 S2 + γ2 S23 + β 2 S13 = 0

(.2)

noting that in addition to the Duffing and 1 : 3 modal interaction nonlinearity, we also
include the dispersive coupling nonlinearity with coefficients α1 and α2 ./ We ignore the
effects of the frequency shift and noise contributions of the dispersive coupling terms,
since they are not needed to capture the behavior of current interest.

88

The external drive clearly acts on all modes, but is assumed here to act only on the
primary mode. This is justified when the drive acts dominantly on the primary mode
and/or it is nonresonant with the secondary mode. The effects of the drive and the form
of nonlinear modal coupling depend on the system details, including any symmetries,
and must be considered on a case-by-case basis.
Approximate techniques, both analytical and/or numerical, are required to investigate the above set of coupled differential equations. It is convenient to consider the
dynamics of the slowly varying complex amplitudes, which describe the modal amplitude and phase dynamics. To this end we express S1 = s1 exp(iΩt) + c.c., S2 =
s2 exp(inΩt) + c.c., S˙ 1 = iΩs1 exp(iΩt) + c.c., and S˙ 2 = iΩs2 exp(iΩt) + c.c.. Differentiating the above again and utilizing the constraints from S˙1 and S˙2 , we obtain S¨1 =

−Ω2 S1 + 2iΩs˙ 1 exp(iΩt) and S¨2 = −9Ω2 S2 + 6iΩs˙ 2 exp(i3Ωt). Substitution of these into
Equations (.1) and (.2) and integrating over a timescale of the slowest period of oscillation,
T = 2π/Ω, we obtain complex amplitude equations:



α1
i
3
iF
3β 1 s1∗2 s2 + i
γ1 |s1 |2 s1 −
s˙ 1 = −Γ1 − i δω1 − |s2 |2  s1 +
Ω
2Ω
2Ω
4Ω
s˙ 2 =

−Γ2 − i δω2 −

α2
| s |2
6Ω 1

β i
3
s2 + 2 s31 + i
γ | s |2 s
6Ω
6Ω 2 2 2

(.3)

(.4)

with detuning parameters defined as

(Ω2 − ω12 )
δω1 =
2Ω
2
(3Ω − ω22 )
δω2 =
2Ω

(.5)
(.6)

These equations are the normal form for this resonance and describe its dynamics, which
are known to be very complicated, as demonstrated by the Argonne device.

89

A.1.2

A Special Case

In order to streamline the analysis, we assume that the second mode is linear, γ2 = 0, and
ignore diffusive coupling, α1 = α2 = 0, so that the second mode complex amplitude is
governed by
β i
s˙ 2 = [−Γ2 − iδω2 ] s2 + 2 s31 .
6Ω

(.7)

The behavior of interest is captured by this simplified model with the benefit of a less
complicated analysis, but at the expense of some accuracy when comparing the model to
experiments. As shown in the main text, for example, the experimental results suggest
that the second mode Duffing nonlinearity does come into play. Fortunately, the effect
on the response is a predictable deviation from the system with a linear second mode.
Another convenient simplification is to assume that ω1 Γ1 << ω2 Γ2 , which implies that
the second mode dynamics decays much faster than the first mode. In this case the second mode adiabatically tracks the first mode. This assumption provides a single degree
of freedom model for the slow mode that complicated by coupling to the fast mode. This
model provides the same steady-state amplitude structure as the model given in Equations (.3) and (.4), but it fails to capture Hopf bifurcations in the system. With these assumptions, we can solve for the steady state response of the second mode response in
terms of the first mode, by assuming s˙ 2 = 0, yielding
s2 =

β2
s31
iΓ
+
δω
( 2
2)

(.8)

where we define parameter β 2 as
β2 =

β2
.
2
6Ω(Γ2 + δω22 )

(.9)

A model for the slow mode dynamics in this case can be obtained by substitution of
Equation (.8) into Equation (.3), resulting in a dynamic model for s1 . The steady state
response of the first mode can then obtained from s˙ 1 = 0, and an equation for steady
90

state involving only s1 is obtained by employing Equation (.8). This results in a fifth
order equation for s1 ,

[−Γ1 − iδω1 ] s1 + β(−Γ2 + iδω2 )|s1 |4 s1 + i

iF
3
γ1 |s1 |2 s1 −
=0
2Ω
4Ω

(.10)

where
β=

nβ 1 β 2
2Ω

(.11)

Re-organizing the terms, we obtain

(−Γ1 − Γ2 β|s1 |4 ) − iδω1 s1 + i

3
iF
γ1 + βδω2 |s1 |2 |s1 |2 s1 −
= 0.
2Ω
4Ω

(.12)

This equation shows that the coupling term has the effect of renormalizing the effective
damping and cubic nonlinearity of the first mode, resulting in higher order nonlinearities
in stiffness and damping.
An expression for the mode 1 steady state real amplitude is found by letting a21 = |s1 |2 ,
and solving for a1 ,
a21

2
−Γ1 − Γ2 βa41 +

2
3
F2
2
4
γ1 a1 + βδω2 a1
=0
−δω1 +
−
2Ω
16Ω2

(.13)

which is a fifth order polynomial in a21 .
A.1.3

The Single Mode Duffing Case

We can also reduce the above to the single mode duffing case by letting β = 0. Letting
z1 = a21 , the equation for the amplitude becomes
2
F2
3
−
=0
γ z
z1 (−Γ1 ) + −δω1 +
2Ω 1 1
16Ω2
2

(.14)

which can be explicitly solved for Ω(z1 ), providing two branches of solutions:
Ω2

=

2z1 (−2Γ21 + 3γ1 z1 + ω12 ) ±

( F 2 − 16Γ21 z1 (−Γ21 + 3γ1 z1 + ω12 ))z1
2z1

91

(.15)

which map out the Duffing frequency response [39]. The phase is obtained from
φ1 = tan−1

−8Γ1 Ω

−8Ωδω1 + 3γ1 z1

(.16)

This single mode response is useful for characterization since one can obtain the mode 1
parameters by operating the first mode away from the internal resonance condition.

A.2

Determining Stability and Bifurcation Criteria

In this section, we describe the stability and types of bifurcations of the steady state responses in the four-dimensional dynamical system representing the slowly varying amplitudes and phases of the modes in the open loop response, Equations (.3) and (.4). This
analysis can be found in standard texts such as in [21], and was used to generate the
bifurcation curves shown in Figure 4.3.
The stability of the solutions is determined by the eigenvalues of the Jacobian, J,
computed for the four dimensional system with states u, u T = [ r1 , r2 , φ1 , φ2 ], where
ai = ri /2, the real amplitudes. The Jacobian describes the linear dynamics for small per˙ = Jδu, where
¯ specifically, δu
turbations about a given steady state operating point u,
¯ The characteristic equation of the Jacobian is Det(λI − J ) = 0, which is a
δu = u − u.
fourth order polynomial and has the following form
p ( λ ) = λ4 + p3 λ3 + p2 λ2 + p1 λ + p0

(.17)

The characteristic polynomial can be written in terms of the trace, tr ( J ), and determinant,
det( J ), as:
p0 = det( J )
p1 = −
p2 =

1
3

(.18)

3
tr ( J 3 ) − tr ( J 2 )tr ( J ) + tr ( J )3
2

1
tr ( J )2 − tr ( J 2 )
2

(.19)
(.20)
(.21)

p3 = −tr ( J )
92

When Re[λ j ] < 0 ∀ j, u¯ is locally stable. When there is one or more Re[λ j ] > 0, u¯ is unstable. The transition between stable and unstable conditions correspond to bifurcations of
¯
u.
The modal coupling gives rise to additional bifurcations, beyond the simple saddle
node bifurcations (SNB) that occur in the usual Duffing model. As one sweeps through a
SNB, the response jumps from the stable branch and will be attracted to another response,
generally far away in phase space. SNBs correspond to a zero eigenvalue and are thus
determined from the characteristic polynomial in Equation .17 using the condition:
p0 = det( J ) = 0

(.22)

Since this is SNB in a four dimensional phase space, we need to consider the stability
of the bifurcating response, since only SNBs from stable response branches will be observable in the experiments, although both types contribute to the overall bifurcation diagram. For a stable branch to exist the polynomial will need to satisfy the Routh-Hurwitz
criterion, which states that all pi ≥ 0, in addition to
∆3 = p1 ( p2 p3 − p1 ) − p23 p0 ≥ 0

(.23)

We are also interested in Hopf bifurcations (HB) since these occur in the open loop
system as it sweeps through the internal resonance region. For a super(sub)-critical HB
the steady state response branch will (de)stabilize and create a(n) (un)stable limit cycle.
The conditions for a HB can be determined by expanding the assumed form of the characteristic equation, (λ + Ω2H )(λ − λ1 )(λ − λ2 ) = 0, where Ω H is the Hopf frequency,
and where λ H = 0 + iΩ H , yielding conditions Ω H = p1 /p3 and ∆3 = 0. As for the SN
bifurcations, only HBs from stable responses will be visible in the experiment.

93

Mode 1 (71kHz)
Mode 2 (222kHz)
(!h)
(a)

(b)

Figure .1 (a) The mechanical device along with the two mechanical modes and the
eigenfrequencies as computed in COMSOL. (b) The layout with the red region being the
mechanical device, the dark green indicating the drive and sense electrodes for mode 1,
and the light green electrodes allow for modal tuning and coupling bias.

A.3

An Alternative 1 : 3 Internally Resonant Device Design

Figure .1 shows the layout of a flexural-flexural in-plane device with dynamic features
similar to those of the ANL device, namely, a 1:3 internal resonance and nonlinear modal
coupling. The linear mode shapes shown are for the mechanical structure without electrostatic transduction or tuning effects. These modes have distinct symmetries and are,
of course, uncoupled. The light green electrodes can apply different types of asymmetric
lateral bias forces, which result in in changes in the linear modes and thus in the nonlinear
coupling, as described above. In contrast to the bias and coupling in the flexural-torsional
device, in the proposed design thes effects are tunable. This is just one of several possible
designs that allow for controllable parameters in internal resonance. COMSOL simulations verify the existence of the tunable nonlinear coupling.

94

BIBLIOGRAPHY

95

BIBLIOGRAPHY
[1] Dynamics-enabled frequency sources. Technical report, DARPA, August 2009.
[2] Dario Antonio, Damián H Zanette, and Daniel López. Frequency stabilization in
nonlinear micromechanical oscillators. Nature Communications, 3:806, 2012.
[3] A. Chenakin. Frequency synthesis: Current solutions and new trends. Microwave
Journal, 50(5):256, 2007.
[4] Y.H. Chuang, S.H. Lee, R.H. Yen, S.L. Jang, J.F. Lee, and M.H. Juang. A wide locking
range and low voltage CMOS direct injection-locked frequency divider. Microwave
and Wireless Components Letters, IEEE, 16(5):299–301, 2006.
[5] AN Cleland and ML Roukes. Noise processes in nanomechanical resonators. Journal
of Applied Physics, 92(5):2758–2769, 2002.
[6] P. Collet and J.P. Eckmann.
Birkhäuser Boston, 1980.

Iterated maps on the interval as dynamical systems.

[7] MRM Crespo da Silva and CC Glynn. Nonlinear flexural-flexural-torsional dynamics of inextensional beams. i. equations of motion. Journal of Structural Mechanics,
6(4):437–448, 1978.
[8] MRM Crespo Da Silva and CL Zaretzky. Nonlinear flexural-flexural-torsional interactions in beams including the effect of torsional dynamics. i: Primary resonance.
Nonlinear Dynamics, 5(1):3–23, 1994.
[9] S. Daneshgar, O. De Feo, and M.P. Kennedy. Observations concerning the locking
range in a complementary differential lc injection-locked frequency divider: part i:
qualitative analysis. IEEE Transactions on Circuits and Systems Part I: Regular Papers,
57(1):179–188, 2010.
[10] P. Danzl and J. Moehlis. Weakly coupled parametrically forced oscillator networks:
existence, stability, and symmetry of solutions. Nonlinear Dynamics, 59(4):661–680,
2010.
[11] B.E. DeMartini, J.F. Rhoads, K.L. Turner, S.W. Shaw, and J. Moehlis. Linear and nonlinear tuning of parametrically excited mems oscillators. Journal of Microelectromechanical Systems, 16(2):310–318, 2007.
[12] Alper Demir. Phase noise and timing jitter in oscillators with colored-noise sources.
Circuits and Systems I: Fundamental Theory and Applications, IEEE Transactions on,
49(12):1782–1791, 2002.
[13] Alper Demir, Amit Mehrotra, and Jaijeet Roychowdhury. Phase noise in oscillators:
a unifying theory and numerical methods for characterization. Circuits and Systems
I: Fundamental Theory and Applications, IEEE Transactions on, 47(5):655–674, 2000.
96

[14] S. Dou, B.S. Strachan, S.W. Shaw, and J.S. Jensen. Structural optimization for nonlinear dynamic response. Philosophical Transactions of the Royal Society of London. Series
A, 2015.
[15] Suguang Dou and Jakob Søndergaard Jensen. Optimization of nonlinear structural
resonance using the incremental harmonic balance method. Journal of Sound and
Vibration, 334:239 – 254, 2015.
[16] Michael M. Driscoll. Phase noise performance of analog frequency dividers. IEEE
Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 37:295–301, 1990.
[17] MI Dykman and MA Krivoglaz. Classical theory of nonlinear oscillators interacting
with a medium. physica status solidi (b), 48(2):497–512, 1971.
[18] MI Dykman, R Mannella, Peter VE McClintock, SM Soskin, and NG Stocks. Noiseinduced narrowing of peaks in the power spectra of underdamped nonlinear oscillators. Physical Review A, 42(12):7041, 1990.
[19] M.J. Feigenbaum. Universal behavior in nonlinear systems. Physica D: Nonlinear
Phenomena, 7(1-3):16–39, 1983.
[20] J. Guckenheimer and P. Holmes. Nonlinear oscillations, dynamical systems, and bifurcations of vector fields, volume 42. Springer-Verlag New York, 1983.
[21] John Guckenheimer and Philip Holmes. Nonlinear oscillations, dynamical systems, and
bifurcations of vector fields. New York Springer Verlag, 1983.
[22] A. Hajimiri. Noise in phase-locked loops [invited]. Proc. SSMSD, 1:1–6, 2001.
[23] Robert G Harrison. A broad-band frequency divider using microwave varactors.
Microwave Theory and Techniques, IEEE Transactions on, 25(12):1055–1059, 1977.
[24] Zainabolhoda Heshmati, Ian C Hunter, and Roger D Pollard. Microwave parametric frequency dividers with conversion gain. Microwave Theory and Techniques, IEEE
Transactions on, 55(10):2059–2064, 2007.
[25] Ville Kaajakari et al. Practical mems: Design of microsystems, accelerometers, gyroscopes, rf mems, optical mems, and microfluidic systems. Las Vegas, NV: Small Gear
Publishing, 2009.
[26] Eyal Kenig, MC Cross, LG Villanueva, RB Karabalin, MH Matheny, Ron Lifshitz,
and ML Roukes. Optimal operating points of oscillators using nonlinear resonators.
Physical Review E, 86(5):056207, 2012.
[27] Christopher L Lee and Noel C Perkins. Nonlinear oscillations of suspended cables
containing a two-to-one internal resonance. Nonlinear Dynamics, 3(6):465–490, 1992.
[28] Thomas H Lee and Ali Hajimiri. Oscillator phase noise: a tutorial. Solid-State Circuits,
IEEE Journal of, 35(3):326–336, 2000.
97

[29] Li Li, Yujin Hu, and Xuelin Wang. Design sensitivity and hessian matrix of generalized eigenproblems. Mechanical Systems and Signal Processing, 43(1-2):272–294, 2014.
[30] Konstantin Konstantinovic Licharev. Dynamics of Josephson junctions and circuits. CRC
Press, 1986.
[31] R. Lifshitz and M. C. Cross. Response of parametrically driven nonlinear coupled oscillators with application to micromechanical and nanomechanical resonator arrays.
Physical Review B, 67, 2003.
[32] R. Lifshitz and M.C. Cross. Nonlinear dynamics of nanomechanical and micromechanical resonators. Review of Nonlinear Dynamics and Complexity, 1:1–52, 2008.
[33] Ron Lifshitz and Michael L Roukes. Thermoelastic damping in micro-and nanomechanical systems. Physical review B, 61(8):5600, 2000.
[34] Robert Lutwak, William J Riley Jr, and Kenneth D Lyon. Mems analog frequency
divider, August 12 2003. US Patent 6,605,849.
[35] Robert J Matthys. Crystal oscillator circuits. New York, Wiley-Interscience, 1983, 244
p., 1, 1983.
[36] Nick Miller. Noise in nonlinear micro-resonators, 2012.
[37] N.J. Miller and S.W. Shaw. Frequency sweeping with concurrent parametric amplification. ASME, 2008.
[38] DT Mook, RH Plaut, and N HaQuang. The influence of an internal resonance on
non-linear structural vibrations under subharmonic resonance conditions. Journal of
Sound and Vibration, 102(4):473–492, 1985.
[39] A. H. Nayfeh. Nonlinear Interactions: Analytical, Computational, and Experimental Methods. Wiley Series in Nonlinear Science. Wiley-Interscience, New York, 2000.
[40] A. H. Nayfeh and D. T. Mook. Nonlinear Oscillations. Wiley-Interscience, 1979.
[41] Ali H Nayfeh and Dean T Mook. Nonlinear oscillations. Wiley. com, 2008.
[42] TA Nayfeh, W. Asrar, and AH Nayfeh. Three-mode interactions in harmonically excited systems with quadratic nonlinearities. Nonlinear Dynamics, 3(5):385–410, 1992.
[43] CT-C Nguyen. Mems technology for timing and frequency control. Ultrasonics, Ferroelectrics and Frequency Control, IEEE Transactions on, 54(2):251–270, 2007.
[44] Sarah H Nitzan, Valentina Zega, Mo Li, Chae H Ahn, Alberto Corigliano, Thomas W
Kenny, and David A Horsley. Self-induced parametric amplification arising from
nonlinear elastic coupling in a micromechanical resonating disk gyroscope. Scientific
reports, 5, 2015.

98

[45] Niels L. Pedersen. Design of cantilever probes for atomic force microscopy (afm).
Engineering Optimization, 32(3):373–392, 2000.
[46] KR Qalandar, BS Strachan, B Gibson, M Sharma, A Ma, SW Shaw, and KL Turner.
Frequency division using a micromechanical resonance cascade. Applied Physics Letters, 105(24):244103, 2014.
[47] J. F. Rhoads, S. W. Shaw, K. L. Turner, J. Moehlis, B. E. DeMartini, and W. Zhang.
Generalized parametric resonance in electrostatically actuated microelectromechanical oscillators. Journal of Sound and Vibration, 296(4-5):797–829, 2006.
[48] Jeffrey F. Rhoads, Steven W. Shaw, and Kimberly L. Turner. Nonlinear dynamics and
its applications in micro- and nanoresonators. Journal of Dynamic Systems, Measurement, and Control, 132(3):034001, 2010.
[49] Hannes Risken. Fokker-Planck Equation. Springer, 1984.
[50] Michael Rosenblum and Arkady Pikovsky. Synchronization: from pendulum clocks
to chaotic lasers and chemical oscillators. Contemporary Physics, 44(5):401–416, 2003.
[51] Kaushik Sengupta, TK Bhattacharyya, and Hossein Hashemi. A nonlinear transient
analysis of regenerative frequency dividers. Circuits and Systems I: Regular Papers,
IEEE Transactions on, 54(12):2646–2660, 2007.
[52] Jinsiang Shaw, Steven W Shaw, and Alan G Haddow. On the response of the nonlinear vibration absorber. International Journal of Non-Linear Mechanics, 24(4):281–293,
1989.
[53] B. Scott Strachan, Steven W. Shaw, and Oleg Kogan. Subharmonic resonance cascades in a class of coupled resonators. Journal of Computational and Nonlinear Dynamics, 8(4):041015, 2013.
[54] B Scott Strachan, Steven W Shaw, and Oleg Kogan. Subharmonic resonance cascades
in a class of coupled resonators. Journal of Computational and Nonlinear Dynamics,
8(4):041015, 2013.
[55] Rouslan L Stratonovich. Topics in the theory of random noise, volume 2. CRC Press,
1967.
[56] Krister Svanberg. The method of moving asymptotes—a new method for structural
optimization. International Journal for Numerical Methods in Engineering, 24(2):359–
373, 1987.
[57] Simon M Sze and Kwok K Ng. Physics of semiconductor devices. Wiley. com, 2006.
[58] M. Tiebout. A CMOS direct injection-locked oscillator topology as high-frequency
low-power frequency divider. Solid-State Circuits, IEEE Journal of, 39(7):1170–1174,
2004.

99

[59] Daniel A Tortorelli and Panagiotis Michaleris. Design sensitivity analysis: overview
and review. Inverse problems in Engineering, 1(1):71–105, 1994.
[60] C Touzé, O Thomas, and Antoine Chaigne. Hardening/softening behaviour in nonlinear oscillations of structural systems using non-linear normal modes. Journal of
Sound and Vibration, 273(1):77–101, 2004.
[61] C. Touzé, M. Vidrascu, and D. Chapelle. Direct finite element computation of nonlinear modal coupling coefficients for reduced-order shell models. Computational
Mechanics, pages 1–14, 2014.
[62] K.L. Turner, P.G. Hartwell, and N.C. MacDonald. Multi-dimensional mems motion
characterization using laser vibrometry. In IEEE Transducers Digest of Technical Papers,
pages 1144–1147, 1999.
[63] K.L. Turner, S.A. Miller, P.G. Hartwell, N.C. MacDonald, S.H. Strogatz, and S.G.
Adams. Five parametric resonances in a microelectromechanical system. Nature,
396(6707):149–152, 1998.
[64] M Urteaga, M Seo, J Hacker, Z Griffith, A Young, R Pierson, P Rowell, A Skalare,
and MJW Rodwell. Inp hbt integrated circuit technology for terahertz frequencies.
In Compound Semiconductor Integrated Circuit Symposium (CSICS), 2010 IEEE, pages
1–4. IEEE, 2010.
[65] C Van der Avoort, R Van der Hout, JJM Bontemps, PG Steeneken, K Le Phan, RHB
Fey, J Hulshof, and JTM van Beek. Amplitude saturation of mems resonators explained by autoparametric resonance. Journal of Micromechanics and Microengineering,
20(10):105012, 2010.
[66] John R Vig. Quartz crystal resonators and oscillators for frequency control and timing applications. a tutorial. NASA STI/Recon Technical Report N, 95:19519, 1994.
[67] L Guillermo Villanueva, Rassul B Karabalin, Matthew H Matheny, Eyal Kenig,
Michael C Cross, and Michael L Roukes. A nanoscale parametric feedback oscillator.
Nano letters, 11(11):5054–5059, 2011.
[68] Ashwin Vyas, Dimitrios Peroulis, and Anil K Bajaj. Dynamics of a nonlinear microresonator based on resonantly interacting flexural-torsional modes. Nonlinear Dynamics, 54(1-2):31–52, 2008.
[69] Ashwin Vyas, Dimitrios Peroulis, and Anil K Bajaj. A microresonator design based
on nonlinear 1: 2 internal resonance in flexural structural modes. Microelectromechanical Systems, Journal of, 18(3):744–762, 2009.
[70] Carl O Weiss and Ramon Vilaseca. Dynamics of lasers. NASA STI/Recon Technical
Report A, 92:39875, 1991.

100

[71] B Yurke, DS Greywall, AN Pargellis, and PA Busch. Theory of amplifier-noise evasion in an oscillator employing a nonlinear resonator. Physical Review A, 51(5):4211,
1995.
[72] N.J. Miller Z. Yie, K.L. Turner and S.W. Shaw. Sensitivity enhancement using parametric amplification in a resonant sensing array. Hilton Head, 2010: A Solid State
Sensor, Actuator, and Microsystems Workshop, 2010.
[73] Xiangdong Zhang, Xuesong Zhou, B Aliener, and AS Daryoush. A study of subharmonic injection locking for local oscillators. Microwave and Guided Wave Letters, IEEE,
2(3):97–99, 1992.

101