of 5
Nonlinear dynamics and chaos in two coupled nanomechanical resonators
R. B. Karabalin, M. C. Cross, and M. L. Roukes
Department of Physics and Kavli Nanoscience Institute, California Institute of Technology, Pasadena, California 91125, USA

Received 10 February 2009; published 13 April 2009

Two elastically coupled nanomechanical resonators driven independently near their resonance frequencies
show intricate nonlinear dynamics. The dynamics provide a scheme for realizing a nanomechanical system
with tunable frequency and nonlinear properties. For large vibration amplitudes, the system develops sponta-
neous oscillations of amplitude modulation that also show period-doubling transitions and chaos. The complex
nonlinear dynamics are quantitatively predicted by a simple theoretical model.
DOI:
10.1103/PhysRevB.79.165309
PACS number

s

: 85.85.

j, 05.45.

a, 62.25.

g
I. INTRODUCTION
Resonant nanoelectromechanical systems

NEMS

1
are
attracting interest in a broad variety of research areas and for
many possible applications due to their remarkable combina-
tion of properties: small mass, high operating frequency,
large quality factor, and easily accessible nonlinearity.
2
Fur-
ther development of NEMS applications, such as superior
mass,
3
force,
4
and charge
5
sensors, or reaching the quantum
limit of detection in mechanical systems,
6
requires address-
ing several important challenges. For example, the nonlinear-
ity of the devices has to be either minimized or utilized for
improving performance.
7
In addition, the large-scale integra-
tion of nanodevices demands a detailed understanding of the
behavior of coupled devices in NEMS arrays.
8
11
In this paper we demonstrate the complex nonlinear be-
havior of a pair of coupled nanomechanical devices that can
be quantitatively understood from the basic physics of the
devices. We show that the linear and weakly nonlinear re-
sponse of one oscillation can be modified by driving the
second oscillation and that the linear-response range of the
first oscillation can be significantly extended. When both os-
cillations are driven into their strongly nonlinear range, more
complicated frequency-sweep response curves are found,
corresponding to the well-known bistability of driven anhar-
monic “Duffing” resonators, but now with switching be-
tween a variety of different stable states of the coupled pair
not seen previously. Spontaneous amplitude modulation os-
cillations may develop, with frequencies characteristic of the
dissipation rates rather than of the intrinsic frequencies or
their sums and differences. These amplitude modulations
show period-doubling bifurcations and chaos, giving a clear
demonstration of these complex dynamics in the context of
nanomechanical systems.
12
The complex dynamics are repro-
duced quantitatively by a simple theoretical model.
II. EXPERIMENTAL SETUP
We study a system of two strongly coupled nonlinear na-
noelectromechanical resonators using a structure of doubly
clamped beams with a shared mechanical ledge shown in
Fig.
1

a

. The devices consist of a stack of three layers of
gallium arsenide

GaAs

: a 100 nm highly
n
-doped layer, a
50 nm insulating layer, and another 50 nm layer that is
highly
p
doped. The piezoelectric property of GaAs results in
a highly efficient integrated actuation mechanism described
in Ref.
13
. A preliminary 120-nm-deep etch step is done to
isolate the actuation electrodes of the two beams, so that the
two beams can be addressed separately while retaining
strong elastic coupling.
Optical interferometry is used for the motion trans-
duction.
14
The GaAs heterostructures are epitaxially grown
on a nearly perfectly lattice-matched sacrificial layer of
Al
0.8
Ga
0.2
As, which can be removed with a high degree of
selectivity to reveal a smooth substrate that forms a mirror.
The top surface of the NEMS device forms a second mirror.
Coupled together, the mirrors produce a primitive Fabry-
Pérot interference cavity. The reflection coefficient for light
that is normally incident on a vacuum-GaAs interface is
about 33% and because of the relatively low light intensity, it
is sufficient to consider a single reflection from each surface.
An infrared laser diode

904 nm

followed by an attenuator
is used as the light source. Its frequency is below the band
gap of GaAs

830 nm

, in order to minimize heating and
photocurrent generation. The beam is focused onto the nano-
0
5
10
15
20
25
3
0
0
50
100
150
200
|
A
II
|
(nm)
Frequency Shift (kHz)
experiment
theory
3m
(a)
(c)
(b)
-5 0 5 10 15 20
0
1
2
3
4
|A|
(ω−ω
0
)
/
γ
up
down
FIG. 1.

Color online

a

SEM image of the system

beam
dimensions: 6

m

500 nm

200 nm

.

b

Up

dark, blue

and
down

light, yellow

frequency sweeps of the amplitude

A

of the
response for a single drive of strength 4.3 times critical as a func-
tion of the frequency relative to the linear resonance frequency

0
and scaled by the width

. The amplitude is plotted in units of the
maximum amplitude at the critical drive strength.

c

Frequency
shift of the weakly driven first mode as a function of the displace-
ment amplitude

A
II

of the strongly driven second mode: points and
solid line—experiment; dashed line—small amplitude theory.
PHYSICAL REVIEW B
79
, 165309

2009

1098-0121/2009/79

16

/165309

5

©2009 The American Physical Society
165309-1
mechanical device with a lens of numerical aperture 0.15,
giving a spot size of approximately 10

m. The laser beam
is adjusted so that both beams are in the illuminated spot.
After reflection from the top surface of the mechanical reso-
nator and the bottom mirror, the beam intensity is modulated
proportional to the mechanical displacement. A high-
bandwidth low-noise photodetector is used to measure the
intensity of the light and thus the mechanical displacement of
the resonators.
III. THEORETICAL MODEL
We model the behavior of the two strongly interacting
nonlinear resonators by a system of coupled equations of
motion for the beam displacements
x
1
,
x
2
in their fundamen-
tal modes
x
̈
1
+

1
x
̇
1
+

1
2
x
1
+

1
x
1
3
+
D

x
1
x
2

=
g
D
1

t

,

1

x
̈
2
+

2
x
̇
2
+

2
2
x
2
+

2
x
2
3
+
D

x
2
x
1

=
g
D
2

t

.

2

As well as the usual terms describing the resonant frequen-
cies, damping, and Duffing nonlinearity, we include a linear
coupling term in the displacements of strength
D
. The terms
on the right-hand side are the external drives applied to the
two beams, which are controlled independently. The linear
terms in the equations, ignoring for now the drive and dissi-
pation, give two modes with frequencies

I
and

II
and cor-
responding eigenvectors
e
I
,
e
II
.
15
The frequency difference of
the modes results from both the intrinsic frequency differ-
ence

1
-

2
and the coupling
D
.
To investigate the behavior of the system, the two beams
are connected to two different sources with independent fre-
quencies and amplitudes. The monitored output variables are
the amplitudes and phases of a linear combination of the
mechanical displacements of the two beams at or near the
two drive frequencies.
16
By applying various combinations
of small amplitude signals to the two beams at frequencies
near the mode resonances, the linear coupling parameters can
be determined. For the particular device shown in Fig.
1

a

,
the mode frequencies are determined to be

I
/
2

= 16.79 MHz and

I
/
2

= 17.25 MHz and the eigenvectors
e
I
=

0.854 , 0.521

and
e
II
=

−0.521 , 0.854

. The frequencies
of the individual resonators determined from inverting the
mode equations are

1
/
2

= 16.71 MHz and

2
/
2

= 16.92 MHz. The frequency separation of 200 kHz is con-
sistent with the fabrication tolerance. The coupling strength

D
/
2

= 2.63 MHz is consistent with the strength of the
elastic coupling found by finite element simulations of simi-
lar devices.
17
Transforming the measured width of the modes
back to the original equations determines the values of

1
,

2
. The nonlinear parameters

1
,

2
are deduced from the
expression for the geometric nonlinearity
7
using the beam
thickness known from the fabrication and lengths calculated
from the beam frequencies

1
,

2
and the material
constants.
18
The response of the system driven near resonance and for
small dissipation and driving can be calculated from Eqs.

1

and

2

using the standard methods of secular perturbation
theory.
2
This approach has previously used for the case of
parametrically driven nanomechanical devices,
9
11
where
hysteretic switches between different stable states in fre-
quency sweeps were also predicted. Briefly, we introduce the
slowly varying complex mode amplitudes
A
I
,
A
II
and forces
F
I
,
F
II
using
x

I

=Re

A

I

e
i


I

t

and
g
D

I

=Re

F

I


t

e
i


I

t


where

I

stands for either I or II

, substitute into the equa-
tions of motion, and retain only the near resonant terms. This
reduces the equations of motion to
2
i

I
A
̇
I
+
i

I

I
A
I
+

I

A
I

2
A
I
+
I

A
II

2
A
I
=
F
I

t

,

3

and a corresponding equation for
A
II
. The mode nonlinearity
parameters


I

and

I

are calculated from

1
,

2
and the
eigenvectors
e
I
,
e
II
so that all the parameters in Eq.

3

are
known from linear measurements and the beam geometry
and material constants.
IV. RESULTS
We first look at the case where each mode responds at the
drive frequency

D

I

which is set near the resonant fre-
quency


I

so that
A

I

e
i


D

I



I


t
, with

A
I

2
=

F
I

2

2

I


D
I

I


I

A
I

2
I

A
II

2

2
+

I
2

I
2
,

4

etc. For a single drive

e.g.,
A
II
=
F
II
=0

, so that the cross-
mode nonlinear coupling proportional to

I

is not involved,
this expression reproduces the regular Duffing response
curve.
19
Prominent features are the shift of the frequency of
the maximum response to larger values for positive


non-
linear spring stiffening

and bistability and hysteresis that
develop above a critical drive strength. An experimental ex-
ample of upward and downward frequency sweeps for a
drive strength 4.3 times critical is shown in Fig.
1

b

.
We first look at the linear and weakly nonlinear responses
of the first mode when the second mode is driven into its
strongly nonlinear regime. Retaining the leading-order effect
of the first mode on the intensity of the second mode, we can
write the numerator in Eq.

4

as

2

I


D
I

I

I

A
II

0


2

̄
I

A
I

2

2
+

I
2

I
2
,

5

where
A
II

0

is the solution for the second mode in the absence
of mode I given by the equation for
A
II
corresponding to Eq.

4

assuming zero
A
I
. The parameter

̄
I
is an effective Duf-
fing nonlinearity coefficient and is given by

̄
I
=

I
−2
I
II

II


A
II

0


2
/


D
II
.

6

The expression

5

predicts two important effects: the fre-
quency tuning of the first mode proportional to the square of
the amplitude of the second mode

upwards for positive
I

and the change in the effective nonlinear coefficient for the
motion of the first mode depending on the excitation strength
of the second mode through the last term in

̄
I
Eq.

6

.
To test the frequency tuning, we excite the second mode
at a drive level approximately 4.3 times the critical value so
that the spectral response is the strongly nonlinear Duffing
curve

see Fig.
1

b


. As the actuation frequency of the sec-
KARABALIN, CROSS, AND ROUKES
PHYSICAL REVIEW B
79
, 165309

2009

165309-2
ond mode is steadily increased in small steps, its vibration
amplitude rises over a wide frequency range until it drops to
the lower amplitude state beyond the maximum. The evolu-
tion of the spectral response of the first mode is monitored at
a driving level approximately four times lower than the criti-
cal value for this mode using a network analyzer. The depen-
dence of the first mode frequency shift on the vibration am-
plitude of the second mode is shown in Fig.
1

c

. The
experimental results for frequency tuning closely follow the
predicted parabolic dependence
I

A
II

2
for amplitudes up to
about 20 nm.
If the actuation level of the first mode is increased above
the onset of nonlinearity, while the second mode drive is kept
at a much higher level then the effective nonlinear coefficient

̄
I
is decreased. The decrease is largest when the second
mode is driven on the portion of the Duffing response curve
where the intensity is increasing linearly with the drive fre-
quency

A
II

2

2

II


D
II

II

/

II
. This gives the minimum
effective nonlinear coefficient

̄
I,min
=

I

1−
I
II

I

II
.

7

This result indicates that if the coupling is strong enough, so
that
I
II

I

II
, the minimum value of the effective nonlin-
ear coefficient is negative. In this case the resonance curve
tilts to the left, as opposed to the usual case for a doubly
clamped beam where the peak leans to the right. It also
means that the nonlinear coefficient vanishes for some drive
strength and frequency of the second mode.
Some experimental results for the first mode driven at
twice the critical strength illustrating this effect are shown in
Fig.
2
. The plots show the shape of the resonance peak for
three values of the drive frequency of the second mode. In
panel

a

, the amplitude
A
II
of the second mode is low and
the first mode spectral response has the regular nonlinear
Duffing shape leaning to the right. For larger
A
II
as in

b

, the
first mode resonance peak shape assumes a form close to a
Lorentzian with little nonlinearity apparent. For even larger
A
II
as in

c

, the sign of the effective Duffing coefficient
becomes slightly negative, causing the spectral response
peak to lean to the left. The quenching of the nonlinearity in

b

could be used to enhance the limited dynamic range
7
of
nanomechanical devices.
If the drive level of the first mode is increased further, a
variety of new effects can be observed. Under these condi-
tions the dynamics of the system is not fully explained by the
steady-state solutions as in Eq.

4

, and so we solve for the
expected behavior by numerically integrating the time-
dependent coupled Eq.

3

and the corresponding equation
for
A
II
. Now the behavior is more complex, as the response
of the first mode becomes large enough to cause transitions
in the second mode response through the nonlinear coupling.
As a result, the spectral response curves acquire peculiar
nontrivial shapes, as shown in Fig.
3
. Numerical simulations
of the equations give good predictions for the complex phe-
nomena.
Figure
3
shows frequency sweeps of the first mode reso-
nance for four increasing drive frequencies of the second
mode.

The second mode response was monitored, but is not
shown in the figure.

The top four plots are the experimental
measurements and the bottom four are the numerical simu-
lations. In panel

a

for experiment and theory, the second
mode response is small, and the sweep of the first mode
shows approximately the hysteretic response of a single Duf-
fing oscillator. The interaction with the second mode is more
important in the succeeding panels. For panels

b

and

c

,on
the upward sweep the amplitude of the second mode is ini-
tially large, which shifts the first mode response to larger
frequencies, and gives a negative effective nonlinearity lead-
ing to a jump of the first mode to a larger amplitude as the
frequency is swept up. This large amplitude is then sufficient
to cause the second mode to transition to a lower amplitude,
so that the remainder of the first mode sweep is close to the
single mode Duffing curve. The downward sweeps in panels

b

and

c

are different depending on whether the second
mode transitions between small and large amplitude states.
In panel

d

, the second mode amplitude is initially small on
the upward sweep, and the first mode follows the single os-
cillator Duffing curve. This continues until the second mode
16.75 16.80
0
2
4
6
8
10
Displacement (nm)
16.75 16.80
Frequency (MHz)
16.75 16.80
(b)
(a)
(c)
FIG. 2.

Color online

Frequency sweeps of the first mode re-
sponse for increasing amplitudes of the second mode: the dark

blue

lines are for upward sweeps; the light

yellow

line in

a

is a
downward sweep, showing the hysteresis and amplitude jumps in
this case. There are no hysteretic jumps for

b

and

c

, and so the
downward sweeps are not shown.
(c)
(a)
(b)
(d)
16.7 16.8 16.9
0
5
10
15
20
25
16.7 16.8 16.9
16.7 16.8 16.9
16.7 16.8 16.
9
0
5
10
15
20
2
5
Frequenc
y(
MHz
)
Displacement (nm)
FIG. 3.

Color online

First mode frequency response as in Fig.
2
for increasing values of the second mode drive frequency but for
stronger driving of the first mode

both modes are driven at about
four times the critical strength

. The top plots are experimental mea-
surements while the bottom ones are theoretical simulations. The
second mode response was monitored but is not shown.
NONLINEAR DYNAMICS AND CHAOS IN TWO COUPLED
...
PHYSICAL REVIEW B
79
, 165309

2009

165309-3
transitions to a large amplitude so that the first mode re-
sponse jumps down and follows a curve close to the one in
panel

c

. The simple nonlinear model explains the complex
behavior of the coupled system both qualitatively and quan-
titatively with a precision of better than 5%.
An obvious difference between the experimental and the-
oretical plots in Fig.
3
is the noisy regions on the theoretical
curves near the up and down transitions. This difference ac-
tually results from the different ways the plots are generated
in theory and experiment, and a more careful investigation
shows consistent and interesting dynamics in both experi-
ment and theory. It turns out that for the drive frequencies
near the transition points, the fast

about 17 MHz

oscillating
response becomes amplitude modulated with a frequency of
about 10–20 kHz.
20
The theoretical simulations show that the
frequency for the amplitude modulation is determined by the
line width, which in our experimental setup is about 8 kHz,
and is not related to sum and difference frequencies of the
two modes. Examples of simulated and measured amplitude
modulation dynamics for three different sets of drive param-
eters are shown in Fig.
4
for a second device with different
parameters to the one used for Figs.
1
3
. The first column
shows examples of numerically calculated phase portraits of
Re
A
II
versus Re
A
I
obtained by solving the time-dependent
Eq.

3

. Phase portraits measured experimentally are shown
in the second column. These are obtained using homodyne
down conversion of the transduced mechanical signals from
both modes, which are then read by independent oscilloscope
channels. We also perform wider frequency band spectrum
analyzer measurements shown in the third column. Since the
nature of the dynamics is sensitive to the precise values of
the system parameters, the drive strength and frequencies in
the experiment are slightly adjusted from the values corre-
sponding to the theoretical plots to produce comparable
phase portraits.
The top row of Fig.
4
shows a relatively simple example
where the motion is periodic but the
A
I
-
A
II
trajectory forms a
small figure-of-eight loop. The dynamics can be roughly un-
derstood in terms of transitions between states with the first
mode at large amplitude and the second mode at small am-
plitude and vice versa. The spectrum of the measured me-
chanical signal shows satellite peaks corresponding to the
anharmonic amplitude modulation.
As the parameters of the system are changed, we have
observed period doubling or quadrupling in the amplitude
modulation, where the phase trajectory takes two or four
revolutions in order to complete the cycle.
21
An example of
period quadrupling is shown in the middle row of Fig.
4
. The
complicated loop structures are visible in both theoretical
and experimental trajectories, and the spectrum reveals am-
plitude modulation peaks at frequencies that correspond to
double and quadruple periods. Period-doubling transitions
are often associated with chaotic dynamics,
21
and indeed for
other parameter values we observe chaos in our coupled na-
-2 -1 0 1 2 3
-4
-2
0
2
4
Re(
A
I
)
-4
-2
0
2
4
-4
-2
0
2
4
6
Re(
A
I
)
14.70
14.75
14.80
0
10
20
30
Signal
V
)
Frequency (MHz)
14.70
14.75
14.80
0
10
20
30
Signal
V
)
Frequency (MHz)
-4 -2 0 2 4 6
-2
-1
0
1
2
3
Re(
A
II
)
Re(
A
I
)
14.70
14.75
14.80
0
10
20
Signal
V
)
Frequency (MHz)
Re(
A
II
)
Re(
A
I
)
Re(
A
I
)
Re(
A
I
)
Re
(
A
II
)
Re(
A
II
)
Re(
A
II
)
Re(
A
II
)
2
0
-1
1
0
2
4
1
2
3
0
1
0
1
3
2
2
1
3
4
FIG. 4. Complex dynamics of two strongly coupled nanomechanical resonators. The three rows correspond to different input parameters

drive frequencies and amplitudes

. The first column shows the theoretical calculation of the phase portrait, the second shows its experi-
mental measurement, and the third column shows the corresponding experimental power spectrum of the optical measurement near one of
the drive frequencies.
KARABALIN, CROSS, AND ROUKES
PHYSICAL REVIEW B
79
, 165309

2009

165309-4
noelectromechanical system and in the theoretical model, as
shown in the bottom row of Fig.
4
. The evidence for the
chaotic dynamics in the experiment is the broad band com-
ponent to the spectrum

evident in the shoulders to the am-
plitude modulation peaks

and a phase portrait trajectory that
does not form a closed loop. The theoretical model shows a
similar phase portrait.
V. CONCLUSIONS
The complex dynamics we observe in coupled nanome-
chanical beam resonators have a number of potential appli-
cations. For example, driving one of the modes can be used
to tune the effective nonlinearity of the other mode. This can
be used to significantly increase the dynamic range of the
resonator by quenching the effective nonlinearity. In a first
approximation, the motion of one mode couples quadrati-
cally to the resonance frequency of the other mode: a phe-
nomenon that has been proposed for quantum nondemolition
measurements in nanomechanical systems.
22
,
23
The full
range of complex dynamics investigated is quantitatively re-
produced by theory, giving us confidence that the nonlinear
behavior of coupled nanomechanical devices can be under-
stood, controlled, and exploited. Applying these methods to
large arrays of nanomechanical resonators that can be fabri-
cated in the near future will enable the quantitative investi-
gation of fundamental issues in the nonlinear dynamics of
large numbers of coupled devices.
ACKNOWLEDGMENTS
Many of the ideas leading to this work were developed in
collaboration with Ron Lifshitz supported by the U.S.-Israel
Binational Science Foundation

BSF

through Grant No.
2004339. We would like to thank Iwijn De Vlaminck and
Gustaaf Borghs from IMEC

Leuven, Belgium

for providing
us with the GaAs material. We thank Matt Matheny for many
useful discussions.
1
M. L. Roukes, Phys. World
14
,25

2001

.
2
For a review, see R. Lifshitz and M. C. Cross,
Reviews of Non-
linear Dynamics and Complexity
, edited by H. G. Shuster

Wiley, Weinheim, 2008

, Chap. 1, p. 52.
3
Y. T. Yang, C. Calegari, X. L. Feng, K. L. Ekinci, and M. L.
Roukes, Nano Lett.
6
, 583

2006

.
4
D. Rugar, R. Budakian, H. J. Mamin, and B. W. Chui, Nature

London

430
, 329

2004

.
5
A. N. Cleland and M. L. Roukes, Nature

London

392
, 160

1998

.
6
M. D. Lahaye, O. Buu, B. Camarota, and K. C. Schwab, Science
304
,74

2004

.
7
H. W. C. Postma, I. Kozinsky, A. Husain, and M. L. Roukes,
Appl. Phys. Lett.
86
, 223105

2005

.
8
E. Buks and M. L. Roukes, Europhys. Lett.
54
, 220

2001

.
9
R. Lifshitz and M. C. Cross, Phys. Rev. B
67
, 134302

2003

.
10
Y. Bromberg, M. C. Cross, and R. Lifshitz, Phys. Rev. E
73
,
016214

2006

.
11
E. Kenig, R. Lifshitz, and M. C. Cross, arXiv:0808.3589

unpub-
lished

.
12
Previous experiments by Scheible
et al.

D. V. Scheible, A. Erbe,
R. H. Blick, and G. Corso, Appl. Phys. Lett.
81
, 1884

2002

,
that suggested chaos in a nanomechanical system demonstrated
a complex behavior of the response at one of the drive frequen-
cies as that parameter was swept, rather than complex dynamics
for fixed system parameters as shown in our work.
13
S. C. Masmanidis, R. B. Karabalin, I. Vlaminck, G. Borghs, M.
R. Freeman, and M. L. Roukes, Science
317
, 780

2007

.
14
D. W. Carr and H. G. Craighead, J. Vac. Sci. Technol. B
15
,
2760

1997

.
15
H. Goldstein, C. Poole, and J. Safko,
Classical Mechanics
, 3rd
ed.

Addison-Wesley, New York, 2001

.
16
The combination measured depends on the geometry of the op-
tical spot relative to the beams, which we determine from the
linear experiments using a variety of drive combinations.
17
R. B. Karabalin, Ph.D. dissertation, Caltech, 2008.
18
The structure of the dynamical system depends on the ratio

I
/

II
which is close to unity for our nearly identical beams. The
absolute magnitude of


approximately 0.007 2

MHz
/
nm

2

determines an overall proportionality constant in the driving
strengths and calibrates the amplitude of the response in nano-
meter.
19
A. H. Nayfeh and D. T. Mook,
Nonlinear Oscillations

Wiley,
New York, 1979

.
20
The theoretical points, which are the end values of a long nu-
merical simulation, depend on the phase of the amplitude modu-
lation at the end of each run, whereas this modulation simply
leads to a drop in the experimentally measured amplitude since
RF power is transferred outside the measurement bandwidth of
the network analyzer.
21
S. H. Strogatz,
Nonlinear Dynamics and Chaos

Addison-
Wesley, New York, 1994

.
22
K. C. Schwab and M. L. Roukes, Phys. Today
58

7

,36

2005

.
23
V. B. Braginsky and F. Y. Khalili,
Quantum Measurement

Cam-
bridge University Press, Cambridge, 1995

.
NONLINEAR DYNAMICS AND CHAOS IN TWO COUPLED
...
PHYSICAL REVIEW B
79
, 165309

2009

165309-5