Modal Behavior of Microcantilevers Arrays with Tunable Electrostatic Coupling

Abstract

We analyse the spectral content and parametric resonant dynamics of an array of elastically and electrostatically coupled interdigitated micro cantilevers assembled into two identical half-arrays. In this uncommon arrangement, within each of the half-arrays, the beams are coupled only elastically. The half-arrays are intercoupled only electrostatically, through fringing fields. First, by using the reduced order (RO) model, we analyse the voltage-dependent evolution of the eigenvalues and the eigenvectors of the equivalent mass-spring system, starting from the small two, three and four beams arrays and up to large beams assemblies. We show that at the coupling voltages below a certain critical value, the shape of the eigenvectors, the frequencies of the veering and of the crossing are influenced by the electrostatic coupling and can be tuned by the voltage. Next, by implementing the assumed modes techniques we explore the parametric resonant behavior of the array. We show that in the case of the sub critical electrostatic coupling the actuating voltages required to excite parametric resonance in the damped system can be lower than in a strongly coupled array. The results of the work may inspire new designs of more efficient resonant sensors.

1. Introduction

The dynamics of large arrays of coupled micro- or nanoelectromechanical (MEMS/NEMS) structural elements have been intensively studied theoretically and experimentally for more than two decades [1] (e.g., see reviews [2,3,4], see [5,6,7,8,9,10] for more recent contributions). Micro- and nanodevices manifest rich behavior and may serve as a convenient platform for investigation of exciting dynamic phenomena, such as parametric resonances (PR) [7,11,12,13,14,15,16,17], intrinsic localized modes [4,18,19], pattern selection [20,21] and synchronization [8], to name just a few. Within the applied physics and engineering communities, the interest in microscale arrays is motivated by their applicability as mass, chemical and biological multi-point sensors [5,9,10,22,23,24,25,26], efficient band pass filters or active metamaterials.

In most of the cases, the MEMS arrays are realized as assemblies of cantilevered or double clamped beams, attached at their ends to a thin flexible plate (an “overhang”) serving as a source of the elastic coupling [4,26]. In these devices, both the onsite (OS) and intersite (IS) stiffness, along with the OS and IS nonlinearities are of the mechanical nature, are fully dictated by design and cannot be varied (stiffness tuning using an integrated piezoelectric actuator was reported in [19,27]; stiffness modulation by laser was explored in [28]). In contrast, electrostatic (ES) coupling makes it possible to change the effective OS and IS stiffness of the array [29,30] and opens new possibilities for a plethora of operational scenarios. Note that while frequency tuning in ES actuated microarrays was intensively investigated, much less attention was paid to the influence of the nonlinear ES forces on the eigenvectors of the arrays. In the architectures based on the beams moving towards each other [12,13,14,20,31], the attractive ES force has a softening effect on the stiffness. As a result, the ES component of the IS stiffness of the equivalent mass-spring system is negative [11].

An alternative design introduced in [16,17] and further studied in [7,21] is based on the beams interacting through fringing ES fields. In these devices, the ES force is restoring and the IS stiffness increases with the voltage. For convenience, the architecture of this structure is briefly reviewed here. The array shown in Figure 1 is composed of N interdigitated cantilevers and is configured as two partially overlapping half-arrays with the voltage difference V applied between them. Each half-array contains $N/2$ beams with the rectangular cross section of the width b and thickness h. The length $L_n$ of the cantilevers can be uniform along the array [17,21] or may vary between $L_{min}$ and $L_{max}$ [7,32]. The beams, designed to deflect in the vertical, z direction, are attached to a compliant overhang of the length $L_o$ and of the thickness h identical to that of the beams. The distance between the centerlines of the neighboring beams within the half-arrays is $2B$, the pitch between the adjacent beams of the array is B. The length of the overlap between the beams of the opposite half-arrays (the electrode length) is $L_e$, and the gap between the beams is g. Unmovable fixed beams positioned at the ends of the array serve as the end electrodes and are provided to enforce the zero-deflection of the endpoints of the array [12,17]. The devices reported in [7,17,32] were fabricated from single crystal silicon using a silicon-on-insulator (SOI) wafer with the device layer of the thickness h.

One of the key, distinguishing, features of the device shown in Figure 1 is the way it combines, in an unusual manner, both the elastic and the ES couplings. Specifically, in absence of the voltage, the cantilevers of the opposite half-arrays are decoupled while the beams of the same half-arrays are coupled only elastically. Application of the voltage results in the appearance of the ES coupling between the two sets of the cantilevers. As we show here, this mirror-type symmetry and well-defined parity may result in interesting modal behavior, when the odd and the even eigenvectors do not mix. Note that an array of similar architecture, which was excited parametrically was considered in [17], the further investigation of the PR dynamics in these arrays was reported in [21]. It was shown that the ES coupling between the half-arrays affects both the natural frequencies and the eigenvectors. However, the influence of the ES coupling on the arrays spectral content, veering and crossing phenomena, were not explored in [21]. The effect of the time-independent (DC) voltage on the natural frequencies in an array incorporating mechanically decoupled electrostatically interacting beams (an architecture different than ours) were studied in [29,30]. While the frequency tuning and veering were demonstrated, the crossing and the voltage-dependent evolution of the eigenvectors was not addressed in these works.

In the first part of the present work, the influence of the ES coupling on the eigenfrequencies and eigenvectors of the array represented as the associated linearized mass-spring chain is explored. In the second part, the resonant parametric dynamics of the array are studied using the RO model. The eigenvectors obtained numerically for the specific values of the ES coupling voltage are used as the base vectors. The analytical model predictions are compared to the numerical results. The main findings of the work are summarized in the conclusions.

2. Formulation

In our development, we follow the general lines of the description in [21,32]. The vibrational dynamics of each beam are described in the framework of the Euler–Bernoulli model, the electrostatic force between the cantilevers is approximated by the Pad’e fit. Implementation of the Galerkin order reduction procedure results is the following system of nonlinear ordinary differential equations in terms of the non-dimensional endpoint deflections $u_n\left(t\right)$ of the beams

$$\begin{array}{c}\hfill {\ddot{u}}_n+\frac{l_n^2}{Q_n}{\dot{u}}_n+l_n^4{u}_n+{\epsilon }^{nl}{u}_n^3-\\ \hfill \eta (u_{n+2}-2u_n+u_{n-2})& -& {\eta }^{nl}\left[{\left(u_{n+2}-u_n\right)}^3-{\left(u_n-u_{n-2}\right)}^3\right]-\hfill \\ \hfill \beta V^2\left[\frac{(u_{n+1}-u_n)}{1+{(u_{n+1}-u_n)}^2}\right.& -& \left.\frac{(u_n-u_{n-1})}{1+{(u_n-u_{n-1})}^2}\right]=0,(n=2\dots N-1),\hfill \end{array}$$

(1)

where $n=1,...,N$ represents the vibrating beams numbers. Note that the first $n=1$ and the last $n=N$ equations corresponding to the first and the last beams are modified to take into account the presence of the fixed electrodes at the two ends of the array. The over-dot $\dot{(•)}=d/dt$ denotes the derivative with respect to the non-dimensional time, $l_n=L_{max}/L_n\ge 1$ is the length ratio parameter, where $L_{max}$ is the length of the longest beam of the array; $Q_n$ is the quality factor associated with the OS damping, and ${\epsilon }^{nl}$ is the OS nonlinearity coefficient. The coefficients $\eta$ and ${\eta }^{nl}$ are associated with the linear and nonlinear closest neighbor (within each half-array) elastic IS coupling [4], respectively. The last line in Equation (1) represent the IS electrostatic force. The geometrical dimensions of the array used in all the calculations in the present work are as in [17,32],

Table 1. Dimensions of the array.

Quantity	Description	Value (μm)
$L_o$	Overhang length	100
$L_{max}$	The longest beam length	500
$L_{min}$	The shortest beam length	350
$b$	Beam width	16
$h$	Thickness of the beam/overhang	5
$L_e$	Electrode length (overlap)	150
$g$	Gap between the beams	5
$B$	Pitch between the adjacent beams	21

. The corresponding non-dimensional parameters are: ${\epsilon }^{nl}=1.606\times {10}^{-4}$, $\eta =0.0548$ ${\eta }^{nl}=\eta /10$; the latter was obtained by re-scaling the ratio ${\eta }^{nl}/\eta$ from [4] (see also [33]) in conjunction with the actual device geometry; $\beta =4.35\times {10}^{-5}$ is the voltage parameter.

Taylor series expansion of the last term of Equation (1) (representing the ES force) up to the cubic order of $u_n\left(t\right)$ yields

$$\begin{array}{c}{\ddot{u}}_n+\frac{l_n^2}{Q_n}{\dot{u}}_n+l_n^4{u}_n+{\epsilon }^{nl}{u}_n^3-\eta \left(u_{n+2}-2u_n+u_{n-2}\right)-\hfill \\ \beta V^2\left(u_{n+1}-2u_n+u_{n-1}\right)-{\eta }^{nl}\left[{\left(u_{n+2}-u_n\right)}^3-{\left(u_n-u_{n-2}\right)}^3\right]+\hfill \\ \beta V^2\left[{\left(u_{n+1}-u_n\right)}^3-{\left(u_n-u_{n-1}\right)}^3\right]=0,(n=2\dots N-1).\hfill \end{array}$$

(2)

Equations (1) and (2) have clear physical interpretation and describe a chain (a lattice) of masses attached to a substrate and coupled by linear and nonlinear elastic and effective ES springs. To further highlight the contribution of various parameters influencing the array response, it is instructive to consider the linearized undamped counterpart of Equation (2)

$$\ddot{u_n}+l_n^4{u}_n-\eta \left(u_{n+2}-2u_n+u_{n-2}\right)-\beta V^2\left(u_{n+1}-2u_n+u_{n-1}\right)=0.$$

(3)

Note that (in contrast to [12,13,14,20,31]) both elastic and electrostatic IS forces are restoring and the associated coefficients $\eta$ and $\beta$ are positive. Since one of our goals is to analyze the influence of the ES coupling on the array dynamics, for clarity and completeness, we briefly present here the main results of the dispersion analysis of the linearized system, Equation (3). Assuming an infinite uniform ($l_n=1$) array, the steady, time-independent voltage $V=V_{\mathrm{dc}}$ and substituting $u_n=e^{ikn-i\omega t}$ (where $\omega$ is the non-dimensional frequency, $\kappa =2\pi B/\lambda$ is the non-dimensional wavenumber and $\lambda$ is the wavelength) into Equation (3) allows us to obtain the dispersion relation of the array [21]

$$\omega =\sqrt{1+2\eta (1-cos\left(2\kappa \right))+2\beta V_{\mathrm{dc}}^2(1-cos\left(\kappa \right))}.$$

(4)

Figure 2 depicts the dispersion curves $\omega =\omega \left(\kappa \right)$ for several values of the coupling voltage $V_{\mathrm{dc}}$. The propagation band is bounded by the lower ${\omega }_L=1$ and the tunable upper ${\omega }_U$ cut-off frequencies. When $V_{\mathrm{dc}}=0$, the coupling is purely elastic and ${\omega }_U={\omega }_U^{MECH}=\sqrt{1+4\eta }$. The corresponding wavenumber is ${\kappa }_U=\pi /2$. Since $\eta$ is dictated by the material properties and the geometry, ${\omega }_U$ cannot be changed without re-design. With the increase in $V_{\mathrm{dc}}$, ${\omega }_U$ increases in accordance with Equation (4) where $\kappa$ is replaced by its upper cut-off value

$${\kappa }_U=\pi -arccos\left(\frac{\beta V_{\mathrm{dc}}^2}{4\eta }\right),\left(\frac{\pi }{2}\le \kappa \le \pi \right).$$

(5)

At certain critical finite value $V_{\mathrm{dc}}=V_{\mathrm{crit}}=2\sqrt{\eta /\beta }=70.98$ V, the wavenumber reaches its maximum ${\kappa }_U=\pi$. The corresponding upper cut-off ${\omega }_U={\omega }_U^{\mathrm{crit}}=\sqrt{1+4\beta V_{\mathrm{crit}}^2}$ is independent of the mechanical coupling parameter $\eta$ (since within each half-array all $u_n$ are identical; see Figure 2). For $V_{\mathrm{dc}}\ge V_{\mathrm{crit}}$, the electrostatic coupling is dominant, and, while $\omega \ge {\omega }_U^{\mathrm{crit}}$ increases with increasing $V_{\mathrm{dc}}$, ${\kappa }_U=\pi$ is voltage-independent.

This is the case, hereafter referred as the “fully coupled” or “supercritical”, which was analyzed in [21]. It was shown that at $V_{\mathrm{dc}}>V_{\mathrm{crit}}$ the increase in the voltage has only minor influence on the shape of the eigenvectors. For this reason, it was possible to use the eigenvectors of a uniform linear mass-spring chain, described by simple analytical expressions, as the base vectors in the assumed modes ansatz. In contrast, when $V_{\mathrm{dc}}\le V_{\mathrm{crit}}$ the ES coupling may have a strong influence on the eigenvectors shape. In this case the use of the simple eigenvectors adopted in [21] as the base vectors in the order reduction procedure can be inappropriate. The present work is focused on the array dynamics in the case of weak, not fully developed, ES coupling, when $V_{\mathrm{dc}}\le V_{\mathrm{crit}}$.

3. Eigenvalues and Eigenvectors of Electrostatically Coupled Arrays

In this section, several cases of small arrays assembled of two, three and four beams are first considered followed by the analysis of the mid size arrays containing 8 and 16 beams. The results for the large 200 beams array are also presented. We show that interesting phenomena, such as eigenvalue crossing and veering, occur during the voltage increase. The analysis of the voltage-dependent eigenvectors is important in order to lie foundations for the appropriate reduced order modeling methodology of the arrays.

3.1. Two-Beams System

Two ideally clamped cantilevers (with additional fixed “electrodes” at each side) of the length $L_1=L_{max}$ and $L_2=L_1-\Delta L$ are coupled only through the electrostatic force, Figure 3a. In this case $\eta =0$, and recalling that $l_1=1$, Equation (3) takes the form

$$\begin{array}{c}{\ddot{u}}_1+\left(1+2\beta V^2\right)u_1-\beta V^2{u}_2=0,\hfill \\ {\ddot{u}}_2-\beta V^2{u}_1+\left(l_2^4+2\beta V^2\right)u_2=0.\hfill \end{array}$$

(6)

We re-scale time $\tau =t\sqrt{1+2\beta V^2}$ and obtain

$$\begin{array}{c}{\ddot{u}}_1+u_1-\gamma u_2=0,\hfill \\ {\ddot{u}}_2-\gamma u_1+{\chi }_2{u}_2=0,\hfill \end{array}$$

(7)

where the overdot is re-defined as $\dot{\left(•\right)}=d/d\tau$. The ES coupling parameter $0\le \gamma <1/2$ and the detuning parameter ${\chi }_2\ge 1$ defined in

Table 2. Non-dimensional parameters used in the development.

Parameter	Two Beams Array	Three+ Beams Array
ES coupling	$\gamma =\beta V^2/(1+2\beta V^2)$	$\gamma =\beta V^2/(1+2\eta +2\beta V^2)$
Detuning	${\chi }_2=(l_2^4+2\beta V^2)/(1+2\beta V^2)$	${\chi }_n=(l_n^4+2\eta +2\beta V^2)/(1+2\eta +2\beta V^2)$
Elastic coupling		$\tilde{\eta }=\eta /(1+2\eta +2\beta V^2)$

are not independent since they are both parameterized by the voltage.

The eigenvalues analysis yields

$$\begin{array}{ccc}\hfill {\lambda }_r& =& \frac{1}{2}\left[1+{\chi }_2+{(-1)}^r\sqrt{4{\gamma }^2+{({\chi }_2-1)}^2}\right],r=1,2,\hfill \end{array}$$

(8)

$$\begin{array}{ccc}\hfill {\lambda }_1& =& 1,{\lambda }_2={\chi }_2,(\mathrm{uncoupled}\mathrm{case}\gamma =0),\hfill \end{array}$$

(9)

$$\begin{array}{ccc}\hfill {\lambda }_1& =& 1-\gamma ,{\lambda }_2=1+\gamma ,(\mathrm{symmetric}\mathrm{case}{\chi }_2=1).\hfill \end{array}$$

(10)

The influence of the ES coupling $\gamma$ and of the detuning ${\chi }_2$ on the eigenvalues is shown in Figure 4a,b, respectively, their influence on the eigenvalues ratio ${\lambda }_2/{\lambda }_1$ are depicted in Figure 4c,d. As expected [34], crossing, namely the situation when ${\lambda }_1={\lambda }_2$, is possible only in the symmetric uncoupled case ${\chi }_2=1,\gamma =0$. Presence of coupling results in the frequency split (veering). In accordance with Table 2 and Equation (8), at large voltages, $\gamma \to 1/2$, ${\chi }_2\to 1$ and ${\lambda }_2/{\lambda }_1\to 3$, suggesting that the influence of the detuning on the array spectral content decreases with the increasing voltage. This feature could be useful in applications as a means to overcome undesirable detuning due to the fabrication-related deviations in the beams geometry.

3.2. Three Beams Array

Three beams system, Figure 3b, is the smallest array where both the mechanical and the electrostatic coupling are present and the frequency crossing could be possible at $\gamma \ne 0$. The non-dimensional equations of motion are

$$\begin{array}{ccc}\hfill {\ddot{u}}_1+\left(1+2\eta +2\beta V^2\right)u_1-\beta V^2{u}_2-\eta u_3& =& 0,\hfill \\ \hfill {\ddot{u}}_2-\beta V^2{u}_1+\left(l_2^4+2\eta +2\beta V^2\right)u_2-\beta V^2{u}_3& =& 0,\hfill \\ \hfill {\ddot{u}}_3-\eta u_1-\beta V^2{u}_2+\left(l_3^4+2\eta +2\beta V^2\right)u_3& =& 0.\hfill \end{array}$$

(11)

By normalizing the time by $\sqrt{1+2\eta +2\beta V^2}$, we obtain the characteristic determinant

$$\begin{array}{c}\left|\begin{array}{ccc}1-\lambda & -\gamma & -\tilde{\eta }\\ -\gamma & {\chi }_2-\lambda & -\gamma \\ -\tilde{\eta }& -\gamma & {\chi }_3-\lambda \end{array}\right|=0,\hfill \end{array}$$

(12)

where the non-dimensional parameters $\tilde{\eta },\gamma ,{\chi }_n,n=1,2,3$ are given in Table 2. The characteristic equation associated with Equation (12) is cubic in terms of $\lambda$, so the eigenvalues ${\lambda }_r,r=1\dots 3$ can be found analytically in a closed form. In a general case, these expressions are too cumbersome to be presented explicitly. Here, the results are summarized in

Table 3. Case study—eigenvalues of the three beams arrays. Non-dimensional parameters are presented in Table 2.

Case	Uniform Array	Symmetric Array
	${\mathbf{\chi }}_{\mathbf{n}}=\mathbf{1}$	${\mathbf{\chi }}_{\mathbf{1}}={\mathbf{\chi }}_{\mathbf{3}}=\mathbf{1},{\mathbf{\chi }}_{\mathbf{2}}>\mathbf{1}$
${\lambda }_1$	$1-\frac{\tilde{\eta }+\sqrt{{\eta }^2+8{\gamma }^2}}{2}$	$\frac{1}{2}\left(1+{\chi }_2-\tilde{\eta }-\sqrt{{\left(\tilde{\eta }+{\chi }_2-1\right)}^2+8{\gamma }^2}\right)$
${\lambda }_2$	$1+\tilde{\eta }$	$1+\tilde{\eta }$
${\lambda }_3$	$1-\frac{\tilde{\eta }-\sqrt{{\eta }^2+8{\gamma }^2}}{2}$	$\frac{1}{2}\left(1+{\chi }_2-\tilde{\eta }+\sqrt{{\left(\tilde{\eta }+{\chi }_2-1\right)}^2+8{\gamma }^2}\right)$
Case	Pure mechanical coupling	Symmetric array, pure ES coupling
	$\gamma =0$	$\eta =0$, ${\chi }_1={\chi }_3=1,{\chi }_2>1$
${\lambda }_1$	$\frac{1}{2}\left({\chi }_3+1+\sqrt{4{\tilde{\eta }}^2+{({\chi }_3-1)}^2}\right)$	$\frac{1}{2}\left(1+{\chi }_2-\sqrt{{\left({\chi }_2-1\right)}^2+8{\gamma }^2}\right)$
${\lambda }_2$	${\chi }_2$	1
${\lambda }_3$	$\frac{1}{2}\left({\chi }_3+1-\sqrt{4{\tilde{\eta }}^2+{({\chi }_3-1)}^2}\right)$	$\frac{1}{2}\left(1+{\chi }_2+\sqrt{{\left({\chi }_2-1\right)}^2+8{\gamma }^2}\right)$

for four particular cases: (i) a uniform array ${\chi }_n=1$; (ii) a symmetric array ${\chi }_1={\chi }_3=1$, ${\chi }_2>1$; (iii) a general array with purely mechanical coupling and $\gamma =0$; (iv) a symmetric array ${\chi }_1={\chi }_3=1$, ${\chi }_2>1$ with purely electrostatic coupling and $\eta =0$.

The influence of the ES coupling on the eigenvalues and the eigenvectors is illustrated in Figure 5. The evolution of the ratios ${\lambda }_r/{\lambda }_1,r=2,3$ with the change of $\gamma$ (as usual, ${\lambda }_r\le {\lambda }_{r+1}$) in the case of a uniform (${\chi }_n=1$) and symmetric non-uniform (${\chi }_1={\chi }_3=1$, ${\chi }_2>1$) array are shown in Figure 5a. Since ${\lambda }_2,{\lambda }_3$ both grow at a different rate with increasing $\gamma$, eigenvalue crossing ${\lambda }_2={\lambda }_3$ may take place at a certain value of $\gamma ={\gamma }_{\times }$. Figure 5a suggests that the crossing values of ${\lambda }_r$ are influenced by the detune and the voltage. In general, the crossing occurs when $\gamma$ follows a certain path on the $\gamma -{\chi }_2$ plane, as illustrated in Figure 5b. For a given $\eta$ and ${\chi }_2$, ${\gamma }_{\times }=\sqrt{\tilde{\eta }\left(\tilde{\eta }-{\chi }_2+1\right)}$. By substituting into this relation the expressions for $\tilde{\eta }$ and ${\chi }_n$ from Table 2 we obtain $V_{\times }=\sqrt[4]{\eta (1+\eta -l_n^4)}/\sqrt{\beta }$. This result suggests that the crossing is possible when $\eta \ge l_n^4-1$. One may conclude therefore that in the arrays with higher initial detuning, higher coupling is required to achieve crossing. In the uniform array where ${\chi }_2=1$, the crossing occurs at $\gamma ={\gamma }_{\times }=\tilde{\eta }\approx 0.045$, or, in terms of the voltage, $V_{\times }=\sqrt{\eta /\beta }$. This value is smaller than the critical voltage $V_{\mathrm{crit}}=2\sqrt{\eta /\beta }$ (hereafter also denoted as $V_{\mathrm{crit}}^{\inf }$), separating the weakly and the fully coupled cases in the infinite array. For the adopted parameters, Table 1, we obtain $V_{\times }\approx 35.5$ V.

The full sets of the eigenvectors of a uniform array for two values of $\gamma$ below and above the crossing ${\gamma }_{\times }$ are shown in Figure 5c. The eigenvectors ${\Psi }^{\left(r\right)},r=1,2,3$ are ordered based on the increasing associated eigenvalues ${\lambda }_r$. For further discussion, we also adopt an alternative ordering based on the number of nodes m. The nodes are the geometric points (not including the array’s ends) where the imaginary line connecting the masses of the array (dashed lines in Figure 5c) crosses zero. Hereafter, to distinguish between the two ordering approaches, the notation ${}^{\left(m\right)}\Psi ,m=0,\dots ,N-1$ is adopted to denote the eigenvectors sorted by m. Figure 5c suggests that the order of ${}^{\left(m\right)}\Psi$ may change, depending on the voltage. At the value of $\gamma =0.043<{\gamma }_{\times }=0.045$ the second eigenvector ${\mathbf{\Psi }}^{\left(2\right)}={}^{\left(2\right)}\mathbf{\Psi }$ is symmetric, while the third ${\mathbf{\Psi }}^{\left(3\right)}={}^{\left(1\right)}\mathbf{\Psi }$ is anti symmetric. At $\gamma =0.047>{\gamma }_{\times }$, as expected, the eigenvectors ${\mathbf{\Psi }}^{\left(2\right)}$ and ${\mathbf{\Psi }}^{\left(3\right)}$ flip over. The situation is further illustrated in Figure 5d,e where the evolution of the eigenvectors ${\mathbf{\Psi }}^{\left(r\right)},{}^{\left(m\right)}\mathbf{\Psi }$ with $\gamma$ is shown. Note that while the eigenvectors are discrete, their visualization as a continuous function is convenient for the illustration of the eigenvectors shape evolvement with $\gamma$. In accordance with Figure 5d,e, in ${\mathbf{\Psi }}^{\left(1\right)}={}^{\left(0\right)}\mathbf{\Psi }$ within the entire range of $\gamma$. In Figure 5d, the discontinuity in ${\mathbf{\Psi }}^{\left(2\right)},{\mathbf{\Psi }}^{\left(3\right)}$ is solely the result of the adopted eigenvectors ordering. The evolution of ${}^{\left(m\right)}\mathbf{\Psi }$ is smooth within the entire range of $\gamma$ including the crossing point, Figure 5e.

3.3. Four Beams Array

The four beams array is the simplest architecture where the beams of each of the half-arrays are coupled elastically while the half-arrays are interacting through the ES forces. In this case the characteristic determinant takes the form

$$\left|\begin{array}{cccc}1-\lambda & -\gamma & -\tilde{\eta }& 0\\ -\gamma & {\chi }_2-\lambda & -\gamma & -\tilde{\eta }\\ -\tilde{\eta }& -\gamma & {\chi }_3-\lambda & -\gamma \\ 0& -\tilde{\eta }& -\gamma & {\chi }_4-\lambda \end{array}\right|=0.$$

(13)

In accordance with Equation (13) the array’s spectral content is defined by five parameters. For simplicity, we choose to focus on (i) symmetric arrays, such that ${\chi }_4={\chi }_1=1$ and ${\chi }_3={\chi }_2$ (see Figure 3c) and (ii) uniform arrays ${\chi }_n=1,n=1\dots 4$. In the symmetric case, the number of parameters is reduced to three. In the case of a uniform array, only two independent parameters $\tilde{\eta }$ and $\gamma$ affect the array’s behavior. The eigenvalues of the symmetric array are

$$\begin{array}{ccc}\hfill & {\lambda }_{1,2}& =\frac{1}{2}\left({\chi }_2+1+\gamma \right)\pm \frac{1}{2}\sqrt{4{(\tilde{\eta }-\gamma )}^2+{\gamma }^2+2\gamma ({\chi }_2-1)+{({\chi }_2-1)}^2},\hfill \\ & {\lambda }_{3,4}& =\frac{1}{2}\left({\chi }_2+1-\gamma \right)\pm \frac{1}{2}\sqrt{4{(\tilde{\eta }+\gamma )}^2+{\gamma }^2-2\gamma ({\chi }_2-1)+{({\chi }_2-1)}^2},\hfill \end{array}$$

(14)

and in the uniform case

$$\begin{array}{ccc}\hfill & {\lambda }_{1,2}& =1+\frac{\gamma }{2}\pm \frac{1}{2}\sqrt{4{(\tilde{\eta }-\gamma )}^2+{\gamma }^2+2},\hfill \\ & {\lambda }_{3,4}& =1-\frac{\gamma }{2}\pm \frac{1}{2}\sqrt{4{(\tilde{\eta }+\gamma )}^2+{\gamma }^2-2}.\hfill \end{array}$$

(15)

Figure 6a shows the evolvement of the eigenvalues ratios ${\lambda }_r/{\lambda }_1,r=1,2,3,4$ of the uniform array with the increase in $\gamma$. In the absence of the ES coupling, $\gamma =0$, two pairs of repeated eigenvalues ${\lambda }_1={\lambda }_2$ and ${\lambda }_3={\lambda }_4$ are obtained. With the increase in $\gamma$, the repeated frequencies split and four different eigenvalues emerge. Dissimilarly to the case of the three beams array, both veering and crossing phenomena are observed here. The crossing between ${\lambda }_3$ and ${\lambda }_4$ occurs at $\gamma ={\gamma }_{\times }=\sqrt{3}\tilde{\eta }$, which corresponds to $V_{\times }=\sqrt{\sqrt{3}\eta /\beta }$ (smaller than $V_{\mathrm{crit}}^{\inf }=2\sqrt{\eta /\beta }$ of the infinite array). 3D rendering of the eigenvalues ratio evolution in the symmetric array ${\chi }_1={\chi }_4=1$ and ${\chi }_2={\chi }_3$ is shown in Figure 6b. The crossing manifests itself as a line in the $\gamma -{\chi }_2$ plane. The full sets of the eigenvectors of a uniform array for $\gamma =0$ and three other values of $\gamma$ in the vicinity of the veering value ${\gamma }_{\asymp }$, as well as slightly below and above the crossing value ${\gamma }_{\times }$ (shown by the vertical dashed lines in Figure 6a), are depicted in Figure 6c. Similarly to the thee-beams case, Figure 5c, the first mode ${\mathbf{\Psi }}^{\left(1\right)}={}^{\left(0\right)}\mathbf{\Psi }$ is symmetric for any $\gamma$. When $\gamma$ passes through the the veering or crossing values, part of the vectors flip over. When $\gamma$ becomes larger than the highest ${\gamma }_{\times }$, the number of nodes m increases monotonically with r (such that ${\mathbf{\Psi }}^{(r-1)}={}^{\left(m\right)}\mathbf{\Psi }$). The evolution of the eigenvectors ${\mathbf{\Psi }}^{\left(r\right)}$ and ${}^{\left(m\right)}\mathbf{\Psi }$ with the increase in the ES coupling is illustrated in Figure 6d,e, respectively. The m-based ordering eliminates the discontinuity in ${\mathbf{\Psi }}^{\left(3\right)},{\mathbf{\Psi }}^{\left(4\right)}$ at $\gamma ={\gamma }_{\times }$. In contrast, in the case of veering, $\gamma ={\gamma }_{\asymp }$, the evolution of ${\mathbf{\Psi }}^{\left(2\right)},{\mathbf{\Psi }}^{\left(3\right)}$ is smooth, Figure 6d, whereas the evolution of ${}^{\left(1\right)}\mathbf{\Psi },{}^{\left(3\right)}\mathbf{\Psi }$ is discontinuous, Figure 6e.

3.4. Large Uniform Arrays

To analyse the behavior of larger arrays, we limit our consideration to the uniform arrays [17,21]. For convenience, the results are presented for smaller eight and sixteen beam arrays. The extension to an even larger array results in qualitatively similar bahavior. In the case of the uniform eight beams array, the characteristic determinant is

$$\begin{array}{c}\left|\begin{array}{cccccccc}1-\lambda & -\gamma & -\tilde{\eta }& 0& 0& 0& 0& 0\\ -\gamma & 1-\lambda & -\gamma & -\tilde{\eta }& 0& 0& 0& 0\\ -\tilde{\eta }& -\gamma & 1-\lambda & -\gamma & -\tilde{\eta }& 0& 0& 0\\ 0& -\tilde{\eta }& -\gamma & 1-\lambda & -\gamma & -\tilde{\eta }& 0& 0\\ 0& 0& -\tilde{\eta }& -\gamma & 1-\lambda & -\gamma & -\tilde{\eta }& 0\\ 0& 0& 0& -\tilde{\eta }& -\gamma & 1-\lambda & -\gamma & -\tilde{\eta }\\ 0& 0& 0& 0& -\tilde{\eta }& -\gamma & 1& -\gamma \\ 0& 0& 0& 0& 0& -\tilde{\eta }& -\gamma & 1-\lambda \end{array}\right|=0,\hfill \end{array}$$

(16)

and contains only two parameters, $\gamma$ and $\tilde{\eta }$.

Figure 7a presents the evolution of the eigenvalues ratios ${\lambda }_r/{\lambda }_1,r=2\dots 8$ with the increase in the ES coupling. Figure 7b presents the “frequency curves”—the full set of ${\lambda }_r/{\lambda }_1$ for several values of $\gamma$. (The lines connecting the markers corresponding to the same $\gamma$ are drawn only for convenience, see Figure 2). In the electrostatically uncoupled array the eigenvalues are repeated pairs; as $\gamma$ increases, the pairs split and the deviation between the eigenvalues of each initially degenerated pair grows, Figure 7b. While the lowest cut-off eigenvalue ${\lambda }_1$ is always equal to unity due to the normalization, an increase in $\gamma$ results in an increase in the upper cut-off ${\lambda }_N$ and, consequently, in the wider propagation band. The evolution of the eigenvectors ${\mathbf{\Psi }}^{\left(r\right)}$ with $\gamma$ is illustrated in Figure 7c. The switching between the eigenvectors is through the veering or crossing mechanism.

Three full sets of the eigenvectors of the 16-beam arrays for three different values of $\gamma$ are depicted in Figure 8. The shape of the first eigenvector ${\mathbf{\Psi }}^{\left(1\right)}={}^{\left(1\right)}\mathbf{\Psi }$ (the first left image in the rows 1, 3, 5 in Figure 8) is qualitatively similar to that of a uniform mass-spring chain with fixed ends. However, since the array is built as two initially decoupled half-arrays, the shape of the higher eigenvectors ${\mathbf{\Psi }}^{\left(r\right)},r=2\dots N$ within the set is affected by the ES coupling. For example, for $\gamma =0$, the shapes of the eigenvectors ${\mathbf{\Psi }}^{\left(r\right)},{\mathbf{\Psi }}^{(r+1)}$ corresponding to the repeated pairs ${\lambda }_r={\lambda }_{r+1},r=1,3,5\dots 15$ are in accordance with the following pattern (the two upper rows in Figure 8): ${\mathbf{\Psi }}^{\left(1\right)}={}^{\left(0\right)}\mathbf{\Psi }$ and ${\mathbf{\Psi }}^{\left(2\right)}={}^{\left(7\right)}\mathbf{\Psi }$; ${\mathbf{\Psi }}^{\left(3\right)}={}^{\left(1\right)}\mathbf{\Psi }$ and ${\mathbf{\Psi }}^{\left(4\right)}={}^{\left(6\right)}\mathbf{\Psi }$; ${\mathbf{\Psi }}^{\left(5\right)}={}^{\left(2\right)}\mathbf{\Psi }$ and ${\mathbf{\Psi }}^{\left(6\right)}={}^{\left(5\right)}\mathbf{\Psi }$; ${\mathbf{\Psi }}^{\left(7\right)}={}^{\left(3\right)}\mathbf{\Psi }$ and ${\mathbf{\Psi }}^{\left(8\right)}={}^{\left(4\right)}\mathbf{\Psi }$. In the more general case of even larger uniform arrays, the pattern is similar, such that for the pairs ${\lambda }_r={\lambda }_{r+1},r=1,3,5\dots N-1$ the corresponding eigenvectors are ${\mathbf{\Psi }}^{\left(r\right)}={}^{\left(\right(r+1)/2-1)}\mathbf{\Psi }$ and ${\mathbf{\Psi }}^{(r+1)}={}^{(N-(r+1)/2)}\mathbf{\Psi }$. As an example, the $\gamma$-dependent evolution of the eigenvector ${}^{\left(15\right)}\mathbf{\Psi }$ position within the r-set is illustrated in Figure 8. Hereafter, the eigenvector distinguished by the the maximal number of the nodal points $m=N-1$ is also denoted as ${\mathbf{\Psi }}^{(\ast )}$), such that ${\mathbf{\Psi }}^{(\ast )}={}^{(N-1)}\mathbf{\Psi }$. It is marked by the red frame in Figure 8. At $\gamma =0$, ${}^{\left(15\right)}\mathbf{\Psi }={\mathbf{\Psi }}^{\left(2\right)}$. When $\gamma$ increases, the number r of the eigenvalue ${\lambda }_r$, that ${}^{\left(15\right)}\mathbf{\Psi }$ is associated with, increases as well. For example, for $\gamma =0.033\le {\gamma }_{\mathrm{crit}}$, ${}^{\left(15\right)}\mathbf{\Psi }={\mathbf{\Psi }}^{\left(6\right)}$. At $\gamma \ge {\gamma }_{\mathrm{crit}}$ (where ${\gamma }_{\mathrm{crit}}={\gamma }_{\times }^{\max }$ is the highest $\gamma$ at which crossing occurs), ${}^{\left(15\right)}\mathbf{\Psi }={\mathbf{\Psi }}^{\left(16\right)}$. In this “supercritical” case, the number of nodes m in the eigenvector ${\mathbf{\Psi }}^{\left(r\right)}$ associated with the eigenvalue ${\lambda }_r$ is independent of $\gamma$ and increases with the mode number r in accordance with the relation $m=r-1$. While the shapes of several lower eigenvectors are qualitatively similar to that of a uniform linear mass-spring chain with fixed ends, the highest eigenvectors are still different, apparently due to the influence of the “boundary conditions” of the array. This behavior is illustrated in the two bottom rows of Figure 8.

Finally, we present the results for the large 200-beam array. Recall that all the calculations are carried out for the array with dimensions reported in [17,21], Table 1. Figure 9 shows the full set of the eigenvalues ratios ${\lambda }_r/{\lambda }_1$ for few values of $\gamma$ below and above the critical value ${\gamma }_{\mathrm{crit}}=0.1415$ (also, here, the markers corresponding to the same $\gamma$ are connected by a line). Note that in accordance with our results, the critical ES coupling parameter for an infinite array ${\gamma }_{\mathrm{crit}}^{\inf }=4\eta /(1+10\eta )$ can serve as a limiting value of ${\gamma }_{\mathrm{crit}}$ for a (sufficiently large) finite array. One observes that for $0<\gamma <{\gamma }_{\mathrm{crit}}$, the specific mode number $r=r_{\mathrm{fold}}$ exists, such that there is a discontinuity (fold) in the slope of the curve. The eigenvector ${\mathbf{\Psi }}^{(\ast )}={}^{(N-1)}\mathbf{\Psi }$ corresponding to $r=r_{\mathrm{fold}}$ is the vector with the maximal number of the nodal points. For $r<r_{\mathrm{fold}}$, the eigenvectors resemble those of the fully coupled array (as in the first five images of the third row in Figure 8). On the other hand, for $r>r_{\mathrm{fold}}$, the eigenvectors are in a sense similar to the independent modes of each of the half-arrays. This behavior is attributed to the fact that the relative displacement between the adjacent elements of the arrays is more pronounced for higher modes, and therefore a higher voltage is required to achieve full coupling.

4. Parametric Resonances and Pattern Switching in an Array with Subcritical ES Coupling

To illustrate the influence of the ES coupling on the array dynamics, the response of a uniform array undergoing parametric excitation is analyzed in this section. The actuating voltage $V\left(t\right)=V_{\mathrm{dc}}+V_{\mathrm{ac}}cos\left(\omega t\right)$ applied between the two halves of the array contains both the steady (dc) and the time-harmonic (ac) components. The emphasis is on the weakly coupled case $V_{\mathrm{eff}}<V_{\mathrm{crit}}$ (where $V_{\mathrm{eff}}=\sqrt{V_{\mathrm{dc}}^2+V_{\mathrm{ac}}^2/2}$ is the effective dc voltage). The analysis for the fully coupled case $V_{\mathrm{dc}}>V_{\mathrm{crit}}$ has been presented in [21]. The eigenvectors obtained in Section 3 are used as the base vectors in the RO model construction.

4.1. Single Degree of Freedom Reduced Order Model

The vector $\mathbf{u}\left(t\right)=\left\{u_n\left(t\right)\right\},n=1\dots N$ of the normalized deflections at the beam’s endpoints is approximated using the expression

$$\mathbf{u}\left(t\right)\approx \mathbf{\Psi }\mathbf{v},$$

(17)

where $\mathbf{v}\left(t\right)=\left\{v_r\left(t\right)\right\},(r=1\dots N)$ is the vector of the generalized (modal) coordinates and $\mathbf{\Psi }=\left[{\mathbf{\Psi }}^{\left(1\right)},{\mathbf{\Psi }}^{\left(2\right)},\dots {\mathbf{\Psi }}^{\left(N\right)}\right]$ is the (unit mass-orthonormal) $N\times N$ modal matrix, whose columns ${\mathbf{\Psi }}^{\left(r\right)}=\left\{{\psi }_n^{\left(r\right)}\right\}$ are the eigenvectors of Equation (3). In contrast to the approach implemented in [21], where the voltage-independent base vectors were used, here, the base vectors ${\mathbf{\Psi }}^{\left(r\right)}$ are affected by the time-independent (dc) component of the voltage $V_{\mathrm{eff}}$. The vectors are found by replacing V by $V_{\mathrm{eff}}$ in Equation (3), and then numerically solving the resulting eigenvalue problem each time for the specific value of $V_{\mathrm{eff}}$. In the framework of the assumed modes procedure [35], Equation (17) is substituted into Equation (2) and multiplied by ${\mathbf{\Psi }}^T$ to yield

$$\ddot{\mathbf{v}}+\mathbf{C}\dot{\mathbf{v}}+\mathbf{K}\mathbf{v}-\mathbf{F}=0.$$

(18)

Here, $\mathbf{K}=\mathbf{I}\Lambda$ is the diagonal stiffness matrix (where $\Lambda =\left\{{\lambda }_r\right\}$ is the vector of eigenvalues) and $\mathbf{C}$ is the damping matrix. All the nonlinearities including the OS nonlinearity, as well as the mechanical and the electrostatic nonlinear IS coupling terms, are incorporated into the vector $\mathbf{F}=\mathbf{F}(\mathbf{v},V(t\left)\right)$.

We consider the vibrations at one particular rth mode of the array and chose the driving signal frequency $\omega$ to be close to twice the natural frequency ${\Omega }_r=\sqrt{{\lambda }_r}$ associated with the rth mode. In this case, the parametric excitation is predominantly due to the term $2V_{\mathrm{dc}}{V}_{\mathrm{ac}}cos\left(\omega t\right)$. As a result, and the equation of motion, Equation (18), takes the form

$$\ddot{v_r}+\frac{1}{Q}\dot{v_r}+\left[{\lambda }_r+P_{r}cos\left(\omega t\right)\right]v_r+\left[G_r+H_{r}cos\left(\omega t\right)\right]v_r^3=0,(r=1..N),$$

(19)

where

$$\begin{array}{ccc}\hfill {\lambda }_r& =& 1+2\eta (1-b_{rr})+2\beta (1-f_{rr})V_{\mathrm{eff}}^2,\hfill \\ \hfill P_r& =& 4\beta (1-f_{rr})V_{\mathrm{ac}}^2\sqrt{{\left(\frac{V_{\mathrm{eff}}}{V_{\mathrm{ac}}}\right)}^2-\frac{1}{2}},\hfill \\ \hfill G_r& =& \epsilon {\kappa }_{rr}^{nl}-{\eta }^{nl}{b}_{rr}^{nl}+\beta V_{\mathrm{eff}}^2{f}_{rr}^{nl},\hfill \\ \hfill H_r& =& \frac{P_r}{2}\left(\frac{f_{rr}^{nl}}{1-f_{rr}}\right).\hfill \end{array}$$

(20)

The summation coefficients appearing in Equation (20)

$$\begin{array}{ccc}\hfill b_{rr}& =& \sum _{n=1}^N{\psi }_n^{\left(r\right)}{\psi }_{n+2}^{\left(r\right)},\hfill \\ \hfill f_{rr}& =& \sum _{n=1}^N{\psi }_n^{\left(r\right)}{\psi }_{n+1}^{\left(r\right)},\hfill \\ \hfill {\kappa }_{rr}^{nl}& =& \sum _{n=1}^N{\left({\psi }_n^{\left(r\right)}\right)}^4,\hfill \\ \hfill b_{rr}^{nl}& =& \sum _{n=1}^N{\psi }_n^{\left(r\right)}\left[{\left({\psi }_{n+2}^{\left(r\right)}-{\psi }_n^{\left(r\right)}\right)}^3-{\left({\psi }_n^{\left(r\right)}-{\psi }_{n-2}^{\left(r\right)}\right)}^3\right],\hfill \\ \hfill f_{rr}^{nl}& =& \sum _{n=1}^N{\psi }_n^{\left(r\right)}\left[{\left({\psi }_{n+1}^{\left(r\right)}-{\psi }_n^{\left(r\right)}\right)}^3-{\left({\psi }_n^{\left(r\right)}-{\psi }_{n-1}^{\left(r\right)}\right)}^3\right]\hfill \end{array}$$

(21)

are calculated for each value of $V_{\mathrm{eff}}$ using the voltage-dependent base vectors ${\mathbf{\Psi }}^{\left(r\right)}$. Note that for the adopted parameters $|H_r/P_r|\ll 1$.

By re-scaling time $\tau =\omega t$, introducing the small parameter ${ϵ}_r=P_r/{\lambda }_r\ll 1$ and the new generalized coordinate $q_r=v_r\sqrt{{ϵ}_r}$, assuming that the linear damping term is of the order of $\mathcal{O}\left({ϵ}_r\right)$, and then neglecting terms of $\mathcal{O}\left({ϵ}_r^2\right)$, Equation (19) is reduced to the canonical nonlinear Mathieu’s (Mathieu–Duffing) equation [36]

$$\frac{d^2{q}_r}{d{\tau }^2}+{ϵ}_r{\mu }_r\frac{d{q}_r}{d\tau }+\left({\delta }_r+{ϵ}_{r}cos\tau \right)q_r+{ϵ}_r{\alpha }_r{q}_r^3=0,$$

(22)

where

$${\delta }_r=\frac{{\lambda }_r}{{\omega }^2},{ϵ}_r=\frac{P_r}{{\omega }^2},{\alpha }_r=\frac{G_r}{{\omega }^2},{\mu }_r=\frac{1}{{ϵ}_r\omega Q}=\frac{\omega }{P_{r}Q}.$$

(23)

We first consider the linearized counterpart of Equation (22) (by setting ${\alpha }_r=0$) and find the transition curves in terms of the physical parameters $\omega$, $V_{\mathrm{eff}}$ and $V_{\mathrm{ac}}$. While carrying out the mapping, the values of $V_{\mathrm{ac}}$ and $V_{\mathrm{eff}}$ are set, $V_{\mathrm{eff}}$ is held constant and the modulation of $P_r$ is achieved only by varying the ac voltage. The dc voltage $V_{\mathrm{dc}}$ is expressed in terms of $V_{\mathrm{ac}}$ and $V_{\mathrm{eff}}$ by requiring that $\sqrt{V_{\mathrm{dc}}^2+V_{\mathrm{ac}}^2/2}$ is equal to the adopted $V_{\mathrm{eff}}$ value. Note that since $V_{\mathrm{eff}}$ is held constant, as the ac voltage increases, the dc voltage decreases. Substituting Equation (23) into the expression for the transition curves [36]

$${\delta }_r=\frac{1}{4}\pm \frac{{ϵ}_r}{2}\sqrt{1-{\mu }_r^2}+\mathcal{O}\left({ϵ}_r^2\right)$$

(24)

yields the expression for the boundaries of the parametric tongues

$${\omega }_{1,2}^2=4{\lambda }_r-\frac{2}{Q^2}\pm 2\sqrt{P_r^2-\frac{4{\lambda }_r}{Q^2}+\frac{1}{Q^4}}$$

(25)

where $P_r=P_r\left(V_{\mathrm{ac}}\right)$ and ${\lambda }_r$ are defined in Equations (20) and (21).

Figure 10a shows the transition curves in terms of $\omega$ and $V_{\mathrm{ac}}$ for $V_{\mathrm{eff}}=40$ V. Each curve corresponds to the case when the arrays vibrate at one of the modes r varying between mode $\#70$ up to mode $\#80$. Figure 10b depicts the same curves for higher modes $\#129$ to $\#140$.

To excite PR in a damped system, the magnitude ${ϵ}_r$ of the time-harmonic term modulation in Equation (22) should be higher than a certain threshold value [36]. It is instructive to illustrate the role of the ES coupling on the actuating threshold. Substituting ${ϵ}_r,{\mu }_r$ from Equation (23) into the threshold condition ${\mu }_r=1$, using the notations Equation (20), while taking into account that at PR $\omega =2{\Omega }_r$ yields the threshold value of the ac voltage

$$V_{\mathrm{ac}}^{\mathrm{th}}=\sqrt{V_{\mathrm{eff}}^2-\sqrt{V_{\mathrm{eff}}^4-\frac{{\lambda }_r}{2Q^2{\beta }^2{(1-f_{rr})}^2}}}.$$

(26)

Here, in accordance with Equation (20), ${\lambda }_r$ depends only on $V_{\mathrm{eff}}$ and not on $V_{\mathrm{dc}}$ and $V_{\mathrm{ac}}$ separately. Figure 11 illustrates the influence of the mode number on $V_{\mathrm{ac}}^{\mathrm{th}}$ for two different values of $V_{\mathrm{eff}}$ below $V_{\mathrm{crit}}$, namely $V_{\mathrm{eff}}=30$ V, $V_{\mathrm{eff}}=40$ V, and one above it $V_{\mathrm{eff}}=80$ V (the fully coupled case). In general, an increase in $V_{\mathrm{eff}}$ results in a decrease in $V_{\mathrm{ac}}^{\mathrm{th}}$. In addition, $V_{\mathrm{ac}}^{\mathrm{th}}$ is generally higher for lower modes of the array. However, the behavior is quite different for the subcritical $V_{\mathrm{eff}}<V_{\mathrm{crit}}$ and supercritical fully coupled $V_{\mathrm{eff}}\ge V_{\mathrm{crit}}$ cases. In the fully coupled case $V_{\mathrm{ac}}^{\mathrm{th}}$ monotonically decreases with the mode number. In contrast, in the subcritical case, this monotonic dependence is preserved only for r below $r_{\mathrm{fold}}$, which is the mode number corresponding to the eigenvector ${\mathbf{\Psi }}^{(\ast )}={}^{(N-1)}\mathbf{\Psi }$ and to the fold in the frequency curve, Figure 9. At $r=r_{\mathrm{fold}}$ the pronounced decrease in $V_{\mathrm{ac}}^{\mathrm{th}}$ occurs, and for $r>r_{\mathrm{fold}}$, the dependence of $V_{\mathrm{ac}}^{\mathrm{th}}$ on r becomes to be irregular. In accordance with Figure 11, for certain r, excitation of PR in the subcritical case can be achieved at lower $V_{\mathrm{eff}}$ and $V_{\mathrm{ac}}$ voltages than in the fully coupled case. This behavior can be attributed to the fact that at $r>r_{\mathrm{fold}}$ the array vibrates at a predominantly staggered mode, when the beams of the two half-arrays vibrate at the anti-phase. In this case, each of the half-arrays is in a sense excited independently, when the beams of the opposite half-array serve as an electrode. It is remarkable that when the eigenvectors ${\mathbf{\Psi }}^{\left(r\right)}=sin[\pi nr/(N+1)]$ of a uniform “fully coupled” mass-spring chain are used as the base vectors in the assumed modes procedure, the decrease in the threshold voltage and the irregular dependence of $V_{\mathrm{ac}}^{\mathrm{th}}$ on the mode number $r>r_{\mathrm{fold}}$ are not captured by the (singe DOF) model.

Looking now for the solution of the nonlinear Mathieu equation, Equation (22), we consider the undamped case where ${\mu }_r=0$. The well-known asymptotic solution to this equation in the vicinity of the primary PR is obtained by setting ${\delta }_r=1/4+{ϵ}_r{\tilde{\delta }}_r$ (where ${\delta }_r,{ϵ}_r$ are defined in Equation (23) and ${\tilde{\delta }}_r$ is the detuning parameter), which yields [36]

$$q_r\left(t\right)=R_{r}cos\left(\frac{t}{2}-\theta \right),$$

(27)

$$R_r=\sqrt{-\frac{4}{3{\alpha }_r}\left(\frac{cos\left(2\theta \right)}{2}+{\tilde{\delta }}_r\right),}$$

(28)

where $R_r$ is the amplitude and ${\alpha }_r$ is defined in Equation (23). Here, $\theta$ takes the values of $0,\pi /2,\pi$ and $3\pi /2$, depending on the branch of the solution. Equation (28) describes curves on the ${\delta }_r$-${ϵ}_r$ plane corresponding to the nontrivial solutions of the undamped nonlinear Mathieu equation. Mapping of the stable and unstable branches (for $G_r<0$) based on Equation (28) allows us to build the frequency-response curves of the array vibrating at a specific mode in terms of the physical parameters

$$\begin{array}{ccc}\hfill v_r^U& =& \sqrt{-\frac{4{\Omega }_r+2P_r-{\omega }^2}{3G_r}},\left(\mathrm{unstable}\right),\hfill \end{array}$$

(29)

$$\begin{array}{ccc}\hfill v_r^S& =& \sqrt{-\frac{4{\Omega }_r-2P_r-{\omega }^2}{3G_r}},\left(\mathrm{stable}\right),\hfill \end{array}$$

(30)

where ${\Omega }_r=\sqrt{{\lambda }_r}$, $P_r$ and $G_r$ are given by Equation (23); $v_r^U$ and $v_r^S$ denote the amplitudes corresponding to the unstable and stable branches, respectively. In the case of $G_r>0$, the solution Equation (29) becomes stable, and the branch Equation (30) becomes unstable. In accordance with Equations (29) and (30), for each mode, the two critical values of the excitation signal frequency

$${\omega }_r^U=\sqrt{4{\Omega }_r+2P_r},{\omega }_r^S=\sqrt{4{\Omega }_r-2P_r},$$

(31)

correspond to the points on the $\omega$ axis, where $v_r^U\left({\omega }_r^U\right)=0$, $v_r^S\left({\omega }_r^S\right)=0$ and where the unstable and stable branches emerge from. The unique stable non-trivial solution exists within the interval ${\omega }_r^U<\omega <{\omega }_r^S$.

The character of the nonlinearity (softening vs. hardening)) of Equation (19) is dictated by the sign of $G_r$, defined in Equation (20) (e.g., see [37]), which, in turn, depends on the mode number r and the voltage $V_{\mathrm{eff}}$ (directly, Equation (20) and, due to the dependence of ${\mathbf{\Psi }}^{\left(r\right)}$, and therefore of the summation coefficients, Equation (21), on $V_{\mathrm{eff}}$). Figure 12 presents the value of $G_r$ as a function of $V_{\mathrm{eff}}$ and of the mode number. At lower modes $G_r>0$, and the system exhibits hardening nonlinearity. As the mode number increases, $G_r$ decreases untill its sign changes to negative and the system exhibits softening nonlinearity. The transition boundary subdividing the regions of positive and negative $G_r$ are shown by a dotted line in Figure 12a. Note that for low $V_{\mathrm{eff}}$ (when compared to $V_{\mathrm{crit}}$), there are isolated regions on the map $G_r(r,V_{\mathrm{eff}})$ where the system can also exhibit hardening nonlinearity at higher mode numbers. A discontinuity in Figure 12a corresponds to $r=r_{\mathrm{fold}}$, which is associated with the eigenvector ${\mathbf{\Psi }}^{\ast }={}^{(N-1)}\mathbf{\Psi }$. For a given $V_{\mathrm{eff}}<V_{\mathrm{crit}}$ and $r<r_{\mathrm{fold}}$ (left on the discontinuity line), as well as for $V_{\mathrm{eff}}\ge V_{\mathrm{crit}}$ (above the discontinuity line), the eigenvectors of the array are similar to those of a uniform mass-spring chain. In contrast, for $V_{\mathrm{eff}}<V_{\mathrm{crit}}$ and $r>r_{\mathrm{fold}}$, the device would behave, to a large extent, as two individual half-arrays. The transition boundary $r_{\mathrm{fold}}=r_{\mathrm{fold}}\left(V_{\mathrm{eff}}\right)$ is shown by the white dashed line in Figure 12a. The behavior of $G_r$ within the region of small mode numbers and small $V_{\mathrm{eff}}$ is detailed in Figure 12b.

Figure 13 illustrates the single mode undamped response of the array at few specific modes $r=40,50,60,70$ and 80. The left column Figure 13a,c presents the weakly coupled case where $V_{\mathrm{eff}}=30\mathrm{V}>V_{\mathrm{crit}}$, while the right column presents the case of the full coupling $V_{\mathrm{eff}}=80$ V. The differences in the (single mode) parametric resonant responses of the array subject to the strong $V_{\mathrm{eff}}\ge V_{\mathrm{crit}}$ and weak $V_{\mathrm{eff}}<V_{\mathrm{crit}}$ ES coupling manifest themselves in several aspects. As expected, for the same mode, higher $V_{\mathrm{eff}}$ results in the stiffening of the array and therefore in the shifting of the frequencies where the PR is exited toward higher values. For example, in the subcritical case and for $r=70$, the width of the PR region is such that ${\omega }_{70}^U=2.131$ and ${\omega }_{70}^S=2.185$ (Figure 13c, purple line). In the supercritical case, Figure 13d, the values are higher and the band is narrower, such that ${\omega }_{70}^U=2.41$ and ${\omega }_{70}^S=2.446$. An additional difference is the overlap between the modes, namely, the range of frequencies where PR can be excited. In the subcritical case, the PR frequency bands corresponding to the modes 60, 70, and 80 are almost fully overlapped, Figure 13a,c; no overlap is observed in the supercritical case, Figure 13b,d.

4.2. Numerical Analysis

Equation (2) completed by small non-zero initial conditions (of the order of ${10}^{-9}$, to make the excitation of the PR in the homogeneous equation possible) was solved numerically using the Verlet algorithm implemented in Matlab. The results are obtained for the mid-size array containing 16 beams. The model allowed us to capture the key features of the large arrays dynamic responses, while being small enough to conveniently represent and analyze the results. Two additional fixed beams were added at each end of the array, $n=-1,0,N+1,N+2$, to account for the “fixed” boundary conditions $u_{-1}=u_0=u_{N+1}=u_{N+2}=0$ (see Figure 3). The actuation voltage was $V\left(t\right)=V_{\mathrm{dc}}+V_{\mathrm{ac}}cos\left(\omega t\right)$. The calculations were performed for $Q=100$, for the cases of the full ES coupling $V_{\mathrm{eff}}=80$ V $>V_{\mathrm{crit}}$ and for $V_{\mathrm{eff}}=30$ V $<V_{\mathrm{crit}}$. In all the cases, the ac voltage was $V_{\mathrm{ac}}=6$ V. In the framework of the adopted excitation scenario, the up-sweep of the frequency was used within the frequency range covering the entire propagation band.

The amplitude–frequency response, in terms of the spatially averaged deflection amplitudes ${\sum }_1^N\sqrt{u_n^2}/N$ is shown in Figure 14a for $V_{\mathrm{eff}}=30$ V, the same response for $V_{\mathrm{eff}}=80$ V is depicted in Figure 14b. In both cases, excitation frequency $\omega$ increases starting from the value below the PR band of the entire array. The response follows the trivial stable branch until the frequency ${\omega }_r^U$, Equation (31), corresponding to the left transition curve of the lowest excited mode is reached. Since, in our case, the actuating voltage $V_{\mathrm{ac}}$ is lower than the parametric threshold $V_{\mathrm{ac}}^{\mathrm{th}}$ for the few lowest modes, these modes are not excited through the PR mechanism. Following the jump at $\omega ={\omega }_r^U$, the response follows the stable branch, Equation (30). The amplitude decreases either down the trivial stable response at ${\omega }_r^S$, or until the excitation frequency reaches the value corresponding to the unstable branch of the next mode in the sequence. The comparison of the responses shown in Figure 14a,b indicates that in the weakly coupled case, the full overlap between the modes is observed within the entire range of frequencies where the PR is excited, starting from the mode $r=6$. In contrast, in the fully coupled case, the lowest modes excited through the PR mechanism are not interacting (modes 5 to 8) and the overlap is observed starting from the mode $r=9$. In addition, consistently with the asymptotic analysis for the single-mode response, the range of the excitation frequencies, where the PR is excited, is wider in the fully coupled case. Moreover, since $V_{\mathrm{ac}}^{\mathrm{th}}$ is generally higher for smaller $V_{\mathrm{eff}}$, the first mode excited in the weakly coupled case is $r=6$ (which for this $V_{\mathrm{eff}}$ also corresponds to the $r_{\mathrm{fold}}$, as in Figure 8, the third row), as opposed to $r=5$ in the fully coupled case. The situation is further illustrated in Figure 14c,d where the standing wave pattern for the cases described above are shown. The separation between the first excited modes 5 to 8 in Figure 14d is clearly observed. In contrast, in Figure 14c the modes are overlapping within the entire frequency range where the PR is excited. Note that since the colors in Figure 14c,d represent the modal amplitudes, the negative deflections are not shown (e.g., see [17]). For this reason, in the fully coupled arrays, namely Figure 14d, the pattern of the modes $r=(N+1)/2+1/2+j$ and $r=(N+1)/2-1/2-j$ look qualitatively similar. In contrast, in the case of the weak coupling, Figure 14c, the first mode excited through the PR is $r=6=r_{\mathrm{fold}}$, which looks similar to the mode $r=1$. As a result, the “wavelength” of the modes always decreases with the increase in the frequency above the value associated with $r_{\mathrm{fold}}$. This behavior is consistent with the results shown in Figure 8.

5. Discussion

The primary goal of this work was to investigate the electrostatically coupled arrays dynamics, and, specifically, the influence of the electrostatic coupling on the eigenvalues and eigenvectors of these structures. We believe however that the results of the work could also be of practical significance. In contrast to the elastic coupling, which is dictated by the array geometry and cannot be changed, the ES coupling can be easily controlled by the voltage. An efficient tuning of the modal behavior of the system, of its natural frequencies and modes, the veering and the crossing points, may open up opportunities for the new operational scenarios of resonant micro- and nano-scale sensors. The sensitivity of resonant sensors is usually lower in the vicinity of veering [38]. An ability to shift the veering frequencies away from the working frequencies of the sensor or replacing veering by crossing may avoid this undesirable sensitivity degradation. The unique ability to tailor the eigenvectors shape, for example, by tuning the veering or crossing frequencies, may allow to change the position of the specific (nth) beam of the array within the vector from the node where the displacement is close to zero up to locations where the modal amplitude is maximal. This feature could be beneficial in so-called “mode localization” based sensors [5,39,40], where the ratio between the modal amplitudes, rather than the frequency shifts, are monitored.

While exploration of possible sensing scenarios, without even mentioning possible design and performance considerations, is out of the scope of the present work, it is instructive to illustrate the influence of the disturbance of the array parameters on its spectral contents and associated modal patterns. Figure 15a shows the evolution of the unperturbed uniform eight-beam array eigenvectors with the increase in the ES coupling parameter $\gamma$. The increase in the mass # 3 by 2% results in a clearly visible change in the modal patterns, Figure 15b. (Relatively large mass perturbation was chosen for the clarity of illustration.) The differences between the perturbed and the unperturbed patterns are shown in Figure 15c. The actual perturbed and unperturbed eigenvectors, obtained for a specific value of the ES coupling, are depicted in Figure 15d. Figure 15e and Figure 15f show the eigenvectors perturbation engendered by the 2% mass increase in the beams # 2 and # 4, respectively. The results suggest that the mass perturbation has a significant influence on the array modal response, especially in the case of higher eigenvectors 5, 6, 7 and 8. Moreover, the same mass added to different beams of the array may result in very different, clearly distinguishable, modal patterns. The change in the location (in terms of $\gamma$) of the veering and the crossing points is also very pronounced. In contrast, the largest relative change in the eigenvalues due to the 2% mass perturbation of different beams is no higher than 0.38%. The relative change in the sum of the eigenvectors is close to the relative added mass. This result is expected since the eigenvalues of a cantilever depend linearly on a small added mass. Note that extremely small, attogram [41] and even zeprogram, masses were shown to be detectable by nanoelectromechanical (NEMS) cantilever sensors. However, eigenvectors, and not only eigenvalues, monitoring in large arrays may open possibilities for identification of the perturbed beams of the array. This feature can be potentially used for multi-sensing without a need to individually address each of the beams. Note that this scenario can be realized in an actual experiment, for example, by identifying and monitoring the resonant frequency associated with the specific mode while sweeping the coupling voltage.

6. Conclusions

In this work, we analyze the eigenfrequencies, eigenvectors, and parametric resonant dynamics of an array of elastically and electrostatically coupled micro beams. One of the unique features of the considered architecture is that it is realized as two identical half-arrays. Within each half-array, the beams are coupled only elastically. Application of the voltage between the halves of the array results in an ES coupling. The main goal of this work was to explore the interplay between the elastic and the ES coupling and their influence on the arrays’ spectral content and the parametric resonant responses. The emphasis is on the cases in which the ES coupling between the two half-arrays is not dominant.

In the first part of the work, the influence of the ES coupling on the eigenfrequencies and the eigenmodes of the array is analyzed. Meanwhile, at zero voltage, the eigenvalues are obtained as repeated degenerate pairs (due to the unique mirror symmetry of the arrays), strengthening of the ES coupling introduces frequency veering and crossing. We show that variation in the ES coupling may allow us to change the frequencies corresponding to the veering and the crossing points, as well as the width (in terms of frequencies) of the propagation band. One of the main conclusions of this work is that the ES coupling affects not only the eigenvalues, but also the eigenvectors of the array. To illustrate the influence of the voltage on the modal shapes, we use two approaches to the eigenvectors ordering. When the eigenvectors are sorted based on the associated ascending eigenvalues, the eigenvectors’ order is unaffected by the voltage. However when the sorting is based on the number of eigenvectors’ nodal points (the points at which the imaginary line connecting the masses of the array crosses zero), the increasing ES coupling may alter, through the crossing or veering mechanism, the order of the eigenvectors within the set. In the considered array, the critical voltage $V_{\mathrm{crit}}$ can be found such that no frequency crossing or veering, and consequently no changes in the (nodal points-based) eigenvector order, takes place above this value. In the “fully coupled” case, the eigenvector ${\mathbf{\Psi }}^{\ast }$ with the largest number of the nodal points corresponds to the largest eigenvalue. The entire array, its eigenvectors, and its eigenvalues behave similarly to those of a uniform mass-spring chain. In contrast, for the voltages below $V_{\mathrm{crit}}$, the number $r=r_{\mathrm{fold}}$ of the eigenvector ${\mathbf{\Psi }}^{\ast }$ within the set depends on the voltage. In this case, the uniform mass-spring chain behavior is observed at the lower modes $r<r_{\mathrm{fold}}$. For $r>r_{\mathrm{fold}}$, the arrays behave, in a sense, as two independent half-arrays.

In the second part of the work, the influence of the weak ES coupling on the parametric resonant behavior of the array is illustrated using the RO model, and also numerically. We show that the weak and the strong couplings may result in qualitatively different behavior. For example, in the case of a full coupling, the threshold voltage $V_{\mathrm{ac}}^{\mathrm{th}}$ necessary to excite the PR in a damped system decreases monotonically with the mode number. In contrast, in the case of a weak coupling, $V_{\mathrm{ac}}^{\mathrm{th}}$ decreases abruptly at $r=r_{\mathrm{fold}}$. In the case of a higher coupling voltage, lower modes can be excited through the PR mechanism, and the entire band of frequencies at which the PR is excited is wider. We also show the importance of the use of not just admissible but actual eigenvectors of the associated linearized problem as the base vectors in the RO model.