Actuation of Higher Harmonics in Large Arrays of Micromechanical Cantilevers for Expanded Resonant Peak Separation

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 [3–11]. 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 [15–19], sensing within arrays of weakly coupled resonators is based on modal shape changes [20–25]. 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 (f_L) and the upper (f_U) cutoff values. While most of the previous studies were focused on arrays of identical or almost identical beams [23,26,33–36], 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 f_L and f_U 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.

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, L_o ≈ 100 μm, pitch B ≈ 50 μm, L₁ = L_max ≈ 500 μm, and L_N = L_min ≈ 350 μm.

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. 2–4. 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 f_i, 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 L₁ = 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.

where A = bh and I = bh³/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/m³ 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 L_n decreases with increasing n, the frequency shift increases with the mode number within each propagation band.

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,41–43]. 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, w_n 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 l_n = L₁/L_n 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 $ŷn$ along each of the beams is normalized by the length L_n of the corresponding beam, we have $y∈[0,1]$⁠.

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].

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.

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)$⁠.

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.

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].

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.

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

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 Q₁ = 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 L₂₅ = 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).

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⁻⁴ Pa (≈ 10⁻⁶ 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.

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 L₁ and L₂₅ 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 L₁ and L₂₅ 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 L₂₅, 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 L₂₅. 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 L_o. 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.

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.

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.

• 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).