# Self-modulation instability in high power ferromagnetic resonance of BiYIG nanodisks

Igor Ngouagnia Yemeli, Salvatore Perna, Diane Gouéré, Amel Kolli, Soraya Sangiao, José María de Teresa, Manuel Muñoz, Abdelmadjid Anane, Massimiliano d'Aquino, Hugo Merbouche, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Igor Ngouagnia Yemeli, Salvatore Perna, Diane Gouéré, Amel Kolli, Soraya Sangiao, et al. Self-modulation instability in high power ferromagnetic resonance of BiYIG nanodisks. 2025. hal-04983769

HAL Id: hal-04983769

https://hal.science/hal-04983769v1

Preprint submitted on 9 Mar 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Self-modulation instability in high power ferromagnetic resonance of BiYIG nanodisks

I. Ngouagnia Yemeli,¹ S. Perna,²,* D. Gouéré,³ A. Kolli,¹ S. Sangiao,⁴,⁵ J. M. De Teresa,⁴

M. Muñoz,⁶ A. Anane,³ M. d'Aquino,² H. Merbouche,¹ C. Serpico,² and G. de Loubens¹,†

¹SPEC, CEA, CNRS, Université Paris-Saclay, Gif-sur-Yvette, France

²Department of Electrical Engineering and ICT, University of Naples Federico II, Naples, Italy

³Laboratoire Albert Fert, CNRS, Thales, Université Paris-Saclay, Palaiseau, France

⁴Instituto de Nanociencia y Materiales de Aragón (INMA), CSIC-Universidad de Zaragoza, Zaragoza, Spain

⁵Laboratorio de Microscopías Avanzadas (LMA), Universidad de Zaragoza, Zaragoza, Spain

⁶Instituto de Tecnologías Físicas y de la Información (CSIC), Madrid, Spain

We study the high power ferromagnetic resonance (FMR) of perpendicularly magnetized BiYIG nanodisks where the uniaxial anisotropy almost compensates the shape anisotropy. We observe a strong saturation of the averaged magnetization upon moderately increasing the amplitude of the rf field and a broadening of the FMR line towards lower and higher magnetic field. Full micromagnetic simulations reveal that a self-modulation of the dynamic magnetization is responsible of this behavior. To get more insight into this unstable dynamics, it is analysed in terms of normal modes. The number of modes involved is found to rapidly increase above the critical threshold. Still, a normal modes model taking into account only a few of them and their mutual nonlinear couplings allows us to qualitatively reproduce the observed phenomenon. The normal modes analysis and micromagnetic simulations also predict a Suhl-like instability at larger excitation power when it is slowly increased from low values, and bistability. Using two-tone spectroscopy, we directly measure the self-modulation spectrum and provide experimental evidence of bistable dynamics. These findings open some perspectives of using high dimensional dynamics in magnetic nanostructures for unconventional information processing.

Nonlinear magnetization dynamics [1–4] exhibit various phenomena, which have been extensively studied for several decades. They range from spin-wave (SW) instabilities [5–9], auto-oscillations and chaos [10–13], to magnetic solitons [14–16] and the condensation of magnons [17–19]. Beyond the fundamental interest of such phenomena and their counterparts in other fields of physics, perspectives of using nonlinear magnetization dynamics for unconventional information processing with magnonic and spintronic devices have recently emerged, such as for symbolic dynamics [20], random number generation [21], and neuromorphic computing [22].

Ferromagnetic resonance (FMR) is well adapted to investigate nonlinear phenomena. There, a uniform rf field is used to resonantly drive the dynamical system in the strong out-of-equilibrium regime [23]. In large ferromagnetic samples, such as millimeter sized yttrium iron garnet (YIG) spheres, a rich phase diagram emerges due to interactions between quasi degenerate SW modes [24]. In out-of-plane magnetized thin films, the main FMR mode lies at the bottom of the SW spectrum, and by decreasing the sample dimensions, the density of normal modes can be lowered, with signatures of this discretization in the nonlinear regime [25, 26]. In nanostructures, the strong geometrical confinement completely lifts the degeneracy between quantized SW modes, so that energy and momentum conservation laws are more difficult to fulfill for nonlinear scattering processes [27]. As a result, SW instability thresholds are pushed back to higher power, and spatially coherent magnetization precession can be driven at very large amplitude in submicron disks, accompanied by a strong shift of the resonance line and bistable foldover due to the large nonlinear frequency shift [28].

This nonlinear frequency shift can be controlled by adjusting the total effective anisotropy of the thin film, in particular,

its sign can be tuned [29]. In the situation where it is close to zero, large angles of coherent precession can be achieved without significant shift of the resonance conditions [30], which can be useful for the realization of certain phenomena. In fact, spin-orbit torque induced emission of coherent propagating SWs [31], SW amplification [32] and magnon condensation [19] were recently demonstrated based on hybrid BiYIG / platinum bilayers, which was only possible thanks to the uniaxial anisotropy compensating the shape anisotropy in the thin bismuth doped YIG layer [33, 34].

In this context, we investigate the nonlinear FMR of out-of-plane magnetized BiYIG submicron disks with nearly vanishing effective anisotropy. From the above arguments, one would expect that large amplitude single mode oscillation without nonlinear frequency shift could be achieved. Instead, we observe a strong saturation of the FMR line when the amplitude of the rf field is increased, together with the appearance of a frequency modulated spectrum and of bistable dynamics. Based on micromagnetic simulations and a normal modes model, we ascribe these surprising results to the onset of self-modulation of the precession profile, where dipolar interactions between quantized SW modes play a major role. Our analysis shows that above a threshold power there is a chaotic attractor that is created, which, for certain microwave power, coexists with an oscillation synchronous to the rf field. For higher rf field power, there is a Suhl-like mechanism where such oscillation becomes unstable and chaotic motion is the only possible observable. Moreover, the number of modes involved in the chaotic dynamics depends exponentially on the rf field amplitude. These phenomena affect the FMR response which looks very different from the typical one of a uniform magnet with uniaxial anisotropy.

The studied sample is composed of a few BiYIG disks with
<figure><img src="image_3.png" /><figcaption>FIG. 1. (a) Out-of-plane MRFM spectroscopy at 5 GHz of a BiYIG nanodisk of diameter 700 nm as a function of rf power, varied from -16 dBm to +9 dBm by steps of 1 dBm. (b) Longitudinal reduction of magnetization extracted from micromagnetic simulations as a function of rf field amplitude, varied from 0.02 mT to 1 mT by logarithmic steps.</figcaption></figure>

diameters 500, 700 and 1000 nm patterned from a 30 nm Bi<sub>1</sub>Y<sub>2</sub>Fe<sub>5</sub>O<sub>12</sub> thin film grown by pulsed laser deposition on a (111)-oriented substituted-GGG substrate, optimized to exhibit close to zero effective magnetization, μ<sub>0</sub>M<sub>eff</sub> = -10 mT, and low damping, α = 7·10<sup>-4</sup> [34]. A 4.5 μm wide gold microwave antenna is patterned on top of them to supply a spatially uniform, linearly polarized in-plane rf field. The BiYIG disks are saturated by a magnetic field H<sub>z</sub> applied along the normal z to the plane. We use a magnetic resonance force microscope (MRFM) to perform the FMR spectroscopy of the individual disks [35]. The MRFM probe, made of a cobalt nanosphere grown at the apex of a soft cantilever [36], sensitively monitors the averaged longitudinal component of magnetization M<sub>z</sub> of the BiYIG disk under rf excitation, which is pulse modulated at the cantilever frequency, f<sub>c</sub> ≈ 13 kHz.

The broadband MRFM spectroscopy performed in the linear regime confirms the good dynamical quality of the BiYIG disks, which all have a damping close to 10<sup>-3</sup>. In the following, we fixed the frequency of the rf field to 5 GHz. Figure 1(a) presents the evolution of the main FMR line of the 700 nm nanodisk, as a function of the rf power P, varied from -16 dBm to +9 dBm. At low P, it has a lorentzian shape with a linewidth of 0.3 mT. When P is increased, the FMR line first grows in amplitude and slightly shifts toward higher field, which is expected from the small and negative effective magnetization of the BiYIG. However, we observe a qualitative change of its aspect beyond P ≈ 0 dBm (green-yellow spectra in Fig. 1(a)): it strongly saturates in amplitude and broadens in both high and low field directions. A similar behavior is observed for all BiYIG disks with diameters ranging from 500 to 1000 nm [37].

To understand these experimental results, we performed micromagnetic simulations [38] using magnetic parameters extracted from linear MRFM spectroscopy. The 700 nm diameter disk was discretized using a 64×64×1 rectangular mesh. Figure 1(b) shows the evolution of ΔM<sub>z</sub>/M<sub>s</sub> = 1 − M<sub>z</sub>/M<sub>s</sub> = 1 − m<sub>z</sub> versus H<sub>z</sub>, swept downwards by steps of 0.1 mT, for different rf field amplitudes spanning from 0.02 mT to 1 mT.

<figure><img src="image_4.png" /><figcaption>FIG. 2. Micromagnetic simulations. (a,b) Time series of m<sub>x</sub>, (c,d) corresponding power spectra and (e,f) mode amplitudes vs. time at μ<sub>0</sub>H<sub>z</sub> = 0.169 T, respectively below (μ<sub>0</sub>h<sub>rf</sub> = 0.4 mT) and beyond (μ<sub>0</sub>h<sub>rf</sub> = 0.6 mT) the critical threshold. (g) Phase diagram in the control plane (H<sub>z</sub>, h<sub>rf</sub>). μ<sub>0</sub>H<sub>z</sub> is varied from 0.1645 to 0.172 T by steps of 0.5 mT and μ<sub>0</sub>h<sub>rf</sub> from 0.1 to 3 mT by steps of 0.1 mT, the rf frequency being fixed to 5 GHz. Red dots correspond to the critical curve obtained by directly applying the rf field to the equilibrium state, blue upwards (downwards) triangles to the one obtained by increasing (decreasing) slowly h<sub>rf</sub> from zero (maximum) amplitude. The dashed curve is the instability threshold deduced from P-mode analysis. (h) Dependence on h<sub>rf</sub> of the number of modes with amplitudes |a<sub>k</sub>| > 10<sup>-3</sup> involved in the dynamics at four different H<sub>z</sub>.</figcaption></figure>

For each $H_z$, the spatially averaged cartesian components of the magnetization are recorded for 100 ns under an harmonic excitation at 5 GHz of constant amplitude $h_{rf}$. The plotted quantity is obtained by averaging $m_z$ over the second half of the simulated time window. A good qualitative agreement with the MRFM spectroscopy data is observed. In particular, the breakdown of the FMR line is reproduced beyond $\mu_0 h_{rf} \approx 0.4$ mT, and $\Delta M_z/M_s$ saturates around 3 to 4%.

To determine the phase diagram of our dynamical system in the control plane (H<sub>z</sub>, h<sub>rf</sub>), we performed additional simulations using the MaGICo micromagnetic code [39]. For each H<sub>z</sub> value, the equilibrium configuration is excited with an rf field of frequency 5 GHz and amplitude h<sub>rf</sub> and the dynamics
is integrated for 1100 ns. The last 100 ns are used to evaluate
the type of regime. The time series of the transverse magne-
tization $m_x = M_x/M_s$ are Fourier transformed to obtain their
power spectrum. At low $h_{rf}$, $m_x$ oscillates at constant ampli-
tude at the forced frequency of 5 GHz, see Fig.2(a), resulting
in a single-peak (SP) spectrum such as shown in Fig.2(c). At
sufficiently high power, however, it does not display a time-
harmonic response anymore. Instead, one observes a low fre-
quency modulation of its amplitude, see Fig.2(b), and a multi-
peaks continuous (MP-C) spectrum such as in Fig.2(d). The rf
field amplitude responsible for the transition between SP and
MP-C dynamical states is plotted as a function of $H_z$ using
red dots in Fig.2(g), which defines a critical curve. Simula-
tions have also been run with $h_{rf}$ being slowly varied. First, $h_{rf}$
is increased starting from a SP state until the oscillations be-
come MP-C (upwards triangles). Then, $h_{rf}$ is decreased start-
ing from a MP-C state until the oscillations become SP (down-
wards triangles). Interestingly, there is a region in the control
plane ($H_z, h_{rf}$) where there is coexistence between oscillations
with MP-C and SP spectra. We also observe that the critical
curve estimated by directly applying $h_{rf}$ (red dots) quantita-
tively agrees with the curve obtained when $h_{rf}$ is slowly de-
creased from a MP-C state (downward triangles).

To analyse in more details the underlying dynamics, we
project snapshots of the magnetization field from micromag-
netic simulations onto the basis of normal modes [40]. The
latter are first obtained numerically by linearizing the Landau-
Lifshitz (LL) equation, which yields a generalized eigenvalue
problem that can be solved [41], and ordered by ascending
frequency [37]. The fundamental mode is the uniform or Kit-
tel mode. Other normal modes are quantized modes along the
azimuthal and radial directions. Only the purely radial modes
couple to the uniform rf field excitation [42]. The profiles
of five of them, 1, 6, 15, 30 and 51, are shown as insets in
Fig.2(e). The mode amplitudes evaluated in this way are plot-
ted as a function of time in Figs.2(e,f). Before the threshold,
only modes 1 and 6 have significant amplitude. The one of the
uniform mode, excited close to resonance by the rf field, is by
far the largest and is constant in the steady state. Above the
threshold, it strongly oscillates and the amplitudes of other ra-
dial modes start to be significant and oscillating, as shown in
Fig.2(f). This oscillation of modes is responsible for the one
of the precession profile visualized in movies extracted from
micromagnetic simulations, in which a strong spatio-temporal
modulation of $m_z$ in the disk is observed [37].

We can also estimate how many modes are involved in the
dynamics from the micromagnetic simulations. Fig.2(h) dis-
plays the number of normal modes with significant amplitude
($|a_k| > 10^{-3}$) as a function of $h_{rf}$ for several values of $H_z$. In-
terestingly there are two types of MP-C states: one where the
number of modes involved is around 10 and the other where
it becomes several times that. When $\mu_0 H_z = 0.169$ T, it is less
than 5 until $\mu_0 h_{rf}$ reaches 1.9 mT, where it increases to 10. A
more dramatic behavior is observed at $\mu_0 H_z = 0.165$ T, where
the number of modes involved explodes from 2 to 30 above
1.5 mT.

<figure><img src="image_5.png" /></figure>

FIG. 3. Time series of $m_x$ (insets: mode amplitudes vs. time and power spectra) calculated at $\mu_0 H_z = 0.172$ T and $h_{rf} = 3$ mT using the NMM including (a) two (b) three and (c) four modes. (d) Phase diagram in the control plane ($H_z, h_{rf}$) comparing the critical curve extracted from full micromagnetics (red dots) to the ones using the NMM with three, four and five modes (blue, orange and green dots, respectively).

Micromagnetic simulations show that several SW modes
participate to the dynamics in the nonlinear regime. Hence, a
normal modes model (NMM) taking into account the relevant
modes should be appropriate to reproduce the self-modulation
instability. Because it only considers a few degrees of free-
dom, it is computationally efficient and it can bring useful
physical insight. Shortly, the LL equation is projected onto
the normal modes basis, resulting in a system of nonlinear or-
dinary differential equations for the mode amplitudes, whose
coefficients can be pre-computed from the normal modes de-
termined previously [43, 44]. We first consider a NMM with
only modes 1 and 6, which have the largest amplitudes in mi-
cromagnetics. Integrating the two coupled equations yields
a stable dynamics at $\mu_0 h_{rf} = 3$ mT, see Fig.3(a); it becomes
unstable only at very large amplitude ($> 4$ mT). By adding
mode 15 into the NMM, the self-modulation is already present
at 3 mT, see Fig.3(b). As can be seen in Fig.3(c), it remains
upon including in the dynamics mode 30, although the spec-
trum changes. Using the NMM with three modes, one can
construct the critical line in the ($H_z, h_{rf}$) plane when the rf field
is pulsed from equilibrium and compare it to the one obtained
from full micromagnetic simulations. A qualitative agreement
is observed in Fig.3(d), but the threshold is overestimated by
the NMM. The inclusion of more modes improves the situa-
tion, still, quantitative agreement could not be obtained with
only five modes. Integrating the NMM when $h_{rf}$ is slowly var-
ied also shows the presence of a coexistence region between
SP and MP-C states [37]. The curve for increasing $h_{rf}$ is quite
close to the one computed micromagnetically, however, the
curves for decreasing $h_{rf}$ differ quantitatively. This is ascribed
to the fact that in the MP-C dynamical state the number of
modes involved is higher than that considered in the NMM.
<figure><img src="image_6.png" /><figcaption>FIG. 4. (a) Single-tone MRFM spectroscopy at 5 GHz of a BiYIG nanodisk of diameter 500 nm for three different rf powers (pulse modulated). (b), (c) and (d) Two-tone MRFM spectroscopy as a function of the up-swept field $H_z$ and the frequency $f_2$ of the low power rf excitation (pulse modulated). The power of the continuous rf excitation at $f_1 = 5$ GHz is fixed respectively to 1, 4 and 10 dBm. (f) MRFM spectrum as a function of $f_2$ extracted at fixed $H_z$ from the spectroscopy data, indicated by the vertical dashed line in (e).</figcaption></figure>

We also performed a P-mode analysis, where the amplitudes of normal modes are calculated at equilibrium in a rotating frame synchronous with the rf field angular frequency $\omega$ and the P-modes equations are derived from the hypothesis $a_k = \bar{a}_k e^{j\omega t}$, which allows to define a Suhl-like instability threshold [37]. It is found to be in good agreement with the threshold estimated from simulations when slowly increasing $h_{rf}$, see dashed line and upwards triangles in Fig.2(g).

To provide some experimental evidence of the scenario deduced from micromagnetics and normal modes analysis, we performed two-tone MRFM spectroscopy [28]. In these measurements, a continuous wave rf excitation of fixed frequency ($f_1 = 5$ GHz) and high power $P_1$ is applied to the 500 nm BiYIG disk to drive its dynamics in the nonlinear regime, while a second rf tone of low power $P_2$ and varying frequency $f_2$ is pulse modulated at the cantilever frequency to probe the dynamical state driven by the first rf tone. At $P_1 = 1$ dBm, the single-tone MRFM spectroscopy (blue spectrum in Fig.4(a)) shows no sign of instability, and the two-tone spectroscopy

<figure><img src="image_7.png" /><figcaption>FIG. 5. (a) MRFM cantilever frequency as a function of $H_z$ under continuous rf excitation at $f_1 = 5$ GHz. (b) and (c) Two-tone MRFM spectroscopy performed at $\mu_0 H_z = 0.1685$ T, coming respectively from low and high field (see red and blue dots in (a)).</figcaption></figure>

map of Fig.4(b) is a signature of a coherent P-mode dynamics, with two symmetric peaks around $f_1$ reflecting the nutation of the magnetization vector around its large amplitude steady state trajectory [28]. At $P_1 = 4$ dBm, the single-tone MRFM spectrum plotted in green in panel (a) strongly distorts and saturates for $0.1695 < \mu_0 H_z < 0.17$ T, and in this region, the two-tone spectroscopy displays a bifurcation of the two nutation peaks followed by a broad multi-peaks spectrum (Fig.4(c)), which reflects the transition to a modulation instability. Finally at $P_1 = 10$ dBm, the instability threshold is crossed for a broad range of applied field (red spectrum in panel (a), $0.1685 < \mu_0 H_z < 0.1705$ T). As a consequence, the two-tone spectroscopy map is multi-peaks almost everywhere, see Fig.4(d). A $f_2$-cut of the data at $\mu_0 H_z = 0.17$ T displaying a modulation frequency of 48 MHz, in the same range as found in simulations, is shown in Fig.4(d).

To experimentally access the bistability between SP and MP-C dynamical states, we monitor the frequency of the MRFM cantilever under continuous rf excitation when $H_z$ is swept up and down through the FMR line, as shown in Fig.5(a). An hysteretic behavior is observed between 0.1675 and 0.169 T. The two-tone MRFM spectroscopy performed at $\mu_0 H_z = 0.1685$ T coming from low and high field, shown in Fig.5(b) and (c), displays distinct signatures, assigned to coherent and frequency modulated dynamics, respectively [37].

Finally, we briefly discuss the origin of the dynamic instability observed when driving these BiYIG nanodisks into the nonlinear regime of FMR. Due to finite size effects and the non-ellipsoidal disk geometry, the fundamental mode excited by the uniform rf field has a non-uniform amplitude, despite having a uniform precession phase (mode 1 in Fig.2(e)). When its amplitude grows, it usually cannot maintain its profile, instead it spatially extends or contracts depending on the sign of the nonlinear frequency shift [28, 45-47]. When the latter is close to zero, it can keep its non-uniform profile in the nonlinear regime. The strong non-uniform dipolar field that it generates in the magnetic body makes it possible to excite other SW modes, which couple weakly to the uniform rf excitation, at a substantial level (particularly mode 6, see Fig.2(e,f)). There is thus an energy from long to shorter wavelength modes, which can lead to the modulation instability described above.
To conclude, we would like to emphasize that two types of phenomena have been unveiled by our theoretical analysis: the coexistence of oscillation synchronous with rf field (SP) and chaotic attractor (MP-C) at moderate power, and a Suhl-like instability at higher power, where the SP state becomes unstable. We have also found that the number of modes involved in the chaotic dynamics depends exponentially on power, which limits the accuracy of the normal modes model, but at the same time opens opportunities to study turbulence in confined systems. Moreover, detailed analysis of this chaotic regime is under investigation, as well as the possibility to use chaotic modes motion and multi-stability of dynamical regimes to perform computations.

This work was supported by the Agence Nationale de la Recherche (ANR) under grant no. ANR-18-CE24-0021 (Maestro), ANR-20-CE24-0012 (Marin), the PEPR SPIN project ANR-22-EXSP-0005 (SpinCom), and by the Horizon2020 Research Framework Programme of the European Commission under grant no. 899646 (k-NET). It was also supported by a public grant overseen by the ANR as part of the "Investissements d'Avenir" program (Labex NanoSaclay, reference: ANR-10-LABX-0035). S.S. and J.M.D.T thank the Gobierno de Aragón through the grant numbers E13_23R and Ayuda CEX2023-001286-S funded by MICIU/AEI/10.13039/501100011033 and European Union "NextGenerationEU"/PRTR.

<footnote>* salvatore.perna@unina.it</footnote>

<footnote>† gregoire.deloubens@cea.fr</footnote>

<footnote>[1] P. E. Wigen, ed., *Nonlinear Phenomena and Chaos in Magnetic Materials* (World ScientificWorld Scientific, Singapore, 1994).</footnote>

<footnote>[2] V. S. L'vov, *Wave Turbulence Under Parametric Excitation: Applications to Magnets* (Springer-Verlag, 1994).</footnote>

<footnote>[3] A. G. Gurevich and G. A. Melkov, *Magnetization Oscillations and Waves* (CRC Press, 1996).</footnote>

<footnote>[4] I. D. Mayergoyz, G. Bertotti, and C. Serpico, *Nonlinear magnetization dynamics in nanosystems* (Elsevier, 2009).</footnote>

<footnote>[5] P. W. Anderson and H. Suhl, Instability in the motion of ferromagnets at high microwave power levels, *Phys. Rev. 100*, 1788 (1955).</footnote>

<footnote>[6] H. Suhl, The theory of ferromagnetic resonance at high signal powers, *J. Phys. Chem. Solids* 1, 209 (1957).</footnote>

<footnote>[7] G. Bertotti, I. D. Mayergoyz, and C. Serpico, Spin-wave instabilities in large-scale nonlinear magnetization dynamics, *Phys. Rev. Lett.* 87, 217203 (2001).</footnote>

<footnote>[8] G. de Loubens, V. V. Naletov, and O. Klein, Reduction of the spin-wave damping induced by nonlinear effects, *Phys. Rev. B* 71, 180411 (2005).</footnote>

<footnote>[9] T. Gerrits, P. Krivosik, M. L. Schneider, C. E. Patton, and T. J. Silva, Direct Detection of Nonlinear Ferromagnetic Resonance in Thin Films by the Magneto-Optical Kerr Effect, *Phys. Rev. Lett.* 98, 207602 (2007).</footnote>

<footnote>[10] G. Gibson and C. Jeffries, Observation of period doubling and chaos in spin-wave instabilities in yttrium iron garnet, *Phys. Rev. A* 29, 811 (1984).</footnote>

<footnote>[11] S. M. Rezende and F. M. de Aguiar, Spin-wave instabilities, auto-oscillations, and chaos in yttrium-iron-garnet, *Proc. IEEE* 78, 893 (1990).</footnote>

<footnote>[12] S. Petit-Watelot, J.-V. Kim, A. Ruotolo, R. M. Otxoa, K. Bouzehouane, J. Grollier, A. Vansteenkiste, B. V. de Wiele, V. Cros, and T. Devolder, Commensurability and chaos in magnetic vortex oscillations, *Nature Phys.* 8, 682 (2012).</footnote>

<footnote>[13] A. B. Ustinov, A. V. Kondrashov, I. Tatsenko, A. A. Nikitin, and M. P. Kostylev, Progressive development of spin wave chaos in active-ring oscillators, *Phys. Rev. B* 104, L140410 (2021).</footnote>

<footnote>[14] A. Kosevich, B. Ivanov, and A. Kovalev, Magnetic solitons, *Phys. Rep.* 194, 117 (1990).</footnote>

<footnote>[15] M. Wu, B. A. Kalinikos, L. D. Carr, and C. E. Patton, Observation of spin-wave soliton fractals in magnetic film active feedback rings, *Phys. Rev. Lett.* 96, 187202 (2006).</footnote>

<footnote>[16] S. M. Mohseni, S. R. Sani, J. Persson, T. N. A. Nguyen, S. Chung, Y. Pogoryelov, P. K. Muduli, E. Iacocca, A. Eklund, R. K. Dumas, S. Bonetti, A. Deac, M. A. Hoefer, and J. Akerman, Spin torque-generated magnetic droplet solitons, *Science* 339, 1295 (2013).</footnote>

<footnote>[17] S. O. Demokritov, V. E. Demidov, O. Dzyapko, G. A. Melkov, A. A. Serga, B. Hillebrands, and A. N. Slavin, Bose-Einstein condensation of quasi-equilibrium magnons at room temperature under pumping, *Nature* 443, 430 (2006).</footnote>

<footnote>[18] M. Schneider, T. Brächer, D. Breitbach, V. Lauer, P. Pirro, D. A. Bozhko, H. Y. Musiienko-Shmarova, B. Heinz, Q. Wang, T. Meyer, et al., Bose-Einstein condensation of quasiparticles by rapid cooling, *Nature Nanotechnol.* 15, 457 (2020).</footnote>

<footnote>[19] B. Divinskiy, H. Merbouche, V. E. Demidov, K. O. Nikolaev, L. Soumah, D. Gouéré, R. Lebrun, V. Cros, J. B. Youssef, P. Bortolotti, A. Anane, and S. O. Demokritov, Evidence for spin current driven bose-einstein condensation of magnons, *Nat. Commun.* 12, 10.1038/s41467-021-26790-y (2021).</footnote>

<footnote>[20] M.-W. Yoo, D. Rontani, J. Létang, S. Petit-Watelot, T. Devolder, M. Sciamanna, K. Bouzehouane, V. Cros, and J.-V. Kim, Pattern generation and symbolic dynamics in a nanocontact vortex oscillator, *Nat. Commun.* 11, 1 (2020).</footnote>

<footnote>[21] N.-T. Phan, N. Prasad, A. Hakam, A. S. E. Valli, L. Anghel, L. Benetti, A. Madhavan, A. S. Jenkins, R. Ferreira, M. D. Stiles, U. Ebels, and P. Talatchian, Unbiased random bitstream generation using injection-locked spin-torque nano-oscillators, *Phys. Rev. Appl.* 21, 034063 (2024).</footnote>

<footnote>[22] L. Körber, C. Heins, T. Hula, J.-V. Kim, S. Thlang, H. Schultheiss, J. Fassbender, and K. Schultheiss, Pattern recognition in reciprocal space with a magnon-scattering reservoir, *Nat. Commun.* 14, 10.1038/s41467-023-39452-y (2023).</footnote>

<footnote>[23] G. Bertotti, C. Serpico, and I. D. Mayergoyz, Nonlinear magnetization dynamics under circularly polarized field, *Phys. Rev. Lett.* 86, 724 (2001).</footnote>

<footnote>[24] P. Bryant, C. Jeffries, and K. Nakamura, Spin-wave nonlinear dynamics in an yttrium iron garnet sphere, *Phys. Rev. Lett.* 60, 1185 (1988).</footnote>

<footnote>[25] R. D. McMichael and P. E. Wigen, High-power ferromagnetic resonance without a degenerate spin-wave manifold, *Phys. Rev. Lett.* 64, 64 (1990).</footnote>

<footnote>[26] V. V. Naletov, G. de Loubens, V. Charbois, O. Klein, V. S. Tiberkevich, and A. N. Slavin, Ferromagnetic resonance spectroscopy of parametric magnons excited by a four-wave process, *Phys. Rev. B* 75, 140405 (2007).</footnote>

<footnote>[27] G. Melkov, D. Slobidianiuk, V. Tiberkevich, G. de Loubens, O. Klein, and A. Slavin, Nonlinear ferromagnetic resonance in nanostructures having discrete spectrum of spin-wave modes, *IEEE Magn. Lett.* 4, 4000504 (2013).</footnote>

<footnote>[28] Y. Li, V. V. Naletov, O. Klein, J. L. Prieto, M. Muñoz, V. Cros, P. Bortolotti, A. Anane, C. Serpico, and G. de Loubens, Nutation spectroscopy of a nanomagnet driven into deeply nonlinear</footnote>
ferromagnetic resonance, Phys. Rev. X **9**, 041036 (2019).

[29] B. Lührmann, M. Ye, H. Dötsch, and A. Gerspach, Nonlinearities in the ferrimagnetic resonance in epitaxial garnet films, J. Magn. Magn. Mater. **96**, 237 (1991).

[30] K. Gnatzig, H. Dötsch, M. Ye, and A. Brockmeyer, Ferrimagnetic resonance in garnet films at large precession angles, J. Appl. Phys. **62**, 4839 (1987).

[31] M. Evelt, L. Soumah, A. Rinkevich, S. Demokritov, A. Anane, V. Cros, J. Ben Youssef, G. de Loubens, O. Klein, P. Bortolotti, and V. Demidov, Emission of coherent propagating magnons by insulator-based spin-orbit-torque oscillators, Phys. Rev. Appl. **10**, 041002 (2018).

[32] H. Merbouche, B. Divinskiy, D. Gouéré, R. Lebrun, A. El Kanj, V. Cros, P. Bortolotti, A. Anane, S. O. Demokritov, and V. E. Demidov, True amplification of spin waves in magnonic nanowaveguides, Nat. Commun. **15**, 10.1038/s41467-024-45783-1 (2024).

[33] L. Soumah, N. Beaulieu, L. Qassym, C. Carrétéro, E. Jacquet, R. Lebourgeois, J. B. Youssef, P. Bortolotti, V. Cros, and A. Anane, Ultra-low damping insulating magnetic thin films get perpendicular, Nat. Commun. **9**, 10.1038/s41467-018-05732-1 (2018).

[34] D. Gouéré, H. Merbouche, A. El Kanj, F. Kohl, C. Carrétéro, I. Boventer, R. Lebrun, P. Bortolotti, V. Cros, J. Ben Youssef, and A. Anane, Temperature-independent ferromagnetic resonance shift in bi-doped yig garnets through magnetic anisotropy tuning, Phys. Rev. Mater. **6**, 114402 (2022).

[35] O. Klein, G. de Loubens, V. V. Naletov, F. Boust, T. Guillet, H. Hurdequint, A. Leksikov, A. N. Slavin, V. S. Tiberkevich, and N. Vukadinovic, Ferromagnetic resonance force spectroscopy of individual submicron-size samples, Phys. Rev. B **78**, 144410 (2008).

[36] S. Sangiao, C. Magén, D. Mofakhami, G. de Loubens, and J. M. D. Teresa, Magnetic properties of optimized cobalt nanospheres grown by focused electron beam induced deposition (febid) on cantilever tips, Beilstein J. Nanotechnol. **8**, 2106 (2017).

[37] See Supplemental Material at [URL will be inserted by publisher] for supplementary figures and movies, and details on the normal modes model and P-mode analysis.

[38] (), http://micromagnetics.org/magnum.fd/.

[39] (), http://wpage.unina.it/mdaquino/index_file/MaGICo.html.

[40] M. Massouras, S. Perna, M. d'Aquino, C. Serpico, and J.-V. Kim, Mode-resolved micromagnetics study of parametric spin wave excitation in thin-film disks (2024).

[41] M. d'Aquino, C. Serpico, G. Miano, and C. Forestiere, A novel formulation for the numerical computation of magnetization modes in complex micromagnetic systems, J. Comput. Phys. **228**, 6130 (2009).

[42] V. V. Naletov, G. de Loubens, G. Albuquerque, S. Borlenghi, V. Cros, G. Faini, J. Grollier, H. Hurdequint, N. Locatelli, B. Pigeau, A. N. Slavin, V. S. Tiberkevich, C. Ulysse, T. Valet, and O. Klein, Identification and selection rules of the spin-wave eigenmodes in a normally magnetized nanopillar, Phys. Rev. B **84**, 224423 (2011).

[43] S. Perna, F. Bruckner, C. Serpico, D. Suess, and M. d'Aquino, Computational micromagnetics based on normal modes: Bridging the gap between macrospin and full spatial discretization, J. Magn. Magn. Mater. **546**, 168683 (2022).

[44] S. Perna, F. Bruckner, C. Serpico, D. Suess, and M. d'Aquino, Normal modes description of nonlinear ferromagnetic resonance for magnetic nanodots, AIP Adv. **12**, 035244 (2022).

[45] Y. V. Gulyaev, P. E. Zil'berman, A. G. Temiryazev, and M. P. Tikhomirova, Principal mode of the nonlinear spin-wave resonance in perpendicular magnetized ferrite films, Phys. Solid State **42**, 1094 (2000).

[46] A. Slavin and V. Tiberkevich, Spin wave mode excited by spin-polarized current in a magnetic nanocontact is a standing self-localized wave bullet, Phys. Rev. Lett. **95**, 237201 (2005).

[47] V. E. Demidov, M. Buchmeier, K. Rott, P. Krzysteczko, J. Münchenberger, G. Reiss, and S. O. Demokritov, Nonlinear hybridization of the fundamental eigenmodes of microscopic ferromagnetic ellipses, Phys. Rev. Lett. **104**, 217203 (2010).
# Supplementary Material: “Self-modulation instability in high power ferromagnetic resonance of BiYIG nanodisks”

I. Ngouagnia Yemeli,¹ S. Perna,²,* D. Gouéré,³ A. Kolli,¹ S. Sangiao,⁴,⁵ J. M. De Teresa,⁴

M. Muñoz,⁶ A. Anane,³ M. d'Aquino,² H. Merbouche,¹ C. Serpico,² and G. de Loubens¹,†

¹SPEC, CEA, CNRS, Université Paris-Saclay, Gif-sur-Yvette, France

²Department of Electrical Engineering and ICT, University of Naples Federico II, Naples, Italy

³Laboratoire Albert Fert, CNRS, Thales, Université Paris-Saclay, Palaiseau, France

⁴Instituto de Nanociencia y Materiales de Aragón (INMA), CSIC-Universidad de Zaragoza, Zaragoza, Spain

⁵Laboratorio de Microscopías Avanzadas (LMA), Universidad de Zaragoza, Zaragoza, Spain

⁶Instituto de Tecnologías Físicas y de la Información (CSIC), Madrid, Spain

## CONTENTS

I. MRFM spectroscopy vs. disk diameter 2

II. Supplementary movies 2

III. Normal modes model 3

A. Calculation of normal modes 3

B. Projection of LLG equation on normal modes basis 3

C. Reduced normal modes model 3

D. P-mode analysis 4

IV. Experimental evidence of bistable dynamics 6

References 8

<footnote>* salvatore.perna@unina.it</footnote>

<footnote>† gregoire.deloubens@cea.fr</footnote>
<figure><img src="image_8.png" /><figcaption>**Supplementary Figure 1** | Dependence on disk diameter. MRFM spectroscopy performed at 5 GHz as a function of the excitation power for BiYIG nanodisks of diameters 1000 nm (a), 700 nm (b), and 500 nm (c). The data plotted in (b) are the same as in Fig.1(a) of the main text.</figcaption></figure>

# I. MRFM SPECTROSCOPY VS. DISK DIAMETER

For each diameter among 500, 700 and 1000 nm, there are two BiYIG disks available under the microwave antenna. We performed MRFM spectroscopy on all of them, and observed a similar saturation and broadening of the main FMR line as the rf amplitude is increased. Suppl. Fig.1 shows these measurements for one particular disk of each diameter. Upper and lower panels show the same data plotted as 2D colormaps and 1D spectra vs. power, respectively.

# II. SUPPLEMENTARY MOVIES

The supplementary movies show the dynamics calculated using micromagnetic simulations at $\mu_0H_z = 0.169$ T and $\mu_0h_{rf} = 1$ mT (rf frequency equal to 5 GHz) in the 700 nm disk as a function of time. The total simulated time is equal to 130 ns, the last 30 ns being displayed in the movies.

The upper panel of SupMov1.gif shows the spatial average over the disk volume of the magnetization components along x, y and z, in blue, green and red lines, respectively. The transverse magnetization, defined as $m_t = \sqrt{m_x^2 + m_y^2}$ is also shown in cyan. The components $m_x$ and $m_y$ oscillate in quadrature at 5 GHz, and exhibit an amplitude modulation of period around 10 ns (i.e., a modulation frequency of about 100 MHz), which can also be clearly seen on the $m_z$ component. The lower left and right panels show the evolution as a function of time of the 1D spatial profiles of the transverse ($m_t$) and longitudinal ($m_z$) components vs. the reduced radial coordinate $r/R$, respectively ($R = 350$ nm is the disk radius). One clearly observes spatio-temporal oscillations of the $m_t$ spatial profile, accompanied by a strong modulation of $m_z$, which exhibits large and sudden reductions at the disk center. There, $m_z$ can even become negative at $t \approx 12$ ns, i.e., the angle of precession $\theta$ becomes larger than $\pi/2$ ($m_z = -0.5$ corresponds to $\theta = 2\pi/3$).

The upper panel of SupMov1bis.gif is the same as SupMov1.gif. The three lower panels display, from left to right, the time evolution of the 2D spatial profiles of the precession phase, $\phi = \arctan(m_y/m_x)$, of the transverse component, $m_t$, and of the longitudinal one, $m_z$. In the lower left panel, the hue color scale codes the [0, $2\pi$] interval, while in the other two panels, the gray contrast codes the [0, 1] interval (white: 0, black:1). One can observe that the phase is not constantly uniform in the disk, instead it shows periodic modulation of its spatial profile, as do the $m_t$ and $m_z$ components.
<figure><img src="image_9.png" /><figcaption>**Supplementary Figure 2** | Micromagnetic normal modes. (a) Normal mode spectrum of the 700 nm BiYIG disk calculated at $\mu_0H_z = 0.169$ T. (b) Absolute value of the spatial average of the x-component of normal modes vs. mode number. (c) Mode profiles of 12 chosen normal modes. The red/blue color represents the phase of the precession, while the contrast represents its amplitude. Red labels emphasize radial modes 1, 6, 15, 30, 51 having a nonzero average in (b).</figcaption></figure>

### III. NORMAL MODES MODEL

#### A. Calculation of normal modes

The normal modes are calculated by MaGICo [1] which solves the generalized eigenvalue problem for normal modes according to the formulation presented in [2]. Suppl. Fig.2 displays the normal modes spectrum (a), the absolute value of the spatial average of the x-component (b) of normal modes, and the norm of the real part of normal modes (c). Among the first 60 normal modes, 5 of them have significative average along the x-direction: 1, 6, 15, 30, 51, see Suppl. Fig.2(b). Their respective frequencies are outlined with red filled dots in Suppl. Fig.2(a). They correspond to radial modes with azimuthal invariance, as shown in Suppl. Fig.2(c).

#### B. Projection of LLG equation on normal modes basis

When the micromagnetic equilibrium is slightly perturbed, the magnetization dynamics can be decomposed in harmonic and independent oscillations of normal modes at their respective natural frequency. The set of normal modes can also be used as a basis to project the nonlinear magnetization dynamics described by the LLG equation. In this way, the following system of coupled nonlinear equations is derived where the unknowns are the amplitudes of normal modes:

$$ \dot{a}_k = \frac{j\omega_k}{(1 + \alpha^2)} \left[ b_k + \sum_h b_{k,h}a_h + \sum_{h,i} c_{k,h,i}a_h a_i + \frac{1}{2} \sum_{h,i,j} d_{k,h,i,j}a_h a_i a_j + \dots \right]. \quad (1) $$

This is the so-called normal modes model (NMM) [3]. In several cases of interest for applications, only few modes and few terms in the NMM are required for a quantitative analysis of the micromagnetic system under investigation. In this way, a hierarchy of approximated models can be defined, whose description is intermediate between the macrospin (uniform mode) approximation and the full micromagnetic description.

#### C. Reduced normal modes model

Here, we consider normal modes that have a significative spatial average along the rf field direction (x). Such modes correspond to index numbers 1, 6, 15, 30, 51,... We take only the first 5 of them, see Suppl. Fig.2. Not all the terms of the NMM
<figure><img src="image_10.png" /><figcaption>**Supplementary Figure 3 | P-mode analysis.** Comparison of the real and imaginary parts of modes 1 and 6 extracted vs. rf field amplitude from micromagnetic simulations (continuous lines) and P-mode analysis (dashed lines) at (a), (b) $\mu_0H_z = 0.168$ T and (c), (d) $\mu_0H_z = 0.171$. (e) Dependence on $h_{rf}$ of the real part of P-mode eigenvalues at $\mu_0H_z = 0.17$ T. (f) Instability threshold vs. $H_z$ from P-mode analysis.</figcaption></figure>

are considered in the numerical analysis with normal modes. Indeed, the spectrum of the averaged magnetization components is centered around the rf field frequency of 5 GHz. The NMM is then written in a rotating frame at the angular frequency of the rf field, $\omega_{rf}$, where the k-th amplitude of normal mode is expressed as $a_k(t) = \bar{a}_k(t)e^{j\omega_{rf}t}$ and $\bar{a}_k(t)$ is a slow time-varying function. In this way, we get the following reduced NMM:

$$ \dot{\bar{a}}_k = \frac{j\omega_k}{(1+\alpha^2)} \left[ \bar{b}_k + \left(b_{k,k} - \frac{\omega_{rf}}{\omega_k}\right) \bar{a}_k + \sum_{h,i} \left( c_{k,h,-i}^{f,cw} \bar{a}_h \bar{a}_i^* + c_{k,-h,i}^{f,cw} \bar{a}_h^* \bar{a}_i + c_{k,h,i}^{f,ccw} \bar{a}_h \bar{a}_i \right) + \frac{1}{2} \sum_{h,i,j} \left( d_{k,h,i,-j}^0 \bar{a}_h \bar{a}_i \bar{a}_j^* + d_{k,h,-i,j}^0 \bar{a}_h \bar{a}_i^* \bar{a}_j + d_{k,-h,i,j}^0 \bar{a}_h^* \bar{a}_i \bar{a}_j \right) \right], \quad (2) $$

where all the terms $b_{k,h}$ with $h \neq k$ have been neglected. The terms $c_{k,-h,i}^{f,cw}$ are the coefficients of the three-magnon scattering that are due to the rf field. The superscript cw means that the clockwise rotating component of the rf field is considered for its evaluation. In the same way, ccw stands for counter-clockwise (the rf field being linearly polarized, it can be decomposed into its cw and ccw components). The terms $d_{k,h,i,j}^0$ are the coefficients of the four-magnon scattering. The zero superscript indicates that they are derived from conservative torque terms of the LLG. The other terms derived from non-conservative torques have been neglected because of their relatively small value.

## D. P-mode analysis

For small rf power, the driven magnetization oscillations in micromagnetic simulations are time harmonic. The spectrum of such oscillations shows a single peak (SP) at the same frequency of the driving field, which therefore will be termed SP solution.
<figure><img src="image_11.png" /><figcaption>**Supplementary Figure 4** | Phase diagram in the control plane ($H_z$, $h_{rf}$) calculated using the NMM including five modes (1, 6, 15, 30, 51). Instability curves with upwards (downwards) triangles have been obtained by increasing (decreasing) slowly $h_{rf}$ from zero (maximum) amplitude. The dashed curve is the instability threshold deduced from **P**-mode analysis (see Suppl. Fig.3(f)).</figcaption></figure>

Then, we look for solutions of eq.2 of the following type $a_k(t) = \bar{a}_k e^{j\omega_{rf}t}$, where this time $\bar{a}_k$ is a phasor not depending on time. In this way, the reduced NMM becomes a set of algebraic nonlinear equations whose solution gives the phasors of the normal modes. Such a set will be called in the following **P**-mode equation.

For a fixed initial solution and fields values ($H_z$, $h_{rf}$), the solutions of the **P**-mode equation are computed using the Newton-Raphson method. Starting from zero rf field where all modes have zero amplitude, numerical continuation is used and a sequence of solutions ($h_{rf,n}$, $\{\bar{a}_1, \bar{a}_2, \dots\}_n$) for $n = 0, 1, \dots$ is calculated at fixed $H_z$. These solutions, which in time domain represent modes with time harmonic oscillations and fixed phases, are then compared with the amplitudes of normal modes computed by full micromagnetic simulations under the same excitation conditions. For each couple ($H_z$, $h_{rf,n}$), taking the equilibrium as initial condition, a simulation of 1100 ns is performed and the snapshot of the final magnetization field is saved. The duration of the simulation is such that the SP solution is reached for values of $h_{rf}$ where it exists. These snapshots are then projected onto the basis of normal modes, so that their complex amplitudes can be computed. The projection technique is discussed in the section Methods of ref.[3]. The comparison of the real and imaginary parts of normal modes computed using the **P**-mode equation and micromagnetic simulations is shown as a function of $h_{rf}$ for two values of $H_z$ in Fig.3(a)-(d). Starting from $h_{rf} = 0$ mT, the agreement is quantitative up to a certain value ($h_{rf} = 1$ mT), after which there is no more agreement. In fact, the analysis of the spectrum obtained from micromagnetic simulations beyond this value shows that it is not SP anymore, but multi-peaks continuous (MP-C).

This mismatch is investigated further with the reduced NMM. First, eq.2 is linearized around a **P**-mode solution calculated by solving the **P**-mode equation with the continuation method, where the rf field is slowly increased. Then, for each $h_{rf}$ in the range 0 – 4 mT, the eigenvalues of the linear dynamical system obtained are computed. The critical field for the occurrence of the instability process is defined as the lowest $h_{rf}$ where the real part of one of the eigenvalues becomes positive. This procedure is visualized in Suppl. Fig.3(e)-(f). It corresponds to a generalization of the second-order Suhl's instability in the framework of normal modes. One main difference with the Suhl's approach is that the **P**-mode equation permits to consider large-angle steady oscillations before the instability.

The instability curve defined in this way looks qualitatively similar but not in quantitative agreement with the curve found from micromagnetic analysis shown in Fig.2(g) (no rate). So we investigated the dynamics with continuous time simulations: $H_z$ being fixed, $h_{rf}$ is first slowly increased to the point where the SP solution becomes an MP-C solution, and then decreased until the SP solution is recovered. This continuous time analysis has been done both with the micromagnetic code MaGICo and with the reduced NMM. In Fig.2(g) of the main text, the results obtained with MaGICo are shown, while those obtained with the reduced NMM are shown in Suppl. Fig.4. They reveal several interesting facts. First, the instability curve found from micromagnetic simulations where the equilibrium was taken as initial condition is very close to the curve corresponding to the transition MP-C→SP computed from continuous time micromagnetics. Second, the curve corresponding to the transition SP→MP-C, that was a prediction of the **P**-mode analysis, has also been found with micromagnetic simulations. Third, the SP→MP-C curve computed by the reduced NMM reproduces well that computed by micromagnetic simulations. Finally, the reduced NMM captures the coexistence between the SP and MP-C solutions. The region of coexistence between these two solutions differs from the one obtained from full micromagnetics because the NMM with only 5 modes describes the dynamics of the MP-C solution only qualitatively.
## IV. EXPERIMENTAL EVIDENCE OF BISTABLE DYNAMICS

To further demonstrate the bistable dynamics driven by nonlinear FMR, we provide additional experiments to the ones presented in Fig.5 of the main text. Suppl. Fig.5 presents two-tone MRFM spectroscopy measurements performed on a BiYIG nanodisk of diameter 500 nm and their comparison to micromagnetic simulations. The evidence of bistability is achieved by sweeping up and down the applied magnetic field $H_z$ under continuous rf excitation at $f_1 = 5$ GHz and $P_1 = 10$ dBm, and by supplying a second rf tone of low power $P_2 = -11$ dBm and varying frequency $f_2$ to probe the spectrum of the dynamics driven by the main rf tone. To better evidence a possible hysteretic behavior, it is important to sweep $H_z$ at constant $f_2$ rather than the opposite (which was the case in the two-tone MRFM spectroscopy maps shown in Fig.4 of the main text).

Suppl. Fig.5(a) and (c) display the two-tone spectroscopy maps acquired by sweeping $H_z$ upwards and downwards, respectively. The comparison of these two panels demonstrates some hysteresis. For the upward sweep, a multi-peaks spectrum is measured from 0.169 T to 0.1703 T, meaning that in this region, the instability threshold is reached. In contrast, for the downward sweep, the multi-peaks dynamical state is measured from 0.1703 T all the way down to 0.168 T. In other words, for $0.168 < \mu_0 H_z < 0.169$ T, the self-modulation instability is only observed upon decreasing the applied field. Additional evidence of this bistable dynamics is provided by monitoring the frequency $f_c$ of the MRFM cantilever using a phase-locked loop while sweeping up and down the field, which is done simultaneously to the two-tone MRFM spectroscopy, and presented in Suppl. Fig.5(b). One observes that the upward and downward traces of the $f_c$ signal do not superpose for $0.168 < \mu_0 H_z < 0.169$ T, an indication that the static change of longitudinal magnetization induced by the continuous rf field is hysteretic in this region [4]. We note that the specific shape of the $f_c$ signal in Fig.5(b) is due to the cantilever oscillations above the disk [5].

Finally, we performed full micromagnetic simulations which reproduce these experiments by slowly sweeping $H_z$ up and down under continuous rf excitation. The longitudinal reduction of magnetization is shown in Suppl. Fig.5(e) for upward and downward sweeps, while the corresponding power spectra are displayed in panels (d) and (f), respectively. Similarly to the MRFM spectroscopy data, a clear bistable dynamics is observed for $0.1682 < \mu_0 H_z < 0.1695$ T between the SP and MP-C states, respectively achieved in this region for increasing and decreasing $H_z$. We also note the appearance of a second narrower hysteretic region for $0.1706 < \mu_0 H_z < 0.171$ T where these two dynamic states can be excited depending on the field sweep direction.
<figure><img src="image_12.png" /><figcaption>**Supplementary Figure 5** | Experimental evidence of bistable dynamics. (a) and (c) Two-tone MRFM spectroscopy as a function of the upwards and downwards swept field $H_z$, respectively, and the frequency $f_2$ of the low power rf excitation (pulse modulated at −11 dBm). The power of the continuous rf excitation at 5 GHz is fixed to 10 dBm. (b) Simultaneously measured MRFM cantilever frequency vs. $H_z$. (e) Longitudinal reduction of magnetization vs. $H_z$ extracted from micromagnetic simulations ($\mu_0h_{rf}$ = 1 mT at 5 GHz). (d) and (f) Corresponding spectral analysis as a function of the upwards and downwards swept $H_z$, respectively. Red and blue dotted vertical lines indicate limits of the instability region on all panels of the figure.</figcaption></figure>
REFERENCES

[1] http://wpage.unina.it/mdaquino/index_file/MaGICo.html.
[2] M. d'Aquino, C. Serpico, G. Miano, and C. Forestiere, A novel formulation for the numerical computation of magnetization modes in complex micromagnetic systems, J. Comput. Phys. 228, 6130 (2009).
[3] S. Perna, F. Bruckner, C. Serpico, D. Suess, and M. d'Aquino, Computational micromagnetics based on normal modes: Bridging the gap between macrospin and full spatial discretization, J. Magn. Magn. Mater. 546, 168683 (2022).
[4] Y. Li, V. V. Naletov, O. Klein, J. L. Prieto, M. Muñoz, V. Cros, P. Bortolotti, A. Anane, C. Serpico, and G. de Loubens, Nutation spectroscopy of a nanomagnet driven into deeply nonlinear ferromagnetic resonance, Phys. Rev. X 9, 041036 (2019).
[5] I. Ngouagnia Yemeli, Magnetization dynamics of nanostructures in strongly out-of-equilibrium regimes, Theses, Université Paris-Saclay (2022).