Magnonic frequency comb through nonlinear magnon-skyrmion scattering
Zhenyu Wang, H. Y. Yuan, Yunshan Cao, Z.-X. Li, Rembert A. Duine, Peng Yan
MMagnonic frequency comb through nonlinear magnon-skyrmion scattering
Zhenyu Wang , H. Y. Yuan , ∗ Yunshan Cao , Z.-X. Li , Rembert A. Duine , and Peng Yan † School of Electronic Science and Engineering and State Key Laboratory of Electronic Thin Films and Integrated Devices,University of Electronic Science and Technology of China, Chengdu 610054, China and Institute for Theoretical Physics, Utrecht University, 3584 CC Utrecht, The Netherland
An optical frequency comb consists of a set of discrete and equally spaced frequencies and has found wideapplications in the synthesis over broad spectral frequencies of electromagnetic wave and precise optical fre-quency metrology. Despite the analogies between magnons and photons in many aspects, the analogue of opticalfrequency comb in magnonic system has not been reported. Here, we theoretically study the magnon-skyrmioninteraction and find that magnonic frequency comb (MFC) can be generated above a threshold of driving ampli-tude, where the nonlinear scattering process involving three magnons prevails. The mode-spacing of the MFCis equal to the breathing-mode frequency of skyrmion and is thus tunable by either electric or magnetic means.The theoretical prediction is verified by micromagnetic simulations and the essential physics can be generalizedto a large class of magnetic solitons. Our findings open a new pathway to observe the frequency comb structurein magnonic devices, that may inspire the study of fundamental nonlinear physics in spintronic platform in thefuture.
Introduction. —Frequency comb is a spectrum consisting ofa series of discrete and evenly spaced spectral lines. It wasoriginally proposed in optical systems through a pulse traingenerated by a mode-locked laser [1] and later realized in mi-croresonators utilizing the Kerr nonlinearity [2]. The opticalfrequency comb enables the mutual conversion of electromag-netic waves in optical and microwave frequencies and furthermanifests its important role in the high-precision frequencymetrology, such as the optical atomic clock, optical ruler, ex-oplanets search, and molecular spectra sampling [3–6]. Thesuccess of photonic frequency comb also inspires the physi-cists to search for alternates of frequency comb generators.Recently, phononic combs have been theoretical predicted[7] and demonstrated in a microscopic extensional mode res-onators through a three-wave mixing process [8].Magnon is the quantum of collective spin excitations in or-dered magnets and it resembles photon as bosonic quasiparti-cle. The magnon spintronics rises for the long coherent time,low-power consumption and convenient electrical manipula-tion of magnons as information carriers. Many photonic phe-nomenon such as Bose-Einstein condensation, squeezed state,antibunching and multi-particle entanglement have been re-ported in magnonic systems [9–14], and cavity photon andmagnon can be even strongly entangled to achieve coherentinformation transfer [15–22]. Such interplay extends the hori-zon of traditional spintronics and further makes the magnonicsystem a promising candidate for quantum information pro-cessing and quantum computing. However, despite theseanalogies, the analogue of frequency comb in magnonic sys-tem is yet to be realized. The major challenge is that the usualcoupling of a ferromagnetic resonance mode with the mi-crowave falls into the linear regime. Under strong microwavedriving, a ferromagnet, however, will be parametrically ex-cited instead of forming a frequency comb due to the weaknonlinearity in uniform ferromagnets.Recently, it has been shown that magnetic textures can sig-nificantly enhance the nonlinearity in magnetic systems, suchas the three-magnon splitting and confluence processes [23– n r n ω ω ω= + ω − ω − ω ω r ω r ω r ω r ω r ω S p i n w a v e a m p lit ud e ω ω n ω ω r ω Incidentspin waveBreathingskyrmionScatteredspin wave
FIG. 1: Schematic illustration of nonlinear magnon-skyrmion scat-tering and the resulting MFC in a magnetic film. The mode spacingin the frequncy comb is determined by the frequency of skyrmionbreathing mode ω r . a r X i v : . [ c ond - m a t . m e s - h a ll ] F e b Hamiltonian formalism. —Let us consider a chiral magneticsystem in xy plane, described by the following Hamiltonian, H = A ( ∇ m ) − Km z + D [ m z ∇ · m − ( m · ∇ ) m z ] , (1)where m is the normalized magnetization, A is the exchangesti ff ness, K is the easy-axis anisotropy coe ffi cient, and D isthe strength of Dzyaloshinskii-Moriya interaction (DMI). Theinterplay of the exchange interaction, anisotropy and DMIwill stabilize the magnetic skyrmion. To investigate the in-teractions of magnons excited from the ground state, wemake a small perturbation around the skyrmion profile, i.e., m ( x , y ) = m ( x , y ) + δ m ( x , y ), where m ( x , y ) is the staticskyrmion profile while the fluctuation δ m ( x , y ) can be fur-ther quantized as the magnon creation and annihilation op-erators through Holstein-Primako ff transformation [26]. Bysubstituting the ansatz back into Eq. (1), we can rewrite theHamiltonian as H = H (0) + H (2) + H (3) + H (4) + ... , where H (0) , H (2) , H (3) , H (4) are respectively the ground state energy,two, three and four magnon scattering processes [27]. Due tothe lacking of analytical skyrmion profile and the complexityto consider the multi-magnon process around this nonlinearstructure, we do not intend to fully solve this spectrum. In-stead, we observe the main contribution to the three-magnonprocess is H (3) ∝ a k a † k a † k + a k a k a † k + h . c . , where a k ( a † k ) ismagnon annihilation (creation) operator around the skyrmionprofile. Usually, this three-magnon process is weak for di-lute magnon gas excited in a ferromagnetic state. However, alarger magnetization gradient near skyrmion wall can enhancethis nonlinearity [23]. If one mode such as a k is driven witha strong power microwave, the three-magnon process may befurther amplified by the amplitude of the driving (cid:104) a k (cid:105) . Onthe other hand, a strong driving may cause significant defor-mation and fluctuation of the skyrmion profile, such that itsinternal mode including breathing and gyrotropic modes canbe excited. Then the breathing mode, for instance, can hybridwith the driving through this three-magnon process and gen-erate the sum and di ff erence frequency signals, while thesesecondary signals can further hybrid with the breathing modeto generate higher-order frequency modes. This chain reactionwill cause a cascade of the magnonic excitation in the system,and finally result in the generation of MFC with the frequencyspacing around the skyrmion breathing mode, as illustrated inFig. 1.To capture the essential physics of the three-magnon pro-cess, we propose a phenomenological Hamiltonian, H = ω k a † k a k + ω r a † r a r + ω p a † p a p + ω q a † q a q + g p ( a k a r a † p + h . c . ) + g q ( a k a † r a † q + h . c . ) + h ( a k e i ω t + a † k e − i ω t ) , (2)where ω k is the incident spin-wave mode, ω r is the frequencyof breathing mode, p and q are respectively the sum-frequencymode with frequency ω p = ω k + ω r and di ff erence-frequencymode with frequency ω q = ω k − ω r . g p and g q are respectivelythe strength of the three-magnon confluence and splitting. Ingeneral, they should depend on the mode overlap of the three types of magnons. Here we treat them as phenomenologicalparameters. The last term is a microwave driving with ampli-tude h and frequency ω .We can eliminate the time dependence term in (2) by trans-forming into a rotating frame V = exp[ − i ω t ( a † k a k + a † p a p + a † q a q )], and the result is, H = ∆ k a † k a k + ω r a † r a r + ∆ p a † p a p + ∆ q a † q a q + g p ( a k a r a † p + h . c . ) + g q ( a k a † r a † q + h . c . ) + h ( a k + a † k ) , (3)where the detunings are defined as ∆ ν = ω ν − ω ( ν = k , p , q ).Based on Hamiltonian (3), we can readily derive the Heisen-berg dynamic equation of the system as, i da k dt = ( ∆ k − i α k ω k ) a k + g q a r a q + g p a † r a p + h , (4a) i da r dt = ( ω r − i α r ω r ) a r + g q a k a † q + g p a † k a p , (4b) i da p dt = ( ∆ p − i α p ω p ) a p + g p a k a r , (4c) i da q dt = ( ∆ q − i α q ω q ) a q + g q a k a † r , (4d)where α ν ( ν = k , r , p , q ) is the damping rate of the corre-sponding magnon mode. In the strong driving limit, the meanfield approximation can be taken such that (cid:104) u i u j (cid:105) = (cid:104) u i (cid:105)(cid:104) u j (cid:105) ,where u = ( a k , a r , a p , a q , a † k , a † r , a † p , a † q ) T . Then the Heisen-berg equations in the steady condition ( d (cid:104) u (cid:105) / dt =
0) al-ways have a trivial solution: (cid:104) a k (cid:105) = − h / ( ∆ k − i αω k ) and (cid:104) a r (cid:105) = (cid:104) a p (cid:105) = (cid:104) a q (cid:105) =
0. This corresponds to the linear re-sponse of the system, i.e., only one magnon mode matchingthe driving frequency ω is strongly excited. When the driv-ing power increases above a threshold, this trivial solution be-comes unstable and the system will enter into the nonlinearregime.To calculate the driving threshold, we expand the averageof first moments around this trivial solution as a ν = (cid:104) a ν (cid:105) + δ a ν and rewrite the dynamic equations as d δ u / dt = M · δ u , where M is the drift matrix [27]. Based on Routh-Hurwitz criteria[28], all eigenvalues of M should have negative real parts toguarantee the stability of the system. This condition impliesthat at least one eigenvalue of M is purely imaginary at thecritical point, which determines the threshold field h c as, h c = α ω (cid:113) ω + ω ω r − ω ω r − ω r g √ ω − ω r , (5)where we assume g p = g q = g , and all magnon modes haveidentical damping rates α k = α r = α p = α q = α (cid:28) ω (cid:29) ω r , then the threshold can be approx-imated as h c ≈ α ω / √ g . The existence of such a thresholdis the first key prediction of our theoretical formalism.Above the threshold h c , the sum-frequency mode ( ω k + ω r )and di ff erence-frequency mode ( ω k − ω r ) will be excited si-multaneously. By solving Eq. (4), we can analytically derive -2-1012 -2-1012 x m δ -1
1 m µ m µ
2 (GHz) ω π FF T a m p lit ud e ( a r b . un it s ) , ( n m ) ss xy ( n m ) s R , ( n m ) ss xy ( n m ) s R (ns) t (ns) t (a) (b)(c) (d) s x s y z m − + z m δ with SKno SKwith SKno SK FIG. 2: (a) Schematic illustration of a magnetic thin film hosting aN´eel-type skyrmion. (b) The internal spectrum of the magnetic filmwith (solid line) and without skyrmion (dashed line). (c)(d) The timeevolution of skyrmion position and skyrmion radius under the drivingwith frequency ω/ π = . the magnon population as, |(cid:104) a k (cid:105)| = h c αω , |(cid:104) a r (cid:105)| = √ h c ( h − h c ) αω , (6a) |(cid:104) a p (cid:105)| = √ h c ( h − h c ) √ α ( ω + ω r ) , |(cid:104) a q (cid:105)| = √ h c ( h − h c ) √ α ( ω − ω r ) . (6b)We immediately see that the magnon number of the breath-ing mode |(cid:104) a r (cid:105)| will keep accumulating as we increase thedriving field and even surpass the driving magnons. In reality,this trend will be suppressed due to the four-magnon scatter-ing ( H (4) term), resembling the stability of magnon densityin a ferromagnetic state under parallel pumping [29]. Withinmean-field formalism, the four magnon interaction will mod-ify the eigenfrequency ω r as ω r + κ (cid:104) a † r a r (cid:105) , where the coef-ficient κ represents the interaction strength. Now the steadypopulation is approximately to be |(cid:104) a k (cid:105)| ∼ h / ( αω k ) , |(cid:104) a r (cid:105)| ∼ (cid:15)ω r /κ , where the dimensionless parameter (cid:15) is close to zeroas κ → [30]. Numerical verification. —We consider a ferromagnetic thinfilm shown in Fig. 2(a) with dimensions 1000 × × .It can host a N´eel-type skyrmion at the film center, if the inter-facial symmetry is broken by growing a heavy metal layer inadjacent to the magnetic film. The magnetization dynamics is simulated by numerically solving the Landau-Lifshitz-Gilbertequation, ∂ m ∂ t = − γ m × h e ff + α m × ∂ m ∂ t , (7)where γ is gyromagnetic ratio, h e ff = − δ H /δ m is the e ff ectivefield including exchange, anisotropy, dipolar, and driving ACfields, and α is Gilbert damping. The magnetic parameters ofCo [31] are used, i.e., exchange sti ff ness A = . × − J / m,saturation magnetization M s = . × A / m, anisotropy co-e ffi cient K = × J / m , DMI strength caused by interfa-cial symmetry breaking D = . / m , magnetic damping α = − . The mesh size is 2 × × .We first characterize the intrinsic magnon excitation spec-trum of the skyrmion by applying a sinc-function field in thefilm plane, i.e., h ( t ) = h sinc( ω H t ) ˆ x with amplitude h = ff frequency ω H / π =
100 GHz. The inter-nal spectrum of the magnetic system is obtained by carry-ing a standard fast Fourier transform (FFT) for each cell andthen averaging over the whole film. Three main peaks cen-tering at 37.7 GHz, 30.7 GHz, and 8 GHz are found, asshown in Fig. 2(b), while, in the absence of skyrmion, onlythe 37.7 GHz mode appears. Through this comparison, 37.7GHz mode is ascertained to be the ferromagnetic resonance(FMR) mode, whose frequency can be analytically calculatedas ω c / π = γ (2 K − µ M s ) / (2 π M s ) = . h ( t ) = h sin( ω t ) ˆ x with ω/ π =
80 GHzin a narrow rectangular area of the film [black bar in Fig.3(a)]. The magnons excited near the source will propagatetoward the skyrmion located at the film center and interactwith it. To analyze the magnon excitation in the skyrmionregion, we make FFT of the time dependent magnetizationand record the response at various fields. Figure 3(b) showsthe response distribution in frequency space as we tune thedriving amplitude h , and three phases can be classified: (i)Below 100 mT, mainly the magnon mode matching the driv-ing frequency ω/ π =
80 GHz and the FMR mode ( ω c ) areexcited. (ii) In the window from 100 to 130 mT, two mainside-peaks at 80 + =
88 GHz and 80 − =
72 GHz are ex-cited. (iii) Above 130 mT, clear frequency comb with mode-spacing of the breathing-mode frequency 8 GHz emerges, asshown in Fig. 3(c). These features are consistent with the the-oretical predictions. Figure 3(d) shows the quantitative com-parison between theory and simulations for the amplitudes offour main modes at 80 GHz (black dots), 8 GHz (green dots),72 GHz (blue dots), and 88 GHz (red dots), respectively. Thephenomenological model thus well captures the trend of thesefour modes, while the main deviation is on the amplitude ofbreathing mode. This is attributed to the simplified three-mode interaction in the phenomenological formalism, and the ( ) ˆsin h t x ω ω ω c ω r ω ( GH z ) ω π (mT) h µ FF T a m p lit ud e ( a r b . un it s )
2 (GHz) ω π (mT) h µ FF T a m p lit ud e ( a r b . un it s ) r ω r ω c ω c ω c ω ω ω ω ω ω ω
10 mT100 mT200 mT
80 GHz ana. sim. ana. sim. 8 GHz72 GHz88 GHz0.1 − + x m FFT amplitude (arb. units) max (a) (b)(c) (d) FIG. 3: (a) Illustration of magnetic film subjected to a microwavesource applied at the boundary (thick black line). (b) The completeresponse of the system as we tune the driving amplitude ( h ). Thedriving frequency is fixed at 80 GHz. The color codes the amplitudeof the excitation. (c) Response of the system at three representativefields µ h =
10 mT (upper panel), 100 mT (middle panel) and 200mT (lower panel), respectively. (d) The amplitude of the four mainpeaks at ω/ π = , , and 80 ± α = . , ω k / π =
80 GHz , ω r / π = , and g = . × − GHz. breathing mode may be further consumed by generating thehigher-order modes inside the MFC.To study the robustness of the magnonic comb against fluc-tuations of the driving frequency, we vary the microwave fre-quency ω / π ranging from 60 up to 200 GHz. Figure 4(a)shows that the magnonic comb is mostly visible in a windowfrom 67 to 82 GHz, while the mode spacing of the comb isalways equal to the frequency of breathing mode (8 GHz)regardless of the driving frequency. To understand this be-havior, we notice that the window is around the double fre-quency of the FMR mode (2 ω c / π =
75 GHz). This suggeststhat the spins near the skyrmion may be more e ff ectively ex-cited though subharmonic process. Such consideration natu-rally implies that the strength of three magnon-process ( g p , g q )will be frequency dependent. Here we assume that it followsa Gaussian profile as g p ( ω ) = g q ( ω ) = g e − ( ω − ω c ) /σ with g being determined by the field threshold shown in Fig. 3(d).Figure 4(b) shows that this consideration can recover the ex-citation window of frequency comb very well. Meanwhile,the amplitude of four main modes modes ω (black dots), ω r (green dots), ω − ω r (blue dots), ω + ω r (red dots) as a func-
60 70 80 90 10002x10 ( GH z ) ω π
2 (GHz) ω π
2 (GHz) ω π
60 70 80 90 100050100150200 FF T a m p lit ud e ( a r b . un it s ) ana. sim. ω ana. sim. r ω r ω ω− r ω ω+ (a) (b) ω ω r ω c ω FFT amplitude (arb. units) max FIG. 4: (a) Response of the system as a function of the frequency ofthe microwave. The driving amplitude is fixed at 150 mT. (b) Theamplitudes of the four main modes ω , ω ± ω r , ω r as a function ofdriving amplitude. The theoretical curves are calculated resemblingthat in Fig. 3 with a driving frequency dependent coupling g . g = . × − GHz , σ = tion of the driving frequency can be described quantitatively.The breathing mode ω r (green lines) is higher than numericalvalues, which, similar to Fig. 3, is due to the neglect of higher-order modes in our phenomenological formalism. How to de-rive the frequency window from a microscopic theory basedon the overlap of various magnon modes is still an open ques-tion.Note that as we increase the driving frequency furtheraround 3 ω c , the MFC can also appear. However, to excitethese modes, a much stronger microwave is required due to thesmallness of the coupling between high-frequency magnonand skyrmion. As a comparison, we simulate the responseof a ferromagnetic domain to the strong driving at 80 GHzand find that only the FMR mode ( ω c ) and double frequencymode (160 GHz) are excited instead of a frequency comb [27].This highlights the essentiality of spin textures to enhance thethree-magnon process. Discussion and conclusion. —To verify our theoretical pre-dictions in experiments, one can locally inject a spin wave ina chiral magnetic thin film hosting skyrmions, and then mea-sure the spectrum of the excited spin waves near the skyrmionthrough the frequency-resolved technique such as Brillouinlight scattering microscope [32]. This allows the real-timeidentification of the frequency modes in the magnonic comb.Moreover, since the breathing-mode frequency of a skyrmionis inversely proportional to the square of skyrmion size ( R s ),i.e., ω r ∝ R − s [33], we estimate that the mode spacing ina magnonic comb can range from GHz regime for a 10-nmskyrmion down to kilohertz regime for a 1000-nm skyrmion.This complements the working frequency of the optical fre-quency comb with the mode-spacing from 50 MHz to 500GHz [3], that may be very useful for spin-wave calibration.In conclusion, we have shown that a strong microwavedriving can induce a hybridization of the driving mode withthe breathing mode of a magnetic skyrmion, generating amagnonic frequency comb. The threshold of the driving am-plitude and the excitation window of the driving frequencyare well modelled by our Hamiltonian formalism involvingthe sum-frequency and di ff erence-frequency process of threemagnons. We envision that the essential physics should alsohold for other magnetic solitons with low-frequency internalmodes, such as antiskyrmions, domain walls, and vortices.Note that the experimental study of the interaction of spin-wave and domain wall was recently reported [34]. Terahertzfrequency comb induced by nonlinear magnon-solition inter-action in antiferromagnet is also an interesting issue. Our re-sults open a novel avenue to study the frequency comb physicscombining the advantages of magnons and general magnetictextures.This work was funded by the National Natural ScienceFoundation of China (Grants No. 12074057, No. 11604041,No. 11704060, and No. 11904048). Z.W. acknowledges fi-nancial support from the China Postdoctoral Science Founda-tion (Grant No. 2019M653063). R.A.D. has received fund-ing from the European Research Council (ERC) under theEuropean Unions Horizon 2020 research and innovation pro-gramme (Grant agreement No. 725509).Z.W. and H.Y.Y. contributed equally to this work. ∗ Corresponding author: [email protected] † Corresponding author: [email protected][1] Th. Udem, R. Holzwarth, and T. W. H¨ansch, Optical frequencymetrology, Nature , 233 (2002).[2] P. Del’Haye, A. Schliesser, O. Arcizet, T. Wilken, R.Holzwarth, and T. J. Kippenberg, Optical frequency comb gen-eration from a monolithic micoresonator, Nature , 1214(2007).[3] T. Fortier and E. Baumann, 20 years of developments in opticalfrequency comb technology and applications, Commun. Phys. , 153 (2019).[4] A. Pasquazi, M. Peccianti, L. Razzari, D. J. Moss, S. Coen, M.Erkintalo, Y. K. Chembo, T. Hansson, S. Wabnitz, P. Del (cid:48) Haye,X. Xue, A. M. Weiner, and R. Morandotti, Micro-combs: Anovel generation of optical sources, Phys. Rep. , 1 (2018).[5] M. G. Suh, X. Yi, Y. H. Lai, S. Leifer, I. S. Grudinin, G. Va-sisht, E. C. Martin, M. P. Fitzgerald, G. Doppmann, J. Wang,D. Mawet, S. B. Papp, S. A. Diddams, C. Beichman, and K.Vahala, Searching for exoplanets using a microresonator astro-comb, Nat. Photon. , 25 (2019).[6] A. Dutt, C. Joshi, X. Ji, J. Cardenas, Y. Okawachi, K. Luke, A.L. Gaeta, and M. Lipson, On-chip dual-comb source for spec-troscopy, Sci. Adv. , e1701858 (2018).[7] L. S. Cao, D. X. Qi, R. W. Peng, M. Wang, and P. Schmelcher,Phononic Frequency comb through nonlinear resonance, Phys.Rev. Lett. , 075505 (2014).[8] A. Ganesan, C. Do, and A. Seshia, Phononic Frequency combvia intrinsic three-wave mixing, Phys. Rev. Lett. , 033903(2017).[9] T. Nikumi, M. Oshikawa, A. Oosawa, and H. Tanaka, Bose-Einstein Condensation of Dilute Magnons in TlCuCl , Phys.Rev. Lett. , 5868 (2000).[10] S. A. Bender, R. A. Duine, and Y. Tserkovnyak, Elec-tronic Pumping of Quasiequilibrium Bose-Einstein-Condensed Magnons, Phys. Rev. Lett. , 246601 (2012).[11] B. Flebus, S. A. Bender, Y. Tserkovnyak, and R. A. Duine,Two-Fluid Theory for Spin Superfluidity in Magnetic Insula-tors, Phys. Rev. Lett. , 117201 (2016).[12] J. Zhao, A. V. Bragas, D. J. Lockwood, and R. Merlin, MagnonSqueezing in an Antiferromagnet: Reducing the Spin Noise be-low the Standard Quantum Limit, Phys. Rev. Lett. , 107203(2004).[13] H. Y. Yuan and R. A. Duine, Magnon antibunching in a nano-magnet, Phys. Rev. B , 100402(R) (2020).[14] A. Kamra, E. Thingstad, G. Rastelli, R. A. Duine, A. Brataas,W. Belzig, and A. Sudbø, Antiferromagnetic magnons as highlysqueezed Fock states underlying quantum correlations, Phys.Rev. B , 174407 (2019).[15] H. Huebl, C. W. Zollitsch, J. Lotze, F. Hocke, M. Greifenstein,A. Marx, R. Gross, and S. T. B. Goennenwein, High Coopera-tivity in Coupled Microwave Resonator Ferrimagnetic InsulatorHybrids, Phys. Rev. Lett. , 127003 (2013).[16] M. Goryachev, W.G. Farr, D. L. Creedon, Y. Fan, M.Kostylev, and M.E. Tobar, High-Cooperativity Cavity QEDwith Magnons at Microwave Frequencies, Phys. Rev. Applied , 054002 (2014).[17] X. Zhang, C.-L. Zou, N. Zhu, F. Marquardt, L. Jiang, and H.X. Tang, Magnon dark modes and gradient memory, Nat. Com-mun. , 8914 (2015).[18] B. Yao, Y. S. Gui, J. W. Rao, S. Kaur, X. S. Chen, W. Lu, Y.Xiao, H. Guo, K.-P. Marzlin, and C.-M. Hu, Cooperative po-lariton dynamics in feedback-coupled cavities, Nat. Commun. , 1437 (2017).[19] T. Wolz, A. Stehli, A. Schneider, I. Boventer, R. Macˆedo, A. V.Ustinov, M. Kl¨aui, and M. Weides, Introducing coherent timecontrol to cavity magnon-polariton modes, Commun. Phys. ,3 (2020).[20] D. Lachance-Quirion, S. P. Wolski, Y. Tabuchi, S. Kono, K. Us-ami, and Y. Nakamura, Entanglement-based single-shot detec-tion of a single magnon with a superconducting qubit, Science , 425 (2020).[21] H. Y. Yuan, Shasha Zheng, Z. Ficek, Q. Y. He, and M.-H. Yung,Enhancement of magnon-magnon entanglement inside a cavity,Phys. Rev. B , 014419 (2020).[22] H. Y. Yuan, P. Yan, S. Zheng, Q. Y. He, K. Xia, and M.-H.Yung, Steady Bell State Generation via Magnon-Photon Cou-pling, Phys. Rev. Lett. , 053602 (2020).[23] D. N. Aristov and P. G. Matveeva, Stability of a skyrmion andinteraction of magnons, Phys. Rev. B , 214425 (2016).[24] B. Zhang, Z. Wang, Y. Cao, P. Yan, and X. R. Wang, Eaves-dropping on spin waves inside the domain-wall nanochannelvia three-magnon processes, Phys. Rev. B , 094421 (2018).[25] L. K¨ober, K. Schultheiss, T. Hula, R. Verba, J. Fassbender,A. K´akay, and H. Schultheiss, Nonlocal Stimulation of Three-Magnon Splitting in a Magnetic Vortex, Phys. Rev. Lett. ,207203 (2020).[26] T. Holstein and H. Primako ff , Field Dependence of the IntrinsicDomain Magnetization of a Ferromagnet, Phys. Rev. , 1098(1940).[27] See the Supplementary at http: // link.aps.org / supplemental / fora detailed derivation of the three-magnon interaction, the ex-plicit form of the drift matrix, the exciplicit form of a self-consistent equations, generation of magnonic frequency combunder microwave driving with higher frequencies, absence offrequency comb in a uniform ferromagnetic state and spin waveprofiles and skyrmion motion in the MFC, which includes Refs.[26, 28].[28] E. X. DeJesus and C. Kaufman, Routh-Hurwitz criterion in the examination of eigenvalues of a system of nonlinear ordinarydi ff erential equations, Phys. Rev. A , 5288 (1987).[29] Z. Haghshenasfard and M. G. Cottam, Quantum statistics andsqueezing for a microwave-driven interacting magnon system,J. Phys.: Condens. Matter , 045803 (2017).[30] A. Vansteenkiste, J. Leliaert, M. Dvornik, M. Helsen, F. Garcia-Sanchez, and F. B. V. Waeyenberge, The design and verificationof MuMax3, AIP Adv. , 107133 (2014).[31] J. Sampaio, V. Cros, S. Rohart, A. Thiaville, and A. Fert, Nucle-ation, stability and current-induced motion of isolated magneticskyrmions in nanostructures, Nat. Nanotech. , 839 (2013). [32] T. Sebastian, K. Schultheiss, B. Obry, B. Hillebrands, and H.Schultheiss, Micro-focused Brillouin light scattering: imagingspin waves at the nanoscale, Front. Phys. , 35 (2015).[33] V. P. Kravchuk, D. D. Sheka, U. K. R¨oßler, J. van den Brink,and Y. Gaididei, Spin eigenmodes of magnetic skyrmions andthe problem of the e ff ective skyrmion mass, Phys. Rev. B ,064403 (2018).[34] J. Han, P. Zhang, J. T. Hou, S. A. Siddiqui, and L. Liu, Mutualcontrol of coherent spin waves and magnetic domain walls in amagnonic device, Science366