A multi-sensing scheme based on nonlinear coupled micromachined resonators

A new multi-sensing scheme via nonlinear weakly coupled resonators is introduced in this paper, which can simultaneously detect two different physical stimuli by monitoring the dynamic response around the first two lowest modes. The system consists of a mechanically coupled bridge resonator and cantilever resonator. The eigenvalue problem is solved to identify the right geometry for the resonators to optimize their resonance frequencies based on mode localization in order to provide outstanding sensitivity. A nonlinear equivalent model is developed using the Euler–Bernoulli beam theory while accounting for the geometric and electrostatic nonlinearities. The sensor's dynamics are explored using a reduced-order model based on two-mode Galerkin discretization, which reveals the richness of the response. To demonstrate the proposed sensing scheme, the dynamic response of the weakly coupled resonator is investigated by tuning the stiffness and mass of the bridge and cantilever resonators, respectively. With its simple and scalable design, the proposed system shows great potential for intelligent multi-sensing detection in many applications.

The resonant MEMS sensors have gained massive success thanks to their ultra-small physical scale and ultra-high sensitivity. Resonant MEMS sensors rely on the ambient environment's influence on the resonators' resonance frequencies in different ways.
The parameters related to the resonance frequencies include the resonator's mass, geometric shape, axial stress level, and stiffness. Besides, the corresponding equipment aiming at measuring resonant frequency shift has proven to reach a high resolution, which further enhances the MEMS sensors' sensitivity [14,15]. In addition, the structure's nonlinearity features appear in the frequency response [16,17], especially jump or bifurcation, further enhancing the detection process [18].
Research aiming to improve the MEMS sensors sensitivity has been conducted extensively in the past few years. Most works assessed the approach to miniaturize further the structure size from the microscale to the nanoscale [19][20][21]. A grand scheme of these sensors is the MEMS-based sensor that employs the phenomenon of mode localization. The latter introduces linear modal interaction behaviours, especially mode crossing and mode veering [22], which show great potential for sensitivity improvement. These modal interactions typically happen as the two modal frequencies of the system approach each other. Some research focuses on utilizing the mode coupling phenomenon to design high-performance sensors through mechanically and/or electrostatically coupled resonators [23,24]. Zhang et al. [25] reported an acceleration sensing scheme based on electrostatically weakly coupled resonators and reached high amplitude ratio, which is 302 times higher than ratio used in resonance frequency theory. Other types of displacement sensors [26], pressure sensors [27], and ultrasmall mass sensors [28] are also credible examples.
Compared to sensors based on linear behaviours, other designs closely related to nonlinear phenomena received significant attention, such as utilizing bistability [29] and detecting bifurcation points [30]. Although a sizeable nonlinear behaviour usually accompanies instability, the benefit of sensitivity and quality overwhelms the disadvantage. Electrostatically actuated MEMS resonators have been designed to detect forces [31], masses [32], and accelerations [30,33]. The issue of stabilizing the response in the nonlinear regime has been widely studied. Kacem et al. [34] researched a dynamic stabilization technique for electrostatically actuated nanoresonators theoretically and experimentally. Their theory focused on actuating both primary and superharmonic resonances simultaneously. Besides, using the internal resonance induced between two bending modes is proven to be a possible solution for stabilizing mechanical oscillators [35]. Also, bifurcation phenomena have been utilized to enhance the sensitivity of MEMS sensors. Kumar et al. [36] operated a piezoelectric cantilever beam near the saddle-node bifurcation at its first natural frequency and detected the sudden jump in amplitude with the mass perturbation on resonator. In addition, the nonlinear designed MEMS resonator shows its potential in multi-sensing applications.
In the past few years, the topic of multiple sensing has received attention due to the increasing need for multiple parameters monitoring applications, especially in medical treatment, aerospace, and industry 4.0 [12,13,17]. The conventional solutions focus on the combination of different kinds of sensors [16]. For instance, Clifton et al. [37] developed multifunctional integrated sensors (MFISES) chip, combining numerous MEMS sensors into an integrated chip and benefits from the ultra-small size and high sensitivity of MEMS devices. Their MFISES device integrates 8 MEMS sensors on a single 2 mm 9 2 mm die which obtains high sensing performance and low power consumption. However, the direct employment of multiple devices inevitably increases the system complexity and production cost. Some research studies investigated the realization the multiple sensing on a single device, e.g., Jaber et al. [13] demonstrated a resonant gas sensor actuated simultaneously near the first and second vibration modes, where monitoring the frequency shifts of these two modes due to physical stimuli, the detection of both environmental temperature and gas concentration is achieved. This paper proposes a new sensing structure using the nonlinear response of electrostatically actuated coupled resonators. The sensor mainly consists of a mechanically coupled cantilever and bridge resonator. Different sensing targets are set for two resonators for multi-sensing purposes: the micro-gravimetric sensing for the cantilever and the stiffness sensing for the bridge. To investigate the rich nonlinear dynamic behaviour of the structure, a nonlinear theoretical model is developed using the Euler-Bernoulli beam theory while accounting for the geometric (i.e., dominated by cubic nonlinearity originated from the midplane stretching effect) and electric (i.e., dominated by quadratic nonlinearity) nonlinearities, respectively. The reduced-order model based on the two-mode Galerkin discretization and shooting technique are used to analyse the system. The influence of the resonators' geometric size, damping, actuation voltage, and actuation scheme is detailed and discussed theoretically. The sensor's multi-sensing ability is demonstrated by tracking the bifurcation frequency shift and peak frequency shift in the modes' shapes induced by mass and stiffness perturbations. The numerical results prove that this new kind of sensing structure could monitor concurrently two different kinds of perturbations with high efficiency through the nonlinear behaviour of the structure.
The present paper's structure is organized as follows. The geometric size of the multi-sensing scheme and the related problem formulation are introduced in Sect. 2. The theoretical results of the system's nonlinear behaviour, the sensing performance, and the parameter analysis are included in Sect. 3. In the end, the main conclusions of the research are summarized in Sect. 4.

Structure description and model
As displayed in Fig. 1, the multi-sensing scheme comprises a weakly coupled resonator, including cantilever and bridge resonators, mechanically coupled by a thin beam. The coupling strength could be controlled by changing the coupling position, the moment of inertia, and the coupling beam's length [38]. The detailed geometric parameters and physical properties of the coupled resonator system, considered to be fabricated from silicon, are listed in Table 1. Both cantilever and bridge resonators are driven electrostatically by a DC polarization voltage V DC , and the AC harmonic actuation V AC is given on one (or both) of them depending on the actuation theory. The system geometry is optimized to ensure that the first and second modes are similar.
The two resonators would be used to operate different functions. The cantilever resonator would detect the mass perturbation while the bridge experiences the stiffness perturbation. The red component on the tip of the cantilever, as shown in Fig. 1, represents the mass perturbation set on the cantilever, which could be generated from the special coating in the practical design. The stiffness of the bridge resonator is another sensing target, which links to many factors like the thickness, stress and stretch, and the environmental temperature.
The sensor structure ( Fig. 1) is composed of two microbeams (a bridge resonator and a cantilever resonator) with lengths L 1 and L 2 (L 1 [L 2 ), respectively, width b and thickness h. The coupling beam is functionalized to provide weak coupling, and its coupling strength significantly affects the sensor's sensitivity. Hence, quantifying and controlling the coupling beam's dimensions is crucial. Most of the sensors use an overhang to connect microbeams [30]. However, this coupling approach is challenging to model and practically control its coupling length [39]. Hence, the coupling beam with length L c , width b c , and same thickness with resonators is chosen and located at distance X c . Under this condition, the structure of sensors could be modelled as two Euler-Bernoulli beams coupled with a rotational spring k r [40]. The torsional stiffness could be represented as: where G denotes the material's shear modulus (69.3 Gpa for silicon), and b is a coefficient depending on the coupling beam's width b c and thickness h. By using Hamilton's principle and the equation of Euler- Bernoulli beams with distributed elements [41], the equations of motion governing the transverse deflectionsw 1 (for the bridge resonator) andw 2 (for the cantilever resonator) are: EIw 0000 In Eqs. (2) and (3), the primes and dots denote the partial differentiation of transverse deflections e w i ði ¼ 1; 2Þ with respect to the beam position e x and the time e t, respectively; A and I are the area and the moment of inertia of the rectangular cross section; d is the Dirac delta function; and e c 1 and e c 2 denote the viscous damping. The parameters and corresponding values are defined in Table 2. Note that the AC voltage is imposed on only one resonator excitation at a time (i.e., if V AC1 = 0 means that V Ac2 = 0 V and vice versa) to investigate the corresponding dynamic response of the proposed structure, while the DC polarization actuation is imposed to both resonators to enhance the associated quadratic nonlinearity.
After non-dimensionalizing and discretizing Eqs. (2) and (3) using the Galerkin procedure, the governing non-dimensional equations are written as follows: where the dk and dm represent the linear stiffness variation on bridge resonator due to the external stimulus and the mass perturbation on cantilever resonator, respectively. The detailed derivation process and definition of parameters appearing in Eqs. (4) and (5) are given in Appendix. Through Eqs. (4) and (5), the Jacobian matrix of the system is: where Solving the eigenvalue problem through the Jacobian matrix gives the two lowest global resonant frequencies of the coupled system [42]: To verify the theoretical model and obtain the vibration mode shapes of the coupled resonators, a 3D multi-physics finite-element model in the commercial software COMSOL is used [43]. Figure 2 depicts the variation of the first two global natural frequencies with respect to the length of the cantilever by solving Eq. (6), while applying 0 V DC load to the system, as well as FEM simulations from COMSOL. The results show good matching between both methods validating the adopted analytical modal. Figure 2b, c shows the associated first two lowest mode shapes of the structure. Figure 2a demonstrates that when the length of the cantilever is L 2 ¼ 277lm, the two curves are close to each other, and the two lowest natural frequencies are around 54 kHz. For the rest of the analysis, the cantilever length will be fixed at We also analyse the effect of DC load on the lowest two natural frequencies by applying the same DC load, V DC , to both resonators. Figure 3 depicts the variation of the first two global natural frequencies with respect to the DC load by solving Eq. (6). It shows that both the two lowest modes' resonance frequency would decrease with increasing DC polarization actuation due to the dominance of the quadratic nonlinearity (i.e., originated from the electrostatic force). Figure 4 reports a parametric study of the preceding eigenvalue problem. The different colours of the figure represent the effect of stiffness on the bridge resonator, while the x-axis represents the effect of mass variation on the cantilever resonator. The light and dark colours denote the 1 st and the 2 nd lowest global natural frequencies, respectively. The results prove that both parameters could influence the resonance frequencies.

Results and discussion
This section discusses the nonlinear behaviour of the weakly coupled system. Continuous analyses using the shooting technique [38] are conducted to verify the effect of the different physical parameters. The referred parameters are given in Table 1. Figure 5 shows the bifurcation points of the frequency response of global mode shapes as actuating the bridge resonator with V DC1 ¼ 40V and V AC1 ¼ 7:3V, and the cantilever resonator with V DC2 ¼ 10V. The solid blue and green dotted lines represent stable and unstable branches, respectively. Three characteristic points (notes as red in Fig. 5) should be highlighted, the peak point in the W2 appears at 53.835 kHz, and two corresponding bifurcation points appear at Fig. 3 Variation of the two lowest natural frequencies of the coupled system with respect to the DC load. The two resonators are subjected to the same V DC Fig. 4 The two lowest natural frequencies versus mass perturbation on the cantilever and stiffness perturbation on the bridge. Light-coloured lines and dark-coloured lines represent 1st and the 2nd natural frequencies, respectively 69.966 kHz and 57.327 kHz. The hardening behaviour in W1 generates a saddle-node bifurcation point at 69.966 kHz, and the peak in the W2 at 53.835 kHz is suitable for sensing. Near these two points, the tiny variation of perturbations could lead to an obvious change of transverse deflections, which show great potential to provide high sensitivity. The following parts will provide detailed discussions of the effects of different parameters on these two phenomena.

Continuous analysis for the characteristic points
The phase portraits and Poincaré sections for all three characteristic points are shown in Fig. 6. The red curves represent the peak point of W2 at 53.835 kHz, the blue curves represent the higher amplitude bifurcation point at 69.966 kHz, and the green curves represent the lower amplitude bifurcation point at 57.327 kHz, respectively. The phase portraits, which are generated from 200 cycles of steady-state time response, demonstrate elliptical orbits, and all  Poincaré sections converge to a single point, proving that the response is periodic motion of period-1. It could be concluded that before the peak point of W2 at 53.835 kHz, both W1 and W2 phase portraits enlarge with increased frequency. Between 53.835 kHz (W2 peak) and 69.966 kHz (higher amplitude bifurcation point), W1 keeps enlarging while W2 starts to shrink. In the last stable branch larger than 57.327 kHz, the W1 and W2 phase portraits shrink with increased frequency. It should be pointed that the nonlinear behaviour determines the maximum transverse deflection of W1 while the maximum transverse deflection of W2 is further related to its natural frequency.

AC actuation level on bridge
In Fig. 7, the DC polarization actuation is set as V DC1 ¼ 40V and V DC2 ¼ 10V on bridge and cantilever, respectively, and only AC actuation of different amplitudes V AC1 is provided on the bridge resonator. Figure 7 indicates the complete response around the two lowest modes under different levels of AC actuation. (The amplitudes W1 and W2 are nondimensional.) In Fig. 7a, b, the five curves are given the actuation V AC1 of 1.0 V, 1.5 V, 2.0 V, 2.5 V, and 3.0 V. It can be noted that under low AC actuation levels, specifically from 1:0V (purple) to 2:0V (green), the frequency response is linear. When V AC1 increases to 2:5V (yellow), hardening behaviour first appears, where the dotted line denotes the unstable branch. Additionally, the bifurcation jump frequency in W1 and the peak frequency in W2 obtain a similar value and lead to a potential 1:1 interaction (yellow curve). Three specific peaks are observed in the response: (i) the main peak appears near 52.8 kHz. It shifts to the right when the actuation V AC1 increases and finally turns to the bifurcation jump. (ii) The small peak at 53.3 kHz, linked to the linear mode localization, remains at the same frequency. (iii) The small peak at 26.4 kHz is due to the order two superharmonic behaviour.
When the actuation is increased from 4 to 7V, as shown in Fig. 7c, d, the bifurcation jump frequency in W1 continues to increase, and the peak frequency in W2 remains constant. When V AC1 increases further from 8 to 10V, as Fig. 7e, f show, the 1:1 interaction disappears, and the response leads to a new unstable softening branch with maximum value of 0.56. To obtain the best performance in sensing through bifurcation jump and eliminate the risk of the unwanted unstable branch, the AC actuation of V AC1 ¼ 7:3V is a convenient choice.

AC actuation level on cantilever
In this section, the system dynamics are simulated for different actuating schemes. More specifically, the DC actuation remains the same while the AC actuation is switched on the cantilever resonator. Compared to the response in Sect. 3.1.2, the linear behaviours of Figs. 7a, b and 8a, b in the two parts are quite similar: both actuating modes return two peaks, including the central peak at 53.3 kHz and the superharmonic peak at 26.4 kHz. Under high AC actuation, the frequency response curve around the second mode (Fig. 8c, d) shows a nonlinear softening behaviour as being dominated by quadratic nonlinearity. Besides, the level of nonlinearity is also lower than the bridge actuating result, which is insufficient to generate a noticeable bifurcation jump. By comparing the maximum amplitude of two frequency responses under high AC actuation (Figs. 7c, d and 8c, d), it could be found that the peak amplitude of W1 under bridge actuation (shown as red in Fig. 7c) is near 0.37, which is more than two times of the W1 (0.15) under cantilever actuation (shown as red in Fig. 8c). Hence, the AC bridge actuation is a better AC excitation method than the AC cantilever actuation because of the higher amplitude of the nonlinear bifurcation jump.

Effects of damping coefficient
Investigating the effect of damping on the sensor's performance is vital for safe sensor operation and calibration. As particularly for gas sensors which represent the main direct application of the proposed system, the damping could be influenced by the gas type, gas concentration, temperature, and other parameters (especially since we are aiming for multi-gas detection); hence, it is crucial to study the effect of damping on the system dynamics.
The effects of damping coefficients (f 1 and f 2 ) on the resonator response are depicted in Fig. 9. The value of damping influences the AC actuation needed to lead to nonlinear behaviour and the amplitude of the two resonators. Note that the values of damping ratios were chosen arbitrarily, but having the same order of magnitude as previous research studies [38,44], to show the damping effect on the numerical simulations. Under ultra-low damping conditions (f 1 ¼ f 2 ¼ 2:222 Â 10 À3 ) in Fig. 9a, b, the nonlinear jump appears in low actuation V AC1 ¼ 1V, while the softening behaviour first appears when V AC1 ¼ 3V.
When f 1 ¼ f 2 ¼ 1:111 Â 10 À2 , the AC actuation value of the first bifurcation jump and softening branch reaches 3 V and 9 V (Fig. 9c, d). Under high damping (f 1 ¼ f 2 ¼ 2:222 Â 10 À2 ), 9 V AC actuation is not sufficient for the system to exhibit all the nonlinear features, the nonlinear jump first appears at 5 V, and the softening behaviour does not exist at such conditions (Fig. 9e, f).

Effects of resonators' length
We investigate the effect of the resonators' lengths L 1 and L 2 on the system response. The AC actuation of V AC1 equals to 7V in both cases, as shown in Fig. 10. It can be noted that the bridge length influences only the bifurcation jump frequency of the first natural frequency response W1 and does not influence the peak frequency of the second natural frequency response W2. In contrast, the cantilever length only affects the peak frequency of the second natural frequency response W2 and does not change the bifurcation jump frequency of the first natural frequency response W1.

Effect of the resonator thickness
The thickness directly changes the stiffness of the two resonators: the thicker the resonator, the stiffer it gets and requires further AC actuation. The results are depicted in Fig. 11 with AC actuation of 7 V. Figure 11a shows that under the same AC actuation of 7 V, the response may exhibit all possibilities: the

Sensing scheme
The purpose of this research is simultaneously detecting two different physical stimuli by monitoring the dynamic response around the first two lowest modes of the single coupled structure. Hence, we desire to find that stiffness variation and mass variation change the response of the whole system in different ways, as Fig. 12 shows exactly. In Fig. 12a, b, the non-dimensional stiffness perturbation from -0.1 to 0.1 is given to the bridge resonator. The results show that bifurcation jumps in the first natural frequency response W1 in Fig. 12a changes due to stiffness variations on bridge resonator, while the peak value of the second natural frequency response W2 in Fig. 12b keeps constant. When it comes to mass perturbation on the cantilever resonator, the nonlinear jump remains in same largely in Fig. 12c. In contrast, the cantilever resonator's peak frequency and peak values in Fig. 12d are varying. Hence, it is possible to detect both the stiffness and mass perturbation applied on the bridge and cantilever resonators, respectively, by monitoring the nonlinear jump and peak frequency of the first two natural frequency responses. Thus, multi-sensing can be robustly performed.
Besides, the coupling between different modes of the structure would be an exciting topic to investigate. In previous research, the two lowest out-of-plane modes of vibration of the proposed structure (occurring at 178.22 kHz and 189.54 kHz) have a frequency ratio of 3.4 with the first two lowest in-plane modes, which suggests no internal resonance activation for the present case study among these modes. However, the activation of the internal resonance between in-plane and out-of-plane modes could happen for optimized structures as proven in the literature [45,46]. This aspect would be addressed in future research investigating the sensing sensitivity as activating such nonlinear modal coupling.

Conclusions
In this paper, the nonlinear dynamics of two mechanically coupled micromachined resonators (cantilever and bridge resonators) were numerically investigated for potential multi-sensing applications. The concept is based on the simultaneous tracking of the resonance frequencies of the first and second lowest vibration modes. Stiffness and mass perturbations of the bridge and cantilever resonators, respectively, were found to have independent influence on the two vibration modes, demonstrating the promising potential of multi-sensing on a single device. Nonlinear behaviour, including bifurcation jumps and peaks as sensing targets, improves the accuracy and sensitivity of the sensor. The numerical model of the coupled system with geometric and electrostatic nonlinear terms is developed, demonstrating the multi-sensing feasibility. The continuous simulation of the structure is obtained, revealing the full nonlinear dynamics of the coupled system and the effect of different parameters, which is vital for the sensor's design. It is worth mentioning that the value of AC actuation directly relates to the system's nonlinearity. Medium levels of AC actuation linked to the nonlinear bifurcation jumps are suitable to gain the sensor's best performance and weaken the risk of unstable softening branches due to excessive driving input. The perturbation simulations demonstrate the response variation under two perturbations. The results indicate that stiffness and mass perturbations of the bridge and cantilever resonators independently influence the first two vibration modes, proving the promising performance of the multi-gas sensing concept.
The research introduces the methodology of the nonlinear coupled resonator in performing multiparameters detection, which shows its potential for mixture gas analysis, vehicle propulsion system monitoring, and any industrial setting that requires multiple sensing. Future work will focus on the practical sensor system design and fabrication, specifically, the multi-gas sensor for binary gas mixture analysis with the matched testing circuit and platform. Additionally, this research chooses the first and second global modes to perform the sensing strategy. Instead of the loworder modes, the scheme actuating in higher-order modes may increase the response amplitude and improve the sensitivity, which is also a potential research direction. The above would include the investigation of potential activation of internal resonance among in-plane modes and out-of-plane modes of vibrations.
; t ¼t s ; By substituting the new introduced parameters in (7) into Eqs. (2) and (3), the non-dimensional equations of motion are presented as: with the following corresponding boundary conditions: After multiplying both sides of Eqs. (8) and (9) by the 1 À w ð Þ 2 , and applying the Galerkin method [41], it yields: 1 À w 1 ð Þ 2 w 0000 1 þ 1 À w 1 ð Þ 2 o 2 w 1 ot 2 þ 1 À w 1 ð Þ 2 c 1 ow 1 ot The solutions of Eqs. (12) and (13)  By substituting Eqs. (14) and (15) into Eqs. (12) and (13), multiplying by u 1;i , u 2;i , and integrating from x ¼ 0 to 1, x ¼ 0 to R L , it yields: where the dk and dm represents the linear stiffness variation on bridge resonator due to the external stimulus and the mass perturbation on cantilever resonator, respectively. The first and second local mode shapes for the bridge and cantilever are: where the coefficient K takes a value ensuring R R L 0 u 2 2 dx ¼ 1.