Weakly nonlinear analysis on synchronization and oscillation quenching of coupled mechanical oscillators

Various oscillatory phenomena occur in the world. Because some are associated with abnormal states (e.g. epilepsy), it is important to establish ways to terminate oscillations by external stimuli. However, despite the prior development of techniques for stabilizing unstable oscillations, relatively few studies address the transition from oscillatory to resting state in nonlinear dynamics. This study mainly analyzes the oscillation-quenching of metronomes on a platform as an example of such transitions. To facilitate the analysis, we describe the impulsive force (escapement mechanism) of a metronome by a fifth-order polynomial. By performing both averaging approximation and numerical simulation, we obtain a phase diagram for synchronization and oscillation quenching. We find that quenching occurs when the feedback to the oscillator increases, which will help explore the general principle regarding the state transition from oscillatory to resting state. We also numerically investigate the bifurcation of out-of-phase synchronization and beat-like solution. Despite the simplicity, our model successfully reproduces essential phenomena in interacting mechanical clocks, such as the bistability of in-phase and anti-phase synchrony and oscillation quenching occurring for a large mass ratio between the oscillator and the platform. We believe that our simple model will contribute to future analyses of other dynamics of mechanical clocks.

Oscillatory phenomena are widely observed in nature and society.There are both desirable and undesirable rhythms.For example, some rhythms in the human body, such as the heartbeat and circadian clock, are essential for homeostasis, while others appear with diseases, such as abnormal oscillations of neuronal action potentials in Parkinson's disease [1] or epilepsy [2].Thus, it is expected that stabilizing the necessary oscillations and removing the abnormal oscillations will contribute to healthy biological rhythms and the treatment of diseases.Mathematically, an oscillatory state corresponds to the limit cycle of a dynamical system [3].Threfore, developing techniques to stabilize or destabilize limit cycles by external perturbations is important.
Several methods for stabilizing unstable limit cycles were established in the 1990s [4,5].On the other hand, regarding the methods to destabilize the stable limit cycles, although various state-transition methods between multiple limit cycles have been extensively studied [6][7][8][9][10], few studies have explored the general principle for state transition from stable limit cycles to stable fixed points [11][12][13][14].Establishing methods for transferring the system state from a stable limit cycle to a stable fixed point is important for annihilating undesirable rhythms [13,14].
In the present paper, we consider the metronome as an example of a bistable system with a stable limit cycle and stable fixed point.We focus on the oscillation-stopping phenomenon (oscillation quenching) of the metronome on a platform in which the metronome stops vibrating after oscillating for a while.Below, we explain the details of our study referring to the previous studies on metronome dynamics.
A metronome is a mechanical device in which a needle with a weight continuously oscillates.Its characteristic structure is that the combination of a spring and gear provides torque in the same direction as the motion of the needle when the needle reaches a certain position [15].This structure, called an escapement mechanism, enables the metronome to oscillate by counteracting the damping force caused by friction.To investigate the various dynamical behaviors of metronomes and their mechanical analogues, e.g., pendulum clocks, numerous experimental studies have been performed [16][17][18][19][20][21][22][23][24][25][26].One of the most famous dynamics of coupled metronomes is synchronization, in which the timing of the metronome oscillation is aligned when multiple metronomes are placed on a common platform.Synchronization was first discovered by Huygens in an experiment using pendulum clocks [17], and many experiments have since been conducted in various settings [16,[18][19][20][24][25][26].In particular, it is known that both in-phase and antiphase synchronizations, or only one of them, can be observed depending on the experimental situation [20].Another unique behavior is oscillation quenching, which has been observed in metronomes on a movable platform [26] or the pendulum clocks suspended on a movable board [17].However, although a lot of previous studies focused on the synchronization of metronomes, few investigated the oscillation quenching.
Numerous modeling studies have also been performed to analyze the dynamics of metronomes and pendulum clocks [15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30].As summarized in Ref [30], several challenges exist with these studies.The first is the modeling of the escapement mechanism.To describe the escapement mechanism, it is considered appropriate to use functions that provide torque in the same direction as pendulum motion.In the previous studies, the van der Pol-type function [16,21,28], the piecewise linear function [18,19,22,23,25,29], Dirac's delta function [26,30], and the function that instantaneously changes the angular velocity [17] or the kinetic energy [24] of the pendulum at specific positions have been used.Another challenge is that the motion equation of a metronome generally becomes a nonlinear ordinary differential equation (ODE) that cannot be solved explicitly.To analyze the nonlinear motion equation, an averaging approximation [16,26,30] and Poincaré map [17,24] have been applied.
Recently, Goldsztein et al. analyzed a mathematical model of two metronomes on a movable platform and obtained a phase diagram for in-phase and anti-phase synchronization [30].In particular, they succeeded in explaining some of the past experiments by considering the metronome's nonlinearity caused by its pendulum structure [30]; that is, they expanded the usual linear small-angle approximation (sin θ ≃ θ) to include the nonlinear term (sin θ ≃ θ + cθ 3 with sufficiently small c).Seeking a better agreement with the experimental results, they also created a more realistic model by assuming Coulomb friction as the damping of the platform [26].Although these studies are elaborate and sophisticated in terms of both modeling methods and analytical techniques, there are several open questions.In their first study [30], the analytical method (i.e., averaging approximation) was applied only when the amplitude of the metronome was larger than a certain threshold value, preventing the analysis of oscillation quenching.In the second study [26], although this issue was resolved (i.e., the averaged system was valid even when the amplitude was small), a stability analysis was not performed because the averaged system contained discontinuous functions.The behavior of the equation of motion was tested only by numerical simulations, and a phase diagram was not created.
We are particularly concerned with oscillation quenching because oscillation quenching can be considered as an example of a state transition from a stable limit cycle (oscillating state) to a stable fixed point (resting state) when the oscillators receive the feedback resulting from their motion via the platform.Thus, this study aims to treat both the synchronization and oscillation quenching in a unified manner, that is, using the same mathematical model.One of the reasons why the analysis of oscillation quenching was difficult in the previous studies [26,30] is that the authors attempted to make the model realistic by modeling the escapement mechanism with a delta function and considering both the nonlinearity of the pendulum structure and the Coulomb friction acting on the platform.Thus, to facilitate the analysis, we model the metronome as a linear spring pendulum, neglect the damping of the platform, and simulate the escapement with several smooth functions, particularly a polynomial of order five.Owing to these simplifications, we analytically and numerically obtain a phase diagram for both synchronization and oscillation quenching, which has not been obtained in the previous studies on metronomes.
The remainder of this paper is organized as follows.First, we consider the case of a single metronome on a movable platform.To analyze the equation of motion using an averaging approximation, we treat the entire system as a weakly nonlinear oscillator by assuming that both the escapement mechanism and damping are sufficiently small.We use several functions to represent the escapement mechanism and confirm that the averaged system reproduces the bistability and oscillation quenching of a real metronome.We then expand our model to the case of two identical metronomes on a movable platform, where we adopt a fifth-order polynomial as the escapement mechanism to make the analysis easier.Assuming that the mass ratio of the metronome to the platform is sufficiently small, we perform an averaging approximation and a linear stability analysis to obtain a phase diagram for the in-phase synchronization, anti-phase synchronization, and oscillation quenching.We verify the analysis by numerically integrating the equations of motion before the averaging approximation and plotting the results on the same phase diagram.Finally, we provide a summary and discussion.

Model
Figure 1 (a) illustrates our model of one metronome on a movable platform.In this model, a point mass m is connected by two springs with spring constant k/2 to a platform of mass M .The platform has one degree of freedom and moves freely.The variables X(t) and x(t) are the positions of the platform relative to the floor and the mass relative to the platform, respectively.
We assume the following two forces that act from the platform to the mass point: the damping force proportional to the velocity of the mass, ẋ := dx/dt, and the active force describing the escapement mechanism given as a function of x and ẋ.By neglecting the mass of the springs, the air resistance of the mass and platform, and the damping of the platform due to friction with the floor, we obtain the following motion equation: where γ is the damping coefficient, δ represents the magnitude of the escapement mechanism, and f (x, ẋ) is a dimensionless function that shows the nature of the escapement mechanism.
We denote the position of the center of mass as x c .Then, Since we assume that there is no external force acting on the entire system, d 2 xc dt 2 = 0 follows.Thus, Substituting Eq. (3) into Eq.(1), we have We then nondimensionalize Eq. (4).By introducing a small dimensionless parameter ε and the following quantities: and renaming τ → t and x → x, we transform Eq. ( 4) into the following dimensionless system: Here, we assume x = O(1), ẋ = O(1), α = O(1) and g(x, ẋ) = O(1).We show the typical dynamics of Eq. ( 6) in Fig. 1 (b).

Analysis
We perform the averaging approximation to the system (6).We rewrite Eq. ( 6) as We transform the variables x(t) and y(t) into the new variables r(t) with r(t) ≥ 0, θ(t) , and ϕ(t) that satisfy where Then, Eq. ( 7) is transformed into ṙ = −ε sin ϕ (αr sin ϕ + g(r cos ϕ, −r sin ϕ)) , (10a) See Supplementary Information for the derivation of Eq. (10).Since r(t) = O(1) follows from the assumptions x(t) = O(1) and y(t) = O(1), Eq. (10) suggests that the time-scale of r(t) and θ(t) are O(ε −1 ) and thus much larger than the time-scale of metronome's oscillation period 2π (i.e., the time-scale of x(t) in Eq. (8a)).Therefore, we can safely replace ṙ and θ with their time average over 2π.Namely, we approximate the right-hand sides of Eq. ( 10) with the time average as below: where Note that we regard r and θ as constants when we calculate the integrals in Eqs.(11a) and (11b): this is because the change of these variables during the integral interval [0, 2π] is O(ε) and thus can be negligible in the first-order approximation.This approximation method is widely known as averaging method [31,32] (or Krylov-Bogoliubov averaging method [33]), and the above discussion is mathematically justified using near identity transformation [34,35].

Features of escapement mechanism
We first discuss appropriate functions to model the escapement mechanism.Considering a real metronome, it is natural to assume that g(x, y) takes non-zero values if and only if x and y have the same sign.This is because, in the real metronome, the repulsive force works only when the pendulum position is the right from the center and moves to the right, or the pendulum position is the left from the center and moves to the left .Thus, we set g(x, y) to match this assumption.Below, we describe the three cases where g(x, y) is the piecewise linear function, the rational function with numerator of degree 3 and denominator of degree 4, and the 5th order polynomial, respectively.In Supplementary Information, we describe another case where g(x, y) is the rational function with a linear numerator and quadratic denominator.

Model (i)
We use the following piecewise linear function, which has been previously used [18,25,29], as the escapement mechanism of the metronome: where x 1 , x 2 are the positive constants that satisfies x 1 < x 2 .This is one of the simplest models of the real escapement mechanism.The shape of the function g with y > 0 is shown in Fig. 2 (a).In this case, the averaging equation ( 11) is calculated as follows: (14c) See Supplementary Information for the derivation of Eqs. ( 14) and ( 15).We consider the dynamics of Eq. ( 14), which is closed with respect to r.According to the stability analysis, described in Supplementary Information, we find the following: • Eq. ( 14) has the trivial fixed point r = 0, which is always stable regardless of the value of α.
• The saddle-node bifurcation occurs at • If α < α SN , Eq. ( 14) has the nontrivial two fixed points, one of which is stable and the other is unstable, whose values are, respectively, given by Figure 2 (b) presents the typical flows described by Eq. ( 14), which shows that the saddle-node bifurcation occurs as α changes.The bifurcation diagram for r is shown in Fig. 2 (c).The green cross marks in Fig. 2 (c) show the numerically obtained equilibrium states of Eq. ( 6) when we increase α, whereas the red dots are those when we decrease α.These results are in good agreement with the black lines, or the analytically obtained bifurcation diagram (i.e., Eqs.(17b) and ( 18)), which validates our analytical method with averaging approximation.Figure 2 (c) indicates that for a sufficiently small alpha, our system is bistable: both the oscillatory state (r = r stable ) and the resting state (r = 0) are stable.In addition, as α increases, the stable limit cycle disappears by the saddle-node bifurcation and the oscillatory state transits to the resting state.Therefore, our piecewise linear model reproduces the bistability and the oscillation quenching observed in the real metronome on a movable platform.However, in exchange for the simplicity of the model, the flow becomes non-smooth, which will actually hamper the analysis of synchronization for two metronomes.Motivated by this fact, we also consider smooth functions as g(x, y) below.

Model (ii)
We consider g(x, y) := The shape of the function g with y > 0 is shown in Fig. 2 (d).In this case, the averaging equations ( 11) are calculated as below: See Supplementary Information for the derivation of Eqs.(20a) and (20b).
Here, we consider the dynamics of Eq. (20a).Obviously, the fixed point of Eq. (20a) satisfies the following transcendental equation: where ) is a sigmoid function of R, we see that Eq. ( 21) has the unique solution (R = 0) if α > α SN and three solutions if α < α SN (See Fig. S3 in Supplementary Information).Here, α SN is given as where R * > 0 satisfies Thus, as there exists a one-to-one relationship between r ≥ 0 and R, we find the following: • Eq. (20a) has the trivial fixed point r = 0, which is always stable regardless of the value of α.
• If α < α SN , Eq. (20a) has non-trivial two fixed points, one of which is stable (r = r stable ) and the other is unstable(r = r unstable ), whose values are the solutions of Eq. ( 21) with r stable > r unstable .
Figure 2 (e) shows the typical flows of Eq. (20a) before and after the bifurcation point.The bifurcation diagram for r is shown in Fig. 2 (f).The numerically obtained equilibrium states of Eq. ( 6) (green cross marks and red dots) are in good the analytically obtained bifurcation diagram (black lines), which validates our analysis with averaging approximation.Figure 2 (f) indicates that our model with the rational function (19) reproduces the bistability and the oscillation quenching observed in the real metronome on a movable platform.However, although the flow becomes smooth in this model, we expect that the analysis of the two-oscillator system will be difficult.This is because (i) the averaged system (20a) includes the log function and (ii) the amplitude of the stable limit cycle, which is the solution of Eq. ( 21), cannot be analytically obtained.Thus, we next use the polynomial function of order 5, which partly imitates Model (i) (Eq.( 13)) and (ii) (Eq.( 19)).

Model (iii)
We consider where a, b ∈ R are positive constants.The shape of the function g with y > 0 is shown in Fig. 2 (g).In this case, the averaging equations (11) are calculated as See Supplementary Information for the derivation of Eqs.(25a) and (25b).
We discuss the dynamics of Eq. (25a).Obviously, the fixed point of Eq. (25a) satisfies 6παr − 3ar 3 + 2br 5 = 0, which is solved as Thus, we find the following: • Eq. (25a) has the trivial fixed point r = 0, which is always stable regardless of the value of α.
• The saddle-node bifurcation occurs at • If α < α SN , Eq. (25a) has non-trivial two fixed points, one of which is stable and the other is unstable, whose values are, respectively, given by Figure 2 (h) presents the typical flows of Eq. (25a) before and after the bifurcation point.The bifurcation diagram for r is shown in Fig. 2 (i).As with the previous cases, the numerically obtained equilibrium states of Eq. ( 6) (green cross marks and red dots) are in good the analytically obtained bifurcation diagram (black lines), which validates our analysis with averaging approximation.Figure 2 (i) indicates that model (iii) reproduces the bistability and the oscillation quenching observed in the real metronome on a movable platform.Moreover, the averaged equations (25a) and (25b) are expressed by the polynomial of r, which we expect will facilitate the analysis of the two-oscillator system.

Model
We analyze the synchronization and oscillation quenching of two metronomes on a movable platform, using model (iii) as the escapement mechanism.The model for the two-oscillator system is shown in Fig. 3 (a).Here, a point mass m i (i = 1 or 2) is connected by two springs with spring constant k i /2 and natural length l i to a platform of mass M .The platform has one degree of freedom and moves freely.The variables X(t) and x i (t) are the positions of the platform (more precisely, the position of the platform's center plate that separates the two metronomes) relative to the floor and the mass relative to the platform, respectively.The origin of x i coordinate is set to the position of the point mass m i in the equilibrium (i.e., ẋ1 = ẋ2 = Ẋ = 0).
By assuming the same situation as in the single metronome model, we obtain the following motion equations: for i = 1, 2, where γ i is the damping coefficient, δ i represents the magnitude of the escapement mechanism, and f i (x, ẋ) is a dimensionless function that shows the nature of the escapement.We denote the position of the center of mass as x c .Then, Since we assume that there is no external force acting on the whole system, d 2 dt 2 x c = 0 holds, which implies that Substituting Eq. ( 32) into Eq.( 30), we have We introduce a small dimensionless parameter ε and the following quantities: By renaming τ → t and xi → x i , we transform Eq. ( 33) into the following dimensionless system: where we assume that , and g i (x i , ẋi ) = O(1).
For simplicity, we consider the case where two metronomes are identical.Namely, we set Removing ẍ1 and ẍ2 terms from right-hand sides of Eq. ( 36), we can rewrite Eq. ( 36) as We show the typical dynamics of Eq. (37) in Figs. 3 (b) and (c).By neglecting the O(ε 2 ) term from the right-hand side of Eq. (37), we finally obtain

Averaging approximation
We analyze Eq. ( 38) with averaging approximation.We rewrite Eq. (38) as for i = 1, 2. We transform the variables x i (t) and y i (t) into the new variables r i (t) with r i (t) ≥ 0, θ i (t) , and ϕ i (t) that satisfy where Then, Eq. ( 39) is transformed into See Supplementary Information for the derivation of Eq. ( 42).Since r i (t) = O(1) follows from the assumptions x i (t) = O(1) and y i (t) = O(1), Eq. ( 42) suggests that the timescale of r i (t) and θ i (t) are O(ε −1 ) and thus much larger than the time-scale of metronome's oscillation period 2π (i.e., the time-scale of x i (t) in Eq. (40a)).Therefore, we can safely replace ṙi and θi with with their time average over 2π, respectively, as below: where ḡ1,2 (•) are given by Eq. ( 12).Here, based on the same arguments as in the one-oscillator system, we regard r i and θ i as constants when we calculate the integrals in Eqs.(43).We also use several integral formulae summarized in Supplementary Information for the derivation of Eq. ( 43).Hereafter, we consider the approximately equal sign (≃) in Eq. ( 43) as the equal sign (=).By introducing Eq. ( 43) is rewritten as In the analysis of the two-oscillator system, we model g(x, y) with Eq. ( 24) because this model reproduces the bistability and the oscillation quenching phenomenon for the one-oscillator system.Then, we finally obtain

Analysis of synchronized states
We perform the linear stability analysis for the averaged system (46).Eq. ( 46) has two fixed points and where Equations ( 47) and ( 48) correspond to the in-phase and anti-phase synchronization states, respectively.Note that the condition is necessary for the existence of these fixed points: if β ≥ β SN , these fixed points disappear with their counterparts (i.e., the fixed points , π that pair with Eq. (47) and Eq.(48), respectively) by the suddle-node bifurcation.By performing the linear stability analysis under this parameter condition (50), we find the following: 1.The fixed point (48), corresponding to the anti-phase synchrony, is always asymptotically stable.
2. The fixed point (47), corresponding to the in-phase synchrony, is asymptotically stable if and only if The details of the stability analysis are summarized in Supplementary Information.
Since Eq. ( 53) is the equation of damped oscillation, we see that the fixed points z r = żr = 0 and z c = żc = 0, which correspond to oscillation quenching, are asymptotically stable.Therefore, we expect that the averaged system (46) converges to the oscillation quenching state in the parameter region where the fixed points that correspond to the synchronous states (i.e., Eqs. ( 47) and ( 48)) disappear.
In fact, the parameter space in which oscillation quenching occurs is wider than the inequality (54) if the two metronomes move in in-phase synchronization.Here, we consider the case when two metronomes have nearly identical movements.By assuming x 1 = x 2 = x, the original motion equation ( 36) is transformed into By setting τ := t √ 1 + 2εµ and renaming τ → t, we obtain Note that g(x, ẋ) = g(x, ẋ√ 1 + 2εµ) since g is given by Eq. (24).In this case, we can remove the assumption that µ = O(1).In other words, as long as εµ = O(1), we can treat Eq. ( 56) as a weakly nonlinear oscillator and thus perform averaging approximation.
Referring to Eqs. (25a) and (25b), we see that the averaging approximation of Eq. ( 56) yields where r and θ are given by Eqs. ( 8) and ( 9).Considering the dynamics of Eq. (57a), we find that the saddle-node bifurcation occurs when which implies that oscillation quenching is observed when Note that the condition (59) is looser than the condition (54).Namely, there exist parameter regions in which the oscillation quenching occurs when two metronomes move in in-phase synchronization while anti-phase synchronization is stable, which has been observed in the previous experiments with pendulum clocks [17].Based on the above analyses, we depict the phase diagram, which is shown as the black lines in Fig. 4. The dash-dotted line shows the line β = β SN , which is the boundary of whether oscillation quenching occurs when the two metronomes move in anti-phase synchronization.In the same way, the solid line shows the line β = β SN_in (µ), which is the boundary of whether oscillation quenching occurs when the two metronomes move in in-phase synchronization.The dashed line is given by µ = µ c (β) in inequality (51).In other words, the stability of in-phase synchronization switches at the dashed line.
The solid and dash-dotted lines in Fig. 4 correspond to the saddle-node bifurcation.Since the Jacobian matrix at the fixed point (47) has a zero eigenvalue on the dashed line (see Supplementary Information), this line corresponds to either of the saddle-node, transcritical, or pitchfork bifurcation [36].According to the symmetry of Eq. (46) (i.e., Eq. ( 46) is invariant if we change r 1 → r 2 , r 2 → r 1 and ψ → −ψ) and the fact that this fixed point (47) does not disappear after the bifurcation, the saddle-node and transcritical bifurcations are unlikely.Since the stable fixed point that satisfies r 1 ̸ = r 2 emerges near the bifurcation point (Fig. 5 (c)), we consider that the dashed line in Fig. 4 corresponds to the supercritical pitchfork bifurcation.

Numerical simulation
We verify the averaging approximation by numerically integrating Eq. (37) and investigate the steady state for different values of β and µ, which are summarized in the colored regions in Fig. 4. In the numerical simulation, we fix a = 4, b = 1, ε = 0.01 and use the following 4 initial conditions: 1. Near in-phase synchronization (x 1 (0) = 5.01, x 2 (0) = 5, ẋ1,2 (0) = 0), 2. In-phase synchronization (x 1 (0) = x 2 (0) = 5, ẋ1,2 (0) = 0), 3. Near anti-phase synchronization (x 1 (0) = 5.01, x 2 (0) = −5, ẋ1,2 (0) = 0), Figure 4: The phase diagram of the system (36).The black lines represent the analytical boundary obtained by an averaging approximation.The dash-dotted and solid lines are given by β = β SN and β = β SN_in , respectively.Namely, the dash-dotted line shows the boundary of whether oscillation quenching occurs when the two metronomes move in anti-phase synchronization, whereas the solid line is the boundary of whether oscillation quenching occurs when the two metronomes move in in-phase synchronization.The stability of in-phase synchronization switches at the dashed line, which is given by µ = µ c .Note that oscillation quenching is asymptotically stable in all of the parameter space.The colored regions, including the white area, show the phase diagram obtained by the numerical simulation of Eq. (37) where g is given by Eq. ( 24).In the white area, oscillation quenching occurs if we use the initial condition that is close to either the in-phase or anti-phase synchronization.The anti-phase synchronization (APS) is stable in the blue area.The in-phase synchronization (IPS) is unstable in the green triangle mark area and stable in the red cross mark area.Note that, in the area painted in blue only, oscillation quenching occurs if we use the initial condition that is either exactly the in-phase synchronization or close to the in-phase synchronization, whereas the anti-phase synchronization is stable.We fix a = 4, b = 1, ε = 0.01 in the simulation.
The white area of Fig. 4 represents the parameter region in which oscillation quenching occurs for any of the four initial conditions.In the blue area, both the initial conditions 3 (near anti-phase synchronization) and 4 (anti-phase synchronization) converge to the anti-phase synchronization, which implies that the anti-phase synchronization is asymptotically stable in this area.In the green triangle mark area, the initial condition 2 (in-phase synchronization) converges the in-phase synchronization, whereas the initial condition 1 (near in-phase synchronization) does not, which means that the in-phase synchronization is the unstable equilibrium state of the system.In the red cross mark area, both the initial conditions 1 and 2 converge to the in-phase synchronization, meaning that the in-phase synchronization is asymptotically stable.Note that oscillation quenching occurs for the initial conditions 1 and 2 in the area painted in blue only.Namely, the in-phase synchronized steady state does not exist in this area.
In Fig. 4, we see that the analytically and numerically obtained phase diagrams are in good agreement, which confirms the validity of the averaging approximation in this study.However, the boundary for the stability of in-phase synchronization shows a slight difference between analysis (the dashed line) and numerical simulation (the boundary between the green triangle marks and the red cross marks).To investigate the cause of this gap, we numerically obtained µ c , the value of µ at which the stability of the in-phase synchronization switches, while changing ε and fixing β.The results are shown in Fig. S4 in Supplementary Information, which shows that µ c approaches the analytically obtained value (i.e., µ c (β) in inequality (51)) as ε decreases.Thus, we conclude that the magnitude of the small parameter ε causes the gap between the boundaries of the stability of in-phase synchronization in Fig. 4.
In the parameter region where in-phase synchronization is unstable (i.e., the green triangle mark area in Fig. 4), several equilibrium states are numerically observed when we use the initial condition 1.These results are shown in Fig. 5, in which we select several sets of parameters β, µ and plot the time series of x i (t) in Eq. (37) after a sufficiently long time.In our simulation results, either of the anti-phase synchronization (Fig. 5 (a)), the alternating increase and decrease of amplitudes which is similar to the beat phenomena in the acoustic wave (Fig. 5 (b)), or the near in-phase movement with slightly different amplitudes (Fig. 5 (c)) are observed.In particular, the dynamics in Figs. 5 (b,c) suggest that there exist limit cycles and other fixed points than Eqs.( 47) and (48) in the averaged system (46), which is an important research theme for future analysis.

Discussion
In this study, we investigate the dynamics of a metronome, which is an example of a bistable system with a stable limit cycle and a stable fixed point.In particular, we focus on the oscillation quenching phenomenon observed in metronomes on a movable platform.First, we construct a mathematical model for a single-oscillator system, ignoring the nonlinearity of the pendulum structure and the damping of the platform.By performing the averaging approximation, we consider several functions to represent the escapement mechanism and find that these functions are appropriate for simulating the real metronome's bistability and oscillation quenching.We conclude that the fifth-order polynomial is the most suitable for both the ease of analysis and reproducibility of the real metronome.Subsequently, we expand our model to include a two-oscillator system.By assuming that the mass ratio of the metronome to the platform, the escapement mechanism, and the damping force are sufficiently small, we perform an averaging approximation and reduce the system to a three-dimensional dynamical system of the amplitudes r 1,2 and phase difference ψ.By performing a stability analysis on the averaged system, we obtain the phase diagram for the in-phase synchronization, anti-phase synchronization, and oscillation quenching.We also numerically integrate the original motion equations and confirm agreement between the numerically and analytically obtained phase diagrams.
Our mathematical model appears to oversimplify real metronomes on a movable platform.However, the motion equations of our model reproduce behaviors similar to those observed experimentally, such as in-phase synchronization, anti-phase synchronization, and oscillation quenching.Thus, we believe that our simplification does not impair the essence of reality.
The analytical method used in this study (i.e., averaging approximation) has been used in previous studies such as those by Pantaleone [16] and Goldsztein et al. [26,30] However, oscillation quenching was difficult to analyze in the models adopted in these studies.This is because the former study [16] used the van der Pol-type function to describe the escapement mechanism, which made the resting state unstable, and the latter studies [26,30] used the delta function, which restricted the detailed analysis of oscillation quenching.We consider several other functions (i.e., Model (i)-(iii)) for the escapement mechanism and show that the averaged system becomes bistable.In particular, for a fifth-order polynomial (Model (iii)), we can explicitly derive the amplitude of the stable limit cycle (Eq.( 28)).This is expected to ease the analysis.Thus, we adopt Model (iii) to describe the escapement of the two-oscillator system and successfully obtain the analytical phase diagram including the area where oscillation quenching occurs.
When considering a real metronome, an unimodal function, wherein the value is zero at the origin and infinity, seems appropriate to model the escapement mechanism.Therefore, Model (ii) is consistent with reality, although the analysis becomes difficult.To facilitate the analysis, we choose to model the escapement mechanism with a 5th-order polynomial (Model (iii)).This model captures the dynamics of the real escapement mechanism near the origin, i.e., ax 3 − bx 5 with x > 0 equals zero at x = 0 and reaches its maximum at x = a/b.However, it diverges negatively at infinity, which means that a large unrealistic negative restoring force acts at a location far from the origin.In general, an infinite degree is required to expand a function that converges to zero at infinity into a power series.Thus, it is natural that we cannot correctly reproduce the behavior of the real escapement mechanism at infinity when we model it using a polynomial of a finite degree.
Our study shows that oscillation quenching occurs by a saddle-node bifurcation when the mass ratio µ increases.This is evident from Eq. ( 5) in the case of one metronome (note that the bifurcation parameter α increases as µ increases) and Fig. 4 in the case of two metronomes.As µ corresponds to the magnitude of the feedback that the metronome receives from the platform, our findings indicate that the oscillation can be stopped by the feedback resulting from the motion of the oscillator.
The occurrence of oscillation quenching by a saddle-node bifurcation was already reported in a previous study [26] for the case of a metronome on a fixed platform.However, our study reveals that the same bifurcation is observed even for a metronome on a movable platform.Further, we analytically find that oscillation quenching occurs when µ increases if two metronomes on a movable platform move almost in-phase synchronization, which has not been shown in previous studies [26,30].In previous experiments with two pendulum clocks suspended from a common plate [17], it was observed that, when two clocks start from the initial condition close to in-phase synchronization, oscillation quenching occurs as the mass ratio between the pendulum clock and the entire system increases.This observation agrees with our results, which suggests that the simplification in our modeling and the analysis with the averaging approximation are valid.
Regarding the state-transition method from a limit cycle to a stable fixed point, our study suggests that an oscillation can be stopped by increasing the feedback resulting from its motion (i.e., the feedback from the platform).In our model, the magnitude of the feedback depends on the mass ratio µ; thus, abnormal oscillations can be controlled by increasing µ instead of directly intervening in the metronome.Figure 4 also indicates that any perturbation that switches from anti-phase to in-phase synchronization can evoke oscillation quenching in certain parameter regions (that is, the area highlighted in blue only in Fig. 4).
There are several open questions in this study.(1) In the analysis of the averaged system (46), we do not clarify the existence and stability of equilibrium states other than the two fixed points (r * , r * , 0) and (r * , r * , π).However, the steady state in which the amplitude of each metronome varies periodically is shown in Fig. 5 (b), suggesting that the system (46) has stable limit cycles.Furthermore, previous experiments revealed a phenomenon called "metronome suppression" [26], in which one metronome oscillates with a larger amplitude than the other.This observation can be analyzed by performing a stability analysis of a fixed point that satisfies r 1 ̸ = r 2 , as shown numerically in Fig. 5 (c).Thus, it is important to investigate the stability of the steady states in the system (46) other than synchronization and oscillation quenching.(2) When modeling a metronome, we neglect the nonlinearity caused by the pendulum structure.There are two reasons for this simplification: 1) the analysis becomes easier, and 2) such nonlinearity does not change the dynamics of amplitude r after the averaging approximation.The second reason is explained as follows.
In the model that considers the weak nonlinearity of the pendulum structure [26,30], term εσx 3 with a real constant σ is added to the motion equation as a result of the Taylor expansion of sin x to the third-order term.Namely, Eq. ( 7) would be ẋ = y, (60a) which implies that ṙ = −ε sin ϕ αr sin ϕ + σr 3 cos 3 ϕ + g(r cos ϕ, −r sin ϕ) , Then, the averaging approximation in Eq. ( 61) yields Comparing Eq. ( 11) with Eq. ( 62), we see that the dynamics of r are the same, whereas the dynamics of θ change.
Because this study mainly addresses the state transition from a limit cycle to a stable fixed point and thus focuses on r dynamics, we neglect the nonlinearity of the metronome, which does not change the dynamics of r.However, as θ dynamics are related to the phase of the oscillator, the results of the stability analysis for the two-oscillator system would change if we adopt Eq. ( 60) as the metronome's motion equation.In a future study, we plan to perform the same analysis for a model that considers the nonlinearity of the pendulum structure.We investigate the dynamics of metronomes on a movable platform using an averaging approximation and numerical simulation.To facilitate the analysis, we ignore the nonlinearity caused by the pendulum structure of the metronome and model the escapement mechanism using a fifth-order polynomial.Finally, we obtain a phase diagram and find that oscillation quenching occurs when the mass ratio between the metronome and platform increases, which agrees with previous experimental results [17].We believe that our simple model will contribute to future analyses of other dynamics, such as clustering [18], the chimera states [21,22], and chaotic dynamics [28].In the derivation of Eq. (S30), we use the relationships cos 2 ϕ = 1 1+u 2 and dϕ = du 1+u 2 .According to Eqs. (S29) and (S30), we see that averaged equations (i.e., Eq. ( 11)) are calculated as We consider the dynamics of Eq. (S31a).Obviously, the fixed point of Eq. (S31a) satisfies the following transcendental equation: where R := r 2 ≥ 0. Since log(1 + R) is a concave function of R, we see that Eq. (S32) has the unique solution (R = 0) if α ≥ 1/π and two solutions if α < 1/π (Fig. S1).Thus, as there exists a one-to-one relationship between r ≥ 0 and R, we find the following: • If α > 1/π, Eq. (S31a) has a stable fixed point r = 0.
Figures 5 (a) and (b) show the typical flows of Eq. (S31a) before and after the bifurcation point.We see that, as α decreases, the trivial fixed point r = 0 changes to unstable and non-trivial fixed point emerges due to the transcritical bifurcation.The bifurcation diagram for r is shown in Fig. 5 (c).The green cross marks show the numerically obtained equilibrium states of Eq. ( 6) when we increase α, which agree with the analytically obtained bifurcation diagram (black lines).However, in contrast to Model (i), (ii), and (iii) in the main article, this model is inappropriate to simulate the dynamics of metronome because it cannot simulate the bistability of the resting state and the oscillatory state.We consider that this problem arises because the escapement mechanism expressed as Eq.(S27) is not sufficiently weak near the resting state.In other words, Eq. (S27) implies that g(x, y) ≃ x in the vicinity of x = 0.These are the reasons why we use in the main article the rational function with numerator of degree 3 and denominator of degree 4 such that g(x, y) ≃ 0 holds in the vicinity of x = 0 by a first-order approximation.6) where g is given by Eq. (S27).The solid and dashed lines represent the analytically obtained stable fixed point (the positive solution of Eq. (S32) if α < 1/π and r * = 0 if α ≥ 1/π) and the unstable fixed point (r * = 0 if α < π), respectively.The green cross marks show the equilibrium solutions obtained by numerical simulation of Eq. ( 6).We first use α = 0.02 and then increase α by 0.02 until α reaches to 0.5.The initial condition for the first simulation is x(0) = 2.0, ẋ(0) = 0.For the following simulations, we use the equilibrium state of the previous simulation as the initial condition.We fix ε = 0.01 for all the simulations in panel (c).

Figure S1 :
Figure S1: This figure shows the solutions of the transcendental equation (S32).The black and blue lines represent παR and log(1 + R), respectively.Both the black dot and black circle show the solutions of Eq. (S32).Note that the black dot corresponds to the stable fixed point of Eq. (S31a), while the black circle corresponds to the unstable fixed point.(a) If α equals to or is larger than 1/π, which is the value when the black line is tangent to the blue curve, Eq. (S32) has the unique solution R = 0. We set α = 0.7 to depict this panel.(b) If α < 1/π, Eq. (S32) has another solution other than R = 0. We set α = 0.15 to depict this panel.

Figure
Figure S2: (a, b) The typical flows of Eq. (S31a), in which we plot ṙ as a function of r.The black dot and black circle in each panel represent the stable and unstable fixed points, respectively.As for the parameter values, we set ε = 0.01, α = 0.7 in panel (a), and α = 0.15 in panel (b).(c) The bifurcation diagram for r obtained by both the averaging approximation and the numerical simulation of Eq. (6) where g is given by Eq. (S27).The solid and dashed lines represent the analytically obtained stable fixed point (the positive solution of Eq. (S32) if α < 1/π and r * = 0 if α ≥ 1/π) and the unstable fixed point (r * = 0 if α < π), respectively.The green cross marks show the equilibrium solutions obtained by numerical simulation of Eq. (6).We first use α = 0.02 and then increase α by 0.02 until α reaches to 0.5.The initial condition for the first simulation is x(0) = 2.0, ẋ(0) = 0.For the following simulations, we use the equilibrium state of the previous simulation as the initial condition.We fix ε = 0.01 for all the simulations in panel (c).