A large array of elastically coupled micro cantilevers of variable length is studied experimentally and numerically. Full-scale finite element (FE) modal analysis is implemented to determine the spectral behavior of the array and to extract a global coupling matrix. A compact reduced-order (RO) model is used for numerical investigation of the array's dynamic response. Our model results show that at a given excitation frequency within a propagation band, only a finite number of beams respond. Spectral characteristics of individual cantilevers, inertially excited by an external piezoelectric actuator, were measured in vacuum using laser interferometry. The theoretical and experimental results collectively show that the resonant peaks corresponding to individual beams are clearly separated when operating in vacuum at the third harmonic. Distinct resonant peak separation, coupled with the spatially confined modal response, make higher harmonic operation of tailored, variable-length cantilever arrays well suited for a variety of resonant-based sensing applications.

## Introduction

Dynamics of large arrays of micro- and nano-electromechanical-coupled resonators have received significant research attention over the last two decades [1]. The first works on micromachined arrays were motivated by the development of miniature electromechanical filters [2]. Coupled resonators with nominally identical or slightly detuned resonant frequencies yield a wider bandwidth when compared to a single resonator. This basic feature continues to stimulate research advancements in the area of micromechanical filter design by exploring diverse architectures and operational principles [311]. Due to their intrinsic filtering feature and ability to increase bandwidth, architectures based on weakly coupled resonators have found applications in inertial sensing [12] and frequency-signature speech processing [13]. Micro- and nano-electromechanical arrays also have a great potential to serve as ultrasensitive detectors for chemical or biological analytes [14]. In contrast to single micromechanical sensing structures where frequency changes are monitored [1519], sensing within arrays of weakly coupled resonators is based on modal shape changes [2025]. Eigenmode changes in an array of mechanically coupled, nearly identical microcantilevers, can be two to three orders of magnitude greater than relative changes in resonance frequencies when a mass is added [26]. Mode localization has also been used in accelerometers [27] and light processing applications [28]. Previous reports highlight the influence of various system parameters such as the mass-ratio [24], measurement noise [23], nonideal clamping [29], coupling stiffness [30], anisotropy of the Young's modulus [31], and global and dissipative coupling [32] on the array dynamics. In these arrayed structures, vibrations are possible only within the allowed propagation band [33], defined as the interval of frequencies between the lower (fL) and the upper (fU) cutoff values. While most of the previous studies were focused on arrays of identical or almost identical beams [23,26,3336], relatively few works considered the case of beams with differing resonant frequencies [13,21,37,38].

In this work, we investigate, both numerically and experimentally, the dynamics of a micromechanical cantilever array, elastically coupled through a flexible overhang. The array is composed of 100 beams with linearly varying length. The fL and fU are defined by the frequencies of the longest and the shortest cantilevers, respectively. Our results show that at excitation frequencies within the propagation band, the vibrations are spatially localized within a region spanning only a few of the neighboring beams. We also demonstrate excitation of the cantilevers at their higher harmonics and map the corresponding propagation bands of the array experimentally and numerically. The spatial localization of vibrations, allowing for excitations of distinct array elements, enables new applications of arrayed cantilever systems in a variety of sensing applications. Direct probing of distinctly separated resonant peaks offers fundamentally new functionalities by enabling large resonator arrays to be utilized for mass and inertial sensing. Furthermore, engineered large resonator arrays mimicking a mechanical fourier transform system could lead to the development of complex networks that could stimulate novel classes of frequency sensing systems.

## Device Architecture

The array shown in Fig. 1 contains N = 100 prismatic cantilevers attached to a compliant overhang and designed to deflect in the out-of-plane (z) direction. The cantilevers and the overhang have the same thickness and are fabricated from the device layer of a silicon-on-insulator substrate [36,39]. The device is attached to the silicon substrate by an underlying ≈ 3 μm thick buried silicon dioxide layer. To allow large unobscured vibrations and to prevent stiction, an opening was created within the handle wafer under the beams. We made cantilever arrays with b ≈ 20 μm, h ≈ 5 μm, Lo ≈ 100 μm, pitch B ≈ 50 μm, L1 = Lmax ≈ 500 μm, and LN = Lmin ≈ 350 μm.

## model

### Finite Element Model.

A full-scale three-dimensional (3D) finite element (FE) analysis of the array was carried out using a commercially available package. Within our analysis, two-dimensional, eight-node, rectangular shell elements with quadratic interpolating functions in each direction were used to mesh the overhang region. A two-node, 12 degrees-of-freedom, three-dimensional beam element with rectangular cross section and consistent mass representation was used to mesh the cantilevers. The overhang plate was clamped at three out of its four edges. Several mesh refinements were performed to assure that the average numerical natural frequency error is less than 1%.

The results of the numerical linear modal analysis are shown in Figs. 24. In general, an infinite number of propagation bands can be obtained, each corresponding to higher harmonics of the cantilevers. Since in this work, we measured the first three propagation bands, the FE analysis is carried out for the first 300 (3N) natural frequencies and the corresponding natural modes of the array. Figure 2 shows three-dimensional snapshots of several natural modes corresponding to different frequencies fi, where i = 1, 2 … 300 is the array's mode number. Figure 3 shows normalized endpoint beam deflections at several natural modes. At each of the natural frequencies only a limited number of beams manifest observable modal amplitudes. The position of this spatial band along the array is governed by the natural frequency. For instance, as shown in Fig. 3, at lower and higher frequencies, the spatial bands comprise longer and shorter beams, respectively. The width of the spatial band, namely the number of the beams vibrating at a specific natural frequency, is approximately the same for all the propagation bands. Figure 4 shows the natural frequency dependence on the mode number. The three frequency ranges correspond to the three, nonoverlapping propagation bands. We anticipate that the dynamics become more complex for overlapping bands. Since N = 100, the first, second, and third bands each contain 100 frequencies, and are attributed to cantilevers vibrating at their first, second, and third harmonics, respectively. The lowest cutoff frequency obtained using the FE model is $fL(1)=f1=24.851$ kHz, and the upper cutoff of the first propagation band is $fU(1)=f100=55.172$ kHz. Hereafter $()(j), j=1,2,3$ denotes the propagation band number and the cantilever harmonic number. For L1 = 500 μm and h = 5 μm, the first three harmonics of the ideally clamped beam, calculated using the Euler–Bernoulli theory, are $f(1)=27.691$ kHz, $f(2)=173.550$ kHz, and $f(3)=485.994$ kHz. The first, second, and third propagation bands are associated with the cantilever vibrations at their first, second, and third harmonics, respectively. The frequency cutoff values for the second band are $fL(2)=f101=151.198$ kHz and $fU(2)=f200=345.638$ kHz, and for the third band are $fL(3)=f201=409.230$ kHz and $fU(3)=f300=845.169$ kHz. Figure 4 shows that the propagation band frequency range is significantly larger at higher harmonics. The FE model results show the frequency range of the first, second, and third band as $f100−f1=30.321$ kHz, $f200−f101=194.440$ kHz, and $f300−f201=435.939$ kHz, respectively. Consequently, the shift between the first two frequencies within the third propagation band $f202−f201=11.575$ kHz is greater than $f102−f101=4.166$ kHz, which is, in turn, greater than $f2−f1=635$ Hz.

To understand the reason of this propagation band stretching, consider the two adjacent ideally clamped cantilevers n and n + 1 of lengths Ln and $Ln+1=Ln−ΔL$, respectively. Here, $ΔL=1.515 μ$m is the difference in length between any two adjacent cantilevers. The shift between the natural frequencies of these two cantilevers vibrating at the harmonic j = 1, 2, 3 is [40]
$fn+1(j)−fn(j)=(λ(j))22πEIρA(Ln−ΔL)4−(λ(j))22πEIρALn4≈(λ(j))2πEIρALn4(ΔLLn) (ΔL≪Ln)$
(1)

where A = bh and I = bh3/12 are the corresponding area and the second moment of area of the rectangular beam cross section, and $λ(j)$ is the eigenvalue corresponding to the jth harmonic of the cantilever. E = 169 GPa and ρ = 2300 kg/m3 are taken as the Young's modulus and the density of the silicon cantilever, respectively. Since for an ideally clamped cantilever $λ(1)=1.875, λ(2)=4.694$, and $λ(3)=7.855$, the frequency shift is higher for the higher harmonics. Since Ln decreases with increasing n, the frequency shift increases with the mode number within each propagation band.

### Reduced-Order Model.

Analysis of the array's dynamics using full-scale FE models is time consuming. For this reason, simplified reduced-order (RO) models of arrays are often constructed [33,34,36,4143]. In these models, the array is represented as a mass-spring lattice chain. Corresponding mass, onsite (OS), and intersite (IS) stiffness parameters are obtained using either simplified beam models [33,34,43] or Galerkin decomposition [36,41,42].

Here, we construct the RO model of the array using modal representation for the deflection of each cantilever. Vibrational dynamics for each cantilever are described in the framework of the Euler–Bernoulli theory. Since our work is focused on the linear modal analysis of the array, we neglect geometric and inertial nonlinearities that are associated with large deflections and rotations of the beam [33]. Under these assumptions the dynamics of the nth beam within the array are described by the following nondimensional partial differential equation:
$∂2wn∂t2+c∂wn∂t+ln4∂4wn∂y4=−∂2zB∂t2$
(2)

Here, wn is the nondimensional deflection of the nth beam, y and t are the coordinate along the beam and time, respectively, c is the coefficient of the linear viscous damping (related mainly to thermoelastic and clamping losses), and ln = L1/Ln is the length ratio parameter. The right-hand side in Eq. (2) represents a uniformly distributed force, which appears due to a kinematic inertial excitation by a piezo-electric transducer resulting from the substrate motion with acceleration $∂2zB/∂t2$. Nondimensional quantities used in Eq. (2) are presented in Table 1. Since the coordinate $ŷn$ along each of the beams is normalized by the length Ln of the corresponding beam, we have $y∈[0,1]$.

Since the bands are not overlapping and the corresponding frequencies are separated (Fig. 4) we assume that the motion of each cantilever can be represented by the expression (Einstein's summation convention is not used)
$wn(j)(y,t)≈qn(j)(t) φ(j)(y)$
(3)
where $qn(j)(t)$ denotes the generalized coordinate corresponding to the jth harmonic of the nth beam. The functions $φ(j)(y)$ are the jth linear undamped eigenmodes of an ideally clamped cantilever given by the expression [40]
$φ(j)=C(j)[sin(λ(j)y)−sinh(λ(j)y)−(sin(λ(j))+sinh(λ(j)))(cos(λ(j)y)−cosh(λ(j)y)) cos(λ(j))+cosh(λ(j))]$
(4)

Here, the constants $C(j)$ are chosen such that $φ(j)(1)=1$ and consequently $qn(j)(t)=wn(j)(1,t)$. The mode shapes $φ(j)(y)$ are obtained as a solution of the eigenvalue problem associated with Eq. (2) with $c=0, z¨B=0$ and subjected to the free-end boundary conditions $∂2φ(j)/∂2y=0, ∂3φ(j)/∂3y=0$ at $y=1$ and ideally clamped boundary conditions $φ(j)=0, ∂φ(j)/∂y=0$ at $y=0$. Since at y = 0 the beam is attached to a flexible overhang, the clamping conditions are nonideal. The resulting dynamics give rise to a decrease of the beam's natural frequency and in the mechanical coupling between the cantilevers [33]. In our work, we first use linear undamped eigenmodes of an ideally clamped beam for the development of the OS terms and then account for mechanical coupling by introducing the IS stiffness terms directly into the RO model [36].

Substitution of Eq. (3) into Eq. (2), followed by the common Galerkin procedure, yields a system of N linear ordinary differential equations. Because the arrayed cantilever devices are coupled through the flexible overhang, further modifications to the model are required to account for this IS mechanical interaction. In general, since mechanical coupling is not local, each beam interacts with beams beyond its nearest neighbors, and the coupling matrix is fully populated. By adding elastic coupling, while neglecting the overhang inertia, we obtain
$mq¨n(j)+cmq˙n(j)+ko(j)ln4qn(j)−∑sk̃ns(j)qs(j)=−z¨Ba(j)$
(5)
where
$m=∫01[φ(j)]2 dy ko(j)=∫01[(φ(j))″]2 dy a(j)=∫01φ(j) dy$
(6)

In Eq. (5), m is related to the mass of the beam, the coefficients $ko(j)$ are associated with the linear OS bending beam stiffness, and $k̃ns(j)$ are the IS stiffness coefficients. For the adopted base functions in Eq. (4), one obtains $m=0.25, ko(1)=3.091, ko(2)=121.3809,ko(3)=951.637, a(1)=0.391, a(2)=0.217, a(3)=0.127$. $ko(j)/m=ω(j)=(λ(j))2$ is the nondimensional $jth$ harmonic frequency of the ideally clamped cantilever. In addition, $( )′=∂/∂y$ and $( )˙=∂/∂t$ denote derivatives with respect to the nondimensional coordinate along the beam and nondimensional time, respectively.

Subdividing Eq. (5) by m and further rescaling time ($τ=tω(1)$), yields
$q¨n(j)+ln2Qnq˙n(j)+ln4qn(j)−∑skns(j)qs(j)=−γ(j)z¨B$
(7)

where $Qn=ω(1)ln2/c$ is the OS damping parameter, $γ(j)=a(j)/m$ is the substrate acceleration parameter, and the over-dot is re-defined as $( )˙=∂/∂τ$. In addition, the nondimensional IS stiffness coefficients are redefined as $kns(j)=k̃ns(j)/k0(1)$.

Equation (7) can be conveniently written in the matrix form
$Mq¨(j)+Cq˙(j)+K(j)q(j)=F(j)$
(8)

where $q={qn(j)}T$ is the displacements vector, $M=I$ is the unit mass matrix, $C=[ln2/Qnδns]$ is the diagonal damping matrix, $F(j)={−γ(j)z¨B}T$ is the force vector, δns is the Kronecker's delta, and ${ }T$ denotes matrix transpose. The fully populated stiffness matrix $K(j)=[ln4 δns−kns(j)]$ contains both OS and IS components.

A simpler and widely used [31,44] counterpart of Eq. (7) can be derived based on the assumption of local, nearest neighbor, mechanical interaction
$q¨n(j)+lnQnq˙n(j)+ln2qn(j)−η(j)(qn+1(j)−2qn(j)+qn−1(j))=−γ(j)q¨B$
(9)

where $η(j)=knn(j)/2 (n=N/2)$ is the nondimensional local IS coupling parameter corresponding to the jth propagation band. The local model, Eq. (9), allows a simple physical representation of the array using a mass-spring system [33].

### Extraction of the Stiffness Coefficients From the Finite Element Model.

In order to evaluate the stiffness coefficients $kns(j)$, the results of a full-scale FE modal analysis were used [35,36]. The analysis provided the values of the 3N frequencies and the corresponding 3N eigenvectors. The first set of N values corresponds to the cantilevers vibrating at their first harmonic within the first propagation band. The second and third sets are associated with the second and the third propagation bands, respectively.

Since the propagation bands do not overlap, the same procedure described below is performed for each band separately. Each eigenvector is obtained as a set of nodal displacements in the FE solution. Then, for each eigenvector, subsets $ψ̃(r)(j),$ are built such that individual vectors $ψ̃(r)(j)$ correspond to out-of-plane modal displacements at the free-end of the cantilever. Here, $( )(r)(j)$ denotes the eigenvector number within the jth propagation band, where $r=1…N$. Since these eigenvectors are obtained numerically, using an approximate FE model, Gram–Schmidt orthogonalization was carried out prior to the formation of the modal matrix. Consequently, we define the modal matrix as [40]
$Ψ̃(j)=[ψ̃(1)(j) ψ̃(2)(j) … ψ̃(N)(j)]$
(10)
Normalizing the eigenvectors using the orthogonality with respect to the unit mass matrix (Eq. (7)) yields
$Ψ(j)=Ψ̃(j)Ψ̃(j)TΨ̃(j)$
(11)
The normalized modal matrix is orthogonal with respect to the stiffness matrix $K(j)$
$Ψ(j)TK(j)Ψ(j)=Λ(j)$
(12)
In Eq. (12), $Λ(j)=[λ(r)(j)δrs]$ is the diagonal matrix where $λ(r)(j)=λ(r)(j)FE/λ(1)(j)FE$ are the normalized eigenvalues extracted from the full-scale FE analysis. As a result, the lowest cutoff normalized eigenvalue of each of the propagation bands is equal to unity $λ(1)(j)=1$. Equation (12) allows us to express the stiffness matrix in terms of the eigenvalues and eigenvectors obtained using the FE model in the following form:
$K(j)=(Ψ(j)T)−1Λ(j) Ψ(j)−1$
(13)
To evaluate the role of the nonlocal coupling, it is instructive to consider the structure of the stiffness matrix. Here, we present a few elements around the main diagonal by including the first, last and middle ($N/2=50$) rows of $K(1)$
$K(1)=(1.2152−0.0006−0.0002−0.00010…1.1216−0.0309−0.0136−0.0061−0.0026−0.0011…1.1124−0.0430−0.0191−0.0085−0.0036−0.0015…1.1192−0.0466−0.0206−0.0091−0.0038−0.0016…⋮…1.9916−0.0984−0.0430−0.0187−0.0077−0.0031…⋮…3.9831−0.2188−0.0666−0.0002SYM…4.0820−0.1584−0.0005…4.2783−0.0026…4.9289)$
(14)

where the dominant diagonal elements contain also the onsite stiffness terms $ln2$.

### Results From Reduction-Order Model.

First, the eigenfrequencies of the array were obtained using the undamped homogeneous counterpart of Eq. (8). The coefficients of the stiffness matrix $K(j)$ were then calculated using Eq. (13). Our results show that the nondimensional frequencies obtained using the RO model coincide with the values provided by the full-scale FE analysis. To quantify the influence of the nonlocal coupling in the stiffness matrix, the eigenfrequencies were also calculated using the local model, Eq. (9) (with the coupling coefficients η = 0.0447, 0.0712, and 0.0965 for the first, second, and third propagation bands, respectively), and the error $|fiLOC−fiFE|/fiFE, i=1…3N$ was evaluated for each of the frequencies. The values of η were obtained using the FE model results, by equilibrating the diagonal terms of the fully populated, Eq. (13), and the local, Eq. (9), matrices $Knn(j)=ln2+2η(j), n=N/2$. Our results show that the local model, while providing a good accuracy (with an average error of 0.47%, 0.65%, and 0.5% for the first, second, and third bands), may lead to an error of up to 9% in the frequencies close to the upper cutoff values. In this work, we use Eq. (7) with the fully populated matrix and Eq. (13) for the numerical analysis of the array behavior.

To illustrate the response of the array to an inertial excitation, Eq. (8) was solved numerically using the Runge–Kutta solver and time-series were obtained for each cantilever. The substrate acceleration parameter and the quality factor used in our calculations were γ = 0.01 and Q1 = 1000, respectively. The array was subject to harmonic driving such that $z¨B=sin(Ωτ)$, where Ω is the excitation frequency. During the numerical frequency sweep, Ω was incrementally varied between the lower and the upper frequency cutoff values within the three propagation bands. At each excitation frequency, the vibrational amplitudes were obtained using a time-series output from the differential equation solver.

Figures 5(a) and 5(b) show the spectral responses for the L25 = 463.64 μm beam at the first and second harmonics, respectively. Insets in Fig. 5 show the frequency interval of ≈330 Hz and ≈2100 Hz for the corresponding first and second propagation bands. Figure 6 shows the modal patterns of the array vibrating within the first and the second propagation bands. Our RO model results (Fig. 6) show good agreement with the full-scale FE model prediction (Fig. 4).

## Experiment

### Setup.

The experimental setup is shown schematically in Fig. 7. Chips containing several arrays were indium bonded to a piezo-electric actuator. The piezo and chip assembly was mounted onto a holder and placed into a high vacuum chamber equipped with an optical quality viewport. The chamber was evacuated to ≈ 10−4 Pa (≈ 10−6 mbar), where air damping is negligible. Electrical feed-throughs were used to apply signals to the piezo-electric actuator [45].

An optical interferometric system was used to monitor the motion of the devices. We used a manual XY translation stage to position a laser spot onto a selected cantilever device. Light from a ≈15 mW He–Ne laser was directed through the microscope objective onto the resonator. The reflected light collected by the objective was measured using an AC coupled photodetector (PD). The signal from the PD was fed to the input of a spectrum analyzer. The amplified spectrum analyzer radio frequency output signal was fed to the piezo-electric actuator to inertially excite device motion. In a typical experiment, a 10× objective was used to focus the laser to a spot of approximately 10 μm diameter on a specific beam to be measured, as shown in Fig. 8. The piezo drive frequency was swept over a span that captures the array resonances. To assure fully developed steady-periodic response and to eliminate the influence of the transient effects, the sweep time was chosen to be at least ≈20 s and up to ≈180 s which is significantly longer than the settling time of the cantilevers vibrations. The resulting spectral response was measured from the modulated PD output using the spectrum analyzer.

## Experimental Results

Figure 9 shows the measured spectral response of several cantilevers at frequencies within the first ((a)–(e)) and second ((f)–(j)) propagation bands. Our results show that the separation between neighboring spectral peaks is larger in the second propagation band. The frequency interval between neighboring peaks within the first propagation band for L1 and L25 are ≈620 Hz and ≈500 Hz, respectively (Figs. 9(a) and 9(c)). In contrast, the frequency interval within the second propagation band is much higher, where values for L1 and L25 are ≈3450 Hz and ≈3200 Hz, respectively (Figs. 9(f) and 9(h)). Figure 9 also shows that the spectral band shifts toward higher frequency values with increasing n. The vibrations of the longest cantilever n = 1 are measured at frequencies between ≈26.13 kHz and ≈29.34 kHz (Fig. 9(a)), whereas the lowest and the highest resonances of the shorter (n = 25) beam are ≈ 27.29 kHz and ≈ 34.33 kHz, respectively (Fig. 9(c)). Since increasing n implies shorter beams, the propagation band will shift toward higher frequencies as n increases.

Our experimental results show that when the drive frequency is swept up, the spatial band propagates along the array from the longest toward the shortest beam (Figs. 9(a)9(e) and Figs. 9(f)9(j)). The width of this spatial band, namely the number of beams vibrating at measurable amplitudes at a specific excitation frequency in the range between the upper and the lower frequency cutoff, is smaller than the length of the array. Consequently, the number of peaks in the measured spectral response of a specific cantilever is smaller than the number of peaks corresponding to the entire array. For example, Fig. 9(c) shows 33 spectral peaks for L25, whereas the total number natural frequency in the first propagation band is 100. Moreover, due to the oscillating character of the natural modes of the array (Fig. 3), the peak amplitudes corresponding to the different modes of the array are not equal, Fig. 9. Due to higher stiffness, the shorter cantilevers, when compared to longer beams, have lower vibrational amplitudes.

Figure 10 shows the measured third propagation band for L25. The vibrations of the cantilever are measured within the interval of frequencies spanning ≈166 kHz. The frequency interval between the two lowest peaks is ≈7700 Hz. Figures 10(b) and 10(c) show the influence of the driving voltage on the measured spectra of the cantilever. Our measured linear drive voltage–amplitude dependence (inset of Fig. 10(c)) show the expected Lorentzian spectral shapes, as shown in Fig. 10(c). As expected for a linear system, the center frequency value of the third harmonic was independent, to within a few Hz, of the drive voltage. Within the linear drive regime, we observed frequency shifts of 6.1±2.6 Hz (mean ± error from the Lorentzian fit). The wide natural frequency separation coupled with high quality factors (varying between Q ≈ 7945 and Q ≈ 17953 on Figs. 10(b) and 10(c)) allow easy identification of individual peaks within the measured spectrum of the arrays driven at higher harmonics.

Figure 11 shows the theoretical and experimental frequencies of the array as a function of the mode number. We measured 57, 60, and 43 resonances in the first, second, and third bands, respectively. The higher mode numbers that are associated with the shorter beams require higher excitation amplitudes. At increased drive amplitudes we observed nonlinear device behavior that leads to structural damage of the array elements.

Both the theoretical (Fig. 5) and the experimental (Figs. 9 and 11) results demonstrate stretching of the propagation band at higher harmonics. Figure 11 also shows that the measured frequencies are higher than the theoretically predicted values. We attribute these quantitative discrepancies to the uncertainty in the device geometry and the overhang length Lo. The latter is governed by the side wall verticality of the backside etch and the alignment between the front-to-back lithographic levels. Higher experimental frequencies imply that the overhang length was smaller than the nominal, as-designed, value.

## Conclusions

In this work we explored, both numerically and experimentally, the collective behavior of 100 cantilevers elastically coupled through a flexible overhang. We carried out the FE analysis for the first 300 natural frequencies and natural modes of the array. The numerical results show three distinct types of modes. The first 100 natural modes of the array correspond to cantilevers vibrating at their fundamental harmonics. The modes between 101 and 200 are associated with the second harmonics and modes between 201 and 300 with the third harmonics of the beams. The corresponding natural frequencies form three propagation bands, each bounded from below and from above by the lower and the upper cutoff frequencies, respectively. The natural modes of the array have localized characteristics whereby limited number of beams oscillates at each of the natural frequencies. Our FE results also show stretching of the propagation bands for cantilevers vibrating at higher harmonics.

We also presented a compact reduced-order model of the array that was built using the Galerkin decomposition. The three terms preserved in the Galerkin series provided a description of beam vibrations at the first three harmonics. This was accomplished by using the stiffness matrix that was evaluated from the results of the FE analysis. The RO model also provided the response of the array to harmonic excitation spanning the lower and upper cutoff frequencies for each propagation band. We developed the time history dynamics of each beam by numerically solving, using a Runge–Kutta solver, a system of ordinary differential equations for a mass-spring system.

Our experimental results of the array dynamics showed distinct resonant spectra and were in good agreement with the FE and RO model predictions for the three measured propagation bands. Furthermore, the mode localization feature coupled with the ability to control the location of the spatial propagation band along the array allows for addressing selected cantilevers by changing the excitation frequency. We anticipate that the array dynamics at higher frequencies, corresponding to higher harmonics of the cantilevers, may reduce the flicker and low-frequency environmental noise. The distinct spectral separation combined with the Lorentzian character of the response present new possibilities in frequency-based sensing with coupled micromechanical arrays.

## Acknowledgment

This work was performed in part at the Center for Nanoscale Science and Technology (CNST) at the National Institute of Standards and Technology (NIST) and the Cornell NanoScale Facility, a member of the National Nanotechnology Coordinated Infrastructure (NNCI), which was supported by the National Science Foundation. We thank Richard H. Rand for helpful discussions and Serhan Ardanuc for help with device measurements. This work was partially supported by the National Science Foundation. Sandia National Laboratories is a multi-mission laboratory managed and operated by Sandia Corporation, a wholly owned subsidiary of Lockheed Martin Corporation, for the U.S. Department of Energy's National Nuclear Security Administration. We acknowledge the support of the Mary Shepard B. Upson Visiting Professorship in Engineering, The Sibley School of Mechanical and Aerospace Engineering, Cornell University. Slava Krylov acknowledges support from the Henry and Dinah Krongold Chair of Microelectronics. Christopher B. Wallin acknowledges support under the Cooperative Research Agreement between the University of Maryland and the National Institute of Standards and Technology Center for Nanoscale Science and Technology through the University of Maryland.

## Funding Data

• Center for Nanoscale Science and Technology (Grant No. 70 NANB14H209).

• National Science Foundation (Grant Nos. 1634664 and 1542081).

• U.S. Department of Energy (Contract No. DE-AC04-94AL85000).

## References

References
1.
,
J.
,
Shaw
,
S.
, and
Turner
,
K.
,
2010
, “
Nonlinear Dynamics and Its Applications in Micro- and Nanoresonators
,”
ASME J. Dyn. Syst., Meas. Control
,
132
(
3
), p.
034001
.
2.
Lin
,
L.
,
Nguyen
,
C.-C.
,
Howe
,
R. T.
, and
Pisano
,
A. P.
,
1992
, “
Microelectromechanical Filters for Signal Processing
,” IEEE
Micro Electro Mechanical Systems
(
MEMS'92
), An Investigation of Micro Structures, Sensors, Actuators, Machines and Robots, Travemunde, Germany, Feb. 4–7, pp.
226
231
.
3.
Ho
,
G. K.
,
Abdolvand
,
R.
, and
Ayazi
,
F.
,
2004
, “
Through-Support-Coupled Micromechanical Filter Array
,”
17th IEEE International Conference on Micro Electro Mechanical Systems
(
MEMS
), Maastricht, The Netherlands, Jan. 25–29, pp.
769
772
.
4.
Elka
,
A.
, and
Bucher
,
I.
,
2008
, “
On the Synthesis of Micro-Electromechanical Filters Using Structural Dynamics
,”
J. Micromech. Microeng.
,
18
(
12
), p.
125018
.
5.
Chivukula
,
V. B.
, and
,
J. F.
,
2009
, “
MEMS Bandpass Filters Based on Cyclic Coupling Architectures
,”
ASME
Paper No. DETC2009-87059.
6.
Judge
,
J.
,
Houston
,
B.
,
,
D.
, and
Herdic
,
P.
,
2006
, “
Effects of Disorder in One- and Two-Dimensional Micromechanical Resonator Arrays for Filtering
,”
J. Sound Vib.
,
290
(
3–5
), pp.
1119
1140
.
7.
Small
,
J.
,
Arif
,
M.
,
Fruehling
,
A.
, and
Peroulis
,
D.
,
2013
, “
A Tunable Miniaturized RF MEMS Resonator With Simultaneous High Q (500-735) and Fast Response Speed (10-60)
,”
J. Microelectromech. Syst.
,
22
(
2
), pp.
395
405
.
8.
Lin
,
Y.
,
Li
,
W.-C.
,
Kim
,
B.
,
Lin
,
Y.-W.
,
Ren
,
Z.
, and
Nguyen
,
C. T.-C.
,
2009
, “
Enhancement of Micromechanical Resonator Manufacturing Precision Via Mechanically-Coupled Arraying
,”
IEEE
International Frequency Control Symposium, Joint With the 22nd European Frequency and Time Forum
, Besancon, France, Apr. 20–24, pp.
58
63
.
9.
Kharrat
,
C.
,
Colinet
,
E.
,
Duraffourg
,
L.
,
Hentz
,
S.
,
Andreucci
,
P.
, and
Voda
,
A.
,
2010
, “
Modal Control of Mechanically Coupled Nems Arrays for Tunable RF Filters
,”
IEEE Trans. Ultrason., Ferroelectr., Freq. Control
,
57
(
6
), pp.
1285
1295
.
10.
,
V.
,
Junghare
,
R.
,
Patrikar
,
R.
, and
Kraft
,
M.
,
2016
, “
Mechanically Coupled Ring-Resonator Filter and Array (Analytical and Finite Element Model)
,”
IET Comput. Digital Tech.
,
10
(
5
), pp.
261
267
.
11.
Nguyen
,
C.-C.
,
2007
, “
MEMS Technology for Timing and Frequency Control
,”
IEEE Trans. Ultrason., Ferroelectr., Freq. Control
,
54
(
2
), pp.
251
270
.
12.
Acar
,
C.
, and
Shkel
,
A.
,
2005
, “
An Approach for Increasing Drive-Mode Bandwidth of MEMS Vibratory Gyroscopes
,”
J. Microelectromech. Syst.
,
14
(
3
), pp.
520
528
.
13.
Haronian
,
D.
, and
MacDonald
,
N.
,
1996
, “
A Microelectromechanics-Based Frequency-Signature Sensor
,”
Sens. Actuators, A: Phys.
,
53
(
1–3
), pp.
288
298
.
14.
Boisen
,
A.
, and
Thundat
,
T.
,
2009
, “
Design and Fabrication of Cantilever Array Biosensors
,”
Mater. Today
,
12
(
9
), pp.
32
38
.
15.
Boisen
,
A.
,
Dohn
,
S.
,
Keller
,
S.
,
Schmid
,
S.
, and
Tenje
,
M.
,
2011
, “
Cantilever-like Micromechanical Sensors
,”
Rep. Prog. Phys.
,
74
(
3
), p.
036101
.
16.
Ziegler
,
C.
,
2004
, “
Cantilever-Based Biosensors
,”
Anal. Bioanal. Chem.
,
379
(
7–8
), pp.
946
959
.
17.
Waggoner
,
P.
, and
,
H.
,
2007
, “
Micro- and Nanomechanical Sensors for Environmental, Chemical, and Biological Detection
,”
Lab Chip
,
7
(
10
), pp.
1238
1255
.
18.
Raiteri
,
R.
,
Grattarola
,
M.
,
Butt
,
H.-J.
, and
Skldal
,
P.
,
2001
, “
Micromechanical Cantilever-Based Biosensors
,”
Sens. Actuators, B: Chem.
,
79
(
2–3
), pp.
115
126
.
19.
Lavrik
,
N. V.
,
Sepaniak
,
M. J.
, and
Datskos
,
P. G.
,
2004
, “
Cantilever Transducers as a Platform for Chemical and Biological Sensors
,”
Rev. Sci. Instrum.
,
75
(
7
), pp.
2229
2253
.
20.
,
J. F.
,
DeMartini
,
B. E.
,
Shaw
,
S. W.
, and
Turner
,
K. L.
,
2006
, “
A SISO, Multi-Analyte Sensor Based on a Coupled Microresonator Array
,”
ASME
Paper No. IMECE2006-13693.
21.
DeMartini
,
B.
,
,
J.
,
Zielke
,
M.
,
Owen
,
K.
,
Shaw
,
S.
, and
Turner
,
K.
,
2008
, “
A Single Input-Single Output Coupled Microresonator Array for the Detection and Identification of Multiple Analytes
,”
Appl. Phys. Lett.
,
93
(
5
), p.
054102
.
22.
Yabuno
,
H.
,
Seo
,
Y.
, and
Kuroda
,
M.
,
2013
, “
Self-Excited Coupled Cantilevers for Mass Sensing in Viscous Measurement Environments
,”
Appl. Phys. Lett.
,
103
(
6
), p.
063104
.
23.
Ryan
,
T.
,
Judge
,
J.
,
Vignola
,
J.
, and
Glean
,
A.
,
2012
, “
Noise Sensitivity of a Mass Detection Method Using Vibration Modes of Coupled Microcantilever Arrays
,”
Appl. Phys. Lett.
,
101
(
4
), p.
043104
.
24.
Glean
,
A. A.
,
Vignola
,
J. F.
,
Judge
,
J. A.
, and
Ryan
,
T. J.
,
2013
, “
Impact of Mass Ratio and Bandwidth on Apparent Damping of a Harmonic Oscillator With Subordinate Oscillator Array
,” ASA International Congress on Acoustics (ICA), Montreal, QC, Canada.
25.
Torres
,
F.
,
Uranga
,
A.
, and
Barniol
,
N.
,
2014
, “
Multi-Cantilever Oscillator
,”
Procedia Eng.
,
87
, pp.
32
35
.
26.
Spletzer
,
M.
,
Raman
,
A.
,
Sumali
,
H.
, and
Sullivan
,
J.
,
2008
, “
Highly Sensitive Mass Detection and Identification Using Vibration Localization in Coupled Microcantilever Arrays
,”
Appl. Phys. Lett.
,
92
(
11
), p.
114102
.
27.
Thiruvenkatanathan
,
P.
, and
Seshia
,
A. A.
,
2012
, “
Mode-Localized Displacement Sensing
,”
J. Microeletromech. Syst.
,
21
(
5
), pp.
1016
1018
.
28.
Oguchi
,
T.
,
Hayase
,
M.
, and
Hatsuzawa
,
T.
,
2005
, “
Micromachined Display Device Using Sheet Waveguide and Multicantilevers Driven by Electrostatic Force
,”
IEEE Trans. Ind. Electron.
,
52
(
4
), pp.
984
991
.
29.
Guillon
,
S.
,
Saya
,
D.
,
Mazenq
,
L.
,
Perisanu
,
S.
,
Vincent
,
P.
,
Lazarus
,
A.
,
Thomas
,
O.
, and
Nicu
,
L.
,
2011
, “
Effect of Non-Ideal Clamping Shape on the Resonance Frequencies of Silicon Nanocantilevers
,”
Nanotechnol.
,
22
(
24
), p.
245501
.
30.
Judge
,
J. A.
,
Woods
,
T. J.
, and
Vignola
,
J. F.
,
2009
, “
Considerations for Use of Square-Paddle Resonators for Arrays of Micro- and Nanoscale Devices
,”
ASME
Paper No. DETC2009-87441.
31.
Choubey
,
B.
,
Boyd
,
E.
,
Armstrong
,
I.
, and
Uttamchandani
,
D.
,
2012
, “
Determination of the Anisotropy of Young's Modulus Using a Coupled Microcantilever Array
,”
J. Microeletromech. Syst.
,
21
(
5
), pp.
1252
1260
.
32.
Sabater
,
A.
, and
,
J.
,
2015
, “
Dynamics of Globally and Dissipatively Coupled Resonators
,”
ASME J. Vib. Acoust.
,
137
(
2
), p.
021016
.
33.
Sato
,
M.
,
Hubbard
,
B.
, and
Sievers
,
A.
,
2006
, “
Colloquium: Nonlinear Energy Localization and Its Manipulation in Micromechanical Oscillator Arrays
,”
Rev. Mod. Phys.
,
78
(
1
), p.
137
.
34.
Buks
,
E.
, and
Roukes
,
M.
,
2002
, “
Electrically Tunable Collective Response in a Coupled Micromechanical Array
,”
J. Microeletromech. Syst.
,
11
(
6
), pp.
802
807
.
35.
Krylov
,
S.
,
Lulinsky
,
S.
,
Ilic
,
B. R.
, and
Schneider
,
I.
,
2014
, “
Collective Dynamics of Arrays of Micro Cantilevers Interacting Through Fringing Electrostatic Fields
,”
ASME
Paper No. DETC2014-34904.
36.
Krylov
,
S.
,
Lulinsky
,
S.
,
Ilic
,
B.
, and
Schneider
,
I.
,
2014
, “
Collective Dynamics and Pattern Switching in an Array of Parametrically Excited Micro Cantilevers Interacting Through Fringing Electrostatic Fields
,”
Appl. Phys. Lett.
,
105
(
7
), p.
071909
.
37.
Ono
,
T.
,
Tanno
,
K.
, and
Kawai
,
Y.
,
2014
, “
Synchronized Micromechanical Resonators With a Nonlinear Coupling Element
,”
J. Micromech. Microeng.
,
24
(
2
), p.
025012
.
38.
,
J.
,
Park
,
H.
, and
Zewail
,
A.
,
2011
, “
Nanomusical Systems Visualized and Controlled in 4D Electron Microscopy
,”
Nano Lett.
,
11
(
5
), pp.
2183
2191
.
39.
Linzon
,
Y.
,
Ilic
,
B.
,
Lulinsky
,
S.
, and
Krylov
,
S.
,
2013
, “
Efficient Parametric Excitation of Silicon-on-Insulator Microcantilever Beams by Fringing Electrostatic Fields
,”
J. Appl. Phys.
,
113
(
16
), p.
163508
.
40.
Meirovitch
,
L.
,
2001
,
Fundamentals of Vibrations
,
McGraw-Hill
,
Boston, MA
.
41.
Gutschmidt
,
S.
, and
Gottlieb
,
O.
,
2012
, “
Nonlinear Dynamic Behavior of a Microbeam Array Subject to Parametric Actuation at Low, Medium and Large Dc-Voltages
,”
Nonlinear Dyn.
,
67
(
1
), pp.
1
36
.
42.
Kambali
,
P.
,
Swain
,
G.
, and
Pandey
,
A.
,
2016
, “
Frequency Analysis of Linearly Coupled Modes of MEMS Arrays
,”
ASME J. Vib. Acoust.
,
138
(
2
), p.
021017
.
43.
Lifshitz
,
R.
, and
Cross
,
M.
,
2003
, “
Response of Parametrically Driven Nonlinear Coupled Oscillators With Application to Micromechanical and Nanomechanical Resonator Arrays
,”
Phys. Rev. B
,
67
(
13
), p.
134302
.
44.
Dick
,
A.
,
Balachandran
,
B.
, and
Mote
,
C.
, Jr.
,
2008
, “
Intrinsic Localized Modes in Microresonator Arrays and Their Relationship to Nonlinear Vibration Modes
,”
Nonlinear Dyn.
,
54
(
1–2
), pp.
13
29
.
45.
Blocher
,
D. B.
,
2012
, “
Optically Driven Limit Cycle Oscillations in MEMS
,” Ph.D. thesis, Cornell University, Ithaca, NY.