Coherent energy transfer in coupled nonlinear microelectromechanical resonators (2025)

Design and characterization of the coupled resonators

The proposed coupled micromechanical resonators system (Fig.1a) was constructed with two double-ended-tuning-fork (DETF) resonators connected via a mechanical disc coupler anchored at the base³⁰. The mechanical coupling between resonators is facilitated by the slight deformation of the coupler, due to the stress generated when the two tines of a resonator vibrate in parallel. The coupling is geometry dependent and can be manipulated by designing the parameters of the coupler, as detailed in Supplementary Information Note1. The resonators are capacitively actuated and sensed via the electrodes on either side of the tines. The device is excited only by thermal noise, when the two switches are off. We therefore obtain the intrinsic natural frequencies of the first two modes: \({\omega }_{1}/2{{{\rm{\pi }}}}\approx 122721.6{{{\rm{Hz}}}}\) and \({\omega }_{2}/2{{{\rm{\pi }}}}\approx 122731.2{{{\rm{Hz}}}}\), as shown in Fig.1c. Such a state of two modes having equivalent amplitudes is named as the symmetric state, at where the frequency split \(\Delta \omega /2\pi \approx 9.6{{{\rm{Hz}}}}\) between the two modes is lowest. This split is referred to as the coupling rate (\(\Delta=9.6{{{\rm{Hz}}}}\)), as it is proportional to the coupling factor of the system, which is defined by the ratio between coupling stiffness (\({k}_{c}\)) and stiffness of the resonator (\(k\)). The coupling factor can be experimentally characterized by using the coupling rate^15,31 \(\kappa=\Delta \omega /{\omega }_{1}\approx 7.82\times {10}^{-5}\).

The other critical parameter of the coupled resonators is the quality factor (\(Q\)), which demonstrates the energy dissipation of resonators to the environmental bath. For our device, \(Q \sim\)42000 is characterized by using the −3 dB bandwidth method³² for resonators at both modes, demonstrating a damping rate of \(\gamma \approx 1.46{{{\rm{Hz}}}}\). This calculated damping is applicable in linear operating conditions without modal overlap. Our subsequent theoretical and experimental studies will explore how damping behaves under conditions of significant energy exchange with structural asymmetry and Duffing nonlinearity. Based on the comparisons of \(\kappa > 1/Q\) and \(\Delta > \gamma\), and the case that there is no modal overlap between modes in the linear regime, the coupled resonators system is considered to be operating in the strong coupling regime^33,34,35 with resonators exchanging energy more rapidly than dissipating to the environmental bath.

The experimentally measured noise spectra of the two resonators demonstrates a base voltage noise level of ~300 nV√Hz. With the variation of the tuning voltage, avoided crossing¹³ and loci veering is observed as shown in Fig.1d, e, splitting this diagram into an upper branch as the representative of the higher-frequency mode (\({\omega }_{2}\)) and a lower branch as the representative of the lower-frequency mode (\({\omega }_{1}\)). The energy gets redistributed and confined to a particular mode with \({{{\rm{\delta }}}}\) deviating from the symmetric point where V_Tuning = −5.85 V. The measurements align remarkably well with the theoretical analysis in Supplementary Information Note2 and Fig.S4. Combining energy redistribution and loci veering, it can be remarked that mode localization^36,37 is evident as an intrinsic character of the coupled resonators. This experimental observation of thermal noise spectra offers an alternative visualization of localization, contrasting with the frequency sweeping method as shown in Supplementary Information Fig.S5.

Frequency responses analysis with nonlinearity

When concerning the coupling, cubic nonlinearity term and noise and assuming the resonators have equivalent damping, the dynamic equations of the system can be written as:

where \({x}_{1}\) and \({x}_{2}\) are the displacements of the coupled resonators, \({\omega }_{0}\) the initial natural frequency, \(\delta=\Delta k/k\) the stiffness asymmetry between resonators, \(k\) the symmetric resonator stiffness, \({\omega }_{0}=\sqrt{k/m}\) the initial frequency of the first mode, \(f=F\sin (\omega t+\theta )\) the applied force to Res 1, and \({\xi }_{1}(t)\) and \({\xi }_{2}(t)\) the random noisy forces. Nonlinear springs for the two resonators are defined as \({k}_{n1}=k(1+{\beta }_{1}{x}_{1}^{2})\), and \({k}_{n2}=k(1+{\beta }_{2}{x}_{2}^{2})\), where \({\beta }_{1}\) and \({\beta }_{2}\) are the Duffing nonlinear coefficients.

With drive level rising, nonlinear response is evident characterized by an increased in the natural frequency with drive level, which is defined as stiffness hardening effect or amplitude-frequency (A-f) effect, as shown in Fig.2. The linear and nonlinear frequency responses can be theoretically modelled based on Eqs. (1) and (2) using the multiple scales method³⁸. In the nonlinear regime, the coupled resonators demonstrated Saddle-Node bifurcations¹³ with transition from the stable branch to the unstable branch for both modes, as shown in the upper plots in Fig.2, and Supplementary Information Fig.S10.

In the symmetric scenario with \(\delta \approx 0\), the coupled resonators exhibit similar stiffness hardening effect, implying comparable A-f curvatures, with equivalent nonlinear coefficients for both modes: \({\beta }_{1}={\beta }_{2} \! \sim \! 1.0\times {10}^{10}{m}^{-2}\). The onset of nonlinearity marking by the critical amplitudes³⁹ is evaluated by numerical simulation, as shown in Fig.S9. The nonlinear terms here are extracted using the fitting equation³⁹ \(\omega={\omega }_{0}(1+\frac{3}{8}{x}^{2}\beta )\), while the backbone curves of the coupled resonators demonstrating the A-f effect are numerically modelled using the derivations (S24a) and (S24b) in Supplementary Information Note3.

In the asymmetric scenario, vibrational energy is redistributed due to mode localization³⁷. In such a case, the vibration amplitude and frequency coupling behaviors are further modified. The apparent A-f effect is modulated simultaneously by nonlinearity and degree of asymmetry. The nonlinear coefficients of the two resonators at the two modes can no longer be fitted well using the above fitting equation³⁹, but should be fitted using the derivations from Supplementary Information equations (S24a, b). This can be explained by the presence of the apparent amplitude-dependent damping effect in asymmetric cases. While our theoretical framework explicitly assumes the absence of nonlinear damping terms in the governing equations, the asymmetric coupled resonator system exhibits apparent amplitude-dependent damping characteristics, as shown in Supplementary Fig.S11c, d. These are observable through comparative analysis of force-normalized response amplitudes under conditions of stiffness mismatch. Such apparent nonlinear responses set limitations in observing energy transfer and localization using conventional frequency-domain analysis due to the absence of discernible mode splitting or intermodal energy transfer in frequency-swept responses in nonlinear operation cases as shown in Fig.S12. This characterization necessitates time-resolved transient analysis of ringdown dynamics, where energy transfer can be observed through decay rate extraction.

Nonlinear characterization with transient responses

The amplitude decay of a single resonator is the transient response when switching off the actuation. It is of an exponential nature with the decay time related to the damping rate γ=mω ⁄ Q^40,41. Energy decay in coupled resonators is more complicated as it may show significant beating²⁹ or Rabi oscillations^41,42,43. However, the beating phenomenon is not always observable in coupled resonators in previous studies, as it requires a sufficient low coupling rate to facilitate noticeable energy transfer. Our basic coupling rate of Δ=9.6 Hz provides sufficient energy transfer per cycle, and our ringdown phase-locked-loop measurement setup ensures the possibility of tracking the evolution of the transient responses.

The time-resolved transient signal amplitudes of the two resonators with initial drive AC of 5 mV and 20 mV are shown in Fig.3a, b. A significant beating phenomenon can be observed. The first pulse is with energy being transferred from Res 1 to Res 2 and exactly in the moment when all energy is transferred, the amplitude of Res 1 keeps constant, and Res 2 starts to transfer energy back. Both spectrogram and PLL-tracking results demonstrate frequency coinciding with the natural frequency of the coupled system i.e., \({\omega }_{1}/2\pi\). If the two resonators operate initially nonlinearly, the spectrogram in Fig.3d shows a progressive decreasing frequency from ~122.740 kHz which is the result of the A-f effect, and ultimately converges to the corresponding small-amplitude eigenfrequency \({\omega }_{1,{{{\rm{linear}}}}} \! \sim \! 122.727{{{\rm{kHz}}}}\) which is defined by the structural asymmetry \(\delta\). More frequency and amplitude ringdown information can be found in Supplementary Information Fig.S18.

The phase difference between resonators, as depicted in Fig.3e, f, demonstrates a locked condition with excitation on, yielding a mean value \({\varPhi }_{0} \sim -0.25\) radians, in good agreement with our measurements using the frequency-phase sweep method (Fig.S8). When ringdown starts, the phase difference shows oscillatory variations between \({\varPhi }_{\max } \! \sim \! 0.38\) radians and \({\varPhi }_{\min }\) close to \({\varPhi }_{0}\). It is known that the in-phase mode and out-of-phase mode introduce a phase difference of \(\varPhi=0\) and \(\varPhi=\pi\), respectively. Such oscillatory observations on phase difference indicate that the resonators are oscillating close to the in-phase mode in the steady state and vibrate in a hybrid mode resulting from the superposition of both modes during the ringdown process.

By extracting the envelopes of transient signals, the amplitude ringdown curves can be obtained. As the frequency of Res 1 has been tracked using a PLL, it is easy to draw the amplitude-frequency plot in the ringdown process, as indicated in Fig.3g–j. It is interesting that the resonator initially vibrates nonlinearly, followed by an amplitude and frequency decrease aligning the A-f curves as that in the frequency sweeping section in Fig.2a, c until converging to its intrinsic low-amplitude frequency, which is also verified in Fig.S18. The ringdown trajectory of Res 1 always follows the nonlinear backbone curves exactly, which can be fitted using the fitting equation³⁹ \(\omega={\omega }_{0}(1+\frac{3}{8}{x}^{2}\beta )\), as shown in Fig. Fig.S16. However, for Res 2, the ringdown trajectories do not follow such fitting equation if the coupled resonators are asymmetrically localized, as indicated by Supplementary Information Fig.S17, which is attributed to the apparent amplitude-damping in the asymmetrically coupled resonators.

Energy transfer with Duffing nonlinearity

Another key result is related to the capacity to manipulate energy transfer by adjusting the initial conditions of the coupled resonators, such as asymmetry and nonlinearity. For this purpose, we changed the initial conditions and recorded a series of ringdown responses for an energy balance analysis. When Res 1 was initially actuated linearly with a drive AC of 1 mV while Res 2 was perturbated by \(\delta\), clear beating phenomenon can be observed in ringdown curves as shown in Fig.4a, b. The simulations in Supplementary Information Fig.S13 closely aligns with the experimental results.

The total energy \({E}_{T{otal}}\) contains the potential energy of each resonator proportional to the summation of the linear and nonlinear square displacement \(\frac{1}{2}{k}_{i}{x}_{i}^{2}+\frac{1}{4}{\beta }_{i}{x}_{i}^{4}\) and the coupling energy \({E}_{c{oupling}}=\frac{1}{2}{k}_{c}{\left({x}_{1}-{x}_{2}\right)}^{2}\) stored in the coupler, as analyzed in the Supplementary Information Note5. \({E}_{T{otal}}\) is always decaying exponentially no matter with or without nonlinearity and asymmetry, revealing a constant damping rate of \({\gamma }_{0}=1.51{{{\rm{Hz}}}}\), as shown in Fig.4i. This damping to the environmental bath includes all the dissipation like anchor damping, thermo-elastic damping and air damping⁴⁰. This \({\gamma }_{0}\) is equivalent to that of the coupled resonators with deep asymmetry (V_Tuning = −4.5 V) as shown in Fig.4a, which is extracted based on the fitting equation of \({x}_{1}=X{e}^{-2{{{\rm{\pi }}}}\gamma t}\) and the method of recording the serial damping rates are introduced in Supplementary Information Note5 and Fig.S20. The value of \({\gamma }_{0}\) also matches with the calculation of \(Q={\omega }_{0}/4\pi {\gamma }_{0}\) using the −3 dB bandwidth measurements in the frequency response section.

The alternating energy exchange pattern of the resonators in suggest timely-variable damping rates. Such variable damping rates are evidence of energy transfer as the exchanged energy offsets the dissipation of Res 1 and Res 2. To quantitively describe the energy transfer between resonators, the time-variable energy decay rates \({\gamma }_{i}\) (i = 1, 2) of each resonator are extracted using the slice-fitting method shown in Supplementary Information Note5. We define the energy transfer rate as \({\zeta }_{i}={\gamma }_{i}-{\gamma }_{0}\). The case of \({\zeta }_{i}=0\) means that the energy of the mode is only dissipated to the environmental bath, and there is almost no internal exchange between resonators. The case of \({\zeta }_{i} > 0\) means that this specific resonator is releasing more energy to its coupled companion than to the environmental bath. Vice versa, \({\zeta }_{i} < 0\) means that the resonator is scavenging from its coupled companion more than dissipation to the environment bath. The specific case of \({\zeta }_{i}=- \! 1.5{Hz}\) is the initial steady state with external energy compensating the dissipation to the environment.

As expected, it is Res 1, the directly actuated resonator, that first exhibits the abrupt energy loss, transiting from \({\zeta }_{1}\approx -1.5{Hz}\) to \({\zeta }_{1}\approx 1.5{{{\rm{Hz}}}}\). The oscillation rate at which the energy transfer occurs varies with \(\delta\), as shown in Supplementary Information Fig.S21. With linear initial states, the statistic maximum and minimum values of \(\zeta\) are \(\sim {\gamma }_{0}\) and \(\sim -{\gamma }_{0}\), suggesting that the energy exchanged between resonators is at most twice the energy dissipated due to damping. The periodic oscillation of \(\zeta\) in Fig.4c, d are attributed to the energy coupling and transfer between resonators. There is no apparent amplitude-dependent damping in this case as shown by the energy ringdown curves in Fig.4i where the mean slopes are similar for Res 1, Res 2 and the total energy during the whole process.

A significant difference with nonlinear initial states is that the maximum value of \(\zeta\) for Res 2 is much larger than \({\gamma }_{0}\), as high as >40 Hz as shown in Fig.4h. This implies that the scavenged and lost energy of Res 2 can significantly surpass its intrinsic energy dissipation. The average damping rate of Res 2 in the deep nonlinear configurations is higher than \({\gamma }_{0}\) as shown in Fig.4h, and the average of Res 1 is lower than \({\gamma }_{0}\) correspondingly. The coupling rate in the nonlinear regime varies with time and finally converges to its intrinsic value defined by the asymmetry \(\delta\) if the resonators are initially vibrating nonlinearly, as indicated by Supplementary Information Fig.S21b.

The ringdown responses of the asymmetrically (V_Tuning = −5.20 V) coupled resonators with increasing drive AC are shown in Fig.4e, f, while the corresponding simulations are provided in Supplementary Information Fig.S12. A key observation is that the amplitudes of the resonators have more visually apparent oscillatory behavior in the nonlinear regime, indicating the enhancement of energy exchange. This enhancement can be quantitatively assessed by examining the power flow²⁷, as illustrated in Fig.4j, l. It is obvious that with linear initial configuration, the highest power flow will significantly decrease by a factor of 4, when the tuning voltage changes from V_Tuning = −5.85 V (the symmetric case) to V_Tuning = −5.20 V (the localized case). On the contrary, if the initial configuration is nonlinear with a drive AC of 26 mV, the highest power flow in the localized case increases by a factor of 2. The maximum power flow during the ringdown process as a function of Drive AC level and perturbation voltage is shown in Supplementary Information Fig.S22, clearly demonstrating the comprehensive manipulation of nonlinearity and localization on energy transfer.

With the impact of the apparent amplitude-dependent damping induced by the asymmetry, the damping rates decrease to a specific converged value for Res 1 whereas that of Res 2 increases from the start of the ringdown, as shown in Fig.4g, h. Apparent amplitude-dependent damping is also observable in Fig.4k where the initial decaying of Res 1 is slower while Res 2 is faster. The effective damping rates of the two resonators depend on the amplitude ratios X₁/X₂ as well as the phase difference, which will be modulated by structural asymmetry and nonlinearity as explained in Supplementary Information Note3.

Furthermore, if we continue increasing the actuation level to > 100 mV, 1:1 internal resonance with frequency locking is seen, as shown in Supplementary Information Fig.S6. In this state, both frequency and amplitude are locked^27,28, as shown in Fig. S7. The ringdown responses with such initial frequency locking are demonstrated in Fig. S23. There is a period where the energy of the two resonators is released at a constant rate of ~40 Hz no matter the level of initial asymmetry. This period is named as the internal resonance release regime.

Energy localization in coupled nonlinear resonators

Previous experiments in Fig.S12 have demonstrated that the energy localization behaviours in the asymmetrically coupled nonlinear resonators cannot be revealed clearly using the frequency-domain analysis even though previous research has demonstrated mass sensors operating in this regime⁴⁴. The eigenstate, represented by the amplitude ratio X₁/X₂, is the other intrinsic metric of the coupled resonators¹⁴. Therefore, its transient dynamic processes would provide additional insight into system dynamics and intrinsic nature. The time-domain transient amplitude ratio measurements provide us with the possibility of revealing further insights into the behavior in this regime, both qualitatively and quantitively.

It can be seen from Fig.5a that the time evolution of the amplitude ratio has similar oscillatory characters as that of the energy transfer rates in Fig.4a, b. We integrated all the amplitude ratio ringdown responses of the coupled resonators with linear initial states and positive and negative perturbations ranging from V_Tuning = −4.00 V to V_Tuning = −7.60 V in Fig.5b, c. This is another kind of portrait of the energy localization phenomenon shown in Fig.1d, e.

If the initial operation state is nonlinear, the heightened asymmetry amplifies the oscillatory frequency as well as the energy transfer as shown in Fig.5e, f. At the same time, only one mode, here the 1st mode can be identified as the device demonstrates considerable stiffness hardening. Notably, the time-domain amplitude ratio evolution method can still reveal the energy localization in Duffing nonlinear or even 1:1 internal resonance regime as provided in Fig.S23, which is a significant advantage over the swept frequency response method. For both linear and nonlinear configurations, the coherence time decreases as structural asymmetry and noise increase. Introducing larger perturbations in the coupled resonators results in shorter lifetimes and noisier amplitude ratios during the ringdown process. The linear configurations exhibit a more pronounced sensitivity to perturbations, leading to accelerated decoherence, whereas the nonlinear configuration demonstrates greater robustness against such disturbances. This is in good match with our previous prediction that the best resolution using amplitude ratio as sensor output metric is around the symmetric line³⁷. As demonstrated previously, the apparent amplitude-dependent damping emerges simultaneously with the mode localization effect, as they both present only in cases when resonators are asymmetrically coupled. Such coexistence results in complicated quantitative determination through systematic analysis of the coupling-governed phase term \(\sin (\varPhi )\) under parameter variations as indicated in Supplementary Information Note3.