Persistent Nonlinear Phase-Locking and Nonmonotonic Energy Dissipation in Micromechanical Resonators

Many nonlinear systems are described by eigenmodes with amplitude-dependent frequencies, interacting strongly whenever the frequencies become commensurate at internal resonances. Fast energy exchange via the resonances holds the key to rich dynamical behavior, such as time-varying relaxation rates and signatures of nonergodicity in thermal equilibrium, revealed in the recent experimental and theoretical studies of micro- and nanomechanical resonators. However, a universal yet intuitive physical description for these diverse and sometimes contradictory experimental observations remains elusive. Here we experimentally reveal persistent nonlinear phase-locked states occurring at internal resonances and demonstrate that they are essential for understanding the transient dynamics of nonlinear systems with coupled eigenmodes. The measured dynamics of a fully observable micromechanical resonator system are quantitatively described by the lower-frequency mode entering, maintaining, and exiting a persistent phase-locked period-tripling state generated by the nonlinear driving force exerted by the higher-frequency mode. This model describes the observed phase-locked coherence times, the direction and magnitude of the energy exchange, and the resulting nonmonotonic mode energy evolution. Depending on the initial relative phase, the system selects distinct relaxation pathways, either entering or bypassing the locked state. The described persistent phase locking is not limited to particular frequency fractions or types of nonlinearities and may advance nonlinear resonator systems engineering across physical domains, including photonics as well as nanomechanics.


INTRODUCTION
Physical systems are never fully isolated from their environment, leading to unavoidable energy exchange.Understanding the origins and effects of energy dissipation has been fundamental for the development of science and engineering fields such as quantum information [1,2], cosmology [3], acoustics [4], timing [5], and sensing [6].The linearized description of such open systems' dynamics, a bedrock of modern physics, leads to damped noninteracting eigenmodes, where damping and stochastic Langevin forces result from averaging over weak interactions with the many degrees of freedom of the external 'thermal bath' [7].While the exact solutions for nonlinear systems are generally hard to analytically obtain, their dynamics can often be studied by treating the nonlinearities as perturbations to the linearized system, retaining the eigenmodes and eigenfrequencies, but allowing higher harmonics, amplitude-dependent eigenfrequencies and time-variable energy flows between eigenmodes.In low-loss nonlinear systems, these result in rich and complex dynamical behavior both in [8][9][10] and out of thermal equilibrium [11][12][13].
In recent years, micro and nanomechanical resonators have emerged as an ideal experimental platform to study nonlinear phenomena, due to the multiple accessible modes with high quality factors and controllable nonlinearities [9,[14][15][16][17][18][19].The individual mode nonlinearities and intermodal couplings allow the resonant frequencies to vary in time, following the changing amplitudes.Furthermore, the systems can dynamically tune in and out of the internal resonances [20][21][22], defined as resonances between system eigenmodes, wherein commensurate eigenfrequencies enable resonantly enhanced energy exchanges between two or more coupled modes [11].Recent experiments involving nonlinear resonators made of carbon nanotubes (CNT), graphene membranes, and nano-and micro-mechanical structures have demonstrated interesting and unexpected dynamical behaviors both during thermalization and in thermal equilibrium.Long-time-scale energy fluctuations were observed in the thermalized state of CNT, indicative of non-ergodic behavior [10].Multilayer graphene membrane resonators have shown abrupt transitions between two different energy relaxation regimes, with the mode-coupled state characterized by a faster relaxation rate for the observed vibrational mode [12].In contrast, the lower-frequency vibrational mode of a nonlinear two-mode micromechanical resonator exhibited a slow, near-zero decay rate in the mode-coupled state [11].These observations offer intriguing glimpses of non-trivial nonlinear dynamics, begging for more intuitive and universal physical descriptions of the underlying mechanisms.However, experimental progress has been hindered so far by the inability of the current experimental platforms to fully isolate and simultaneously observe all the interacting modes and to selectively prepare them in arbitrary initial excited states.
Here, we fully observe the thermalization transient dynamics of a prototypical nonlinear two-mode coupled system from an arbitrarily prepared excited state and develop a simple yet accurate quantitative model describing it.We discover that during the free-ringdown relaxation the system can either enter or bypass a persistent phase-locked state at the 1:3 internal resonance.For the given starting energies, the outcome is determined solely by the initial relative phase between the modes.The phase-locked state, arising during the free ringdown, can last up to ≈ 2.5× longer than the longest thermalization time in the system and is characterized by continuous energy transfer from a higher-loss-rate mode to a lower one, decreasing the overall system dissipation and making it a nonlinear function of system energy.
The system transient dynamics can be quantitatively understood from the perspective of the lower frequency mode experiencing a time-dependent modulation force from the higher frequency mode [13], which creates stable degenerate Period Tripling States (PTS) [23][24][25].While the PTS persists, the energy flows between the modes, maintaining the exact 1:3 relationship between their changing, amplitude-dependent frequencies.As a result, the lower frequency mode experiences a negative dissipation rate, exhibiting a striking non-monotonic energy evolution.

II. ENERGY-AND PHASE-DEPENDENT RELAXATION PATHWAYS FOR COUPLED NONLINEAR RESONATORS
Fig. 1(a) depicts two isolated low-loss vibrational eigenmodes of a nonlinear micromechanical resonator used in this study, having frequencies  2 / 1 ≈ 3, and coupled by a cubic nonlinear term [24,25].The lower frequency Mode 1 and the higher frequency Mode 2 have independent dissipation rates to the thermal bath, Γ 1 and Γ 2 , and can exchange energy at a rate Γ  ≫ Γ 1 , Γ 2 .Mode 1 exhibits a pronounced stiffening Duffing nonlinearity, and Mode 2 is nearly linear (see Fig. 1(b)).By independent actuation and detection at two separate frequencies, we can prepare the system at prespecified initial amplitudes (solid dots in Fig. 1(b)) and observe each mode's evolution toward thermal equilibrium, shown as the black and red arrows [see Appendix A and Supplementary Note 2].
We find that, starting from the same mode energies and depending on the initial relative phase, the system follows one of two radically distinct relaxation pathways leading to different system energy and oscillation frequency  1 () during ringdown (Fig. 1(c),(d)).One pathway is characterized by a phase-locked state at the 1:3 internal resonance with a lifetime much longer than the longest thermalization time in the system (solid blue line in Fig. 1(c)).The distinctive feature of this system state is the lack of decay in the amplitude and frequency of Mode 1 supported by a constant rate of energy transfer from Mode 2. The intermodal energy transfer and the distinct dissipation rates of the modes lead to different system energy loss rates along the two pathways shown in Fig. 1(d).
From the perspective of Mode 1, the nonlinear coupling with Mode 2 manifests as an applied force modulated at  2 ≈ 3 1 , thus the equation of motion for Mode 1 can be generally written as [Supplementary Note 8]: (1) where  1 is the modal displacement, and  1 ( 1 2 ) accounts for the mode's nonlinearity.Our Mode 1 is well described by a Duffing oscillator with a nonlinear coefficient  1 :  1 2 ( 1 ) =  1, 2 +  1  1 2 .To obtain the 1:3 resonant coupling, the interaction with Mode 2 is described by a nonlinear force term ∝  1 2  2 , which in Eq. 1 is represented by the effective drive  2 ()cos [ 2 ()] ∝  2 ().If the nonlinear coupling can be treated as a perturbation, under the conditions discussed in Supplementary Note 8, Mode 1's effect on Mode 2 phase  2 () can be neglected, and Mode 2 performs as a ringdown Duffing oscillator with gradually varying amplitude influenced by slow energy exchange with Mode 1.In this limit, the amplitude  2 () and the frequency  2 () =  ̇2() vary slowly compared to the Mode 1 dynamics and the fullycoupled equations of motion can be reduced to Eq. ( 1) for Mode 1 subject to a period-three drive from Mode 2. As discussed below, the amplitude evolution of both modes on the long timescales of their damping rates can be quantitatively described by accounting for the energy exchange between them.
A form of Eq.1 for a Duffing oscillator with constant  2 and  2 has been a subject of recent theoretical interest [23][24][25] due to the existence of stable phase-locked PTS in the system, occurring when the drive force magnitude exceeds a certain nonzero threshold:  2 >  2,threshold [24].As we experimentally demonstrate below, the dynamics of our fully coupled two-mode system can be well described by the PTS model for the lower frequency Mode 1 under the influence of Mode 2 of gradually changing amplitude and near-constant frequency.Furthermore, this model may form a suitable basis for describing the dynamics of a variety of low-loss coupled nonlinear systems [see Supplementary Note 4,7].
Such PTS for Mode 1 are schematically illustrated as solid red dots on the phase diagram shown in Fig. 1(e) and f [see Appendix B].The sizes of the PTS basins of attraction (colorfilled areas) are determined by the magnitude of the modulating force, i.e., the amplitude of Mode 2 [26].For a substantial Mode 2 amplitude, Fig. 1(e), the PTS basins of attractions occupy much of the high-amplitude portion of the phase diagram, and the excited system is more likely to fall into the phase-locked state (solid black line) rather than bypassing it (dashed black line).Once Mode 1 is locked in a PTS, it remains there until the amplitude of decaying Mode 2 becomes so small that the PTS disappears (Fig. 1 [see Supplementary Notes 7].Compared to the experimental conditions, the schematic of the PTS separatrix uses a much smaller detuning term [see Supplementary Notes 3] for clarity of presentation.The experimental phase diagram is presented within Supplementary note 9, Figure S8.

A. Phase-locked state at internal resonance
Our experimental system is a vibrating clamped-clamped micromechanical beam that exhibits a fundamental in-plane flexural mode (Mode 1) that can couple to a higher frequency out-of-plane torsional one (Mode 2) [11,21].After preparing the system in an initial excited state we measure the motion of both modes during subsequent free relaxation [see Supplementary Note 1], conducting two series of experiments using different initial states.
In the first series, we prepare the system initially in the phase-locked state by solely driving Mode 1 at the internal resonance [see Supplementary Note 2].During the subsequent free ringdown, as shown in Fig. 2(a), there is a finite period of time, the "coherence time"   (colored area in Fig. 2), where the energy of Mode 1 remains approximately constant (variation < 16 %) while Mode 2 loses energy much faster than its dissipation rate into the thermal bath.Beyond   , the modes unlock from each other and follow their individual exponential decays toward equilibrium [11,12].
The observed evolution of the amplitudes is well described by a simple model, which assumes that each mode's coupling to the thermal bath is independent and linear, and the full system energy is the sum of the energies of individual modes (i.e., time-averaged interaction energy is negligible).This condition is a good approximation for our system, though it does not generally apply to internal resonance systems.Each mode energy is assumed proportional to the square of the amplitude of the fundamental harmonic oscillation,  1(2) ∝  1(2) 2 .
During   , the conservation of energy equates the extra power lost by Mode 2 to the power received by Mode 1, after accounting for their intrinsic losses.The exchange power  is constant because it exactly offsets the dissipation of Mode 1, maintaining a constant amplitude on its Duffing curve at the locked frequency of  2 /3 during   .We therefore can simultaneously fit both mode amplitudes, calibrate them, and plot them on a common energy scale using only one fitting parameter  [see Supplementary Note 5].
In Fig. 2(b), we plot the experimentally measured system energy as a function of time.For comparison, the dashed line corresponds to the expected energy relaxation trajectory for the system without coupling between the modes.The slower relaxation of the coupled-mode system can be explained by noting that in our system Mode 2 has a higher loss rate (Γ 2 /2 ≈ 3.3 Hz) than Mode 1 (Γ 1 /2 ≈ 1.5 Hz), and thus, the energy flow from a higher loss rate mode to a lower one decreases the overall system dissipation.
The temporal evolution of each mode's phase and frequency are shown in Figures 2(c)-(e).A key observation is that the modes remain phase-locked for the duration of   (Fig. 2(c)) and, as a consequence, the relative phase  1 −  2 /3 remains constant, at around 0 rad for this experiment (Fig. 2  As previously mentioned, when the modes are locked, the frequencies of both Mode 1 and Mode 2 remain essentially constant (relative change < ±0.05 %).In the internal resonance Mode 1 exhibits small oscillations (≈ 10 Hz root mean square) around the smoothly varying 1 3  2 /2 ≈ 66.5 kHz.The amplitude oscillations corresponding to the frequency oscillations apparent in Fig. 2(d), and similarly Fig. 3(b), are small due to the large Duffing term [Supplementary Note 9], and cannot be resolved.Mode 2 shows gradual changes smaller than ≈ 10 Hz (≈ 0.015 %), and a larger frequency "pull" and "rebound" of ≈ ±30 Hz (≈ ±0.05 %) at the unlocking point.Away from the unlocking point, these minor Mode 2 deviations could be attributed to its small softening nonlinearity and off-resonant nonlinear interactions.
Since the frequency of Mode 2 remains essentially constant throughout the internal resonance, we posit that here Mode 2 can be considered an adiabatically changing external nonlinear driving force, locking Mode 1 into one of three PTS [see Supplementary Note 8].As the system loses energy, the amplitude of Mode 2 decreases, shrinking the PTS basin of attraction, and when its amplitude drops below a fixed threshold  2,threshold , the PTS disappears [Fig.1(e), (f)].
We repeat the ringdown measurement multiple times and find the relative phase during locking randomly choose one of the three degenerate states of PTS. Figure 2(f) shows the relative phase of 15 groups of ringdown data with nearly the same coherent time.Their relative phase shows discrete distribution as shown in the PTS diagram shown in Fig. 1(e).

B. Nonmonotonic energy-dependence of dissipation rate and phasedependent relaxation pathways
One important advance in this work is that we have full experimental access with sufficient sensitivity/bandwidth to measure the dynamics of higher-frequency Mode 2, and, therefore, the relative phase, as well as the system energy, which could not be done in the previous studies [11,12,21].Another important advance is the ability to independently control the initial conditions of both Mode 1 and Mode 2. As PTS is largely determined by the dynamics of Mode 2 ( 2 ∝  2 ), including the relative phase as well as the amplitude, the independent initial condition control is important for understanding how the parameters of Mode 2 influence the relaxation of Mode 1, and thus the system's pathway toward equilibrium.For this purpose, in the second series of experiments we prepare the system in the unlocked initial state with starting energies above the internal resonance (black dots on the two Duffing curves in Fig. 3(a)).After this initial state is prepared by sequentially sweeping the two drive signals [See Appendix A], they are turned off simultaneously and the motion of both modes is recorded.
Since Mode 1 and Mode 2 have Duffing characteristics with opposite sign,  1 and  2 reach  1 : 2 = 1: 3 at particular amplitudes  1,0 and  2,0 during the ringdown, as indicated by green dots in Fig. 3(a).Before that, the modes relax independently with their frequencies changing according to their amplitude-eigenfrequency curves (Fig. 3(b)).Surprisingly, when the frequencies become commensurate, phase-locking does not always occur.We repeated the experiment up to 100 times with the same initial amplitudes but arbitrary relative phases.Once the frequencies reach the internal resonance, Mode 1 frequency either locks to Mode 2 (blue dots) or bypasses it (black dots).These results indicate that the relative phase among modes is key for entering the phase-locked state, as shown schematically in Fig. 1(e) and with experimental data in Fig. S8(c).When the modes lock, Mode 1 frequency is observed to oscillate with gradually decreasing amplitude around the Mode 2 frequency, as Mode 1 settles toward the PTS.Fig. 3(c) shows the Mode 1 energy for the phase-locked and unlocked cases.When the modes lock, the energy can flow between them, while the time-averaged interaction energy remains low compared to the individual mode energies -no abrupt frequency shifts are observed at the locking and unlocking times and the modes continue to follow their individual amplitude-frequency relationships.In the locked state, the net system energy continues to decrease through dissipation, while the frequencies of both modes evolve together in the direction that lowers the net energy, following their frequency-amplitude relationships.Here,  1 =  2 /3 increase together and Mode 2 amplitude drops, while Mode 1 amplitude increases with time [Fig.3

(a) blue arrow].
Separate from the oscillations observed in Mode 1 frequency as it settles toward the PTS, small slow variations in the energy of Mode 1 are observed, possibly arising from the nonadiabatic evolution of the pseudopotential.The resulting dynamic energy exchange between the modes is not described by the simple PTS model, requiring further investigation.For the case where locking is bypassed (black dots), a sharp loss rate increase is visible at the bypassing of internal resonance, consistent with the PTS model [see Supplementary Note 3].Qualitatively, for a certain range of the relative phase, energy transfers away from Mode 1 to Mode 2, therefore briefly accelerating the Mode 1 decay.
The measured time-averaged energy decay rate of Mode 1 is shown in Fig. 3(d) for the locked (blue) and unlocked state (black).As it is evident from the data, the coupling among modes qualitatively changes the effective energy decay rate of Mode 1: it exhibits a striking non-monotonic energy evolution with a negative dissipation rate.When the modes do not lock, the energy loss of Mode 1 remains positive during the entire relaxation time, with a transient peak observed when the frequencies cross the internal resonance.This transient rapid energy loss observed during bypassing is distinct from the sustained faster energy loss observed at higher amplitudes in Ref. 12.It indicates that the energy gain and loss shown in here and Ref. [12], respectively, are not due to their possibly different initial conditions.
The energies of both modes are displayed in Fig. 3(e) and their sum in Fig. 3(f).The data clearly demonstrates that there are two distinct energy relaxation pathways for the system and that they are determined solely by the initial relative phase between the two modes [Supplementary Note 9].Furthermore, because of the ability of the system to dynamically shift the energy between multiple modes with different modal dissipation rates, the system energy dissipation becomes a nonlinear function of the net system energy for as long as the modes are locked.

C. Tunable coherence time and phase-lock probability
With the added flexibility of controlling the initial conditions of the two modes independently, we observe different coherence times.We discover that   can be controlled solely by changing the amplitude of Mode 2 at locking,  2,0 , and is nearly independent of the other state variables [see Supplementary Note 5].Fig. 4(a) shows the measured   normalized by the longest thermalization time constant 1/ 1 in the system, for different  2,0 .The blue squares correspond to the coherence times measured when Mode 1 is prepared initially at the internal resonance (data shown in Fig. 2), and the green diamonds show the measured coherence times when Mode 1 and 2 decay from an initially unlocked state (data shown in Fig. 3).Among all our experiments, which differ by their initial states,   is fully determined by  2,0 , independent of any other variables describing the state of our system.Indeed, the   is given by the time needed for  2 to decay form  2,0 to a threshold amplitude  2,ℎℎ , a minimum amplitude required to maintain PTS, which is a property of the system independent of its state.As  2,0 is controllable,  2,ℎℎ is directly measured (black dashed line in Fig. 4(a)), and the energy dissipation of Mode 2 is fully captured based on energy transfer between modes (fits in Fig. 2, 3 using Eq.S7 ),   can be obtained via calculating the ringdown time of Mode 2 from  2,0 to  2,ℎℎ , Eq. S8 [see Supplementary Note 5].The calculated coherence time is shown as the black line in Fig. 4(a) without using any adjustable parameters.
While the tunable   and non-monotonic energy-dependent dissipation rates offer a novel strategy for dissipation engineering, the probabilistic nature of reaching a phase-lock state during thermalization from a random relative phase must be taken into consideration.When the modes are relaxing from an unlocked state, the initial relative phase of the modes can be set independently from other parameters.Fig. 4(b) shows the statistical results for 100 energy decay measurements under random initial relative phase conditions for each  2,0 .The data indicate that a minimum  2,0 =  2,ℎℎ is required to reach the phase-locking condition and that the likelihood of locking increases with  2,0 approaching unity for large amplitudes.These observations are in agreement with the PTS [24,26] based model, where the basins of attraction for the stable states become larger with increasing  2,0 ∝  2 .Figures 4(c) and 4(d) are schematic of the PTS basins of attraction for large and small  2,0 respectively.The purple dashed circle depicts the Mode 1 states with a given amplitude and random initial relative phases.For larger  2,0 , the majority of the points on the ring lie in the color-filled PTS basins of attraction and thus, the chance for Mode 1 to fall into a PTS is significantly larger.Compared to the experimental conditions, the schematic of the PTS separatrix uses a much smaller detuning term [see Supplementary Notes 3] for clarity of presentation.The experimental phase diagram is presented within Supplementary note 9, Figure S8.

IV. CONCLUSION AND OUTLOOK
A wide class of low-loss systems can be described by a set of independent nonlinear eigenmodes exchanging energy when two or more of their frequencies match, forming integer multiple relationships.
In contrast to trivial brief energy exchanges occurring when the dynamically evolving eigenfrequencies intermittently match, here, using a simple, fully observable, two-mode system, we have shown how eigenmodes can phase-lock at internal resonances and remain locked for time periods up to 6x longer than the system relaxation times (i.e., 1/ 2 ).This allows the modes to rapidly exchange large amounts of energy, producing non-monotonic modal energy evolution, making the locked system energy dissipation nonlinear, and resulting in completely different energy states compared to when the phase-locking has been bypassed.While locked, the energy exchange rate depends on the individual mode's eigenfrequency-energy relations and dissipation rates but does not significantly depend on the exact value of the coupling between them.
We have shown that entering or bypassing the phase-locked state depends sensitively on the initial relative phase, and thus the system's final state is a high gain nonlinear amplifier and discriminator of the relative phase and relative frequency.This phase dependence distinguishes the observed behavior from nonlinear damping recently observed in driven two-mode systems [27].
The experimentally measured relaxation dynamics of the two-mode system, including phase locking, is completely described by extending the well-understood model of a single mode system subject to a commensurate nonlinear drive, e.g., the period tripling state of a Duffing oscillator.The PTS model is accurate when the higher frequency mode dynamics are slow compared to the Mode 1 precession frequency in the PTS pseudopotential, i.e. when Mode 2 undergoes gradual energy change and the corresponding slow frequency shift according to its Duffing nonlinearity.This is clearly valid for our experiment and would be valid for many other systems with low natural decay rates.Additionally, for a system with small time-averaged interaction energy relative to the modes' energies, the long timescale energy evolution can be modeled by a single-parameter energy transfer equation (Eq.S6), providing an intuitive picture simplifying the general fully-coupled equation of motion description (Supplementary Note 8).While phase-locking is predicted within the PTS model, more general theoretical approaches [13] are needed to study the complete parameter range where the discovered phase-locking phenomenon may be observed outside the validity range of the PTS model.
Numerical simulation of our PTS model can qualitatively reproduce and intuitively explain the contradictory experiment results of the rapid energy loss [12], gain [11,21], and non-monotonic energy (this work) for coupled states using different effective Duffing coefficients of Mode 2 (Supplemental Note 4).However, we cannot confirm that our model is applicable to the specific experimental systems other than the one studied here.More complex phenomena could also lie at the root of such contradictory results.Within this model the direction of the energy exchange depends on the specifics of each mode's eigenfrequency-energy relationship near the internal resonance and can be the same or opposite to the one observed in the present experiments.For example, Mode 2 with a linear (zero) eigenfrequency- richer, such systems may exhibit extended locked states, strongly influencing their dynamics under some conditions.Our work opens up a new window for exploration and harnessing of such a general class of nonlinear systems, both in and out of thermal equilibrium.Applying them to many-body systems with large numbers of coupled elements could lead to additional interesting applications, such as tailoring the energy flow direction and efficiency [9] and the broadband sensing with large assemblies of coupled oscillators [28].These results are not limited to nonlinear mechanics, and shed light on coupled nonlinear resonator systems across physical domains, including, for example, nonlinearly coupled optical modes in chip-scale photonic resonators.

APPENDIX A: PREPARATION OF THE SYSTEM'S INITIAL STATE
We use two methods for setting the initial condition of Mode 1 and Mode 2 in the ringdown experiment.The first method is the same as the one we used previously in the internal resonance research where we drive Mode 1 solely to the nonlinear regime and match the drive frequency  1 to  2 /3 [11,21].Due to resonantly enhanced energy transfer, energy flows between the two modes periodically, generating a limit cycle [29] [Supplementary Note 2].In such a case, the two modes are prepared in the phase-locked state before ringdown, i.e., Mode 1 is at the PTS initially.By finetuning the turning-off time on the limit cycle, we can control the initial amplitude (energy) of Mode 2.
When the drive to Mode 1 at internal resonance is strong, the strong interaction between Mode 1 and Mode 2 terminates the limit cycle and both modes drop to the ground states.It limits the maximum  2,0 that can be obtained in this way, and, therefore, the longer   , shown as blue squares in Fig. 4(a).To obtain even longer coherence time and explore larger parameter space, we use the second method that drives Mode 1 and Mode 2 via sequential sweeps of two simultaneously applied independent drives.We first sweep the driving frequency of Mode 1 to a value much higher than  2 /3 while keeping the drive of Mode 2 off.After reaching the high Mode 1 frequency, we maintain the Mode 1 drive while turning on the drive at Mode 2 frequency and sweep it down until Mode 2 reaches the desired initial amplitude.As Mode 2 has a lower frequency at high amplitude due to its spring softening effect, the frequency mismatch between the modes is large and interaction is weaker from the internal resonance.We then tune the driving frequency of Mode 1 down to the desired value.Two modes maintain their initial amplitudes until we start the free decay by simultaneously turning off the drives [black dots in Fig. 3(a)].The initially uncoupled modes allow us to study the locking probability and directly observe the two different energy relaxation pathways of the system.

APPENDIX B: PHASE DIAGRAM OF PERIOD-TRIPLING STATES
In this specific experiment, the mode nonlinearity can be well described by the Duffing nonlinearity and therefore, Eq. ( 1) is rewritten as: It is convenient to write Eq. ( 2) in the rotating frame with a reference of  2 /3 where P and Q are two dimensionless quadratures and  = � 8 1,  3 1 with detuning  =  2 /3 −  1, .In the rotating frame, the dimensionless quasienergy can be written as [24]: where  = is the dimensionless driving strength.The equation of motion for Mode 1 is written as: where  = (5) [26] [see Supplementary Note 3].

REFERENCE
Park, MD 20742, USA

S1-Measurement Setup
In this work, we use a well-established doubly clamped beam micromechanical resonator as the experimental platform to study the energy decay processes of two nonlinearly coupled modes.The silicon resonator consists of a beam with lateral comb-drive electrodes for electrical actuation and detection [1,2].In addition to the electrical measurement, we also perform optical measurement on the resonator out-of-plane motion as a separate measurement channel.Figure S1 shows a detailed schematic of the measurement setup.It is noteworthy to mention that the capacitive signal has a larger readout gain for the in-plane displacement since it is measured as the capacitive variance between the beam and the lateral electrode.On the contrary, the optical measurement, which is more sensitive to the out-of-plane displacement, was performed by a laser Doppler vibrometer (LDV) with a helium-neon laser (λ = 633 nm) and maximum vertical detection limit of ±75 nm.
During the experiment, the doubly clamped beam is placed in a vacuum chamber of the pressure of ≈ 0.1 Pa (≈ 10 -3 Torr), at room temperature.Two a.c.drives of independently set frequency are applied on the side electrodes with a d.c.bias of 6 V.The capacitive (via a low-noise current amplifier) and the optical readout were simultaneously monitored with a lock-in amplifier and an oscilloscope.For the ringdown experiment, we switch off the drives and monitor the motion of the beam from the two readout channels of the oscilloscope whose bandwidth is set to be much larger than the relaxation rates of the modes.The experimental setup developed in this study advances the conventional setups shown in previous works [3], giving us additional control and detection ability.On the one hand, it allows us to independently drive the modes of interest by combining electrostatic forces of different frequencies from the lateral electrodes in an open-loop configuration, making it possible to independently set the arbitrary desired initial condition of the modes.On the other hand, the two channels of detection enable us to detect the displacement of in-plane and out-of-plane modes capacitively and optically, respectively, for achieving the best readout gain and measurement sensitivity.With the appropriate calibration, it allows us to accurately measure the total energy of the coupled-modes-system which is not achievable before [2,3].

Figure S1. Measurement setup:
The schematic shows the clamped-clamped beam resonator and the wiring of the capacitive (black) and optical (blue) detection.In addition to the electrical measurement through the lateral comb electrode, we measure the out-plane mode using a laser Doppler vibrometer.

S2 Two modes of interest and preparation of their initial states
We consider two mechanical modes on the doubly clamped beam.They initially have nearly 1:3 commensurate eigenfrequencies.Figures S2 (a  To monitor the energy decay process of the two modes, we need to prepare the two modes to preset initial conditions.We used two methods to perform the preparation of the initial state.In the first set of experiments [Fig. 2 in the main text], we prepare the two modes at internal resonance [black arrow in Fig. S2 (c)] by only driving Mode 1.The drive frequency is swept up to the internal resonance at  2 /3.The strong nonlinear intermodal coupling transfers energy between Mode 1 to Mode 2 back and forth, inducing limit cycles.Figures S3 (a) and (b) present the time-varying amplitude of Mode 1 and Mode 2, respectively, at a constant drive amplitude and a fixed frequency near the internal resonance.The amplitudes of Mode 1 and Mode 2 periodically vary, with short time intervals when energy is exchanged quickly (faster than ringdown times) back and forth between them.We turn the drive to Mode 1 off at a specific time, for example, as labeled by the black dashed line.The initial conditions for the two modes at the beginning of the free ringdown are the corresponding values of amplitude at that time.It is noteworthy that although we only drive Mode 1, the maximum energy store in Mode 2 can exceed Mode 1 energy by up to 10 times at certain points in the cycle shown in Fig. 2(a).In the second set of experiments [Fig.3 in the main text], we separately drive Mode 1 and Mode 2 by two independently controlled drives, by which we can set the initial conditions of the two modes to whatever condition we want.However, later experimental results show that the coherence time only strongly depends on the amplitude of Mode 2 at the beginning of locking ( 2,0 ).

S3-Simulation of Period tripling states
As shown in Eq. ( 1), Mode 1 performs as an anharmonic oscillator with resonance frequency dependent on amplitude, and it is subject to the external excitation from Mode 2 via the nonlinear coupling term ∝  1 2  2 .In this specific experiment, the amplitude-frequency relationship for Mode 1 is well described by the Duffing form, therefore, Eq. ( 1) is rewritten as: where  1 stands for the Duffing coefficient and The product of the cos( 2 ) term and the approximatelyharmonic displacement term  1 2 ∝ cos( 2 /3) 2 provides a driving force with a harmonic component at the frequency  2 /3 =  1, +  to drive Mode 1.
Following Ref [4], it is convenient to write Eq. (S1) in the rotating frame with a reference frequency of  2 /3 where P and Q are two dimensionless quadratures and  = � 8 1,  3 1 . In the rotating frame, the dimensionless quasienergy can be written as: where  = By solving for  ̇= 0 and  ̇= 0 , we obtain 7 solutions.Then we use  ̈< 0 and  ̈< 0 to select the 3 stable states (red dots in separatrix plot), 3 saddle points, and 1 trivial-state.
The separatrix in the P-Q phase diagram is on the ridge of (, ) which separates the phase diagram into 4 regimes.To obtain the separatrix, we set the initial condition P(0), Q(0) at around the saddle points which are on the minimum point of a ridge of (, ), and do gradient ascendant of P, Q , i.e. inverting time in (S4) and numerically integrating, as: The trajectory of P(t)-Q(t) is the separatrix.Note, the present separatrixes in the paper have  = 0.1, while our real system losses are much smaller,  < 0.001.We use bigger  for visual clarity of the diagrams presented in our paper, to avoid the helix of separatrix winding too dense and narrow.The performance of the system is qualitatively the same for different .In order to show the separatrix under different driving strength, we set medium  = 0.0800 for Fig. 1(e) and the red line in Fig. 2(f), and small  = 0.0034 for Fig. 1(f) and the blue line in Fig. 2(f).Large  = 0.25 and medium  = 0.05 are used for Fig. 4 (c) and (d), respectively.
To illustrate the evolution of an oscillator under a slowly time-varying driving strength (mimicking the gradual amplitude decay of Mode 2), we generate a movie (Supplementary Movie 1) where the oscillator is starting from two different initial conditions having the same amplitude but slightly different phases (<0.1 rad).For the Movie,  = 0.1 is set as a constant while  linearly changes from 0.25 to 0 during the ringdown.In the beginning, the oscillator freely decays.The red dot depicts the trajectory of the oscillator with an initial phase that allows it to lock into the PTS, while the blue one shows the trajectory that bypasses the PTS and decays to the trivial state directly.At the time when the gradually decreasing  reaches a threshold  0 ≈ 0.033 , the PTS attractors disappear.The oscillator (red dot) unlocks from the PTS and rings down to the trivial state.The simple PTS model qualitatively presents the dynamics of Mode 1 observed in Figure 2.
The same dynamics are also depicted in Figure S4, presenting the amplitude as a function of time during the decay from the two initial conditions, numerically integrating (S4) with f gradually decreased in time, as noted above.The sharp amplitude drop visible when the resonance is being bypassed (blue) is consistent with the measurement in Fig. 3c.The linear-damping relaxation rate is recovered once the modes are unlocked.Note, here the frequency  2 is maintained constant while  is decreasing, therefore, the locked state (red) shows a constant amplitude.This is consistent with the measurements in Fig. 2a, where Mode 2's frequency is nearly unchanged (Fig. 2d).The dynamics of the PTS under time-varying modulation frequency and force is more fully considered in the next Supplemental section, showing Mode 1 exhibiting a non-monotonic energy change.

Figure S4 Amplitude of the oscillators shown in Supplementary Movie 1.
The red and blue lines correspond to the dots in Supplementary Movie 1 with the corresponding colors.

S4-PTS under dynamically varying parameters
As we discussed in the previous section, Mode 1 can be regarded as an oscillator subject to PTS, governed by the Eq.(S1), where Mode 2 provides the external drive with frequency  2 ≈ 3 1 and force  2 ∝  2 , whenever any fast dynamics of Mode 2 due to the coupling can be neglected.This simple model not only can qualitatively explain the behavior of Mode 1 in this experiment, such as the locking and bypass trajectories [illustrated in Fig. 1 (e), (f); Fig. S4], the locking probability [Fig 4(c), (d)] and its striking energy decay reversal.The model can also provide a common general approach for describing a multitude of different experimentally observed energydependent dissipation rates in systems at internal resonance [1][2][3], accounting for both an increased as well as decreased and even negative dissipation in the locked state.
Figure S4 uses the PTS model to explain how PTS locks an oscillator [Mode 1] and makes it persist without decay for a period of time (coherence time), longer than its intrinsic dissipation time, during which the driving force is gradually decreasing.It describes the behaviors of Mode 1 in Fig. 2(a) where Mode 2 does not show strong nonlinearity, i.e., the driving frequency is nearly constant  2 ≈  2, .In order to explain the non-monotonic dissipation rate showing in Fig. 3(c) by this PTS model, we need to let the driving frequency from Mode 2  2 vary together with the force  2 , so that we can take into account the amplitude-dependent frequency of nonlinear Mode 2. As the reference frequency  2 /3 of the rotating frame is also changing in time in this case, we use the stationary frame to simulate the dynamics [i.e.directly numerically integrating Eq.S1] and extract the amplitude of the oscillation to avoid any additional approximations or ambiguity.However, we note that this is equivalent to slowly changing variables  2 and  2 /3 in the rotating frame, which now rotates at a variable rate.
The equation of motion is shown as Eq.(S1) where  2 decreases from 1.0 × 10 5 rad 2 s -2 nm -1 to 0 rad 2 s -2 nm -1 between t = 0.1 s and t = 0.33 s and /2 increases from 100 Hz to 150 Hz between t = 0.1 s and t = 0.2 s.In Fig. S5(a), the oscillator is locked to the stable PTS between t = 0 s and 0.1 s.As  2 increases, the oscillator exhibits an energy gain, qualitatively similar to the one observed in our experiment shown in Fig. 3(c).With decreasing  2 , it unlocks from PTS at  ≈ 0.2 s where  2 < 0.6 × 10 5 rad 2 s -2 nm -1 and exponentially decays with its intrinsic dissipation rate afterward.The non-monotonic amplitude (energy) clearly indicates that the drive pumps energy into the oscillator.Intuitively, if we treat the period-3 drive as a perturbation, the oscillator (Mode 1) with a positive  1 follows its hardening eigenfrequency-energy dependence, where the amplitude increases with the locked increasing frequency.The sweeping-up frequency of the drive (Mode 2) drags the phase-locked oscillator (Mode 1) to a higher frequency, leading to a larger amplitude (pumping energy into the oscillator, gain effect).This simple simulation qualitatively explains the novel non-monotonic dissipation rate observed in our experiment.More generally, the PTS model can also explain the rapid energy loss effect for oscillators at internal resonance, as has been recently observed experimentally in another system [3].In this case, the nonlinear drive pumps energy out from the oscillator in the locked state.For this numerical experiment, the parameters of the oscillator remain the same, we only reverse the sign of  , making  2 decrease with decreasing  2 , shown as the red lines in Fig. S5(b),(c).As a result, the oscillator decays faster than its intrinsic dissipation rate while locked in the PTS with time-varying driving parameters [between approximately 0.1 s and 0.2 s].The faster decay rate compared to its intrinsic decay rate clearly shows that the time-varying drive pumps energy out of the oscillator locked to the PTS in this case.
It is noteworthy that the rapid energy gain/loss is mainly due to the time-varying detuning which drags the locked mode.The time-varying nonlinear force does not contribute much to the amplitude change.The black line in Fig. S5(a) shows the case where the detuning /2 = 100 Hz is constant, while the driving force  2 decreases following the same trend as before (black lines in Fig. S5(b),(c)).Here the amplitude of the mode does not show rapid changes.It shows a similar trend as Fig. S4 where the mode maintains a nearly constant amplitude during locking ( ≈ 0.1 − 0.2 s) even though the force is decreasing.After unlocking at  ≈ 0.2 s where  2 < 0.6 × 10 5 rad 2 s -2 nm -1 , the mode begins to lose energy with its intrinsic loss rate.
Although the different behaviors at internal resonance shown in Ref. [2,3] (main text Ref. [11,12]) and present work are also observed within our PTS model by simply changing the sign of Mode 2's Duffing coefficient, we do not have enough information on whether this explanation applies to Ref. [3].

S5-Energy transfer between modes and calibration
Since the time-averaged interaction energy is neglected, the conservation of energy during the coherence time equals the energy lost by Mode 2 to the energy gained by Mode 1 after accounting for their intrinsic losses.It allows us to simultaneously fit both mode amplitudes, calibrate them, and plot them on a common energy scale.
After turning off the external driving force, the amplitude of Mode 1 remains constant or even increases for a finite period of time, and the energy of Mode 2 decays much faster than its expected exponential relaxation via linear damping.The equations describing the amplitude of each mode during ringdown for a coupled two-modes system at internal resonance are: where  1 =  1  1 2 and  2 =  2  2 2 are the energies of Mode 1 and Mode 2, which are proportional to their measured amplitude square with proportionality constants of  1 and  2 , respectively, and  is the energy exchange rate between the two modes, which is approximately constant in our system, since  1  ≈ 0 during the coherence time.
By fitting the ringdown amplitude of the two modes with Eq. (S6), we obtain the calibration which is used for unifying the unit of energy of Mode 1 and Mode 2.
By solving Eq. (S6), we obtain the time-dependent amplitude of Mode 1 and Mode 2: where  1,0 and  2,0 are the initial amplitude of Mode 1 and Mode 2 at the beginning of locking, respectively.These equations were used to fit the data of Figure 2(a) and 3(e) with / 2 as the only fitting parameter (black lines in the figure).Without loss of generality, the energy units are chosen such that / 1 = 1.After unlocking, the two modes dissipate energy at the rate given by their intrinsic dissipation rate  1 and  2 , as expected.The blue and red dashed lines in Fig. 2(a) indicate the independent exponential decays of the modes ( As Mode 2 is regarded as the driving force to Mode 1 to support the PTS, we can solve Eq. (S7) to get a coherence time   by considering the ringdown time of mode 2 from its initial amplitude to the threshold of  2,ℎℎ � 0 � ≈ 0.006  as measured in Figure 4(a).With known  2 , / 2 , and  2,0 , the coherence time can be written as:

S6-Uncertainty analysis
There are two groups of data (blue squares and green diamonds) shown in Figure 4(a), distinguished by how their initial condition are prepared.Their error bars are obtained separately for the two data types via slightly different methods.
For the case where Mode 1 is locked to Mode 2 initially (blue squares), the coherence time   is defined as shown in Figure 2(e), where the start of coherence time is well-defined as the time we turn off the drive ( = 0).The end of coherence time is defined as the point of unlocking, and is determined with a larger uncertainty due to the oscillation of frequency.This is the leading uncertainty in   for these data.For quantifying it, we use one period of oscillation at unlocking, shown as positive and negative error bars (y-axis) for each data point.The error bar of  2,0 (x-axis) indicates one standard deviation due to the detection noise of the  2 measurement at the beginning of the ringdown (when the drive is turned off).The measurement noise is constant and the differences in the oscillation frequency are small, thus the data points share approximately the same horizontal and vertical uncertainties, as indicated by the black error bars in Fig. 4(a).
For the case where Mode 1 is not locked to Mode 2 initially (green diamonds), the coherence time is defined as indicated in Figure 3(b), where both the start and the end of locking are uncertain due to the frequency oscillation.The uncertainty of   is defined as an average of the oscillation periods at the start and the end of locking.Regarding the uncertainty of  2,0 , it is a combination of two statistically independent contributions.The first part is from the detection noise of  2 , the same as above.The second part is arising from the uncertainty at the starting point of the coherence period (locking time) due to the oscillation of frequency, as shown in Fig. 2(b).We define the peak-to-valley difference of  2 during the first frequency oscillation period after locking as the second part of the uncertainty.We combine these two parts of uncertainty for  2,0 via room-meansquare and depict as the error bars in the x-axis for the green dots.The oscillation of frequency is initially relatively large and amplitude and frequency vary between data, we characterize the uncertainties separately for each measurement, shown as the red error bars on top of the corresponding groups of green diamond data points shown in Fig. 4(a).

S7-Nonlinear coupled modes in resonance
We conjecture that the model we have presented may provide a useful approximate description for a broad range of systems near integer fraction resonances.Specifically, many such systems may exhibit stable locked states between two nonlinear coupled modes, whenever at least one of them is described by an amplitude-dependent eigenfrequency.The direction of the energy flow needed to maintain the locked state is determined by the modes' eigenfrequency-amplitude dependences (as illustrated in the numerical experiment in Section S4 and further discussed below).The stability of the locked state for different signs of these dependencies, as well as the conditions for applicability of the adiabatic assumption, are interesting subjects for future research, and may be studied further using existing general theoretical approaches [6], however, the following consideration may serve as a useful guide.
In our simplified model, we make a few assumptions.We start by assuming the energy of the system is predominantly described by the sum of the energies of the two modes, neglecting the timeaveraged energy of the interaction term.We further neglect the higher harmonics and describe each mode by its harmonic eigenfrequency being dependent on its energy.Duffing oscillator is just one specific example of such an energy-eigenfrequency relationship.Closely related is the assumption that the modes are constrained to their eigenfrequency-amplitute relationships, largely independent of the motion of the other mode: neglecting the time-averaged interaction energy makes the eigenfrequency-energy relationships for each mode independent of the energy of the other mode.These frequency-energy relationships  , (  );  = 1,2 are the first cornerstone of the proposed model.
Within this description, modes are allowed to exchange energy when their frequencies are commensurate.Energy conservation, after accounting for independent dissipation by each mode, is the second cornerstone of the model, where () is the energy exchange rate (power flow) between the modes, which is time-dependent, generally.(Cf.Eq.S6, where power flow is near-constant for our specific experimental system).The third cornerstone is the persistent phase-locked state on resonance, similar to the period tripling state for a single nonlinearly-driven Duffing oscillator, with the higher frequency mode acting as the nonlinear drive.The final key simplification here is the assumption that the energy, and therefore also the frequency, of the higher frequency mode is changing sufficiently slowly, such that the stable states of period tripling states driven by the higher frequency mode evolves adiabatically compare to the dynamics of lower frequency mode.This specifically assumes the interaction with Mode 1 does not result in significant dynamical changes to Mode 2 frequency and amplitude other than due to the gradual energy gain/loss of Eq. (S9).This assumption is valid under our specific experimental conditions, with only small deviations, i.e., the slight amplitude oscillations observed in the Mode 1 and Mode 2 data in Figure 3c,e, and a small frequency 'pull' near the unlocking point if Figure 2d.Importantly, the Mode 2 phase remains constant or varies very smoothly everywhere except near the unlocking point.The observed Mode 1 phase oscillations are expected for an oscillator settling into its PTS following our model.The exact general conditions for the validity of this assumption remain to be investigated, such as by considering the perturbation from the small Mode 1 motion near the PTS on Mode 2, and the corresponding perturbation force back on Mode 1 using the more general coupled description [6].
The PTS-like locking to an external drive  cos( 2 ) at frequency  2 ≈  1 is not limited to a specific value of  = 3, but may exist for arbitrary integer  ≥ 1, provided a nonlinear drive term is present, of the form  1 −1 cos( 2 ) .Such effective drive term may arise from a nonlinear interaction or from a linear interaction with a nonlinear resonator mode (i.e., linear drive at the -th harmonic).Here  = 1 is a simple linear drive on resonance, resulting in one non-zero state, and  = 2 is the common parametric drive at twice the frequency, resulting in two non-zero states, etc.As has been noted by others [5], the PTS at  = 3 is distinct from the lower .For example, the stable trivial state co-exists with the PTS states over a frequency range.
The existence of such symmetry-breaking states for  ≥ 2 is not limited specifically to Duffing nonlinearity, but rather is the property of any oscillator with an energy-dependent eigenfrequency  1 ( 1 ) near the internal resonance, with the energy  1 =  1, , where  1, is defined by the period-n driving frequency  2 , i.e.  1 � 1, � =  2 /.When the drive is applied by a higher-order coupled mode (Mode 2), we have  1 � 1, � =  2 � 2, �/.Since we are interested in describing systems near internal resonances, we can linearize the   (  );  = 1,2: For a Duffing oscillator with a linear resonance frequency of  0 and a Duffing coefficient , the eigenfrequency-energy relationship is given by  2 =  0 2 + 〈 2 〉 and  = 〈 2 〉, accounting for both kinetic and potential energies, with a stiffness  =  0 2 for a modal mass .To simplify, we assume the , energy, and force are scaled such that  = 1.Assuming 〈 2 〉 ≪  0 2 we have 3 , or, in the form of Eq (S10), It is evident that for a nonlinear mode near a PTS with given   , (  ) and = 0, we have the familiar PTS-like state with a constant-frequency nonlinear drive from Mode 2.Here we note that for Duffing Mode 1 changing a stiffening nonlinearity with a positive detuning to a softening nonlinearity with a negative detuning retains the same PTS but with a mirrored phase diagram in the  2 /3 rotating frame (Figs.1e, f).More generally for modes of arbitrary energy-frequency relationship, the mirrored phase diagram will also present in the vicinity of the internal resonance upon changing . In the rotating frame, near the resonance, the Mode 1 frequency-amplitude dependence is fully defined by and the �  , Eq (S9) and (S13) can be solved to find the energy exchange rate (), the time evolution of the system energies   () and, therefore, the frequencies   �  ()�.Eq. (S13), when combined with Eq. (S9), describes the system evolution in the locked state on the long timescale as it is ringing down, neglecting any dynamics near the locked state.
Understanding the exact fully-coupled dynamics on the short timescale near the locked state and the general conditions for the existence of the persistent locked state is beyond the scope of this study.However, the PTS-like locked state for  2  2 ≠ 0 can be semi-quantitatively understood in the simplified limit, Mode 2 is still treated as the nearly-constant external drive for Mode 1, but now with Mode 2 frequency varying with changes in its energy as a function of time.Similar to previously considering Mode 1 in the rotating frame of a constant frequency  2 /, we now consider Mode 1 in the reference frame rotated at the variable frequency  2 ( 2 )/ defined by Mode 2. Neglect energy loss due to dissipation,  1 +  2 =  1, +  2, , and using Eq (S10) near a common resonance frequency we obtain Eq. (S14) has the same form as Eq.(S12) for the locked state dynamics with linear Mode 2, except for having an effective Mode 1 frequency-energy dependence Therefore, Mode 1 dynamics near the locked state with both modes nonlinear may be quantitatively described similarly to the case with linear Mode 2, but using an effective nonlinearity for Mode 1..However, we admit that a more consistent treatment of the dissipation for each mode, as well as a consistent fully coupled treatment of the system in general, has to be applied for a more accurate description.The analysis provided in this section is only intended to illustrate how the two-mode PTS-like locked state model we developed to describe our experiments might be generalized to include higher integer frequency ratio  as well as an arbitrary eigenfrequency-energy dependencies   (  ) for each of the two nonlinear modes.In conjunction with Supplemental Section 4, this illustrates how, depending on each Mode's nonlinearities, the persistent locked state leads to either an increase or a decrease in the apparent observed decay rate of a given mode.When appropriate nonlinear coupling terms exist, multiple modes and applied drive stimuli can resonantly exchange energy when the frequency matching conditions of the form ∑     = 0 are satisfied for some integer   .Here the external stimuli are mathematically equivalent to infinitely stiff and massive lossless modes.Since the lowest-order nonlinear terms not forbidden by symmetry are often dominant, the coupling is strongest for the lowest number of participating modes and smallest integers.Such coupled nonlinear systems are commonly encountered and studied, including photonic, nanomechanical, cavity optomechanical, RF and other domains, where nonlinear signal generation and transduction are of high interest.These systems are less constrained and general conditions leading to the persistent locked states are unclear.However, the existence of such stable resonances is apparent, considering, for example, the case of Mode 1 being nonlinearly driven by a combination of Modes 2 and 3 subject to  1 =  2 +  3 and having high energy, stiffness and modal mass, such that the backaction of Mode 1 on them is negligible.Perhaps the simplest possible example of such interaction is encountered by subjecting the Mode 2 of our system to a constant linear drive.Initializing Mode 1 with sufficient energy to lock into a PTS leads to a system locked in a period-tripled symmetry-broken steady state.

S8-The fully-coupled model of two nonlinearly coupled oscillators at the 1:3 internal resonance and how it is reduced to a PTS-like model
The free-ringdown system of two-coupled Duffing oscillators at the 1:3 internal resonance can be described by two equations of motion as [3,6]: with an interaction energy term   =    18) and ( 19), agree with the measured energy data as well as the timeaveraged energy evolution fit in Fig. 2a from the simple model.Beyond that, the complete model reproduces the experimentally observed fast phase oscillation as the mode1 rotates in the PTS pseudopotential with precession frequency Ω, also shown in Movie 1.
Another option is using the calibration constant obtained from SI section 5 which is based on the assumption of constant energy flow W. The calibration constant in Eq. (S6) �c1/c2 =  1  2 is calibrated to be ≈ 8.8, used to put the mode energies on the common energy scale in Fig. 2 and 3.By fixing the ratio  Eq.S19 can be rewritten as Experimentally we observe that mode 2 frequency and amplitude change gradually, within the experimental resolution, in all regimes except very near the unlock time.This suggests considering Eq.S20 in the limit where the rate term is small.In this limit, mode 2 oscillates at the frequency dictated by its linear frequency and Duffing nonlinearity backbone, regardless of the dynamics of mode 1 (seen in Fig. 2d, 3b, red dots), i.e. mode 1 dynamics barely affect the frequency and phase of mode 2. Specifically, the phase of mode 2 oscillation becomes fully determined by the Duffing nonlinearity and changes slowly with the decreasing mode 2 amplitude, on the timescale determined by the damping rates of the system.This allows considering Eq 1 for mode 1 with interaction force of slowly changing amplitude and frequency, and applying the insights that have been previously developed for a single mode period tripling state [4].
Within this single-mode PTS description, near one of three stationary lock points, a high-Q mode (e.g.mode 1) undergoes precession in the local pseudopotential at a frequency Ω, which defines the fastest dynamical timescale.The approximate PTS-like description is applicable when the force changes slowly on this timescale, i.e. both the PTS stationary point location and the pseudopotential profile change adiabatically.Qualitatively, the PTS description is self-consistent when the rate term �, can be comparable to its weak intrinsic damping Γ 2 .Therefore we retain this part of the interaction, which leads to the energy flow  between modes in Eq. (S6), and increases the overall mode 2 loss rate, as seen in Figure 2a, red, and 3e, red dots.
To summarize, at the above analytically discussed and experimentally demonstrated limit, the major dynamics of mode 2 (including phase/amplitude) in the rotating frame are on the time scale of ~1/Γ 1,2 .Given the low energy loss rate 1/Γ 1,2 , mode 2 acts as an adiabatically varying nonlinear drive to mode 1. Mode 1 dynamics can be described using the PTS picture with a slowly varying period-3 force at the frequency of  2 ≈ 3 1 .Specifically, from the perspective of mode 1 at this limit, mode 2 provides a period-3 harmonic drive with slowly varying amplitude and frequency, i.e.  2 () ∝  2 ()cos [ 2 ()] where  2 ̇() =  2 .As the result, Eq. (S16) is reduced to Eq. ( 1) where the coupling term becomes  1 2  2 ()cos [ 2 ()] from The PTS model provides a valid description of our system, evident by the slowly evolving  2 () and  2 () (given by the measured Mode 2 dynamics) relative to the fast frequency and the small-amplitude oscillations of the locked Mode 1, and frequency-amplitude relationship defined by Duffing backbones.However, as shown by the numerical modeling using the fullycoupled model in Fig. S7, the locking phenomena can even occur with parameters where the PTS model does not strictly apply.A general approach, such as Ref. [6], can be used to explore the broader parameter range where the discovered phase-locking phenomenon can occur.Here we illustrate the generality of the observed phenomena by numerically simulating the fullycoupled model with coupling strength, Duffing nonlinearities, and damping rates varied within two orders of magnitude.As shown in Fig. S7, the phase-locking state remains there, showing this physical behavior is present across a wide range of multiple system parameters.In Fig. S7(b), (c), (e), (g), and (j), although the condition of PTS locking by a force with a smoothly varying amplitude and frequency already breaks down, evident by the oscillating Mode 2 amplitude, the phase-locked state still exists, evident by the sustained constant time-averaged Mode 1 amplitude and constant relative phase.The description of our system motion can be further simplified from the complete model when an additional condition is valid, namely that Mode 1 also approximately follows its Duffing backbone, i.e. the modification of Mode 1 frequency-amplitude relation by the interaction term is small.This condition is valid in our system since the measured frequencyamplitude relation in Fig. 2 and Fig. 3 are consistent with its Duffing backbone.This condition means that the time-averaged interaction energy is negligible compared to either of the mode energies, including in the locked PTS state.As the result, the extra energy loss from Mode 2 in addition to natural damping equals the extra energy gain of Mode 1, as described by Eq. (S6).Meanwhile, the phase-locked mode amplitudes are fully determined by their common frequency at all times.The time evolution of the amplitudes is determined by the intrinsic energy loss rates decreasing the total energy of the system, while the energy flow between the modes is such that their individual amplitudes remain on their backbones at some common frequency changing in time.In particular, when the mode 2 backbone is linear (near-vertical), mode 1 time-averaged frequency must remain constant to maintain the lock, which results in mode 1 time-averaged amplitude remaining constant during lock as shown in Fig. 2a.Eq. (S6) with a single parameter W achieves good agreement for the energy-time evolution in While this energy transfer view is particularly compelling in its simplicity, it is not a necessary condition for the PTS description of the Mode 1 evolution.The only benefit of this condition in our work is providing a simple and intuitive way to quantitatively fit the slow energy evolution of both modes and the coherent time, while neglecting the faster features of the full PTS model dynamics from Mode 1, such as the small but observable oscillations of Mode 1 phase and amplitude at the precession frequency Ω.
In summary, while we have numerically confirmed that the system can be modeled well with fully-coupled equations, the reduced model advances the physical understanding of the mode interaction process, while providing a good quantitative description.The reduced model makes a connection between fully-coupled modes and the PTS for the lower mode, which intuitively explains the novel dynamics observed in this work: initial-phase-dependent bypassing or phase-locking at 2  /3 relative phase increment; frequency and amplitude evolution, including the non-monotonic lower mode energy evolution; long coherent time and locking probability as a function of initial conditions.

S9-Dynamics of Mode 1 in the PTS pseudopotential and initial-phasedependent trajectory
The dynamics of Mode 1 in the PTS pseudopotentials are shown in Fig. S8.Fig. S8(a) presents in the phase diagram, four groups of experimental data obtained in the locking/bypassing experiment (Fig. 3).Three groups lock to the PTS Basins around distinct stationary points, while one group (black) bypasses the PTS and directly decays to the ground state.Due to the strong Duffing nonlinearity, there is a rapid frequency change with amplitude, squeezing the Basin in the radial direction.Therefore, we replot the same data in Fig. S8 Only initial conditions within the basin of attraction for each PTS stationary point (areas covered by the locking data points) lead to a locked state while other initial conditions (white empty areas) result in bypassing/ground-state, consistent with the illustration in Fig. 1(e) and Fig. 4(c) (d).In the paper, we use the schematic generated with smaller frequency detuning to more clearly illustrate all the features on a single radial scale (no need to zoom in) for facilitating understanding.

S10-Mechanical spectrum of the device
We drive and measure the mechanical spectrum of the resonator up to 1 MHz (≈15x the fundamental frequency) shown as the blue dots in the upper panel of Fig. S9.The gray dots are the measurement noise floor measured at the stationary reference substrate.Besides, we numerically simulate the eigenfrequency of the devices and quantitatively match all the modeled mechanical modes of the system, shown as the lower panel of Fig. S9, eliminating the possibility of hidden modes.The dashed line labels  ×   where   =  2, /3 is the internal resonance frequency and  = 1,2,3, … Except for mode 1 and mode 3 (labeled red), there is no other commensurate frequency matching.As described on the paper Mode 1's linear eigenfrequency is detuned lower than   , but it matches   at a specific amplitude due to its Duffing-nonlinearity frequency shift.
FIG. 1 Energy exchange and relaxation paths of a system consisting of two nonlinear coupled oscillators.(a) Schematic diagram showing the energy exchanged between two nonlinear coupled oscillators with commensurate frequencies (here  1 :  2 = 1: 3) and with the thermal bath.(b) Nonlinearity of Mode 1 allows tuning of its oscillation frequency  1 above, at, or below the internal resonance  2 /3 as a function of Mode 1 amplitude, as shown in a schematic of the individual modes' spectral responses.(c) Oscillation frequency of Mode 1 (blue) and Mode 2 (red) during ringdown.The oscillation frequency of Mode 1 decreases during ringdown due to positive Duffing nonlinearity.At the internal resonance (black circle), Mode 1 enters (solid line) or bypass (dashed line) a phase-locked state where resonantly enhanced energy exchange with Mode 2 FIG. 2 Measured relaxation dynamics of Mode 1 and Mode 2 from a phase-locked state.(a) Energies E 1 (blue) and E 2 (red) for Mode 1 and 2, respectively, are plotted on a common scale, calibrated from the measured amplitudes.The black lines are fits based on energy transfer between modes.The dashed lines are hypothetic uncoupled energy decays calculated with separately measured individual mode relaxation rates.The yellow highlighted area represents the phase-locked time period.(b) Measured system energy E 1 + E 2 (green dots), and the model fit (black line), with the hypothetical uncoupled (purple dashed line) decay for comparison.(c) Measured phase of Mode 1 (blue dots) and Mode 2 (red dots) with reference frequency of  2, /3 and  2, , respectively.We define  2 ( = 0) = 0 as a start phase.(d) Measured frequency.Inset shows the zoom-in of the dashed area, where  1 =  1 −  2, /3 and  2 /3 = ( 2 −  2, )/3.(e) Measured relative phase,  1 −  2 /3 is constant in the phase-locked state.(f) Repeat measurement of relative phase,  1 −  2 /3 for 15 times.They show signature feature of 2/3 relative phase difference for the period-tripling states.

Fig. 2 (
Fig.2(d)shows the measured frequencies of the modes during and after coherence time.As previously mentioned, when the modes are locked, the frequencies of both Mode 1 and Mode 2 remain essentially constant (relative change < ±0.05 %).In the internal resonance Mode 1 exhibits small oscillations (≈ 10 Hz root mean square) around the smoothly varying

FIG. 3
FIG. 3 Measured relaxation dynamics of Mode 1 and Mode 2 from an initially uncoupled high energy state.(a) The modes are prepared at the initial amplitudes marked by black dots by independently driving them along the shown trajectories.With the drives off, following an initial free decay, the modes reach internal resonance at the amplitudes labeled by green dots, and either enter or bypass the phase-locked state, depending on their relative phase.The arrows show the Mode 2 amplitude-frequency trajectory (red arrow) and the Mode 1 trajectory for the unlocked (black arrow) and the locked (blue arrow) states.Both modes' amplitude-frequency trajectories are confined to their respective individual nonlinear eigenfrequency-amplitude relationships.The narrow dip on the amplitude-frequency curve of the driven Mode 1 corresponds to crossing the 1:3 internal resonance at  2 /3, while the drive to Mode 2 is off.(b) Oscillation frequency of locked Mode 1 (blue dots) and Mode 2 (red dots), with comparison to that of the unlocked Mode 1 (black dots).(c) Energy of Mode 1. Blue and black dots correspond to the locked and unlocked cases.The amplitude oscillations corresponding to the frequency oscillations apparent in (b) are small due to the large Duffing term [Supplementary Note 9], and are not visible in the energy data downsampled to reduce noise.The dashed (solid) line is a 0.03 s moving average of the unlocked (locked) data.(d) Energy decay rate of Mode 1 extracted from the averaged locked (blue) and unlocked (black) data in c.The increasing energy of Mode 1 is supplied by Mode 2. (e) Two distinct energy relaxation trajectories with different final energies from the same initial energies for Mode 1 and Mode 2. (f) System

FIG. 4
FIG. 4 Coherence time and the probability of locking.The coherence time and the locking probability increase as a function of the Mode 2 amplitude  2,0 at the initial time of locking and are independent of all other parameters.(a) Measured   in the units of Mode 1 decay time 1/ 1 .Blue squares are experiments preparing Mode 1 and Mode 2 at the phase-locked state initially, while green diamonds are measured when the modes are initially unlocked.The black dashed line  2,ℎℎ ≈ 0.006 V is obtained from the averaged  2,0 at   ≈ 0 s.The black calculated coherence time is obtained via calculated the relaxation time of Mode 2 from  2,0 to  2,ℎℎ without using any adjustable parameters.The one standard deviation (SD) uncertainties for the blue squares are constant for each data point, labeled by the black cross, while the one SD uncertainties for each group of the green diamonds are labeled on the corresponding data groups [see Supplementary Note 6].(b) The probability for locking increases as a function of  2,0 .The results are from 100 measurements with random initial relative phase  1 −  2 /3 for each  2,0 .The uncertainties are one standard deviation of the corresponding Binomial distribution.(c) PTS basins of attraction for strong driving strength.(d) PTS basins of attraction for weak driving strength.The purple dashed rings represent the initial condition of Mode 1 with a fixed amplitude  1,0 slightly above the PTS (red dots) and a random

Γ 1 2𝛿𝛿𝜔𝜔.
Trajectories in the phase diagram, as shown in Figs.1(e), (f), are obtained by numerically integrating (5) from different initial conditions.The separatrices of the phase diagram are further calculated by doing gradient ascendant from the saddle points of (, ) based on Eq.
) and (b) show the measured linear resonance of Mode 1 and Mode 2. As discussed in the previous section, we measure the displacement of the modes both electrically and optically, shown as the blue and orange dots.The Lorentzian fit of the mode shows  1, /2 ≈ 64630.6Hz , Γ 1 /2 ≈ 1.5 Hz , and  2, / 2 ≈ 199881.8Hz , Γ 2 /2 ≈ 3.3 Hz .With increasing the driving force, Mode 1 turns into a nonlinear regime shown as Fig.S2(c) where a spring hardening effect is observed (positive Duffing coefficient).When the a.c.voltage is larger than 0.06 V, the Duffing curve of Mode 1 reaches the internal resonance ( 1 :  2 = 1: 3) during the upward frequency sweep.Due to the strong nonlinear intermodal coupling, it cannot pass the internal resonance with Mode 2 unless the driving voltage is higher than 0.45 V[1,2].For such voltage Mode 1's Duffing curve can be swept to higher frequencies shown as the red line for example.Unlike Mode 1, Mode 2 shows the spring softening effect (negative Duffing coefficient) as shown in Fig.S2(d).The fit based on Duffing nonlinearity (black line) describes well the electrically/optically measured mechanical response spectrum of Mode 2.

Figure S2 .
Figure S2.Two mechanical modes with nearly 1:3 commensurate frequency relationship: (a) Linear, small-amplitude resonance of Mode 1.The blue and orange dots present the measured data from the electrical and optical detection, respectively.The black line is from the Lorentzian fit.(b) Linear, small amplitude resonance of Mode 2. (c) Nonlinear response of Mode 1 for the a.c.component of the diving voltage ranging from 0.01 V to 0.81 V with a 0.05 V increment.(d) Nonlinear resonance of Mode 2 with a hysteresis.The black line is a Duffing model fit.The arrows label the sweeping directions.

Figure S3 .
Figure S3.Preparation of the initial state at the internal resonance.Electrically measured displacement of (a) Mode 1 and (b) Mode 2. A limit cycle under a steady drive is followed by the free decay of the coupled system after the drive is turned off (black dashed line) near Mode 2's maximum amplitude.The orange line labels the vibration amplitude.The electrical capacitive measurement gain for torsional Mode 2 is much smaller than for translational Mode 1.

Figure
Figure S5 Numerically modeled dynamics of oscillators in PTS with time-varying parameters. 1, /2 = 10 kHz , Γ 1 /2 = 0.1 Hz,  1 = 1.0 × 10 4 nm 2 s 2 ⁄ are fixed. 2 ,  change as shown in (b), (c), resulting in the ringdown shown in (a) marked with the corresponding colors.The opposite tuning direction of driving frequency between t = 0.1 s and t = 0.2 s drags the locked mode to a higher/lower oscillation frequency, exhibiting rapid energy gain/loss faster than its intrinsic dissipation rate.

)
Since the PTS (n = 3) for a Duffing oscillator have been extensively studied, it is useful to reduce the motion of any nonlinear mode with   (  ) ≠ 0 sufficiently close to a PTS at   = (  ) to the equivalent Duffing oscillator having the same   , (  ) and   (  ).Here we derive expression (S10) for the Duffing oscillator and explicitly provide the equivalent Duffing PTS parameters  0 ,  and  for a mode with given   , (  ) and   (  ).

((
) , the equivalent Duffing PTS parameters are  =     (  ),  0 = (  ) −     (  ), and   ).The equivalent Duffing oscillator loss rate remains the same as the nonlinear mode loss.Now, we consider again our model where an adiabatically-varying higher-frequency mode drives the lower-frequency mode of energy-dependent eigenfrequency.In a simple case of linear

2 𝑚𝑚 1 𝑚𝑚 2 ≈ 9 .
Fig. S6-1 Fit by complete model.The fitted results are shown as black lines for amplitude (a) and phase (b).The red and blue symbols and the purple line are the experimental data.The complete model, i.e.Eq. (18) and (19), agree with the measured energy data as well as the timeaveraged energy evolution fit in Fig.2afrom the simple model.Beyond that, the complete model reproduces the experimentally observed fast phase oscillation as the mode1 rotates in the PTS pseudopotential with precession frequency Ω, also shown in Movie 1.
the interaction term as a change in the damping rate, 1 is much smaller than Ω.While mode 2 frequency is barely affected by Mode 1, its energy loss rate to mode 1,

Figure S7 .
Figure S7.Simulated energy and relative phase of the two modes using the fully-coupled twomode model with parameters scaled as follows relative to the ones used for Figure S6-2: /10 (a);  ×10 (b);  1 /10 and  2 /10 (c);  1 ×10 and  2 ×10 (d);  1 /10 and  2 /10 (e);  1 ×10 and  2 ×10 (f);  1 /10 (g);  1 ×10 (h);  2 /10 (i);  2 ×10 (j); Within a wide range of the parameters, the persistent locked states are shown.The blue and red lines denote modes 1 and 2's energy, respectively.The black line denotes their relative phase.When the system is initialized with  2 ≫  1 and the relative phase near the expected PTS locked state, an extended phase-locked period of time with near-constant Mode 1 amplitude is observed for each parameter set.
(b), zooming in on the data at R > 0.033 (labeled by the red dashed line in panel (a)), by redefining the radial coordinate as  ′ +  ′ = ( − 0.033)  .The trajectory of ringdown is consistent with our schematic figure in Fig.1(e), clearly demonstrating the locking/bypassing trajectory in the PTS pseudopotentials, and the precession(Ω) and decay() of Mode 1, as is also seen in Movie 1.To illustrate the initial phase dependence of the phase-locking/bypassing, in FigureS8(c), we present the 44 series of data from all those experiments where locking was observed out of the 100 total repeat experiments starting from random, uniformly distributed phases.The conditions correspond to the A2,0 = 0.03 data in Fig.4(b)[the set with locking probability ≈ 50 %].We focus on the data just before entering the basin regions, by further zooming in and setting origin at R = 0.040, i.e.  ′′ +  ′′ = ( − 0.040)  , and truncating the data traces to include only R>0.042 (red dashed line in (b)).The 44 groups of locking data are colored based on which Basin they end up with.The remaining 56 groups of bypassing data, not shown here, concentrate within the white spaces between the discrete colored spirals, uniformly filling the phase space due to the uniform distribution of the initial relative phase.

Figure S8 .
Figure S8.The dynamics of Mode 1 in the PTS pseudopotential with a reference phase of  2 ()/3.(a) Four groups of locking/bypassing data taken at the same experimental conditions as shown in Fig. 3. (b) Same as (a) but rescaling the radial scale by setting origin at R = 0.033 (red dashed circle in (a)) instead of 0. Three groups of data lock and precess around three PTS stationary points, while one group bypasses and rings down to the ground state.(c) 44 groups of experimental data that resulted in locking out of 100 repeated experiments with uniformly-distributed random initial phases.They are colored based on which basin they end up with.The other 56 groups of data (not shown) are all bypassing the lock.The bypassing data concentrate within the white space between the colored spirals.It is comparable to Fig. 4d.For (c), we focus on the data just before entering the basin regions, by further zooming in and setting origin at R = 0.040, i.e.  ′′ +  ′′ = ( − 0.040)  , and truncating the data traces to include only R>0.042 (red dashed line in (b)).

Figure S9 .
Figure S9.Measured mechanical spectrum (blue dots) and simulated mode shape are shown as upper and lower panels.The simulated (up-pointing black arrows) and measured eigenfrequencies match each other, with disagreement <5 % due to fabrication imperfections.Dashed lines are integer multiples of the internal resonance frequency.They indicate there is no other matching frequency at the 1:3 internal resonance.The gray dots are the measurement noise floor.The arrows in the modal shape panel depict in-plane flexural modes, out-of-plane flexural modes, and torsional modes.Mode 8 to 11 are four degenerate beam modes.
energy relationship can lock Mode 1 at a constant frequency and amplitude as shown in Fig 2, while Mode 2 with a hardening eigenfrequency-energy relationship can drag the locked Mode 1 to decay even faster than its natural relaxation rate, opposite to the energy gain effect shown in this work.Beyond the coupled Duffing oscillators, this model can be extended to describe two modes with arbitrary nonlinearities near an internal resonance (Supplementary Note 7).During the preparation of this manuscript, the authors learned that the energy transfer between coupled modes at 1:3 internal resonance is studied by Y. Yan et.al.Complementing the present study of the transient relaxation, Yan et.al.independentlydemonstrate that the same period-tripling approachquantitatively describes the steady-state dynamics of such systems subject to a continuous external drive of the higher frequency mode.Demonstrating the rich dynamics possible in such systems, when the nonlinear second mode is driven sufficiently strongly, it responds with one of two distinct amplitudes.In combination with the period tripling for the lower mode, this leads to a total of 6 distinct discrete timetranslation symmetry breaking states.Parametric 2f drive for mode 2 also results in 6 coexisting states.A single nonlinear mode can be driven parametrically at arbitrary frequency multiples.Similarly, pairwise mode locking can occur at resonances of different integer fractions, other than 1:3 studies here.
Moreover, the proposed description framework may be fruitfully applied not only to driven systems[Yan   et.al.], but also to systems with more than pairwise resonances, i.e., when more than two modes (and/or drives) participate in the resonant energy exchange (Supplementary Note 7).While, generally, much 3 Materials Research Institute, Penn State University, University Park, PA 16802 1,  1 �, together with the loss rate, describes Mode 1 dynamics near the internal resonance when dynamic backaction of Mode 1 on Mode 2 can be neglected.Whenever the locked state  1 ( 1 ) =  2 ( 2 )/ is maintained in time, by taking a time 1  3 2 for 1:3 internal resonance and the values of coupling constant , modal masses  1 and  2 and coefficients  1 and  2 dependent on the choice of the scales of the generalized coordinates  1 and  2 .Takin the experimentally measured voltage values directly as the generalized coordinates, the Duffing coefficients of the two modes obtained from their measured Duffing spectrum are  1 ≈ 2.41• 10 8 s −2 V −2 and  2 ≈ −1.03 • 10 7 s −2 V −2 .=  2   2,  +  2 *  −2,  , and neglecting non-With the separately measured parameters, Γ 1 , Γ 2 , ω 1, , ω 2, ,  1 and  2 , we can fit the results using only two fitting parameters  12 and  21 .For example, data in Fig. 2a can be fit by this complete model shown as Fig. S6-1 where  12 ≈ 3.32 • 10 5 s −2 V −2 , and  21 ≈ 2.98 • 10 6 s −2 V −2 , are fitted by achieving the least-residual between the fit and the amplitude data.For comparing with experiment, we have chosen to use the experimental voltages directly as the generalized coordinates.This choice results in the unequal modal masses  1 and  2 for the two modes, and therefore unequal values of  12 = .For theoretical analysis, it is conventional to choose the generalized coordinates  � 1 and  � 2 such that  � 1 =  � 2 and therefore  � 12 =  � 21 .This can be easily done by rescaling the generalized coordinate for mode 2 as  � 2 =