www.sciencemag.org/content/36
3/6431/eaav7932
/suppl/DC1
Supplementary
Material
s for
Exotic states in a simple network of nanoelectromechanical oscillators
Matthew H.
Matheny
, Jeffrey
Emenheiser
, Warren
Fon
, Airlie
Chapman,
Anastasiya
Salova
,
Martin
Rohden, Jarvis
Li, Mathias
Hudoba de Badyn, Márton Pósfai
, Leonardo Duenas
-Osorio
,
Mehran
Mesbahi
, James P.
Crutchfield
, M. C.
Cross
, Raissa M.
D’Souza
, Michael L.
Roukes
*
*Corresponding author.
Email:
roukes@caltech.edu
Published
8 March
201
9,
Science
363
, eaav7932
(201
9)
DOI:
10.1126/science.
aav7932
This PDF file includes:
Materials and Methods
Supplementary Text
Figs. S1 to S6
Tables S1 to S3
References
Other supplementary material for this manuscript in
cludes:
Movies S1 to S32 (as identified in Tables S1 to S3)
2
Materials and Methods
NEM
S
resonator and oscillator
A ring network of 8 nanomechanical oscillators was constructed. In each oscillator, a
nanoelectromechanical resonator
was
bonded onto a printed circuit board hosting the
control and feedback electronics. A scanning electron micro
graph of the
resonator
(with
false color added for clarification of the device structures)
is
shown in Fig.
S1A
, top. It
is
an ultrathin fully clamped plate with an embedded aluminum nitride piezoelectric film
(
30
,
61
)
. Each plate
had
a resonant frequency s
et to
푓
"
=
2
.
21000
±
.
00002
MHz with a
quality factor
푄
"
=
3950
±
550
. The resonator clamping g
ave
a geometric Duffing
-
type
nonlinearity (
37
) and the capability of linear voltage
-
frequency tuning
(62
). The resonators
were
mounted inside a vacuum chamber and
wirebonded onto a printed circuit board (PCB)
(Fig.
S
1
B
) to form the self
-
sustained oscillator. Fon, et al. (
30
) give
s
the details for the
oscillator circuit. Tuning of resonator frequency (and
therefore
natural oscillator frequency
휔
/
) and nonlinearity
훼
we
re
done with an on
-
board 16
-
bit ultrastable Digital
-
to
-
Analog
-
Converter (DAC) (white box in
Fig S1B
) and a 7
-
bit variable attenuator (blue box in
Fig.
S1B),
respectively.
Network construction
The oscillator
PCB
was
connected through
an
edge connecto
r
to the network (purple box
in
Fig. S1B
). The network
PCB connect
ed
the 8 oscillators into a ring topology so that each
oscillator
was
coupled to their two nearest neighbors
(Fig.
S
1
C
). The network PCB
had
variable attenuators controlling the
strength (an
d to ensure the proper phase
shift through
the oscillator loop
) of
coupling
훽
along
the
edges (yellow box in
Fig. S1D
). Power and
digital control signals
were
directly provided
to the network PCB
and routed to the
oscillator boards
. Even though our contro
ls
were
adjusted via digital commands, signals
passing between the nodes and edges
were
analog signals and never digitized, and our
network evolves according to a continuous time dynamic.
Data acquisition, control, and calibration
Python scripts were us
ed for taking and analyzing data. Python libraries were also used for
numerical simulations.
A personal computer (PC) was connected to a digital controller (Raspberry Pi) (dark blue
box in Fig. S1D) via Ethernet. Python scripts used C libraries within t
he
controller
which
switch
e
d
settings within digital potentiometers, variable attenuators, and DACs within our
network
.
The
electronic components (DACs and attenuators)
control
led
the natural
frequencies
휔
/
, nonlinearities
훼
, and couplings strengths
훽
of the nodes and edges. Every
node and edge
were
comprised of their own components which
could
be individually
addressed by the controller. In addition, each node
was
readout
in an individual channel of
a simultaneous
-
sampling
8
-
channel oscilloscope.
In
Eq. 2
b
(main text), when the
j
th
oscillator
was
isolated such that
푎
/
3
4
,
푎
/
6
4
=
0
, we see
that changes in
훼
,
훽
, and
휔
/
are
associated with changes in
its
frequency. Each node's
3
parameters
were
individually calibrated with measurement of the isolated
oscillator
frequency for a given choice of parameters. The scaling of these parameters
was
either
linear or quadratic in PCB signals. Quadratic and linear fits yield
ed
a simple calibration
for the parameters; these fits eliminate
d
the need for lookup table
s. We refer the reader to
Matheny, et al. (14) for more details on the calibration.
In order to maintain throughput in data acquisition and storage, the
‘
fast
’
oscillator signals
(at
~
2
MHz)
were
captured with full fidelity in the system calibration but subsequently
undersampled at
200
kHz. Note that since the slow time scale
was
휏
9:;<
≈
>
?
@
?
≈
2
ms
in
our experiment
, our sampling frequency
was
~
1000
×
that of the Nyquist limit relative to
the
evolution of the slow time complex amplitude in Eq. 1, main text. We perform the
Hilbert transform on the raw time
-
domain signal,
퐻
[
푉
/
(
푡
)
]
, and extract the angle to acquire
the phase (
63
). The magnitudes
푎
/
are the magnitudes of the transformed signals
,
normalized by the magnitudes when oscillators are isolated
,
I
퐻
J
푉
/
(
푡
)
K
I
/
I
퐻
[
푉
/
(
푡
=
0
)
]
I
(with the system uncoupled at
푡
=
0
)
. The phases are further binned and averaged by a
factor of 40 to smooth the time domain data.
Supplementary
Text
Additional exotic stat
es
When
훽
<
1
and
훼
>
0
.
06
the experimental system exhibit
s additional states of interest
-
several kinds of traveling waves and a noise
-
driven chimera. The following sections
describe these briefly.
The system
exhibit
ed
additional traveling wave states ot
her than those outlined in the main
text. Figure
S2A
displays a traveling wave state centered around the
k
=
2,6 fixed points.
It consists of a smooth, sinusoidal variation in the phase difference, much like the main
text’s 2
-
TW
-
I state. In contrast to the
2
-
TW
-
I state, we
observed
a spatial wavevector of
휋
/
4
.
This state
was
found to be stable in the simulation when the parametric interactions
(
described later in the supplementary text
)
wer
e
added.
Figure
S1B
shows a
pulse in the phase difference
which
propagate
d
around the ring.
We
found a similar state
to
occur in the main
text’s
2nd order phase model.
In Fig
.
S2C
, we
plot a complex traveling wave state. This state
did
not appear in the 2nd order phase model,
and
was
only stable at large coupling and l
ow nonlinearity,
훽
>
0
.
7
,
훼
<
0
.
1
.
We also
f
ound
a chimera
-
like state. Various definitions of
“chimera”
on a finite network
have been proposed
(
45
,
46
).
Here we
show
a state that has two types of clusters: one type
that
was
coherent and stable in time and an
other that
was
sensitive to noise and appears to
drift. It differs from the “weak chimeras” mentioned in the main text.
The state’s experimental data appears in Fig.
S3
. Figure
S3A
shows that the state consists
of two locked pairs of phase difference: two
at high frequency and two at low frequency.
These pairs
were
in an antiphase configuration (see
M
ovie
S15
in the Supplementary
Material). Between these pairs
we
re
noisy oscillators, which regularly jump by odd
4
multiples of
휋
,
as shown in Fig
.
S3
B and
D
.
These jumps occur
red
when the neighboring
phases
we
re
antiparallel. Oscillators across the ring
we
re
usually
휋
out of phase at most
times, with occasional
2
휋
slips between noisy pairs.
This state
ha
s
at least three frequencies (and hence the invariant
subspace must be greater
than 2 dimensional) and antiphase correlations between the stable oscillators. It
show
s
a
pattern
{
a,b,c,d,
-
a,
-
b,
-
c,
-
d}
and ha
s
the symmetry of the isotropy subgroup
푍
T
(
푝
=
1
)
.
This
state
was
only stable at high
er values of
coupling
훽
>
0
.
7
. In this state, a small amount of
noise shifts the absolute
phase
of the noisy oscillators. The noisy oscillators synchronize
d
for short periods of time, between periodic slips at random distance (Fig.
S3D
).
Figure
S3E
shows
the long
-
term deviati
on of three oscillators,
훿
휙
/
=
휙
/
−
Ω
/
푡
. We
s
aw
that there
was
no long
-
term drift of phase for oscillator 3
, but
oscillators 2 and 4 perform
ed
random walks. The walks
we
re
seemingly driven by noise, since they d
id
not appear in
deterministic simulations (not shown). There, the
`
noisy’
oscillators slip
ped
regularly in
the sequence
{
휋
,
−
휋
,
휋
,
−
휋
,
...
}
.
Despite the presence of noise, we call this state a chimera
since there
was
a coherent set of oscillators and an incoh
erent set.
Stability of splay states and pattern formation within inhomogeneous synchronized states
Next
,
we examine the stability of the
splay states. While done in Emenheiser, et al. (32)
,
here we derive the result differently to elucidate the formation of the patterns in the main
text’s discussion of inhomogeneous
synchronization.
Absent frequency disorder, starting from the equations for the magnitudes and phase
differences (see the mai
n text for the starting complex amplitude equation),
푑
푎
/
푑푇
=
푎
̇
/
=
1
−
푎
/
2
−
훽
2
J
푎
/
3
4
sin
Δ
/
−
푎
/
6
4
sin
Δ
/
6
4
K
푑
Δ
/
푑푇
=
Δ
̇
/
=
훼
c
푎
/
3
4
T
−
푎
/
T
d
+
훽
2
f
g
푎
/
푎
/
3
4
−
푎
/
3
4
푎
/
h
cos
Δ
/
+
푎
/
3
T
푎
/
3
4
cos
Δ
/
3
4
−
푎
/
6
4
푎
/
cos
Δ
/
6
4
k
(S1)
with
Δ
/
=
휙
/
3
4
−
휙
/
as in the main
text.
Each splay state has a fixed phase difference
Δ
/
=
Δ
=
T
lm
n
(
with
k={0,1,...,N
-
1}
for
N
oscillators) and a magnitude
푎
/
=
1
. Thus, the stability of each state is given by,
훿
푎
̇
/
=
−
푎
/
2
−
훽
2
sin
Δ
c
훿
푎
/
3
4
−
훿
푎
/
6
4
d
−
훽
2
cos
Δ
c
훿
Δ
/
−
훿
Δ
/
6
4
d
(S2)
5
훿
Δ
̇
/
=
2
훼
c
훿
푎
/
3
4
−
훿
푎
/
d
+
훽
2
J
c
3
훿
푎
/
−
3
훿
푎
/
3
4
+
훿
푎
/
3
T
−
훿
푎
/
6
4
d
cos
Δ
−
c
훿
Δ
/
3
4
−
훿
Δ
/
6
4
d
sin
Δ
K
Assuming variations in the state take the form of plane waves
훿
푎
/
=
훿푎
푒
p>/
푒
qr
,
훿
Δ
/
=
훿
Δ
푒
p>
(
/
3
4
)
푒
qr
with spatial
wavevector
푄
and inserting into
Eq
. S2
,
we get,
휎훿푎
=
훿푎
2
J
−
1
−
훽
sin
Δ
c
푒
p>
−
푒
6
p>
d
K
−
훽
2
cos
Δ
훿
Δ
t
푒
p>
T
−
푒
6
p>
T
u
휎훿
Δ
=
2
훼훿푎
t
푒
p>
T
−
푒
6
p>
T
u
+
훽
2
v
cos
Δ
훿푎
t
3
푒
6
p>
T
−
3
푒
p>
T
+
푒
w
p>
T
−
푒
6
w
p>
T
u
−
sin
ΔδΔ
c
푒
p>
−
푒
6
p>
d
y
(S3)
This gives a stability matrix of the form,
z
휎
+
1
2
+
푖훽
sin
Δ
sin
푄
푖훽
cos
Δ
sin
푄
2
−
4
푖훼
sin
푄
2
−
푖훽
cos
Δ
t
sin
3
푄
2
−
3
sin
푄
2
u
휎
+
푖훽
sin
Δ
sin
푄
|
(S4)
whose determinant yields solutions of the form,
휎
=
−
푖훽
sin
Δ
sin
푄
+
1
4
}
1
±
~
1
+
64
훽
cos
Δ
sin
T
푄
2
t
훼
−
훽
cos
Δ
sin
T
푄
2
u
(S
5
)
For the pattern forming states (such as
the inhomogeneous
states from the main text), from
a quench we expect the mode which
was
established to be the one with the maximum
growth rate
Γ
from the real part of
Eq
. S
5
. Maximizing with respect to the wavevector
yields
푄
=
2
sin
6
4
훼
/
2
훽
.
For a single period (wavevector
푄
=
l
), the instability occurs
when
훼
=
훽
cos
Δ
sin
T
휋
/
8
. The same result was found in Emenheiser, et al
. (32).
A
two
-
period
mode with a wavelength of 4 oscillators
(
푄
=
휋
/
2
) becomes unstable for
훼
>
훽
/
2
.
However, the inhomog
eneous states
we
re not found to be stable in the
simulations for such parameter values (see Fig.
3C,
main text). Thus, for the
inhomogeneous states (developing from the uniform in
-
phase state) on the 8
-
ring we expect
(and we f
ound
) that all the patterns ha
ve a single period. Thus, the only inhomogeneous
state we
f
ound
arising from the in
-
phase state
was
the pattern with a wavelength of 8
oscillators, corresponding to the first Fourier mode.
6
Origin of the decoupled traveling wave states
We establish the connection between a nonlinear coupling term and the decoupled traveling
waves in this section.
Nonlinear coupling between the NEMS can be established via parametric terms driving
neighboring oscillators. In Villanueva, et al
. (
15
)
, it w
as shown that piezoelectric NEMS
parametrically respond to harmonics of the feedback: feedback which tunes the frequency
of the NEMS at twice the resonant frequency excites mechanical motion. Harmonics of the
feedback can appear due to compression in the e
lectronic circuits used for amplification of
the small NEMS signal. We f
ound
that a second harmonic at the 5% level (compared to
the first harmonic) generate
d
the decoupled traveling wave states (2
-
TW
-
I)
observed
in the
experimental oscillator ring
.
Suppose that a term appears in the equation of motion for the displacement
푥
/
of oscillator
푗
,
푥
̈
/
=
⋯
+
휁
푥
/
3
4
T
푥
/
(S6)
where we have ignored the contribution from the
푗
−
1
oscillator (whose contribution can
be added in a similar fashion at the end o
f the derivation). We are interested in the dynamics
of Eq.
S6
at the slow time scale
(
64
,1
5
,1
6
,1
7
)
. Assuming the displacement can be written
as a slowly varying complex envelope on top of a fast waveform
푥
/
=
퐴
/
푒
pr
+
퐴
̅
/
푒
6
pr
,
2
푖휔
퐴
̇
/
푒
p
r
−
퐴
̅
̇
/
푒
6
pr
=
⋯
+
휁
J
I
퐴
/
3
4
T
I
+
퐴
/
3
4
T
푒
T
pr
+
퐴
̅
/
3
4
T
푒
6
T
pr
K
×
J
퐴
/
푒
pr
+
퐴
̅
/
푒
6
pr
K
(S
7
)
where
퐴
/
is the conjugate of the complex amplitude
퐴
/
, and
휔
is the frequency of oscillation.
Gathering the secular terms at
푒
pr
,
writing out
퐴
/
=
푎
/
푒
p
, and dividing through by
푒
p
gives,
2
휔
c
푎
̇
/
+
푖
푎
/
휙
̇
/
d
=
⋯
−
휁
푎
/
3
4
T
푎
/
J
푖
c
1
−
cos
c
2
J
휙
/
3
4
−
휙
/
K
d
d
+
sin
c
2
J
휙
/
3
4
−
휙
/
K
d
K
(S
8
)
which when inserted in the full magnitude
-
phase equation, normalizing the time by
휔
and
putting in the
푗
−
1
terms,
7
푑
푎
/
푑푇
=
1
−
푎
/
2
−
훽
2
J
푎
/
3
4
sin
c
휙
/
3
4
−
휙
/
d
+
푎
/
6
4
sin
c
휙
/
6
4
−
휙
/
d
K
−
휁
푎
/
2
J
푎
/
3
4
T
sin
c
2
J
휙
/
3
4
−
휙
/
K
d
+
푎
/
6
4
T
sin
c
2
J
휙
/
6
4
−
휙
/
K
d
K
푑
휙
/
푑푇
=
훼
푎
/
T
−
훽
t
1
+
휁
2
J
푎
/
3
4
T
+
푎
/
6
4
T
K
u
+
훽
2
푎
/
J
푎
/
3
4
cos
c
휙
/
3
4
−
휙
/
d
+
푎
/
6
4
cos
c
휙
/
6
4
−
휙
/
d
K
−
휁
2
J
푎
/
3
4
T
cos
c
2
J
휙
/
3
4
−
휙
/
K
d
+
푎
/
6
4
T
cos
c
2
J
휙
/
6
4
−
휙
/
K
d
K
(S9)
Note that we have generated explicit higher harmonic phase coupling
through these
additional parametric coupling terms in the original equation of motion for the oscillators
.
For the appropriate choice of t
he sign of
휁
, we found
th
at
these h
igher order couplings
dr
ove
the trajectories away from the fixed points associated with
Δ
=
휋
2
⁄
,
3
휋
2
⁄
. This
create
d
an instability since the trajectories would decay into these fixed points in the
absence of this additional nonlinear parametric coupling.
We performed a numerical simulation of
Eq
. S9
with
휁
=
0
.
05
훽
, and obtain
ed
behavior
similar to that demonstrated by the
experimental
2
-
TW
-
I state in the main text. The results
are plotted in Fig
.
S4
.
These two plots can be qualitatively compared to the exp
erimental
data from Fig. 4
B
in the main text.
Comparing the results quantitatively, t
he frequency
of
the oscillation in the phase differences found in simulation was within 5
% of the
experiment.
Simulations of the extended phase model
Here we show that the states found in experiment (except for the decoupled traveling wave
and 0
-
TW
-
I states)
we
re
also found in simulations of the extended phase model. We have
already shown a simulation of the WC
-
II state in the main text, so will not rep
roduce it
here. Also, the splay states associated with
푘
=
{
3
,
4
,
5
}
are stable within the Kuramoto
model, so we do not show these states either. All the plots in this section
we
re
composed
with
numerical data from simulations of Eq. (3) in the main text.
W
e start with the
푘
=
{
0
,
1
}
splay states (the
푘
=
7
state has the same stability as the
푘
=
1
state) and the inhomogeneous synchronized state. These are shown as phase difference
plots in Fig.
S
5A
-
C
. Incidentally, the only coupling beyond the Kuramoto
-
Sak
aguchi terms
in the phase model (with repulsive coupling) required for stabilizing the
푘
=
1
w
ere
the
triadic interactions. Likewise, only the biharmonic interactions
we
re
required for the in
-
phase state. In Fig
.
S5D
, we show the pattern of phase
differences for the data from
C
.
Next, we show simulations of the decoupled precessing state (2
-
precess) and the first weak
chimera (WC
-
I). We only show the frequency space
-
time plots in Fig
.
S
6
, since the
clustering can be immediately compared to the exp
erimental data in the main text. We do
not present the scale of the frequencies in the plots
; these frequencies
are determined by
8
the
system’s
disorder (as in the 2
-
precess state) or
were
quantitatively incorrect (as in the
WC
-
I state).
Simulations of th
e extended phase model
T
he supplementary movies show experimental data for
states mentioned within the text
and supplementary information. There are two representations in video form: "spin type"
and "firefly type". In the spin type movies, oscillations
are
represented by arrows with
angle corresponding to phase, and length corresponding to magnitude. In the firefly type,
the phase of the oscillation
is
represented by a color on a cyclic color scale, with no
representation of magnitude.
9
Fig. S1
:
Experimental setup
.
A
)
The ring network is constructed from eight individual
piezoelectric plate resonators (SEM micrograph, top
, with false color to contrast functional
layers
) with dimensions
푎푟푒푎
×
푡
ℎ
푖푐푘푛푒푠푠
=
1600
휇
푚
T
×
.
145
휇푚
.
Each resonato
r is
bonded (bottom) and housed in the white vacuum chamber on the oscillator printed circuit
board (PCB).
B
)
On the 'oscillator' board, the feedback to generate self
-
oscillation is
applied. Controls for natural frequency and nonlinearity are shown by blue
and white
boxes.
C
)
Eight oscillator boards are installed on the network board to form the ring
network.
D
)
Full diagram of the system. The PC sends digital commands to the controller
over ethernet (dark blue box) who translates those commands to the di
gital bus (black
cable) which connects to the network board. The data is read out via the yellow SMA cables
and sent to an 8
-
channel oscilloscope. The controls for one edge
are
shown by the yellow
box.
10
Fig. S2
:
Other tr
aveling wave states
.
A)
2
-
TW
-
II. Additional traveling wave state around
the
푘
=
2
,
6
fixed points. It has a spatial wavevector of
휋
/
4
.
B)
Phase pulse
.
C)
Comp
lex
traveling wave state 0
-
TW
-
I.
0
.
865
0
.
875
N
o
d
e
Num
b
er
0
π
2
π
1
2
3
4
5
6
7
8
tim
e
(s)
1
2
3
4
5
6
7
8
N
o
d
e
Num
b
er
1
2
3
4
5
6
7
8
N
o
d
e
Num
b
er
Ι
(
φ
i
;φ
j
)
0.8
2.4
4.0
-600
-200
Δ
f (Hz)
∆
j
2-TW-II
N
o
d
e
Num
b
er
N
o
d
e
Num
b
er
0
π
0
.
795
0
.
800
-400
200
-
π
1
2
3
4
5
6
7
8
N
o
d
e
Num
b
er
1
2
3
4
5
6
7
8
tim
e
(s)
1
2
3
4
5
6
7
8
0.8
2.4
4.0
Ι
(
φ
i
;φ
j
)
Δ
f (Hz)
B
∆
j
A
0
.
625
0
.
630
tim
e
(s)
−
1200
400
0
π
2
π
Δ
f (Hz)
0.8
2.4
4.0
Ι
(
φ
i
;φ
j
)
∆
j
N
o
d
e
Num
b
er
1
2
3
4
5
6
7
8
N
o
d
e
Num
b
er
1
2
3
4
5
6
7
8
N
o
d
e
Num
b
er
1
2
3
4
5
6
7
8
C
Pulse
0-TW-I
11
Fig. S
3
:
Noise driven chimera.
A
) Phase differences, frequencies, and pairwise mut
ual
information of the noisy chimera state. Note that the frequencies of 4 of the oscillators
display aperiodic jumps in frequency. The mutual information shows a partial
풁
ퟐ
symmetry.
B
) We plot the phases of three of the oscillators in the rotating frame where the
phase of oscillator 4 does not precess between slips. This shows that without slips, the
oscillator is at the mean frequency of its neighbors.
C
) Phase differences across the
ring
show slips across the ring occur in steps of
2
휋
.
D
) The histogram of oscillator events shows
that all the phase slips are odd multiples of
흅
. The nonzero mean of the histogram shows
that additional phase precession
occurred
due to these odd
-
휋
slip
s.
E
) The deviation of the
phases from their mean precession (
훿
휙
/
=
휙
/
−
Ω
/
푡
) shows that the noisy oscillators
were
highly unstable compared the locked oscillators and appear driven by noise
12
Fig. S
4
:
Simulation of parametrically coupled ring
.
The
simulation shows that two of
the nearest neighbor phase differences (j=2 and j=4) have opposite phase and oscillate
around
ퟑ흅
ퟐ
(top)
and t
he space
-
time plot of frequencies shows the exact same wavevector
as the experimental data
(bot)
.
550
560
570
580
590
600
n
o
rmalize
d
time
0
3
π
/
2
2
π
1
2
3
4
5
6
7
8
N
o
d
e
Num
b
er
−
0
.
25
0
.
00
Δ
f
∆
j
π
/
2
π