Magnomechanical backaction corrections due to coupling to higher order Walker modes and Kerr nonlinearities

The radiation pressure-like coupling between magnons and phonons in magnets can modify the phonon frequency (magnomechanical spring effect) and decay rate (magnomechanical decay) via dynamical backaction. Such effects have been recently observed by coupling the uniform magnon mode of a magnetic sphere (the Kittel mode) to a microwave cavity. In particular, the ability to evade backaction effects was demonstrated [C.A. Potts et al., arXiv:2211.13766 [quant-ph] (2022)], a requisite for applications such as magnomechanical based thermometry. However, deviations were observed from the predicted magnomechanical decay rate within the standard theoretical model. In this work, we account for these deviations by considering corrections due to (i) magnetic Kerr nonlinearities and (ii) the coupling of phonons to additional magnon modes. Provided that such additional modes couple weakly to the driven cavity, our model yields a correction proportional to the average Kittel magnon mode occupation. We focus our results on magnetic spheres, where we show that the magnetostatic Walker modes couple to the relevant mechanical modes as efficiently as the Kittel mode. Our model yields excellent agreement with the experimental data.

Magnons can also couple to other degrees of freedom, opening the opportunity of probing and manipulating these via their coupling to the hybridized magnonmicrowave polaritons [13]. In particular, magnetoelastic effects [25][26][27][28] couple the magnetization and the mechanical vibrations of a magnetic material, yielding an interaction between magnons and phonons [29]. Such magnomechanical coupling can be either resonant or parametric [30]. The resonant coupling is relevant for specific geometries where certain magnon modes are resonant with the elastic vibrations of the medium, for example, for magnetic spheres with radii ranging from ∼ 10 nm to ∼ 10 µm [30] and in magnetic films [31][32][33]. The second type of coupling, parametric coupling, is relevant for geometries in which the magnon frequency is far detuned from the phonon, such as for micrometer-sized magnetic spheres ical systems have so far focused on the coupling to a single magnon mode, the uniform precession of the magnetization called the Kittel mode. Nevertheless, a magnetic sphere supports a whole set of magnon modes called Walker modes [46,47] which can also couple to a given vibration mode, in principle even stronger than the Kittel mode. Weak inhomogeneities in the microwave field can drive such higher-order magnon modes, modifying the backaction effects. Furthermore, magnon nonlinearities due to crystalline anisotropy [48,49] can also affect dynamical backaction. Nonlinearities appear as selfand cross-Kerr terms in the magnomechanical Hamiltonian, inducing static frequency shift of the modes which have been reported in Refs. [50,51]. Additionally, the magnon Kerr nonlinearity can yield a bistable behavior of the magnons [52,53] and mechanics [51], and induce quantum phase transitions [54], among other phenomena. Nevertheless, no description of how such nonlinearities affect dynamical backaction beyond the static frequency shift has been reported. For instance, in optomechanical systems, effects of the cavity Kerr nonlinearity in an electromechanical system under sideband cooling have been recently observed [55].
In this work, we extend the theory of dynamical backaction in cavity magnomechanical systems to include Kerr nonlinearities and the coupling of a phonon mode to several magnon modes. We consider the framework depicted in Fig. 1, where a microwave cavity mode couples strongly to a magnon mode and weakly to a set of additional magnon modes. Those in turn exhibit nonlinearities and interact via a radiation pressurelike coupling to a single phonon mode. We derive the phonon self-energy, describing the frequency shift and the magnomechanical decay rate, generalizing previous results [34,36]. The overall effect of the coupling to the additional magnon modes is a correction proportional to the average number of Kittel magnons. We evaluate our model for the case of a magnetic sphere, computing numerically the coupling rates between the (magnetic) Walker modes and the mechanical mode probed in Ref. [43]. At low driving powers our model introduces corrections that agree well with the measured data in Ref. [43], explaining the observed shift in the magnomechanical decay rate. At higher driving powers, there are further deviations which are not captured by our model. These are, however, only relevant for driving frequencies detuned from the backaction evasion point.
This paper is organized as follows. In Sec. I we present a brief review of the description of dynamical backaction in cavity magnomechanics for the cases described in the literature, e.g., Refs. [29,34]. In Sec. II we include Kerr nonlinearities and the coupling to weakly driven additional magnon modes, and the phonon self-energy. Since those corrections depend on how strongly the additional magnon modes couple to the phonon mode, we specialize further our model to a magnetic sphere geometry as probed in Ref. [43]. In Sec. III, we briefly review the derivation of the magnomechanical coupling  1. (a) Schematic representation of the cavity magnomechanical system. A magnetic element is loaded in a microwave cavity. The magnon modes couple resonantly with a mode from the cavity and parametrically with the mechanical vibrations of the magnet. (b) Schematics of the model describing the cavity magnomechanical system with several magnon modes coupled to the same phonon mode [see Hamiltonian (7)]. Our model assumes that only one of the magnon modes, denoted bym, couples strongly to the cavity mode. The red arrow indicates the microwave drive.
following the literature [30], and in subsection III A we use the model to numerically evaluate the coupling between Walker modes with frequencies in a small range around the Kittel mode frequency and a relevant mechanical mode of a magnetic sphere. In Sec. IV, we compare our model to the experimental results presented in Ref. [43] and show that our generalized theory quantitatively accounts for the observed magnomechanical decay for a large range of parameters. Finally, in Sec. V we present our conclusions.

I. PHONON SELF-ENERGY AND DYNAMICAL BACKACTION EVASION
The dynamics of a cavity magnomechanical system consisting of a microwave mode (â with frequency ω a ) and a magnon mode (m with frequency ω m ) coupled parametrically to a phonon mode (b with frequency Ω b ) is described by the Hamiltonian [29] H = ω aâ †â + ω mm †m + Ω bb †b + g am âm † +â †m + g 0 (1) The magnon-microwave coupling rate g am is due to a magnetic dipole interaction between the ferromagnetic resonance of the material and the microwave cavity. The parametric magnon-phonon coupling, with the single magnon coupling rate g 0 mb , is due to magnetoelastic effects. The last term in Eq. (1) describes the coherent drive of the microwave cavity at a frequency ω d with an amplitude d = P/ ω d , where P is the drive power and κ e is the decay rate of the cavity to the external drive port.
In the weak magnomechanical coupling limit, both the phonon frequency Ω b and the decay rate Γ b are modified by the coupling to the driven magnons. The respective shifts are given by where Σ[ω] is the phonon self-energy, obtained by analyzing the linearized dynamics of the system [34,36], and reads From now on we refer to δΩ b as the magnomechanical frequency shift and to Γ mag as the magnomechanical decay rate. Here, g mb = g 0 mb m is the enhanced magnomechanical coupling rate, with | m | 2 the steady-state magnon population. The function Ξ[ω] is a modified Kittel mode susceptibility given by which depends on the magnon susceptibility The detuning between the microwave (magnon) mode and the drive is ∆ a(m) = ω d − ω a(m) . γ m is the magnon decay rate and κ the total microwave decay. The value and the sign of the magnomechanical decay rate depend on the drive frequency, which can tune scattering processes that upconvert or downconvert excitations in the system. Depending on the drive, it is possible to make one of such processes more efficient than the other, yielding a positive (cooling) or negative (amplification) magnomechanical decay rate, in a situation akin to what is found in standard optomechanical systems [35,56,57]. Different from optomechanics, in cavity magnomechanical systems magnons and microwaves hybridize, yielding the unique situation where the different scattering processes that contribute to dynamical backaction are associated with mechanical sidebands of hybridized modes [29,43]. The hybrid magnon-microwave modes have frequencies ω ± that are separated by We will refer to the mode with frequency ω + as the upper hybrid mode, and the mode with frequency ω − as the lower hybrid mode. When the microwave drive is set at a frequency between ω + and ω − , the scattering from the blue sideband of the lower hybrid mode can be balanced by the scattering to the red sideband of the upper hybrid mode, yielding dynamical backaction evasion: the magnomechanical decay rate vanishes. Consequently, the drive at which dynamical backaction evasion happens can be obtained from the condition In a system satisfying the two-phonon triple resonance condition ω + − ω − = 2Ω b , and for resonant magnons and microwaves, such drive frequency is exactly at (ω + + ω − )/2 [43]. The ability to tune a magnomechanical system in the dynamical backaction evasion regime was recently demonstrated in Ref. [43], and is a requirement for implementing a magnomechanical-based primary thermometer [36]. Equation (3) only takes into account the interaction of a phonon mode with a single magnon mode. However, multiple magnetostatic modes can couple to a given phonon mode [30], modifying the magnomechanical decay rate. The different scattering processes to and from the additional magnomechanical sidebands can thus change the frequency at which dynamical backaction is evaded. For instance, in the experimental data shown in Ref. [43], the measured magnomechanical decay rate exhibits a shift with respect to the theoretical prediction obtained from the Hamiltonian given in Eq. (1). This shift was taken into account by adding to the magnomechanical decay rate a phenomenological correction proportional to | m | 2 , which depends on the average number of magnons driven by the microwave tone. Such correction can be attributed to the coupling to additional magnon modes, which are weakly driven by their coupling to the microwave cavity, and to magnon nonlinearities.
While nonlinearities in magnetic spheres are generally weak, the microwave drive combined with the strong magnon-microwave coupling can make nonlinear effects prominent. This is the case provided that the power of the drive is strong enough to induce an average number of magnons above a certain threshold [58], with implications for magnetoelastic effects [59]. Even for drive powers below the nonlinear threshold, magnon nonlinearities can affect the hybrid system dynamics. For instance, in a cavity-magnonic system, the Kittel mode self-Kerr nonlinearity was shown to yield considerable cavity and magnon frequency shifts under moderate driving powers [50]. Experimentally, a phonon frequency shift, as well as mechanical bistability, was reported recently [51], which points to the importance of considering such nonlinearities in the description of dynamical backaction effects.
In what follows, we include in the description of dynamical backaction both the coupling to additional magnon modes as well as magnon nonlinearities.

II. INCLUSION OF KERR NONLINEARITY AND COUPLING TO ADDITIONAL MAGNON MODES IN THE PHONON SELF-ENERGY
To derive the correction term to the self-energy, we consider adding to the Hamiltonian in Eq. (1) self-and cross-Kerr nonlinearities, and coupling to N additional magnon modes, each with annihilation operators {m j } and frequencies ω j (j = 1, ..., N ). The total Hamiltonian is thuŝ From now on, we will refer to the magnon modem as the Kittel mode, since this is typically the magnon mode that has the strongest coupling to the cavity, while we refer to the magnon modesm j as additional magnon modes. We specify such additional magnon modes for the case of a magnetic sphere in Sec. III A. We furthermore assume a rotating wave approximation for the magnon-microwave coupling, as is done to obtain Eq. (1), which eliminates any term of the formm (j)â andm † (j)â † . The rotating wave approximation is also assumed for the magnomechanical coupling, as explained in Sec. III.
Compared with the Hamiltonian in Eq. (1), the above equation includes the following terms: the additional magnon modes; the coupling between these and (i) the microwave mode, each with a coupling rate g amj , and (ii) the phonon mode, each with a coupling rate g 0 mjb ; the self-Kerr term for the Kittel mode; the cross-Kerr term between the Kittel and the phonon modes [60]; and the cross-Kerr term between the Kittel mode and the other magnon modes. For a sphere, the values of those nonlinear terms depend on the relative orientation between the crystallographic axis [100] of YIG and the bias field [48,50,61], which was not perfectly aligned in the experiment [43] that we use as the case of study. The values for K m , K cr , and K (j) cr,m that we consider here are effectively smaller than their values in the perfectly aligned case, which we indicate as, e.g., K 0 m . Provided that the external bias field is aligned with the aforementioned magnetic anisotropy axis of the YIG sphere, the Kittel mode self-Kerr nonlinear coefficient is given by K 0 m = 13 K an γ 2 /(16M 2 s V ), where K an = −610 J/m 2 at room temperature and V is the sphere volume [49]. For the experiment in Ref. [43] this corresponds to K 0 m /2π = −5.15 nHz. The magnon-phonon and magnon-magnon cross Kerr nonlinear coefficients depend on the overlap between these modes and the Kittel mode. In general, the magnon-magnon cross-Kerr coefficient is around the same order of magnitude as K 0 m [60], while the magnonphonon cross-Kerr coefficient is ∼ −5 pHz [51]. Figure 1 shows a schematic of the model, including the different coupling terms.
The values of the magnomechanical couplings depend on the geometry of the magnet, which defines the magnon and phonon mode profiles. In the system under study, the coupling between the Kittel mode and a relevant mechanical mode of a sphere, discussed in Sec. III A, is g 0 mb /2π = 4.56 mHz. In Sec. III we discuss in detail the values of g 0 mjb for a magnetic sphere. It is important to point out that, due to better mode overlap, in principle, g 0 mjb can be comparable to or larger than g 0 mb for some modes. The coupling between magnons and microwaves depends on the microwave field at the magnet position. For homogeneous fields, only the coupling to the Kittel mode does not vanish. Nevertheless, small inhomogeneities could yield a small microwave-magnon coupling g amj which would, in turn, drive weakly such magnon modes. For different cavity geometry, such couplings can be strong [7,8,62,63], a framework which we do not consider here.
In correspondence with the experiment [43], we assume that the additional magnon modes are weakly driven via their coupling to the cavity mode, such that we expect a small steady-state amplitude for those modes. Thus we can safely disregard any self-and cross-Kerr nonlinearity of the formm † km km † jm j . The Heisenberg-Langevin equations describing the dynamics of the coupled modes in the rotating frame with the drive frequency arė (8) In the above equations, κ = κ i + κ e denotes the total microwave cavity decay rate, which is composed of the intrinsic cavity decay κ i and the decay into the external port κ e . The additional magnon modes decay rates are indicated by γ mj which, for magnetostatic modes of a sphere, have the same value of the Kittel mode decay [64]. The magnon decay γ m included in our formalism corresponds to the Gilbert damping term included in the Landau-Lifshitz equation to describe magnetic damping. All parameters appearing in Eq. (8) are summarized in Table I, with the values that will be used throughout this paper. The noise terms denoted byξ η (t) with n Th,η = [exp( ω η /k B T ) − 1] −1 the number of thermal excitations of mode η at a temperature T .
The steady state in a mean-field approximation is obtained by taking the expectation values of the operators in Eqs. (8) and ignoring any quantum correlations, i.e. mb ≈ m b . Since we are assuming that the magnon modes {m j } are weakly coupled to the microwaves, g amj g am k g amj g am , we discard any other indirect coupling between the additional magnon modes via the cavity. These approximations yield where we have defined For m j , we have also discarded the term ∝ g 0 mjb . The steady state of the Kittel mode reads where Equation (12) is solved numerically. Depending on the drive power and the detuning, the equation can have two bistable solutions. We will focus our analysis on a detuning range lying in between the hybridized Kittel magnon-microwave modes. In the considered range, the magnomechanical decay rate of Eq. (2) changes its sign, and in such a region, the nonlinear equation for m has only one solution. Furthermore, we can discard the terms proportional to K (j) cr,m , K cr , and g 0 mb to obtain the solutions of Eq. (12).

A. Linearized dynamics
We can now consider the fluctuations around the steady-state values. We writeô = δô + ô , and discard any terms involving more than two fluctuations. The quadratic Hamiltonian describing the dynamics of the fluctuations is given bŷ with the coupling terms included inĤ Int given bŷ The interacting terms appearing in the Hamiltonian of Eq. (7) induce frequency shifts for the fluctuations, which are given bỹ (16) The coupling rates between the fluctuations are enhanced and modified with respect to the bare ones. Their expressions are shown in Table II.
cr,m m mj * gB,j K The phonon self-energy is obtained by solving the linear Heisenberg-Langevin equations describing the coupled dynamics of the fluctuations for the phonon operator. To compute the effects of backaction in the response of the phonon mode to noise, we consider the Fourier transformed operators defined bŷ j , δb ( †) . We skip the algebraic steps, but outline the main differences with respect to the results in Refs. [34,36]. After writing the cavity operator in terms of the magnon operators, we obtain the following equation for the additional magnon modes: whereξ mj represents the noise term modified by the interaction with the cavity, and we have defined the effective magnon susceptibility in correspondence with the previous case in Eq. (4), We notice that the first term on the right-hand side of Eq. (18) includes the indirect coupling between the j-th magnon mode and the Kittel mode via the cavity. A similar term related to the coupling between the additional magnon modes appears in the last line of Eq. (18), and since g amj g am , we discard these contributions. After using Eq. (18) to eliminate the additional magnon modes in favor of the Kittel mode and the phonon fluctuations, we obtain the following set of coupled equations We have included all the noise terms inξ m,b [ω]; all other functions appearing in the equations below are defined in the following. The coupling to the additional magnon modes has the following effects: the introduction of a selfenergy term on the phonon mode, the modification of the coupling between the Kittel mode and the phonon mode, and a modification of the Kittel mode susceptibility and squeezing. We briefly comment on each of these effects. At this intermediate step, the phonon susceptibility is modified by the self-energy term which is defined in analogy with the self-energy term derived in Refs. [34,36] and given in Eq. (3). Such a term represents the direct dynamical backaction of the coupling between the phonon mode and the additional magnon modes. The additional magnon modes modify the couplings between the Kittel mode and the phonon mode. In fact, the effective coupling constants appearing in Eqs. (20), whose explicit forms are given in Appendix A, include two types of modifications. The first is an indirect coupling between the additional magnon modes and the Kittel mode via the cavity. The second is a term proportional to g B(R),j , which in turn (see Table II) is due to the magnon cross-Kerr nonlinearity,m † jm jm †m in the Hamiltonian of Eq. (7). The relevance of these corrections for a given drive frequency is determined by the susceptibilities of the additional magnon modes.
The Kittel mode squeezing term Λ m [ω] reads where the first term is due to the self-Kerr nonlinearity, while the second term is a combination of the magnon cross-Kerr nonlinearity with the indirect coupling between the magnon modes via the microwave cavity. The Kittel mode susceptibility is also modified by the term η m [ω], which is given by The three terms in Eq. (23) describe the effects of a twomode squeezing between the Kittel mode and each of the additional magnon modes. The first term is related to the indirect coupling of the modes via the microwave cavity while the last term is the direct two mode squeezing induced by the magnon cross-Kerr nonlinearity. The second term is a combination of both effects. From Eqs. (20), we eliminate the Kittel mode operator and obtain an equation for the phonon mode operator, whereΥ b [ω] includes all the noise terms driving the phonon fluctuations, Λ b [ω] describes phonon squeezing, and Σ Tot [ω] is the total self-energy. We focus only on the self-energy term, which is given by where the contribution of Kittel mode to the self-energy is given by The modified magnomechanical couplingsG m,R(B) [ω] are given in Appendix A. The Kittel mode susceptibilitỹ which includes both the Kittel magnon squeezing in Λ m [ω] and the effects of two-mode squeezing interactions with the additional magnon modes in η m [ω].

C. Magnomechanical decay rate corrections
We turn our attention now to the effects on the magnomechanical decay rate, in connection with the observed shift reported in Ref. [43]. The total change in the phonon linewidth is given by The contributions of self-and cross-Kerr nonlinear terms are included in the above self-energy. The corresponding frequency shift is different than the static one, observed in Ref. [51], which has already been included in the modified phonon frequencyΩ b defined in Eqs. (16). The self-Kerr nonlinearity changes the behavior of the magnomechanical decay rate with the detuning. This is due to three main factors: the modification of the steady-state number of magnons | m | 2 , the induced static magnon frequency shift, given in Eq. (16), and the generation of squeezing in the magnon fluctuations. This has a consequence for both dynamical backaction evasion, as we will show, and for backaction cooling. The latter has been recently reported in a system where a mechanical oscillator is parametrically coupled to a nonlinear cavity [55]. The magnomechanical frequency shift is given by which has a decomposition similar to the one shown in Eq. (28). The relevance of the different corrections due to the additional magnon modes depend on their frequencies with respect to the drive and their coupling rates to the cavity. Recalling that G j = g 0 mjb m j , we use the steady state from Eq. (10). For the experimental parameters in consideration, given in Table I Within these approximations, the contribution of such Walker modes to the magnomechanical linewidth is given by cr,m gives a contribution to the magnomechanical decay rate which is proportional to | m | 2 . The contribution to the magnomechanical decay rate also depends on the detuning between the drive and the magnon frequency.
While Γ j [ω] quantifies the direct effect of the coupling between the phonon mode and the additional magnon modes, there are also indirect effects due to the coupling between the additional magnon modes and the Kittel mode via the microwave cavity. Those are included in the term 2Im[Σ m [ω]] via the modified coupling ratesG m(b),R(B) [ω] and the modified Kittel mode sus-ceptibilityΞ m [ω]. In general, the corrections included in those terms are proportional to |g ms | 2 , to |G R(B),j | 2 or to |g R(B),j | 2 . Following the same procedure outlined above, it is possible to show that those are all proportional to the steady-state Kittel mode occupation. Their frequencydependent coefficients have a more complicated form due to the explicit dependence of the modified couplings and susceptibilities on frequency. We can describe all the corrections by the expression where Γ 0 mag [ω] is given in Eq. (3), without including additional magnon modes and nonlinearities. Such a correction was used phenomenologically in Ref. [43] to explain the observed discrepancies of the experimental results with Γ 0 mag [ω].

III. PARAMETRIC MAGNOMECHANICAL COUPLING FOR A MAGNETIC SPHERE
In the previous section, we obtained the modifications of the phonon self-energy due to the coupling between the phonon mode and additional magnon modes. Such contributions depend on the relative strength of the magnomechanical coupling to the additional magnon modes with respect to the coupling to the Walker mode, which in turn depends on the geometry of the magnet. Before specifying the latter, we briefly review the main points of the derivation of the magnomechanical coupling Hamiltonian, done in detail in the literature [29,30,65,66], starting from the magnetoelastic energy [25,26] where M x,y,z are the magnetization components, M S is the saturation magnetization, and is the linear strain tensor with u(r, t) the displacement. The magnetoelastic coefficients B 1 and B 2 are materialand temperature-dependent constants. For YIG at room temperature, M S = 140 kA/m [49], B 1 = 3.48 × 10 5 J/m 3 and B 1 = 6.4 × 10 5 J/m 3 [29]. The magnomechanical Hamiltonian is then obtained by quantizing the magnetization and displacement fields. The procedure yields both resonant and parametric interactions. The first is relevant for nanometer-sized magnets [30], while the second is relevant when the magnon mode is not resonant with the phonon mode, which is typically the case for larger magnets. Details are shown in Appendix B.
Here, we focus on the parametric phonon-magnon coupling, which is given by the Hamiltonian where {j = k} indicates that the sum is over all j's not equal to k and without repeating combinations. In the above equation we have separated the coupling terms between one magnon mode and one phonon mode and the terms involving two different magnon modes and a phonon mode. The explicit forms of the couplings are given in Ref. [30] and in Appendix B in Eq. (B8). They depend on overlap integrals involving the magnon mode functions as well as the derivatives of the displacement field.
Focusing now on the coupling to a specific phonon mode, the magnomechanical Hamiltonian is given bŷ where the coupling terms arê In the above equations, we have separated the terms of the Kittel mode, which from now on we do not label, while the other Walker modes are labeled by the index j. As shown in Appendix D, we can absorb the phase of the coupling to the Kittel mode in the phonon field, such that g 0 mb is real after such a transformation.

A. Magnomechanical coupling rates for a sphere
The magnomechanical coupling rates depend on the specific geometry of the sample and the direction of the applied magnetic field, which defines the mode functions and their overlap. We evaluate now the coupling for a YIG sphere, corresponding to the experimental configuration of Ref. [43].
A sphere supports magnetostatic modes called Walker modes [46,47,67], which have frequencies that can be Even though the magnon and phonon mode functions are given in terms of well-known special functions, the coupling constant (see Eq. (B8)) involves a nontrivial combination of derivatives of those. Furthermore, the derivation of the Walker modes involves a transformation to a nonorthogonal coordinate system, which is not easily invertible. While this is not a problem when computing the coupling to the Kittel mode, which is uniform, the exact expressions for coupling rates are not elucidating. Differently from other parametrically coupled systems, it is hard to infer, for example, selection rules. We, therefore, compute the overlap integrals numerically and evaluate how the couplings g 0 mjb compare with the coupling to the Kittel mode g 0 mb . It is also important to notice that the magnomechanical couplings depend on both the intensity of the bias magnetic field and its direction. In fact, the coupling to the Kittel mode can even vanish for specific relative orientation of the magnetic field [29]. We consider the case of a fixed bias field at a direction that maximizes the coupling between the Kittel and the S 122 modes, as depicted in Fig. 2. Figure 3 shows the frequencies ω mj of the Walker modes, the ratios |g 0 mjb |/|g 0 mb | between the magnomechanical coupling rate to the Walker mode (lmν) and to the Kittel mode, and φ j − φ, the relative phase between g 0 mjb and g 0 mb . Results are shown for l up to 4 and for Walker modes lying in a frequency range close to the Kittel mode. Due to better mode overlap, some higher order Walker modes, for example, the (200), couple strongly with the phonon mode in comparison with the coupling to the Kittel mode. In the theoretical analysis of Sec. I, we have not included in the magnomechanical Hamiltonian in Eq. (7) the last two terms of Eq. (D1). Those describe scattering processes between different magnon modes via the phonon mode. For the considered case, |g 0 mmjb |, |g 0 m k mjb | |g 0 mb |, |g 0 mjb |, and those processes can be safely discarded. Nevertheless, it is possible that for some relative orientation between the magnon modes and the phonon mode, set by the external bias field, those processes can have a stronger coupling rate.

IV. EVALUATION OF THE MODEL FOR THE PHONON SELF-ENERGY ON DYNAMICAL BACKACTION EVASION
The self-energy obtained in Eq. (28) includes contributions due to the Kerr nonlinearity and to the couplings to higher-order Walker modes. We focus our analysis now on the effect of such contributions to the magnomechanical decay for detunings in the vicinity of the backaction evasion point. In correspondence with the experiment [43], we consider that the microwave drive frequency is varied between ω d,− and ω d,+ inside the frequency range {ω − , ω + } between the hybrid mode frequencies, given by Eq. (5). The Walker modes contributing appreciably to the phonon self-energy lie between (ω d,− − Ω b ) and (ω d,+ + Ω b ). Since the system is in the resolved sideband regime, any modes outside this frequency range would not allow efficient scattering of phonons, and thus can be neglected. The first condition corresponds to the lower drive frequency corresponding to the blue sideband of a magnon mode, while the second corresponds to driving the red sideband of a magnon mode. For the parameters summarized in Table I, only the mode (4, 3, 0) is in this frequency range. In fact, the (4, 3, 0) mode is degener-  [36]. The magnomechanical coupling to the (4, 3, 0) Walker mode corresponds to that shown in Fig. 3. The driving power is 15 mW. Parameters in correspondence with the experiment [43], given in Table I. ate with the Kittel mode. The frequency configuration is shown in Figs. 4(a,b), and the mode profile of the Walker mode (4, 3, 0) is shown in Figs. 4 (c,d).
To address the effects of nonlinearities and coupling to the higher order Walker mode, we define the dimensionless parameters η c = g amj /g am and η K = K m /K 0 m . The first parameter quantifies the strength of the coupling between the Walker modes and the cavity compared to the coupling between the Kittel mode and the cavity. The second parameter quantifies the strength of the self-Kerr nonlinearity compared to the value shown in Table I K 0 m = −2π×5.15 nHz. We call here η K the dimensionless Kittel magnon self-Kerr nonlinearity. This parameter depends on the alignment between the anisotropy axis of the magnet with the external magnetic field, which has not been taken into account in Ref. [43].  [36]. The magnomechanical coupling to the (4, 3, 0) Walker mode corresponds to that shown in Fig. 3. The driving power is 15 mW. Parameters in correspondence with the experiment [43], given in Table I. function of the detuning. For a fixed Kerr nonlinearity, the additional magnon mode shifts down the magnomechanical decay; that is, the weakly driven Walker mode adds energy to the vibrational mode, yielding a negative contribution to the decay rate.
The magnomechanical frequency shift is also modified by the Kittel mode nonlinearity and by the coupling to the additional magnon mode. This is depicted in Fig. 6, which shows the magnomechanical frequency shift δΩ = −Re[Σ Tot [Ω b ]] for (a) several values of the self-Kerr nonlinearity and (b) several values of the coupling to the additional magnon mode. Whereas the Kerr nonlinearity induced a tilt in the slope of the magnomechanical decay rate, its effect on the magnomechanical frequency shift consists of an extra negative shift. We also notice that the magnomechanical frequency shift does not vanish for a drive at the frequency where the magnomechanical decay vanishes. This is the case because for the parameters considered, the Kittel mode frequency does not match the microwave frequency. For perfectly matching Kittel mode and microwave frequencies, and in the absence of additional magnon modes, both the magnomechanical decay and the magnomechanical frequency shift vanish at the same drive frequency [36].
In Fig. 5 one notices that the drive frequency at which the magnomechanical decay vanishes changes with both η K and η c . For applications where evading backaction is important, it is necessary that such modifications are taken into account. We show in Fig. 7 the drive frequency for backaction evasion (with respect to the upper hybrid mode frequency) as a function of the drive power for (a) several values of the Kittel self-Kerr nonlinearity and (b) several values of the coupling to the additional magnon mode at a fixed η K . For the case without nonlinearities and without coupling to the additional magnon mode, the backaction evasion frequency has a weak dependency on power (not perceptible in the plot). When the corrections are included, a stronger linear dependency of the backaction drive frequency with the power is induced. For the parameters in consideration, the difference can be of the order of ∼ 0.1 MHz at moderate powers of 10 mW.  [36]. Parameters in correspondence with the experiment [43], given in Table I.
In order to quantify the agreement between our model and the measured data in Ref. [43], we study the difference between the theoretical magnomechanical decay Γ mag [Ω b ] rate and the experimental data Γ exp . In Fig. 8, we show the absolute difference |Γ mag [Ω b ]−Γ exp | between theory and experiment as a function of the drive power at the device for different drive frequencies. We should notice that in [43], there is a power loss of ∼ 2.38 dBm between the generator and the device that has already been taken into account for this plot. Our proposed model agrees well with the experimental data, besides the difference at higher powers and drives farther from the upper hybrid mode, as it is evident in Fig. 8(a). In the worst case, the model proposed here improves the discrepancy between data and theory from ∼ 120 Hz (dashed red curve in Fig. 8), to a difference of ∼ 50 Hz (solid red curve in Fig. 8). We also notice a significant overlap at low drive powers and for all the depicted detunings due to measurement errors. Otherwise, we notice good agreement between theory and experiment for drive powers up to ∼ 14 mW. At such powers, the coherent number of magnons generated by the microwave drive | m | 2 , see Eq. (12), is between ≈ 6.0 × 10 13 at a detuning from the upper hybrid mode ∆ + = −11 MHz and ≈ 7.4 × 10 13 at ∆ + = −14 MHz. We should notice that for the parameters considered here, the system is not in a bistable regime.  (4), while the solid lines correspond to the theory developed in this paper. The shaded region corresponds to the experimental errors of the data obtained in Ref. [43]. Theory predictions use parameters in correspondence with the experiment [43], given in Table I. In Figs. 9 we show in (a-c) the magnomechanical decay as a function of the drive frequency detuning from the upper hybrid mode. The powers indicated in those figures are at the device and already take into account the 2.38 dBm losses between the generator and the device. As in the discussion above, we consider the coupling only to the (4, 3, 0) Walker mode, and we choose η K = 0. 25 and FIG. 9. Comparison between the magnomechanical decay rate Γmag[Ω b ] predicted by Eq. (3) (dashed gray line), by Eq. (28) (blue line) and the experimental data measured in Ref. [43] (magenta points). In these plots we have used ηc = 0.3 and ηK = 0.25, which yields a good agreement between our model and the experimental data, especially close to the point of dynamical backaction evasion. (a)-(c): Magnomechanical decay rate as a function for the detuning; (d)-(f) absolute difference between the theory and the experiment as a function of the detuning (we have omitted the error bars in these plots for a better visualization). Theory curves with parameters in correspondence with the experiment [43], given in Table I, and the driving powers indicated are at the device. η c = 0.3, which yields a good agreement between theory and data. In the plots of Figs. 9 (d-f), we show the difference |Γ mag [Ω b ] − Γ exp |. While the correction due to the coupling to the (4, 3, 0) Walker mode improves the agreement between theory and data with respect to the previous theory framework [34], we notice, as shown in Fig. 8, that at higher drive powers, there is a further discrepancy with the experiment for drives away from the dynamical backaction evasion points. We also notice that the errors of the data shown in Fig. 9 (a) do overlap with both the present theory and the one used in Ref. [34], as it is also shown in Fig. 8.
We attribute the discrepancy at higher powers to other nonlinear phenomena that were not considered here. Specifically, the scattering of magnon modes into the spin wave continuum via three-and four-magnon processes are known to generate an instability of the spin waves, the Suhl instabilities [68][69][70]. Above a certain threshold drive power, the amplitude of the spin waves increases at the expense of a reduction in the amplitude of the magnon modes. This phenomenon scales with the drive power and is more prominent for drive frequencies close to the Kittel mode frequency. A power-dependent reduction in the average number of Kittel magnons would imply a power-dependent reduction of the magnomechanical decay rate, which, according to the requirements for the onset of the Suhl instabilities, would occur more prominently for drives at frequencies closer to the Kittel mode. Such effect is compatible with the behavior exhibited by the data, for example, in Fig. 9, where a reduction in the magnomechanical decay happens for drive frequencies that are relatively close to the Kittel mode frequency (see the schematic in Fig. 4). In Ref. [14], it was shown that the nonlinear behavior could occur for spheres of 100 µm diameter at driving powers ∼ 1 mW. Larger spheres, such as the ones considered here, would require a stronger driving power, but the drive frequency and powers at which the discrepancy is noticeable are compatible with the requirements for the onset of Suhl instabilities. A formal evaluation of such effects goes beyond the scope of this work and will be treated elsewhere.

V. CONCLUSIONS
Dynamical backaction effects in magnomechanical systems are a consequence of the radiation pressure like coupling between magnons and phonons [29,34], which can be exploited for applications ranging from generating entangled states to noise-based thermometry [36]. In this paper, we have extended the description of dynamical backaction in cavity magnomechanics by including in the system's dynamics self-and cross-Kerr nonlinearities, and the coupling between the phonon mode and additional magnon modes. While nonlinearities are intrinsic to magnetic systems due to, e.g., magnetic anisotropy [49], magnon modes other than the uniform Kittel mode are always present and can couple to phonons as efficient as (if not more than) the Kittel mode. A nonuniform microwave field can weakly drive such modes, which modifies the backaction-induced decay rate and frequency shift of the phonon mode. Our framework considers a single phonon mode, an assumption that can be readily generalized.
We have obtained the phonon self-energy, including the aforementioned interactions and showed that, provided that the additional magnon modes couple only weakly to the microwave mode, the overall correction to the magnomechanical decay rate is proportional to the average number of Kittel magnons. We have then focused our results on the case of a magnetic sphere, in connection with the experiment performed in Ref. [43]. Our model explains the observed shift in the magnomechanical decay rate close to the dynamical backaction evasion drive frequency. In this context, we have also evaluated the effects of the different corrections. Specifically, we showed that the drive at which the dynamical backaction decay is zero depends linearly on power. This is a consequence of the corrections being proportional to the steady-state number of Kittel magnons, which scales linearly with the drive power.
A small discrepancy with the experimental data is still present at higher drive powers and for detunings far from the upper hybrid mode. We attribute this difference to nonlinear processes, such as scattering into spin waves generating Suhl instabilities [14,[68][69][70]. Furthermore, even higher-order Walker modes can lie in a frequency range close to the microwave drive and, at higher powers, can modify the magnomechanical decay. Preliminary calculations including five more Walker modes have shown an incremental improvement of our model. We should also point out that the experimental setup in Ref. [43] has particularities not included here. For instance, the magnetic sphere is glued on a dielectric post, which modifies the photon, phonon, and magnon mode profiles. This in turn can change the magnomechanical coupling constants as well as the frequency of the Walker modes. A precise evaluation of such effects requires a more refined numerical analysis, for example, using finite difference software and micromagnetic simulations, which goes beyond the scope of our analysis.
While nonlinear effects in cavity magnomechanical sys-tems have been previously computed for the nonlinear dynamics of magnons [14,51], the evaluation of such effects on the response of the mechanical degree of freedom to noise, as computed by the self-energy, is a step forward in the characterization of these systems as platforms for quantum technologies. Our analysis was restricted to evaluate the effects of all the corrections included in the model of Sec. II in the framework of dynamical backaction evasion set by the experiment of Ref. [43]. Nevertheless, the model derived in Sec. II shows that several phenomena play a role in the modification of dynamical backaction, for example, magnon squeezing and twomode squeezing. It would be interesting to investigate scenarios in which those terms can be harnessed to reduce noise for quantum metrology. Furthermore, the inclusion of the additional magnon modes opens new possibilities for cavity magnomechanical systems, such as the manipulation of the mechanics by driving different sidebands of the different magnon modes in a Floquet-like setup [71]. In this case, it would be interesting to go beyond the approximation used here, where only the Kittel mode couples strongly to the microwave cavity. In fact, several experiments have shown fingerprints of a strong coupling between Walker modes of a sphere and microwaves [7,50,63]. As we have numerically shown, Walker modes other than the Kittel mode can couple better to the phonons, which can be harnessed to applications, such as nonreciprocal transport between phonons and microwaves [72].
The second terms in Eqs. (A1) represent the effect of the indirect coupling between the additional magnon modes and the Kittel mode via the cavity. The third terms are proportional to g B(R),j , which in turn (see Table II of the frequency of the magnon modes j and k, where such an interaction describes the creation of a pair of magnons via the annihilation of a phonon; and (iii) parametric phonon-magnon coupling ∝ g 0 mjm k bαm † jm kb +H.c. The off-resonant termsm jmkbα andm † jm † kb † α can be eliminated via a rotating wave approximation.
For the phonon modes, we consider an unpinned sphere and stress-free boundary conditions [76]. There are two families of mechanical modes of a homogeneous sphere: torsional (T) and spherical (S) modes. Torsional modes are purely shear modes, while spherical modes involve both shear and compression. Both families of modes are labeled by three indices {νlm}, where l and m are polar and azimuthal indices −l ≤ m ≤ l while ν is a radial index. We focus here on S modes, whose frequencies are given by [76]