1 Introduction

Micro/nano resonator has been widely implemented as inertial sensors for mass, force and acceleration sensing applications [1,2,3,4]. Coupled resonators have been suggested to improve the sensitivity and functionality of those sensors. There are several advantages of mechanically or electrically coupling them together. First, the mode localisation effect opens up new sensing methods, such as multi-mode sensing and the use of eigenvector (amplitude), to push down the limit by several orders of magnitude [5,6,7,8]. Second, it is possible to drive and readout from only one single resonator to find all eigenvalues (resonance frequencies) of a large array thanks to the collective behaviour [9]. This reduces the actuation complexity and turns the system into a single-input-single-output (SISO) multiplexed sensors array. In addition, coupling facilitates the characterization of large resonators array utilizing an inverse eigenvalue analysis (IEA) method. By perturbing the terminal element of a coupled array, two unique sets of eigenvalues are recorded to extract the system matrix using inverse problem theory. The system matrix is associated with important device properties such as Young’s modulus [10], mass sensitivity [11] and process variability [12]. However, existing IEA techniques require solving complicated eigenvectors or orthogonal polynomials, and hence become complicated when the array size is large.

Herein, we present a simplified IEA technique to circumvent this issue, thereby providing a new solution to system matrix determination. An array of five nearest-neighbour coupled square plate resonators were fabricated via the multi-user MEMS processes (MUMPs) to validate the ideas. Using the new IEA method, important sensor information including sensitivity and coupling ratio can be provided, which are of substantial interest for resonator designers.

The article is organised as follows: Section 2 provides the mathematical foundation of the simplified IEA technique to extract eigenvectors and system matrix, and analyses its complexity and accuracy. Section 3 describes the MEMS square plate resonator used in this research, including the fabrication process, electromechanical properties, experimental setup and result of a prototype sensor. Section 4 studies the collective behaviour of coupled systems using an example of five coupled resonators. Section 5 characterizes the sensor using perturbation analysis and the proposed IEA technique. System matrix and other parameters such as spring constant and coupling ratio are therefore determined. Section 6 gives the conclusion.

2 Inverse eigenvalue analysis

2.1 Coupled systems

An array of n nearest-neighbour coupled resonators can be modelled by the following equation

$$\begin{aligned} {\mathbf {M}} \ddot{{\mathbf {v}}} + {\mathbf {B}} \dot{{\mathbf {v}}} + \mathbf {Kv} = {\mathbf {F}} \end{aligned}$$
(1)

where \({\mathbf {M}},{\mathbf {B}}\) and \({\mathbf {K}}\) are mass, damping and stiffness matrix. \({\mathbf {v}}\) and \({\mathbf {F}}\) are two column vectors representing the displacement and driving force of each element. For low damping, the equation transforms into the classic eigenvalue problem

$$\begin{aligned} {\mathbf {S}} \mathbf {x}_{\mathbf {i}} = \lambda _i \mathbf {x}_{\mathbf {i}} \end{aligned}$$
(2)

where \({\mathbf {S}} = \mathbf {M^{-1/2}} {\mathbf {K}} \mathbf {M^{1/2}}\) is the symmetric tridiagonal (Jacobian) system matrix. At the ith mode, \({\mathbf {S}}\) has eigenvalues \(\lambda _i = \omega _i^2\) and eigenvectors \(\mathbf {x}_{\mathbf {i}}\) that are related to resonance frequencies and amplitudes. Ideally, there are n eigenvalues already available from the n resonance peaks, which can be observed from any element of the array due to collective behaviour. A mass or stiffness perturbation upon one resonator would lead to another set of n unique eigenvalues \(\lambda _i^*\).

2.2 New characterization algorithm

With the two sets of readily available eigenvalues before and after perturbation, the previous Choubey’s algorithm can be used to reconstruct the system matrix  [13]. This method, however, requires numerous multiplication and division operators to derive all eigenvectors. It is computationally expensive and impractical for real-word sensing when the array size n becomes large. The algorithm can be linearised by ignoring high order terms while preserving adequate precision. When a perturbation happens on the pth element, the new system matrix becomes

$$\begin{aligned} \mathbf {S}^*= {\mathbf {S}} + \varDelta {\mathbf {S}} \end{aligned}$$
(3)

where \(\varDelta {\mathbf {S}}\) is a matrix filled with zeros except for the (pp)th entry \(\varDelta S_{p,p}\), which can be calculated using the trace property of the system matrix before and after perturbation

$$\begin{aligned} \varDelta S_{p,p} = \sum _{i=1}^n{\lambda _i^*} - \sum _{i=1}^n{\lambda _i} = \sum _{i=1}^n{\varDelta \lambda _i} \end{aligned}$$
(4)

The eigenvalue problem has been updated to

$$\begin{aligned} ({\mathbf {S}} + \varDelta {\mathbf {S}})(\mathbf {x}_{\mathbf {i}} + \varDelta \mathbf {x}_{\mathbf {i}}) = (\lambda _i + \varDelta \lambda _i)(\mathbf {x}_{\mathbf {i}} + \varDelta \mathbf {x}_{\mathbf {i}}) \end{aligned}$$
(5)

Since \(\mathbf {x}_{\mathbf {i}}^{T} \mathbf {x}_{\mathbf {i}} = 1\) and \(\mathbf {x}_{\mathbf {i}}^{T} {\mathbf {S}} = \mathbf {x}_{\mathbf {i}}^{T} \lambda _i\), ignoring high order terms yields

$$\begin{aligned} \mathbf {x}_{\mathbf {i}}^{T} \varDelta {\mathbf {S}} \mathbf {x}_{\mathbf {i}} = \varDelta \lambda _i \end{aligned}$$
(6)

which leads to a simple solution

$$\begin{aligned} x_{p,i}^2 = \frac{\varDelta \lambda _i}{\varDelta S_{p,p}} \end{aligned}$$
(7)

This can be compared with the previous approach [13]

$$\begin{aligned} x_{p,i}^2 = \frac{\varDelta \lambda _i}{\varDelta S_{p,p}} \prod _{j =1,j \ne i}^n{\frac{\lambda _i - \lambda _j^*}{\lambda _i - \lambda _j}} \end{aligned}$$
(8)
Fig. 1
figure 1

Algorithmic complexity of the new and previous method using MATLAB computation time

2.3 Complexity and accuracy

Figure 1 shows the algorithmic complexity of the new method and the previous Choubey’s approach for different number of resonators using MATLAB computation time analysis. The new method has reduced the complexity from \(O(n^2)\) to O(n), which facilitates signal processing requirement and paves the way for rapid characterization of massive sensors array. With the amplitudes of the pth resonator known for all modes, the (pp)th element of the system matrix is determined by

$$\begin{aligned} S_{p,p} = \left( \begin{array}{ccc} x_{p,1}&\quad \cdots&\quad x_{p,n} \end{array} \right) \left( \begin{array}{ccc} \lambda _1 &{}\quad &{}\quad \\ &{}\quad \ddots &{}\quad \\ &{}\quad &{}\quad \lambda _n \\ \end{array} \right) \left( \begin{array}{c} x_{p,1} \\ \vdots \\ x_{p,n} \end{array} \right) \end{aligned}$$
(9)

When perturbation happens on the terminal element (\(p = 1\) or \(p = n\)), inverse Lanczos’ algorithm can be applied to extract all other entries in the system matrix iteratively [14, 15]. Hence all mass normalised stiffnesses and coupling strengths have been determined by simply recording the two sets of eigenvalues from any element.

To appreciate the accuracy of the new method for system matrix extraction, the mean absolute error of diagonal elements and off-diagonal elements are depicted in Fig. 2 using a simulated array of five coupled resonators. Perturbation applied to the terminal element, hence inverse Lanczos’ algorithm was used to determine the full system matrix. Various perturbation levels and coupling strengths were considered. Results suggest that for small perturbation (\(\varDelta S_{p,p}/S_{p,p} < k_c/k\)), the relative error can be well controlled below 1%. It is worth noting that off-diagonal elements, which correspond to coupling spring constants, are more prone to error, particularly for a weaker coupling. For weak coupling, these spring constants are very small compared with diagonal elements, hence more sensitive to numerical error. In practice, diagonal elements are of more importance as they contain mass and stiffness for individual resonators.

Fig. 2
figure 2

Accuracy of the new method for extracting a diagonal elements and b off-diagonal elements

3 Square plate resonator

3.1 Device architecture and fabrication process

Square plate resonator is a simple structure that utilises the transverse mode vibration of a centre-stemmed square membrane, and has already been applied in the fields of radio frequency and communication  [16]. The design benefits from high quality factor Q, low impedance and convenient eigenvalue readout channel, which can be potentially used as sensors or transducers.

The device was designed and fabricated through the PolyMUMPs programme provided by the MEMSCAP. The fabrication typically involves a three-layer polysilicon surface micromachining process [17]. Two layers of polysilicon, i.e. Poly 1 and Poly 2, were used in our design, as shown in Fig. 3(a). These polysilicon layers are conductive with 10–20 sheet resistance. Poly 1 shaped the driving and readout electrodes as well as interconnections at the bottom. Poly 2 acted as the centre-stemmed square plate membrane. The thickness t of the square plate (\(t = 1.5\,\upmu {\hbox {m}}\)) and the gap distance g between the top plate and bottom electrode (\(g = 0.75\,\upmu {\hbox {m}}\)) were defined by the process. As Poly 0 layer was not used in our design, the device was anchored to the Nitride substrate by filling an anchor hole with Poly 1. The hole was formed by patterning silicon dioxide with a mask layer and reactive ion etching (RIE).

Fig. 3
figure 3

a Cross-section view of the PolyMUMPs layers used for device fabrication and b electrical connection for one resonator

3.2 Electromechanical properties

Fig. 4
figure 4

FEM simulation of static displacement under different DC voltages

The natural frequency \(f_0\) of such square plate membrane with Young’s modulus E and density \(\rho\) is given by [18]

$$\begin{aligned} f_0 = 0.9697 \sqrt{\frac{E}{\rho }} \frac{t}{L^2} \approx \sqrt{\frac{E}{\rho }} \frac{t}{L^2} \end{aligned}$$
(10)

The side length L is designed to be 100 \(\upmu {\hbox {m}}\) to provide around 1.2 MHz resonance. The mass and stiffness of the resonator are calculated by

$$\begin{aligned} m = \rho t L^2 \end{aligned}$$
(11)

and

$$\begin{aligned} k = E \left( \frac{2 \pi }{L} \right) ^2 t^3 \end{aligned}$$
(12)

Hence the designed mass and stiffness for each resonator are 34 ng and 2020 N/m. Figure 3(b) depicts the electrical circuit for the sensor. The device implements a direct drive and readout scheme. A DC bias \(V_{DC}\) superimposed with a small AC signal \(v_{in}\) is applied to electrode 1 and 3. Electrode 2, 4 and the top plate are grounded. Change in resonance magnitude results in an AC current \(i_{in}\) which is linearly dependent on the vibration velocity  [19]. The conductance is calculated by  [16]

$$\begin{aligned} |Y| = \left| \frac{i_{in}}{v_{in}} \right| = \frac{Q \varepsilon ^2 A_{eff}^2 V_{DC}^2}{\sqrt{km} g^4} \end{aligned}$$
(13)

with \(\varepsilon\) being the permittivity and \(A_{eff}\) the effective area where top plate and driving electrodes overlap. The typical conductance of each resonator, assuming \(Q = 10^4\) and \(V_{DC} = 10\,{\hbox {V}}\), is estimated to be 14.5 \(\upmu\)S.

A large DC bias voltage is preferred for strong signal amplitude; however, may cause non-linear and pull-in effect that can potentially destroy the device. Figure 4 shows the static displacement of the membrane under different bias voltages using COMSOL finite element method (FEM) simulation. The pull-in voltage is projected to be 50 V. The spring constant of the resonator can be tuned by tweaking the DC voltage, hence introducing a negative electrostatic stiffness \(k_E\) [13]

$$\begin{aligned} k_E = - \frac{\varepsilon A_{eff}}{g^3} V_{DC}^2 \end{aligned}$$
(14)

that reduces the resonance frequency. Such spring softening effect provides a flexible way of perturbing the resonator as required by the IEA.

3.3 Single resonator

Fig. 5
figure 5

SEM image of a single square plate resonator and its measurement setup

The SEM image of a single resonator is illustrated in Fig. 5. An impedance analyser with frequency sweep functionality (Hioki IM3570) supplies the required test voltages and measures the impedance or conductance with below 0.08% error. Each frequency sweep takes less than 1ms by recording up to 800 frequency points.

Fig. 6
figure 6

Conductance and phase response of a single resonator

Figure 6 shows the conductance and phase response of a prototype sensor by applying \(V_{DC}=10\) V and \(v_{in} = 50\) mV under a low pressure level \(P < 100\) mPa. Fabricated resonators demonstrate a natural frequency between 0.8 and 0.85 MHz, which is approximately 30% less than the nominal frequency 1.2 MHz due to fabrication tolerance. This was likely to be caused by the extra side length when depositing Poly 2 layer onto the uneven silicon dioxide formed during the patterning process, thereby increasing the effective length by approximately 15%.

Damping effect has been characterized for a single resonator. The quality factor of the mechanical resonator can be estimated using the frequency to bandwidth ratio measured from the frequency spectrum

$$\begin{aligned} Q = \frac{f_0}{f_H - f_L} \end{aligned}$$
(15)

where \(f_H\) and \(f_L\) are the high and low cut-off frequencies where the power is half the maximum value. The relation between the calculated quality factor Q and ambient pressure P is plotted in Fig. 7(a). The quality factor varies between 10 and \(1.6 \times 10^4\) when the pressure drops from 3.6 kPa to 27 mPa. Figure 7(b) shows the square root of the electrostatic stiffness as a linear function of \(V_{DC}\). The resonator collapses at a pull-in voltage about 25 V, which causes short circuit between the top plate and bottom electrodes, therefore burning the device. To prevent such irreversible process from happening and ensure good linearity, the recommended DC operating range is within 10 V, which allows a maximum 2.1% perturbation. The lower-than-expected pull-in voltage is attributed to the same over-sizing effect of the side length.

Fig. 7
figure 7

Experimental characterization of a single resonator: a quality factor under different pressures and b square root of electrostatic stiffness as a linear function of DC bias voltage

4 Coupled resonators array

Weakly coupled MEMS resonators demonstrate collective resonance behaviour at a unique set of eigen-modes that are different from a single resonator, and have been extensively explored as sensors with enhanced sensitivity, reduced connections and bandwidth, as well as the capability of multi-sensing. Coupling is first seen in micro cantilevers array, where the natural overhang structure introduces parasitic crosstalk between neighbouring elements.

To understand the coupling mechanism of the square plate resonators, Figure 8 presents the geometry of two adjacent plates. The two neighbouring membranes are coupled by arranging them in the same diagonal orientation and connecting two opposing corners using a mechanical beam that has a width \(W_C\) and a length \(L_C\). The coupling beam is usually small hence providing a weak coupling ratio.

Fig. 8
figure 8

Geometry of coupled square plate resonators

4.1 Eigen-mode analysis

Figure 9 reveals the eigen-modes of five corner coupled square plate resonators using COMSOL eigenfrequency analysis, along with the lumped-parameter spring mass model. Since all resonators are assumed to be identical, particular eigen-modes become too weak to be distinguished (e.g. the third resonator at the second and fourth mode) due to Anderson localisation [5]. In the presence of noise, a decrease in resonance amplitude would deteriorate the signal-to-noise ratio (SNR) and add complexity to subsequent signal processing. Hence for the conventional IEA method that applies a single actuation scheme by driving and reading out from one resonator, identifying all the eigenvalues of a large array can be challenging. However, one can drive and readout multiple resonators simultaneously to improve the power handling capability and increase the SNR.

Fig. 9
figure 9

Eigen-modes of five corner coupled square plate resonators and the spring mass model

4.2 Array configuration

The microscopic image of the designed five coupled array is shown in Fig. 10(a). A 20 \(\upmu {\hbox {m}} \times 10\,\upmu {\hbox {m}}\) mechanical beam with the same Poly 2 material was used to couple the corners of two neighbouring square plates. Table 1 summarises the dimensions of the prototype MEMS resonators array. The interconnection wires have a designed width of 60 \(\upmu {\hbox {m}}\) and the entire sensor takes up a size of 920 \(\upmu {\hbox {m}} \times 330\,\upmu {\hbox {m}}\).

Table 1 Dimensions for the proposed corner coupled square plate MEMS resonators array

To test the device, the ground electrodes of all five resonators were connected together. The first four elements (resonators 1–4) also shared the same DC supply, while the last element (resonator 5) was connected to an independent bias voltage \(V_P\) to enable terminal perturbation. \(v_{in}\) was added to both \(V_{DC}\) and \(V_P\) and the current summation from all elements was recorded. Figure 10(b) shows the five eigen-modes of frequency response. These eigenvalues can be used for subsequent perturbation analysis based on IEA, which will be discussed in the next section.

Fig. 10
figure 10

a Microscopic image of the designed five coupled array and b its frequency response

5 Sensor characterization

5.1 Perturbation analysis

Having determined all eigenvalues, perturbation can be performed by increasing the bias voltage of the terminal element \(V_P\), which would introduce a negative electrostatic stiffness to such element. However, all the eigenvalues of the system would decrease due to collective coupling behaviour, and satisfy the following interweaving condition  [10]:

$$\begin{aligned} \lambda _1^*< \lambda _1< \lambda _2^*< \lambda _2< \lambda _3^*< \lambda _3< \lambda _4^*< \lambda _4< \lambda _5^*< \lambda _5 \end{aligned}$$
(16)

To validate this experimentally, \(V_P\) was changed from 0 to 10 V, while \(V_{DC}\) remained constant at 10 V. The frequency responses before and after perturbation are presented in Fig. 11. The eigen-frequencies shifted from 0.785, 0.822, 0.829, 0.837 and 0.847 MHz to 0.783, 0.819, 0.829, 0.835 and 0.845 MHz. \(\lambda\) and \(\lambda ^*\) can therefore be calculated to extract the system matrix.

Fig. 11
figure 11

Frequency response of five coupled array before and after terminal element perturbation

5.2 System matrix determination

Using the new IEA algorithm from Sect. 2, the system matrix (before perturbation) is determined to be

$$\begin{aligned} {\mathbf {S}} = \left( \begin{array}{ccccc} 2.68 &{}\quad -\,0.13 &{}\quad 0 &{}\quad 0 &{}\quad 0 \\ -0.13 &{}\quad 2.57 &{}\quad -\,0.13 &{}\quad 0 &{}\quad 0 \\ 0 &{}\quad -\,0.13 &{}\quad 2.68 &{}\quad -\,0.05 &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad -\,0.05 &{}\quad 2.76 &{}\quad -\,0.02 \\ 0 &{}\quad 0 &{}\quad 0 &{}\quad -\,0.02 &{}\quad 2.72 \\ \end{array} \right) \times 10^{13} \end{aligned}$$
(17)

with the perturbation amount being \(\varDelta S_{5,5} = -\,6 \times 10^{11}\), or 2.2% in terms of percentage. The result can be compared with the system matrix extracted from the previous Choubey’s deterministic approach

$$\begin{aligned} {\mathbf {S}} = \left( \begin{array}{ccccc} 2.71 &{}\quad -\,0.13 &{}\quad 0 &{}\quad 0 &{}\quad 0 \\ -0.13 &{}\quad 2.57 &{}\quad -\,0.13 &{}\quad 0 &{}\quad 0 \\ 0 &{}\quad -\,0.13 &{}\quad 2.67 &{}\quad -\,0.05 &{}\quad 0 \\ 0 &{}\quad 0 &{}\quad -\,0.05 &{}\quad 2.75 &{}\quad -\,0.02 \\ 0 &{}\quad 0 &{}\quad 0 &{}\quad -\,0.02 &{}\quad 2.72 \\ \end{array} \right) \times 10^{13} \end{aligned}$$
(18)

The average discrepancy for diagonal and off-diagonal elements are 0.3 and 1.5% respectively, matching well with the prediction in Fig. 2.

5.3 Spring constant and coupling ratio

Sensor parameters such as mass, stiffness and coupling ratio can be derived from the extracted system matrix. The mass normalised spring constants of resonator 1–5 are \((2.55, 2.31, 2.51, 2.69, 2.70) \times 10^{13}\), which correspond to natural frequencies 0.804, 0.765, 0.797, 0.826 and 0.827 MHz. Hence the actual frequency variation of the PolyMUMPs fabrication is determined to be 3.2%.

The system matrix suggests a weak coupling varying from 0.8 to 4.9%, as these coupling springs are very small hence more prone to process error. Since the sensitivity of coupled resonators is highly dependent on the coupling ratio [20], it is important to optimise the design before using them as sensors. To understand how the geometry of these mechanical beams affects coupling ratio, resonators with different coupling beam lengths and widths can be characterized with the IEA method. This can be done using a simple array of two coupled resonators in COMSOL. The coupling ratios for different beam geometries are derived and presented in Fig. 12 using the same perturbation analysis and IEA method. It can be seen that coupling strength is three orders of magnitude more sensitive to \(W_C\) than \(L_C\). By changing the width alone, the coupling ratio is tunable from 0.003 to 20%.

Fig. 12
figure 12

Coupling ratios for various beam geometries

6 Conclusion

We introduce a new inverse eigenvalue technique to investigate coupled micro/nano systems. With reduced algorithmic complexity and less than 1% system matrix error, the method has been experimentally verified using an example of five coupled square plate MEMS resonators with eigenfrequencies around 0.85 MHz. By obtaining two sets of eigenvalues before and after terminal perturbation, the IEA has been performed to derive the system matrix and hence important sensor and process information. The method offers improved accuracy and simplified readout, which can be used to actuate and characterize large array of coupled resonators, including ultrasound transducers and multi-function inertial sensors. As a proof of concept, only five resonators are examined in this paper. Theoretically, the maximum number of resonators in a coupled array is limited by the mode liaison effect, i.e. \(n < \kappa Q \approx 100\)  [20]. Hence we have fabricated larger arrays of twenty to fifty coupled resonators, which will be tested and used as sensors in our future work.