Synchronization by Magnetostriction

We show how to utilize magnetostriction to synchronize two mechanical vibration modes in a cavity magnomechanical system. The dispersive magnetostrictive interaction provides necessary nonlinearity required for achieving synchronization. Strong phase correlation between two mechanical oscillators can be established, leading to the synchronization robust against thermal noise. We develop a theoretical framework to analyze the synchronization by solving the constraint conditions of steady-state limit cycles. We determine that the strong cavity-magnon linear coupling can enhance and regulate the synchronization, which offers a new path to modulate synchronization. The work reveals a new mechanism for achieving and modulating synchronization and indicates that cavity magnomechanical systems can be an ideal platform to explore rich synchronization phenomena.

Here, we show that the recently developed cavity magnomechanical (CMM) system [33][34][35][36] can exactly be such a candidate system.In the CMM system, magnons, quanta of collective spin excitations, in a ferrimagnetic yttrium-iron-garnet (YIG) sphere couple to vibration phonons via the magnetostrictive interaction, which is a dispersive interaction [37,38] and thus provides necessary nonlinearity for achieving synchronization in the system.Such nonlinearity also plays an essential role in preparing macroscopic quantum states [34,[39][40][41][42] and designing novel quantum technologies [43][44][45][46][47][48][49][50][51].In addition, magnons further couple to microwave cavity photons via the magnetic-dipole interaction.Due to the high spin density of YIG, the strong cavity-magnon coupling can be easily achieved, leading to cavity polaritons [52][53][54].Such a coupling is adjustable by changing the position of the YIG sphere in the microwave cavity.The intrinsic nonlinearity and tunable strong coupling of the CMM system make it an ideal platform to explore synchronization.
Specifically, we show that it is possible to achieve robust synchronization of two mechanical vibration modes protected by strong phase correlation under feasible parameters even at room temperature.The synchronization in the CMM system can be analytically decomposed by mapping the constraint conditions of steadystate limit cycles into the parameter space, which provides us a simple way to understand the complicated dynamics of the synchronization.We find that the strong cavity-magnon coupling provides a new degree of freedom, which plays an important and active role in enhancing and modulating the synchronization.This represents a new path to the modulation of synchronization and fundamentally differs from the synchronization mechanism in other systems, e.g., optomechanical systems [24][25][26][27][28][29][30][31].

II. THE MODEL
We consider a typical cavity magnomechanical system [33][34][35][36], as depicted in Fig. 1.It consists of a microwave cavity and a macroscopic YIG sphere placed inside the cavity, which supports a magnon (spin wave) mode and a series of mechanical vibration modes, among which we focus on two mechanical modes and study the synchronization between them.The Hamiltonian of the system reads: where â, m and bj (ω a , ω m , and ω j ) are the annihilation operators (resonance frequencies) of the cavity, magnon and j-th mechanical modes, respectively, satisfying [ Ô, Ô † ] = 1 ( Ô = â, m, bj ).The magnon frequency can be adjusted by altering the bias magnetic field H 0 via ω m = γ 0 H 0 , with the gyromagnetic ratio γ 0 /2π = 28 GHz/T.g j denotes the bare coupling rate between the magnon and the j-th mechanical mode and g ma is the cavity-magnon coupling rate, which can be (much) stronger than the cavity and magnon dissipation rates κ a and κ m [52][53][54].To enhance the magnetostrictive interaction, the magnon mode is driven by a microwave field with frequency ω 0 and amplitude B 0 , and the corresponding Rabi frequency is Ω = ( √ 5/4)γ 0 √ N B 0 [34], where N = ρV is the total number of spins, ρ = 4.22×10 27 m −3 is the spin density of YIG, and V is the volume of the sphere.
In the frame rotating at the driving frequency ω 0 , and by adding dissipative and input noise terms, we obtain the following quantum Langevin equations (QLEs): where ∆ a = ω a − ω 0 and ∆ m = ω m − ω 0 .κ a , κ m and γ j (â in , min and bin j ) are the decay rates (input noise operators) of the cavity, magnon and j-th mechanical modes, respectively.The input noises are assumed Gaussian and white noises, of which the correlation functions are ⟨ Ôin (t) Ôin † (t ′ )⟩ = ( NO + 1)δ(t − t ′ ) and ⟨ Ôin † (t) Ôin (t ′ )⟩ = NO δ(t − t ′ ), with Ô = â, m, bj , and NO = [exp(ℏω O /k B T ) − 1] −1 (O = a, m, j) being the mean thermal excitation number of the corresponding mode, k B the Boltzmann constant and T the bath temperature.

III. PHASE NOISE ANALYSIS
To study synchronization at a finite temperature, thermal noises of the system must be included, as the mean thermal occupation NO ≫ 1 at a high temperature, e.g., room temperature.We therefore apply stochastic Langevin equations (operators Ô are replaced with complex variables O) [31,55] to describe the system dynamics and simulate them numerically up to the long-time limit.The stochastic Langevin equations associated with Eq.
(2) are given by [56,57]: The operators â, m and bj in the QLEs are replaced with c-number complex variables a, m and b j , and the input noise operators are replaced with classical complex random noises with modified correlation functions: because the c-numbers lose the commutation relation.
After repeatedly calculating the stochastic Langevin equations N times (N should be large), the disturbance of the noises to the synchronization can be characterized by the phase-space probability distribution of the considered phases and the phase correlation (particularly, the phase difference), which are defined as: where N θ j(−) (θ) is the number of the results satisfying and the superscript i denotes the i-th stochastic trajectory in the simulation, θ j is the phase of the slowly varying complex amplitude of the j-th oscillator (see Appendix A).The ensembleaveraged quantities and their quantum fluctuations can be estimated by We simulate the stochastic Langevin equations for a time interval γ 1 t = 19, i.e., the same time interval used to obtain the phase diagram in the following section, and repeat the calculations 10 4 times.The noise analysis corresponds to the following four cases, and each case corresponds to a specific stable phase difference (or a set of parameters).The state of the system can be well described by the four cases, i.e., the two oscillators are: i) synchronized; ii) anti-synchronized; iii) at the critical point between synchronization and anti-synchronization; iv) antisynchronized with low energy (see the phase diagram for the detailed description).The corresponding parameters are respectively: i) log 10 (Ω/Ω 0 ) = −0.4,g ma /ω 1 = 0.8; ii) log 10 (Ω/Ω 0 ) = −0.8,g ma /ω 1 = 0.8; iii) log 10 (Ω/Ω 0 ) = −0.5168,g ma /ω 1 = 0.7 and iv) log 10 (Ω/Ω 0 ) = −1, g ma /ω 1 = 0.5.We use experimentally feasible parameters in getting Fig. 2 [33][34][35][36]: × 50 mHz and Ω 0 = 7 × 10 14 Hz (corresponding to the drive magnetic field B 0 = 3.8 × 10 −5 T and power P = 8.3 mW [34]).
The solid lines in Fig. 2(a-d) show the phase probability distribution of the mechanical oscillator 1.We find that when the effects of thermal noises are taken into account, the phase of the mechanical mode is progressively diffused with phase variance ⟨δθ 2  1 ⟩ = 4.8 × 10 −4 corresponding to case i, ⟨δθ 2  1 ⟩ = 2.2 × 10 −2 in case ii, ⟨δθ 2  1 ⟩ = 3.5 × 10 −4 in case iii, and ⟨δθ 2 1 ⟩ = 1.3 × 10 −2 in case iv.The circles show that the distribution of the phase difference is obviously narrowed, which can be described quantitatively by defining a compression ratio: η equals to 1 for two uncorrelated oscillators due to ⟩ in this case.We obtain η = 1.19 × 10 −4 , 1.16 × 10 −4 , 1.75 × 10 −5 , and 5.34 × 10 −4 for the four cases, which reveal the emergence of strong phase correlations between the two oscillators.In Fig. 2(e), we show the 10 4 random results of θ 1 and θ − for the case i.The horizontal axis represent the i-th random result (i = 1, 2, 3, ..., 10 4 ).We can determine that synchronization in the CMM system is robust to thermal noises with the protection of the phase correlation.In addition, we emphasize that the mixture of two (or more) limit cycles never emerges in our simulation results (the considered parameters are not limited to the four points shown), which means that the adjacent attractors are too distant in the phase space.Therefore, the noises can not support the mechanical modes jumping from one stable solution to another.
These results prove that a strong phase correlation is established between two vibrational modes, leading to the synchronization between them very robust against thermal noises.It thus suggests that synchronization can be well studied in the noiseless case, which is considered in the following sections.

IV. SYNCHRONIZATION PHASE TRANSITION
For the system under study, it does not exhibit chaotic behavior, which occurs only under an extremely strong driving field.Therefore, we can characterize synchronization in terms of the phase difference [55]: , where θ j is the phase of the j-th mechanical oscillator, and P = −1, 0 and 1 correspond to the π-phase, non-and zero-phase synchronization, respectively.
The synchronization phase diagram is shown in Fig. 3 by taking time average of P(t) for a sufficiently long time interval ensuring stable values [58], i.e., t ∈ [9/γ 1 , 19/γ 1 ].Clearly, the π-and zero-phase synchronizations are present in a large parameter regime, and a prominent phase transition of synchronization appears for a sufficiently strong (small) coupling g ma (mechanical frequency difference ∆ω).The synchronization is quite robust and even close to the phase-transition boundary, thermal noises and random initial conditions have negligible impact on the synchronization.Our theoretical analysis indicates that the system is bistable or even multistable.However, Fig. 3 displays only one of the steady states of the limit cycles.This is because we do not traverse all possible initial states, but only assume the system is in a thermal initial state (see Sec. V.).It is worth noting that, the cavity-magnon coupling g ma (a controllable parameter that can be tuned in a wide range) can effectively modulate the synchronization phase, c.f. Fig. 3(a).For a moderate value of g ma , the two mechanical modes can be π-phase synchronized even for a large frequency difference ∆ω > 0.1ω 1 and at a low driving power of 83 µW, as shown in Fig. 3(b).This reveals a distinct advantage of the CMM system for re- alizing and modulating synchronization compared with other systems.

V. MECHANISM OF SYNCHRONIZATION AND MULTISTABLE SYNCHRONIZED LIMIT CYCLES
In order to explain the complex limit cycle dynamics of the system, we take the slowly varying amplitude (SVA) equations approach [31,59], and study the long-time dynamics of the two mechanical oscillators in the frame rotating at a fast reference frequency ω, i.e., b j (t) = β s j + B j e −iωt (see appendix A), where β s j are the equilibrium positions, B j are slowly varying complex amplitudes, and ω = (ω 1 + ω 2 )/2.Substituting it into the noiseless Langevin equations, we obtain the formal solutions of the cavity and magnon modes, which can be expressed as the sum of a series of sidebands at the frequencies of nω, with n being an integer.Substituting these solutions into the equations of the oscillators, we obtain the following amplitude equations (see appendix A for the detailed derivation): where g = g 2 1 + g 2 2 , and the dimensionless function , with B = j g j B j .M n is the amplitude of the n-th mechanical sideband, which can be determined via iterative computation or other numerical methods.Equation ( 6) indicates that the backaction of the cavity-magnon system on the dynamics of the mechanical oscillators is fully manifested in the F -function, which renormalizes the frequencies and dissipations, and more importantly, provides an effective coupling between the two oscillators.By rewriting Eq. ( 6) in terms of the modulus I j and phase θ j of the complex amplitude B j = I j e iθj , we obtain the following Kuramoto-like equations (KLEs) [3]: where Γ j = g 2 j F i /g − γ j , F = F r + iF i and the phase difference θ − = θ 1 − θ 2 .The KLEs provide us a powerful tool to describe the self-sustained mechanical oscillations.They can be further simplified to stationary equations by setting the derivatives to zero, which describe two synchronized oscillators as two amplitude-stable limit cycles will be of a constant phase difference.The stationary Ffunction F s can be written as a function of the stationary modulus I s j , i.e., where we define R s = I s 1 /I s 2 for convenience.Substituting Eqs.(8) into Eq.( 7) yields the following state constraint equation on R s and θ s − : which determines the behavior of the stationary synchronization of the two oscillators.Note that the above constraint equation depends only on the mechanical system but not on the cavity-magnon system [60].Hence, the constraint of the synchronization is essentially an intrinsic property of the two oscillators.The solutions of Eq. ( 9), manifested as the identical red lines in Fig. 5(a)-(f), are thus the necessary conditions for the synchronization, which are satisfied by all allowed synchronization states under the given parameters, while any other states outside the red lines are actually the unstable states of the limit cycles.In particular, the perfect zero-phase synchronization θ s − = 0 requires R s = g1γ2 g2γ1 .By contrast, the perfect π-phase synchronization is unattainable for the conventional parameters, as it requires R s = − g1γ2 g2γ1 .For a given θ s − , the solution of R s is symmetric, and θ s − has a single maximum at R s = g1 g2 γ2 γ1 , which yields an optimal P for the π-phase synchronization, i.e., . Apparently, ∆ω ≫ γ 1,2 is the basic condition for the occurrence of the π-phase synchronization.
Utilizing the analytical expression of Eq. ( 8), we now discuss in detail the mechanism of synchronization phase transition.We still study the synchronization dynamics in the long-time limit.As shown in Fig. 4, the steadystate range F s (red and purple lines obtained by solving Eqs.(8)(9)) and the dynamic range of F (lines L 1 -L 4 obtained by solving Eq. ( A9)) is plotted.The intersection point of F s and F indicates that the limit cycles have a steady-state solution.The points A, B, C and D (also marked in the inset) represent the value of the steadystate solution F s at the corresponding parameter of the phase diagram.These points are obtained by first simulating the noiseless Langevin equations with the thermal initial states for a sufficiently long time ensuring stable values (i.e., θ s − and R s ), then substituting them into the F s (θ s − , R s ) function.Clearly, our theoretical approach fits perfectly with the numerical simulation, and more importantly, it reveals the bistability of the limit cycles, i.e., the points A ′ -D ′ .Here, the points A ′ -D ′ are obtained by simulating the noiseless Langevin equations under appropriate initial conditions.To be specific, starting with the final state of A and using the parameters of B (or increasing the drive power), we can obtain a new phase B ′ .Similarly, starting with the final state of B and using the parameters of A (or lowering the drive power), we achieve a new phase A ′ .Therefore, the phase transition is essentially induced by the initial thermal states, which will be different if under different (appropriate) initial conditions.

VI. MODULATION OF SYNCHRONIZATION
According to the value of g ma , we classify the states of the mechanical oscillators (the asymptotic steady states marked in gray are excluded) in the phase diagram of  9), and each blue point corresponds to a pixel in Fig. 3(a).The gray points represent the non-synchronization states, of which the phase difference θ− are nonstationary, but their time average θ s − are around π/2.The parameters are the same as in Fig. 3(a).
Fig. 3(a) and plot their characteristic variables θ s − and R s in Fig. 5(a)-(f).The blue scatter points are "hitched" by the solution of the constraint equation, as excepted, but they do not completely smear the red lines, confirming that the constraint equation is only a necessary condition.As g ma increases, the blue points tend to distribute to both poles, implying that the system has a distinct feature of zero-or π-phase synchronization.This indicates that the synchronization can be enhanced and modulated by adjusting the cavity-magnon coupling rate.The synchronization properties beyond the constraint equation are reflected in the aforementioned F -function, which is entirely determined by the cavity-magnon system.This suggests dividing the whole system into two parts, as sketched in Fig. 6(a): the mechanical system that constrains the range of the legal synchronization states, corresponding to the steady-state modulation; and the cavity-magnon system that selects which kind of the synchronization state that can be finally obtained, corresponding to the nonlinear modulation.These two types of modulation are mutually independent, which greatly simplifies the procedures for synchronizing the oscillators to a given target phase.Specifically, the procedures are summarized as follows: i) Solving R s from the constraint equation Eq. ( 9) with a given θ s − ; ii) Substituting θ s − and R s into Eqs.(8) to obtain the conditions that F should fulfil, denoted as F s ; iii) Adjusting relevant parameters of the cavity-magnon system to satisfy F = F s , corresponding to the KLEs having steady-state solutions and thus the occurrence of the target synchronization state.
Due to the nonlinearity of the F -function, the last procedure is more easily realized by checking the intersection points after plotting F and F s in the parametric space, as shown in Fig. 6(b).As we are interested in prominent synchronization phenomena, the phase differ-   (c) Steady-state modulation by controlling the mechanical system (i.e., via controlling F s ).The red (purple) lines denote the zero-phase (π-phase) synchronization regime with γ1/2π = 100 Hz, γ2/2π = 150 Hz (solid lines), and γ1/2π = 100 Hz, γ2/2π = 60 Hz (dotted lines).Inset shows the impact of γ2 on the synchronization, where g0/2π = 60 mHz.In (b) and (c), we take Ω = 0.1Ω0 and gma = 0.1ω1.The other parameters are the same as in Fig. 3.
ence is restricted to a small range satisfying |P| > 0.995 and the nonlinear modulation is realized by controlling F via changing Ω and g ma .When the driving power and the coupling rate are small (green triangles), there will be no solutions.As the power increases, the curve is shifted along the F s r axis, leading to the appearance of multistable synchronized limit cycles (blue circles).Comparing with the blue circles (with Rabi frequency Ω 0 and cavity-magnon coupling rate 0.1ω 1 ), the orange dots (with Rabi frequency 0.3Ω 0 and cavity-magnon coupling rate ω 1 ) can produce both the zero-and π-phase synchronizations even for relatively small couplings, e.g., g 1 /2π = 16 mHz and g 2 /2π = 18 mHz (dashed lines), which can be easily achieved in the CMM experiments [33][34][35][36].This benefits from a new mechanism of synchronization: as g ma increases, the curve is rotated around the original point, which can sweep over a much wider area in the parametric space.Inset shows how the coupling g ma modulates the synchronization.Clearly, increasing the cavity-magnon coupling can significantly enhance the synchronization.In Fig. 6(c), we explore the steady-state modulation via controlling F s realized by altering γ j .The results indicate that the zero-phase synchronization can be enhanced by reducing the dissipation rates, as more clearly shown in the inset.

VII. CONCLUSIONS
We present a new mechanism of synchronizing mechanical oscillators in a CMM system exploiting the nonlinear magnetostriction.We find that a strong phase correlation can be established between two mechanical oscillators, leading to their synchronization which is robust against thermal noise.We also develop a theoretical framework to analyze the synchronization and determine the active role the cavity-magnon coupling plays in enhancing and modulating the synchronization.All of these indicate that the highly controllable and tunable CMM system can be a promising new platform for studying and modulating synchronization.The work can be extended straightforwardly to study synchronization between two or multi YIG spheres.It can also be applied to other systems that share a similar Hamiltonian as the CMM system, e.g., synchronizing two mechanical oscillators in exciton-optomechanics systems [61][62][63].Synchronized mechanical oscillators can be exploited to achieve the synchronization between two optical cavities, e.g., by means of an opto-magnomechanical configuration [38,64], and between two atomic ensembles by further coupling each cavity to an atomic ensemble [65].This provides possibility to distribute synchronization or quantum states in a complex quantum network [13].where κ a = κ in + κ ex is the total cavity decay rate, with κ in being the intrinsic cavity decay rate, and ∆ j = ω mj − ω 0 .The above QLEs can be well approximated by a set of coupled noiseless Langevin equations [57], considering the synchronization is very robust against thermal noise: ȧ = −(i∆ a + κ a )a − ig ma (m 1 + m 2 ) + Ω a , ṁj = −(i∆ j + κ j )m j − ig ma a − ig j m j (b * j + b j ), ḃj = −(iω j + γ j )b j − ig j |m j | 2 . (B3) In Fig. 7, we plot the synchronization phase diagrams by using Eq.(B3) and taking time average of the phase difference P(t) = cos[θ 1 (t) − θ 2 (t)] for a sufficiently long time interval ensuring stable values, i.e., t ∈ [9/γ 1 , 19/γ 1 ].In this system of two YIG spheres, the effective coupling between the two mechanical modes is much more indirect (via the mediation of two magnon modes and a common cavity), compared with the case studied in the main text.This makes it more difficult to synchronize two mechanical modes of two YIG spheres, reflected by the fact that the parameter regime for achieving synchronization is much smaller than the single-sphere case.To achieve zero-phase (π-phase) synchronization, a much stronger cavity-magnon coupling rate is needed, about g ma ≃ 2ω 1 .Such a strong coupling can, however, be easily obtained in cavity magnonic experiments, thanks to the high spin density of YIG.This indicates the advantage of the system: the achievable very strong cavity-magnon coupling can effectively enhance and modulate the synchronization of two mechanical modes of either one YIG sphere or two YIG spheres.

FIG. 1 .
FIG. 1. Schematic diagram of the cavity magnomechanical system used for achieving synchronization of two mechanical modes.It consists of a microwave cavity mode a, a magnon mode m, and two mechanical vibration modes b1,2 (with resonance frequencies of ω1,2).

FIG. 2 .
FIG. 2. (a), (b), (c) and (d): Phase difference probability distribution θ− (circle) and phase probability distribution of the oscillator 1 θ1 (triangle) obtained by 10 4 times calculations of the stochastic equations with different values of Ω and gma (from left to right and from top to bottom, corresponding to the cases i, ii, iii and iv, respectively).The solid and dashed lines represent the results that fit the corresponding original data with a Gaussian distribution.(e): 10 4 statistical results of θ− (blue) and θ1 (red) for the case i.

FIG. 3 .
FIG. 3. Synchronization phase diagram as a function of (a) Rabi frequency Ω and coupling gma; (b) Rabi frequency Ω and mechanical frequency difference ∆ω = ω2 − ω1.We take ∆ω = 0.01ω1 in (a), gma = 0.5ω1 in (b), and ∆a = ∆m ≃ −ω1 in both plots.The gray areas denote the stable regime when the QLEs only have asymptotic steady state solutions, and they can be determined by the Lyapunov stability criterion.The other parameters are the same as in Fig. 2.

FIG. 5 .
FIG. 5. Steady-state distribution of R s and θ s− with different values of gma.The red line is the solution of Eq. (9), and each blue point corresponds to a pixel in Fig.3(a).The gray points represent the non-synchronization states, of which the phase difference θ− are nonstationary, but their time average θ s − are around π/2.The parameters are the same as in Fig.3(a).

ACKNOWLEDGMENTS
This work has been supported by National Key Research and Development Program of China (Grant No. 2022YFA1405200) and National Natural Science Foundation of China (Grant Nos.92265202, 11704205, 12304389 and 12074206).We also acknowledge the support of the European Union Horizon 2020 Programme for Research and Innovation through Project No. 732894 (FET Proactive HOT) and the Scientific Research Foundation of NEU (Grant No. 01270021920501*115).