Dynamical phase-field model of cavity electromagnonic systems

One main goal of the cavity electromagnonics is to realize strong and dynamically tunable coupling between magnons (quanta of spin waves) and cavity photons (quanta of confined electromagnetic waves)^1,2,3, with application potential in quantum storage^4,5, quantum transduction^6,7 and quantum sensing⁸. The strong coupling between the Kittel mode magnon (spatially uniform precession of magnetization) and the cavity photon was theoretically predicted by Soykal and Flatté^9,10 and experimentally observed in hybrid systems that involve yttrium iron garnet (YIG) bulk crystals ^11,12, permalloy thin-film stripe¹³, YIG film^14,15 mounted on a coplanar microwave resonator, or YIG bulk crystal sphere(s)/slab inside a three-dimensional (3D) microwave cavity^{4,16,17,18,19,20,21,22,23,24,25}. One key feature of such strong coupling is the mode frequency splitting with an avoided crossing in the frequency spectrum, which indicates the hybridization of magnon and photon into a new quasiparticle called magnon polariton^14,19,26. In the time-___domain, the energy of magnon polaritons is constantly exchanged between the magnon and the photon system with 100% conversion efficiency.

To realize practical quantum operation such as mode swapping and storage²⁷, it is necessary to dynamically control the exchange process between the two hybrid modes of magnon polaritons upon the completion of transferring a single quantum of excitation². For example, Floquet engineering²⁸—which herein refers to the simultaneous application of a periodic driving magnetic field—has been successfully implemented to in situ control the transition between the two hybrid modes and even induce further splitting of each mode into two energy levels associated with different Floquet modes²¹, analogous to the Autler-Townes splitting in atomic physics.

In addition to the studies on magnon-photon resonance, tripartite coupling among the photons, magnons, and phonons have also been demonstrated experimentally in a cavity electromagnonic system^{15,22,23,24,25}, which can also be called a cavity electro-magno-mechanical system in this case. For example, Zhang et al. reported a resonant coupling among the Kittel mode magnons, cavity photons, and high-overtone bulk acoustic phonons—all having the same frequency of a few gigahertz (GHz)—in a Gd₃Ga₅O₁₂(GGG, substrate)/YIG(film, 200-nm-thick) mounted on a split-ring resonator¹⁵. In a 0.25-mm-diameter YIG sphere placed in a 3D photon cavity, Zhang et al. ²² demonstrated a coherent coupling between a GHz magnon polariton (with a frequency ω₊ or ω_-) and a megahertz (MHz) acoustic phonon (frequency: ω_p) by parametrically driving the cavity with a strong microwave signal at a frequency ω_d, with ω_d-ω_-=-ω_p or ω_d-ω₊=ω_p.

The main objective of this article is to report a 3D dynamical phase-field model that enables simulating and predicting the coupled dynamics of photons, magnons, and acoustic phonons in a cavity electromagnonic system comprised of a magnon/phonon resonator placed in a bulk 3D photon cavity, which is one of the most used structures in experiments^{4,16,17,18,19,20,21,22,23,24,25}. In contrast to the fact that Hamiltonian-based theoretical analyses (e.g. ^21,29) need to take the mode coupling strength as the input and are therefore not predictive, the present dynamic phase-field model allows for predicting the spatiotemporal evolution of coupled modes in 3D photon cavity and magnon resonators of arbitrary size and geometry under various operation conditions, using only the fundamental materials parameters as the input. Therefore, it can be used to guide the design of cavity structure and control conditions for realizing desirable quantum operation.

The dynamically evolving physical parameters in a cavity electromagnonic system include the magnetic-field component (H^EM) of the electromagnetic (EM) wave in the microwave cavity, the magnetization (M) and elastic strain (ε) of the magnon resonator (e.g., YIG). The propagation of H^EM is governed by the Maxwell’s equations, while dynamical evolution of M and ε, which represent the magnon and phonon subsystem, are usually described by the Landau-Lifshitz-Gilbert (LLG) equation and elastodynamic equation, respectively. Crucially, the dynamics of M is modulated by the H^EM via the Zeeman torque, while the M and ε are coupled via the magnetoelastic interaction³⁰. Therefore, a complete, direct numerical simulation of the dynamical processes in a cavity electromagnonic system requires the simultaneous solution of the coupled LLG, elastodynamic, and Maxwell’s equations.

Thus far, there are only a few advanced computational models that include coupled dynamics of M and EM wave^31,32,33,34 but excludes either the exchange coupling field (i.e., macrospin approximation) in the LLG equation^31,33,34 or the displacement current in the Maxwell’s equations³². Recently, models that include coupled dynamics of ε, M, and EM wave have also appeared^{35,36,37,38,39}, but these models are limited to 1D^36,37,38 or 2D³⁹ system or employ the Newton’s equation³⁵ as a simplification of the elastodynamic equation. Furthermore, these models^{35,36,37,38,39} have not yet been applied to a cavity electromagnonic system. The present dynamical phase-field model addresses the coupled dynamics of ε, M, and EM wave in a 3D cavity electromagnonic system by solving the coupled LLG, elastodynamic, and Maxwell’s equations (see “Methods”). All numerical solvers are accelerated by graphics processing unit (GPU) to increase the computation throughput. As examples of application, we use the dynamical phase-field model to simulate the dynamics of excitation and control of magnon polariton modes in a cavity electromagnonic system comprised of a YIG magnon resonator placed in a 3D photon cavity. Typical coherent gate operations including Rabi oscillation and Ramsey interference are computationally demonstrated. Furthermore, we design a cavity electro-magno-mechanical system, which contains a bilayer YIG/SiN membrane placed in a 3D photon cavity and permits a resonant interaction between the magnon polaritons and the acoustic phonons. We then use the dynamical phase-field model to simulate the coupled mode dynamics under such triple phonon-magnon-photon resonance condition.

Simulation system set-up

Figure 1a schematically shows the cavity electromagonic system. YIG, which has been widely used in hybrid magnonic systems^{16,17,18,22,40,41,42,43,44} due to its ultralow magnetic damping, is used as the magnon resonator. The 3D microwave cavity has a dimension of 45×9×21 mm³, which supports the TE₁₀₁ mode of the standing EM waves with a frequency ω_c/2π = 7.875 GHz. To excite the TE₁₀₁ cavity mode, a point charge current pulse J^c(t) in the form of a Gaussian function \(t{e}^{-{t}^{2}/2{\sigma }_{0}^{2}}\) is applied along the y axis (i.e., only the \({J}_{y}^{{\rm{c}}}\) component is nonzero) at the position (22.5 mm, 1.5 mm, 10.5 mm) of the cavity, where σ₀ is a free parameter that controls the pulse duration and chosen to be 70 ps so that the frequency window of the pulse covers the ω_c. The simulated H^EM has a vortex-like distribution in the xz plane, as shown in Fig. 1a. The YIG resonator is placed at 1.5 mm below the top surface center of the cavity, where the magnitude of H^EM is relatively large. At the initial equilibrium state, the magnetization in the YIG m⁰ is along +z ([001]) due to a bias magnetic field \({{\bf{H}}}^{{\rm{bias}}}=(\mathrm{0,0},{H}_{z}^{{\rm{bias}}})\) applied along the same direction. Note that energy dissipation of both the cavity photon (arising from the imaginary component of the relative dielectric permittivity tensor ε_r) and magnon (arising from the effective magnetic damping) are both set to be zero to study the magnon-photon coupling under the most ideal situation.

a Hybrid magnon-photonic system that contains a YIG cube (not to scale) inside a 3D microwave cavity. The bias magnetic field H^bias is applied along the +z direction. A Gaussian-shaped charge current pulse J^c(t) is injected into the cavity to excite the cavity mode of the standing EM wave. The vectors indicate the direction of the local microwave magnetic field H^EM (TE₁₀₁ mode), where the vector length is proportional to magnitude of the H^EM. b Dynamics of the on-resonance Kittel mode magnon, represented by Δm_x= m_x(t)-m_x(t = 0) and cavity photon, represented by the microwave electric field component \({E}_{y}^{{\rm{EM}}}\) at the detection point which is at 1.5 mm above the bottom surface center of the cavity. Trend lines showing the evolution of the amplitudes of the two modes are added. The values of the simulated \({E}_{y}^{{\rm{EM}}}\) are multiplied by \(\sqrt{{\varepsilon }_{r}}\) to show the electric fields in the cavity of the original size. c, Spatial distribution of local microwave electric field E^EM and the polar plot of the magnitude of the precessing magnetization component Δm_x and Δm_y at (from left to right) Δt = 2.11 ns, 4.4 ns, and 5.71 ns, respectively. The circles indicate |Δm | = 0.01 (innermost), 0.02, 0.03, and 0.04 (outermost), respectively. d Frequency spectrum of the Δm_x(t) in (c). e Simulated mode splitting spectra of the magnon polaritons as a function of ΔH^bias under different YIG sizes of 0.4×0.4×0.4 mm³ (red circles), 1 × 1 × 1 mm³ (green circles) and 2 × 2 × 2 mm³(blue circles), and their analytical fitting curves (lines). ΔH^bias=H^bias- \({{\bf{H}}}_{0}^{{\rm{bias}}}\), where \({{\bf{H}}}_{0}^{{\rm{bias}}}=0.2315\) MA/m (~0.291 T) is the bias magnetic field that ensures magnon-photon on-resonance (ω_m = ω_c).

Full size image

The entire system is discretized into three-dimensional (3D) computational cells with a cell size ∆x = ∆y = ∆z = 2 nm. Numerically, a nm-scale computational cell is necessary to ensure spatial uniformity of the magnetization (i.e., the formation of Kittel mode magnon) via the H^exch (proportional to |∇m | ²) between neighboring spins. Moreover, a basis for any micromagnetic simulation is that the cell size needs to be smaller than the exchange length \({l}_{{\rm{ex}}}=\sqrt{{A}_{\text{ex}}/(0.5{\mu }_{0}{M}_{\text{s}}^{2})}\), which is about 16.3 nm for YIG with an exchange coupling coefficient A_ex of 3.26 pJ/m and M_s = 140 kA/m⁴⁵. However, discretization of a 45×9×21 mm³ system using a cell size of ∆x = ∆y = ∆z = 2 nm would lead to a total of about 10²¹ cells, which is computationally unaffordable. To address the issue, the ε_r is tuned to scale down the EM wavelength and hence the size of the microwave cavity. For example, in the case of 1 × 1 × 1 mm³ YIG cube, we set all three diagonal components of the ε_r to be 6.25 × 10¹⁰ for both the YIG cube and the microwave cavity, thus the EM wavelength is scaled down by 2.5 × 10⁵ \((=\sqrt{{\varepsilon }_{{\rm{r}}}})\) times. Accordingly, the size of the microwave cavity can be reduced from 45 × 9 × 21 mm³ to 180×36×84 nm³ (i.e., 68,040 cells) without changing the spatial profile and the frequency of the TE₁₀₁ mode EM wave. Meanwhile, the size of the 1 × 1 × 1 mm³ YIG cube should be scaled down to 4 × 4 × 4 nm³ to maintain a constant volume ratio of the YIG cube to microwave cavity. Although such size down-scaling makes it not possible to simulate the high-order magnon modes (spatially non-uniform precession of local magnetization) that may occur in a mm-scale YIG, it would not influence the present work on the interaction between Kittel mode magnon and cavity photon. Importantly, although the larger ε_r leads to a smaller E^EM, the magnitude of H^EM, which interacts with the magnetization, remains constant (see Eqs. (6)–(7) in “Methods”). As a result, the simulated coupled magnon-photon dynamics remains the same as that in the original mm-scale system. Furthermore, the use of a larger ε_r allows using a larger time step which significantly reduces the computation time in long-term dynamics simulation (see Methods).

Magnon-photon coupling

To demonstrate the validity and high numerical accuracy of our computational model, we first simulate the formation of the commonly observed k = 0 mode magnon polariton (k is wavenumber), which features the hybridization between the k = 0 (Kittel) mode magnon \(\hat{m}\) and k ≈ 0 mode cavity photon \(\hat{c}\). As illustrated in Fig. 1a, the H^EM, which is perpendicular to the initial equilibrium magnetization m⁰ (see Fig. 1a), drives the magnetization precession. Due to the exchange coupling, all local magnetization vectors m in the YIG precess in phase, resulting in the excitation of the desirable Kittel mode magnon. Since H^EM is largely uniform around the YIG cube (i.e., wavenumber k ≈ 0) and the magnon-photon interaction time is sufficiently long in the present 3D cavity, the k = 0 mode magnon polariton should form if the angular frequency of the Kittel mode magnon ω_m, or the ferromagnetic resonance (FMR) frequency, can be magnetically tuned to match the angular frequency of the cavity ω_c. Specifically, for an initial equilibrium magnetization along [001], one has \(\frac{{\omega }_{{\rm{m}}}}{2{\rm{\pi }}}=\gamma ({H}_{z}^{{\rm{bias}}}-\frac{1}{3}{M}_{s}+\frac{{2K}_{1}}{{\mu }_{0}{M}_{s}})\)³⁶, where γ = 27.86 GHz/T is the gyromagnetic ratio and K₁ = 620 J/m³ is the magnetocrystalline anisotropy coefficient of the YIG. Accordingly, H^bias = (0, 0, 0.291 T) is applied to have ω_m/2π=ω_c/2π = 7.875 GHz.

Figure 1b shows the dynamics of the Kittel magnon mode and the photon, where Δt = 0 refer to the moment at 20 ns after the injection of the Gaussian-shaped current pulse J^c(t) at t = 0 ns. Typical behavior of coherent beating oscillation similar to a two-level system^37,46 is observed. Specifically, the peak amplitudes of the two modes (indicated by the trend lines) show a Rabi-like oscillation¹⁷ with a period of ~6.9 ns (frequency ~145 MHz), suggesting a back-and-forth energy transfer between the YIG and the microwave cavity. Note that we focus on the peak amplitudes rather than the instantaneous values of \(\Delta {m}_{x}\) and \({E}_{y}^{{\rm{EM}}}\), because the energy of the magnon mode is swapped instantaneously between \(\Delta {m}_{x}\) and \(\Delta {m}_{y}\) while the energy of the cavity photon mode is swapped instantaneously between \({E}_{y}^{{\rm{EM}}}\) and \({H}_{x}^{{\rm{EM}}}\). For clearer illustration, Fig. 1c shows the simulated magnon state and spatial distribution of the radiation electric field E^EM(t) at a few representative moments. As shown in the left panel, when the E^EM reaches its peak amplitude, the magnetization m aligns almost along its initial direction [001], indicating an almost zero free energy change in the YIG. As the energy is being transferred from the cavity to the YIG, the amplitude of E^EM in the cavity decreases while the amplitude of the precessing magnetization (or) \(|\Delta {\bf{m}}|=\sqrt{{(\Delta {m}_{x})}^{2}+{(\Delta {m}_{y})}^{2}}\) increases, as shown in the middle panel. After a half period of the energy transfer (3.45 ns = 6.9 ns/2), almost all the EM wave energy is absorbed by the YIG, which is indicated by the negligibly small E^EM in the cavity and relatively large |Δm | , as shown in the right panel of Fig. 1c. Figure 1d shows the frequency spectrum of the temporal waveform Δm_x(t) in Fig. 1b, which reveals two peak frequencies at 7.8 GHz and 7.945 GHz, respectively. The two peak frequencies are symmetric with respect to the ω_m/2π=ω_c/2π = 7.875 GHz with a frequency gap of 145 MHz, indicating the formation of magnon polariton with two different hybrid modes \({\hat{d}}_{+}\) and \({\hat{{d}}}_{-}\). The frequency gap (denoted as δ_D) is consistent with the frequency of the Rabi-like oscillation and defines a magnon-photon coupling strength g_cm = δ_D/2 = 2π×72.5 MHz. It is worth remarking that the magnon subsystem of the YIG is also coupled to the phonon subsystem, because the precessing m generates dynamical strain via the magnetoelastic feedback and the dynamical strain in turn modulates the dynamics of m via the H^mel (see Part 1 of the "Methods"). Moreover, the stiffness damping coefficient β in the elastodynamic equation creates an additional channel for energy dissipation. However, in the system shown in Fig. 1a, the energy exchange between magnon and phonon subsystems is negligible because the magnitude of the dynamical strain is negligibly small ( ~ 10⁻⁷) due to the relatively small magnetoelastic coupling constant of the YIG.

Figure 1e shows the numerically simulated mode frequencies (indicated by hollow circles) as functions of the bias magnetic field in three different hybrid systems where the sizes of the YIG cube are 0.4 × 0.4 × 0.4 mm³, 1 × 1 × 1 mm³ and 2 × 2 × 2 mm³, respectively, and the cavity size remains to be 45 × 9 × 21 mm³. The goal of these simulations is to computationally verify the theoretical relation of \({g}_{{\rm{cm}}}={g}_{0}\sqrt{N}\) (ref. ¹⁶), where g₀ is the coupling strength of a single Bohr magneton to the cavity; N is the total spin number in the YIG and increases linearly with its size. The presence of avoided crossings in all three systems indicate the formation of magnon polaritons. The corresponding magnon-photon coupling strength g_cm can be extracted from the frequency gap under on-resonance \(\hat{m}\) and \(\hat{c}\) modes (where ΔH^bias = 0), which are 2π×18.3 MHz, 2π × 72.5 MHz, and 2π × 182.2 MHz, respectively. One can evaluate that the g_cm is largely proportional to the square root of the cube size and hence the \(\sqrt{N}\). The g_cm in the case of 2 × 2 × 2 mm³ is smaller than the value of 2π × 204.6 MHz obtained from linear extrapolation due to the <100% spatial mode profile overlapping, which is consistent with experimental observation¹⁷. Based on the extracted g_cm, we fit the simulations results via the expression \({\omega }_{\pm }=\frac{1}{2}({\omega }_{{\rm{m}}}+{\omega }_{{\rm{c}}})\pm \frac{1}{2}\sqrt{{({\omega }_{{\rm{m}}}-{\omega }_{{\rm{c}}})}^{2}+4{g}_{{\rm{cm}}}^{2}}\), which describe the angular frequencies of the two hybrid modes \(({\hat{d}}_{+}{\rm{and}}\,{\hat{d}}_{-})\) in a two-level system⁴⁷. The excellent fitting indicates the validity of our model set-up and high numerical accuracy of our dynamical phase-field model.

Floquet-induced magnonic Autler-Townes splitting

Based on the cavity electromagnonic system with a YIG resonator of 1 × 1 × 1 mm³, a periodic dynamical magnetic field h_D(t) is applied along the same axis (z) with the H^bias to implement Floquet engineering. The h_D(t)=|h_D|sin(ω_Dt) is applied uniformly on the YIG resonator, where |h_D| and ω_D are the amplitude and the angular frequency of the Floquet drive, respectively. Figure 2a shows the frequency spectra of the Δm_x(t) simulated under a fixed |h_D| of 2000 A/m but different ω_D/2π varying from 0 to 400 MHz. A static bias magnetic field of 0.2315 MA/m (~0.291 T) is applied to have the magnon mode on resonance with the cavity photon and generate two magnon polariton modes \({\hat{d}}_{+}\) and \({\hat{d}}_{-}\), which have frequencies of 7.945 GHz and 7.8 GHz, respectively. Each polariton mode has several sidebands created by the Floquet drive. The frequencies of these sidebands either increase or decrease as the ω_D increases. When ω_D is equal to the (cavity and magnon) mode splitting frequency ( = 2g_cm = 2π×145 MHz), the first inner sideband of the magnon polariton mode \({\hat{d}}_{+}({\hat{d}}_{-})\) resonantly interact with the other magnon polariton mode \({\hat{d}}_{-}({\hat{d}}_{+})\). As a result, the two energy levels corresponding to the \({\hat{d}}_{+}\) and \({\hat{d}}_{-}\) modes split into four energy levels associated with different Floquet modes (as illustrated in the inset), where the frequency gap between the two newly split energy levels is denoted as Δω_AT, with Δω_AT/2π = 34.8 MHz. Such magnonic Autler-Townes splitting²¹ and the onset of avoided crossing indicate the realization of strong coupling between the two hybrid modes of the magnon polaritons by Floquet drive.

a Frequency spectrum of the magnon polariton as a function of Floquet driving frequency f_D, obtained from (a), dynamical phase-field simulations and (b), Hamiltonian-based theoretical calculations of the absorption spectrum where square root for each data point is take to make the higher order sideband visible. The inset shows the energy level diagram. Evolutions of the magnon (represented by Δm_x, upper panel) and the photon (represented by E_y, lower panel) at the detection point under (c), ω_D/2π = 300 MHz, and (d), ω_D/2π = 145 MHz. The magnitude of the Floquet driving field |h_D | =2000 A/m in (a–d).

Full size image

To gain further insights on the simulated spectrum, we analytically calculate the absorption spectrum of the magnon mode based on Floquet theory described in ref. ²¹, using the numerically simulated on-resonance frequency ω₀ (=ω_m = ω_c) and the magnon-photon coupling strength g_cm as the inputs. As shown in Fig. 2b, the main structure of the calculated spectrum (shown in Fig. 2b) reproduces the simulation results. The calculated spectrum can be understood by writing the Floquet Hamiltonian as following (details of derivation are in ref. ²¹)

Here \({\omega }_{\pm }={\omega }_{0}\pm {g}_{{\rm{cm}}}\), \(\Omega =\gamma |{{\bf{h}}}_{{\rm{D}}}|\), and ℏ is the reduced Planck’s constant. A list of symbols for various modes and related quantities is provided in Supplementary Note 1. The Floquet drive creates a series of sidebands of the \({\hat{d}}_{+}\) modes at frequencies \({\omega }_{\pm }\pm {n\omega }_{{\rm{D}}}\), where n is the sideband order. The last term on the right-hand side of Eq. (1) describes the interactions between different sidebands, where J_n is the nth Bessel function of the first kind. From Eq. (1), one can determine that the coupling strength between the first inner sideband of \({\hat{d}}_{-}\) and the \({\hat{d}}_{+}\) mode (i.e., the magnonic Autler-Townes splitting Δω_AT) is approximated as \(\Delta {{\omega }}_{{\rm{AT}}}\approx {2g}_{{\rm{cm}}}{J}_{1}\left(\Omega /{\omega }_{{\rm{D}}}\right)\). Plugging in the numbers yields a theoretically predicted value of \(\Delta {{\omega }}_{{\rm{AT}}}/2{\rm{\pi }} \sim 34.1\,{\rm{MHz}}\), which is in good agreement with the simulated value of 34.8 MHz. Moreover, since the sum in Eq. (1) only involves odd terms, some of the sidebands (e.g., the first inner sideband of \({\hat{d}}_{+}\) and \({\hat{d}}_{-}\) modes) are not directly coupled, which is revealed by the crossing in both the simulated and calculated spectrum. Detailed discussion on this point can be found in ref. ²¹.

As \(|{{\bf{h}}}_{{\rm{D}}}|\) increases, Δω_AT varies in an oscillatory fashion but always remain nonzero (see Supplementary Fig. S1a in Supplementary Note 2). This trend cannot be explained by the analytical approximation \(\Delta {{\omega }}_{{\rm{AT}}}\,\approx\, {2g}_{{\rm{cm}}}{J}_{1}\left(\Omega /{\omega }_{{\rm{D}}}\right)\), which is only valid when \(|{{\bf{h}}}_{{\rm{D}}}|\) is relatively small. At large \(|{{\bf{h}}}_{{\rm{D}}}|\), the \(|{{\bf{h}}}_{{\rm{D}}}|\)-dependent Δω_AT can be better quantified by the following analytical expression (see detailed derivation in Supplementary Note 3),

which predicts a repeated occurrence of \(\Delta {{\omega }}_{{\rm{AT}}}={2g}_{{\rm{cm}}}\) at \({J}_{0}=0\) and the presence of a Floquet ultrastrong coupling regime²¹ where \(\Delta {\omega }_{{\rm{AT}}} > {2g}_{{\rm{cm}}}\) and \({J}_{0} < 0\). Both features are shown in the frequency spectrum of the magnon polaritons (Supplementary Fig. S1a) obtained by dynamical phase-field simulations.

Our dynamical phase-field simulation results in Fig. 2c further shows the temporal profile of the Δm_x(t) and \({E}_{y}^{{\rm{EM}}}(t)\) for ω_D/2π = 300 MHz. According to the frequency spectra in Fig. 2a, the magnon polariton is still dominated by the intrinsic hybrid modes \({\hat{d}}_{+}\) and \({\hat{d}}_{-}\) with no magnonic Autler-Townes splitting. Correspondingly, the amplitudes of both the ∆m_x(t) and \({E}_{y}^{{\rm{EM}}}(t)\) display a Rabi-like oscillation with a period of 6.9 ns, which is the same as in Fig. 1b. By comparison, for ω_D/2π = 145 MHz where the magnonic Autler-Townes splitting is prominent, the corresponding temporal profiles of ∆m_x(t) and \({E}_{y}^{{\rm{EM}}}(t)\), as shown in Fig. 2d, are clearly composed of components of more than two frequencies. Specifically, there are four major frequency components at \({\omega }_{+}+0.5\Delta {\omega }_{{\rm{AT}}},{\omega }_{+}-0.5\Delta {\omega }_{{\rm{AT}}},{\omega }_{-}+0.5\Delta {\omega }_{{\rm{AT}}}\), and \({\omega }_{-}-0.5\Delta {\omega }_{{\rm{AT}}}\), respectively, corresponding to the four split energy levels as shown in Fig. 2a. Despite the more complex temporal profile, the evolution of the peak amplitudes of ∆m_x(t) and \({E}_{y}^{{\rm{EM}}}(t)\) are still complementary, indicating that the back-and-forth energy exchange still occurs between the Kittel magnon mode \(\hat{m}\) and the cavity photon mode \(\hat{c}\). The beam-splitter type coupling between the \({\hat{d}}_{+}({\hat{d}}_{-})\) and the first inner sideband of the \({\hat{d}}_{-}({\hat{d}}_{+})\) mode (Fig. 2a) can also be interpreted as the energy exchange (i.e., Rabi-like oscillation) between the two energy levels (i.e., the hybridized modes \({\hat{d}}_{+}\) and \({\hat{d}}_{-}\)) through the frequency matching provided by the Floquet drive, i.e., \({\omega }_{+}={\omega }_{-}+{\omega }_{{\rm{D}}}\). Similar Floquet-driven Rabi-like oscillation between two hybridized modes has also been demonstrated experimentally in a two-level photonic system²⁷.

Dynamical control of the energy exchange rate between two magnon polariton modes

The Rabi-like oscillation between the two magnon polariton modes \(({\hat{d}}_{+}{\space}{\rm{and}}\,{\hat{d}}_{-})\) corresponds to a rotation along the real axis of a Bloch sphere, as illustrated in Fig. 3a inset. Here, we model this process in the time ___domain and computationally demonstrate the dynamical control of the energy exchange rate between the \({\hat{d}}_{+}\) and \({\hat{d}}_{-}\) modes (namely, Rabi flopping frequency) by dynamically varying the amplitude |h_D| of the Floque drive. To this end, the system is initialized by pumping a 15-cycle sinusoidal charge current pulse J^c(t)=J₀sin(ωt) at ω = 2π×7.945 GHz along the y axis (i.e., only the \({J}_{y}^{{\rm{c}}}\) component is nonzero) to populate (excite) the \({\hat{d}}_{+}\) mode. A continuous Floquet drive h_D(t) at the cavity-magnon mode splitting frequency (ω_D/2π = δ_D/2π = 145 MHz) is then applied to the system. Figure 3a shows the temporal evolution of the magnetization amplitude of the \({\hat{d}}_{+}\) mode magnon polariton \(({\rm{denoted}}\; {\rm{as}}\,{|\Delta {\bf{m}}|}_{{\hat{d}}_{+}})\) with |h_D | =5000 A/m. Here \({\left|\Delta {\bf{m}}\right|}_{{\hat{d}}_{+}}(t)\) is obtained by first extracting the temporal evolution of \(\Delta {m}_{x}(t)\) and \(\Delta {m}_{y}(t)\) by performing inverse Fourier transform for 7.945-GHz (±50 MHz) peak in their frequency spectra and then calculating its magnetization amplitude via \(\sqrt{{(\Delta {m}_{x})}^{2}+{(\Delta {m}_{y})}^{2}}\). The \({\left|\Delta {\bf{m}}\right|}_{{\hat{d}}_{+}}(t)\) displays a Rabi-like oscillation with a period of 11.6 ns, corresponding to a flopping frequency of 86.2 MHz. We note that the \({\hat{d}}_{-}\) mode magnon polariton (at 7.8 GHz) was also excited after the initial current pulse injection, and the dynamics of \({\left|\Delta {\bf{m}}\right|}_{{\hat{d}}_{-}}(t)\) complements the \({\left|\Delta {\bf{m}}\right|}_{{\hat{d}}_{+}}(t)\), as shown by Supplementary Fig. S2 in Supplementary Note 4, suggesting a dynamical energy exchange between the \({\hat{d}}_{+}\) and \({\hat{d}}_{-}\) modes.

a Dynamics of the magnetization amplitude of the \({\hat{d}}_{+}\) mode magnon polariton under a continuous Floquet driving h_D(t) with |h_D | =5000 A/m and ω_D/2π = 145 MHz. Δt = 0 is the moment when the application of h_D(t) begins after the current pulse J^c(t) injection is complete. The \({\hat{d}}_{+}\) and \({\hat{d}}_{-}\) modes swap at the frequency of Δω_AT/2π along the real axis of the Bloch sphere (inset). b Dynamics of the magnetization amplitude of the \({\hat{d}}_{+}\) mode magnon polariton under different |h_D| but the same frequency of ω_D/2π = 145 MHz.

Full size image

Figure 3b further summarizes the \({\left|\Delta {\bf{m}}\right|}_{{\hat{d}}_{+}}(t)\) simulated under different |h_D| varying from 1000 to 5000 A/m. As shown, the period of Rabi-like oscillation decreases as |h_D| increases, leading to an increase in the corresponding Rabi flopping frequency (f_Rabi). By analytically solving the Heisenberg equation of the \({\hat{d}}_{+}\) and \({\hat{d}}_{-}\) modes under the rotating wave approximation (RWA), we obtain the mode square \({|{\hat{d}}_{+}\left(t\right)|}^{2}={\cos }^{2}\left(\frac{\Omega }{4}t\right)\) (see Supplementary Note 5), yielding \({f}_{{\rm{Rabi}}}=\Omega /2=\gamma \left|{{\bf{h}}}_{D}\right|/2\), which is only valid when ω_D=2g_cm = 2π×145 MHz. It is worth noting that Fig. 3a shows a nonzero offset, which differs from the RWA-based prediction and may be due to the existence of other magnon modes. As shown in the inset of Fig. 3b, the analytically calculated f_Rabi agrees well with the values extracted from the simulated \({\left|\Delta {\bf{m}}\right|}_{{\hat{d}}_{+}}(t)\) with only small deviation at larger \(\left|{{\bf{h}}}_{D}\right|\) values, where the system’s behavior deviates from the RWA.

As \(\left|{{\bf{h}}}_{D}\right|\) further increases, it is no longer possible to numerically extract the \({\left|\Delta {\bf{m}}\right|}_{{\hat{d}}_{+}}(t)\) and hence the f_Rabi because multiple magnon modes coexist at 7.945 GHz. However, it is expected that that the variation trend of f_Rabi with \(\left|{{\bf{h}}}_{D}\right|\) would deviates significantly from the RWA-predicted linear relation. Specifically, it is reasonable to speculate that the \({f}_{{\rm{Rabi}}}{\approx }\Delta {{\omega }}_{{\rm{AT}}}/2{\rm{\pi }}\) because the present Rabi oscillation is based on beam-splitter type coupling between the \({\hat{d}}_{+}({\hat{d}}_{-})\) and the first inner sideband of the \({\hat{d}}_{-}({\hat{d}}_{+})\) mode, as discussed above. In this regard, the \(\left|{{\bf{h}}}_{D}\right|\)-dependent f_Rabi should follow the \(\left|{{\bf{h}}}_{D}\right|\)-dependent \(\Delta {{\omega }}_{{\rm{AT}}}\), which is oscillatory at large \(\left|{{\bf{h}}}_{D}\right|\) as shown in Supplementary Fig. S1a.

Dynamical control of the relative phase between two magnon polariton modes

Experimentally, it has also been shown that driving the transition between two hybridized modes of a two-level photonic system with detuned pulses enables a dynamical control over the relative phase of the two hybrid modes²⁷. By analogy to the protocols described in ref. ²⁷, here we computationally demonstrate the dynamical control of the relative phase between \({\hat{d}}_{+}\) and \({\hat{d}}_{-}\) mode magnon polariton modes (namely, magnonic Ramsey interference). We first excite the system to the \({\hat{d}}_{+}\) mode by injecting a 15-cycle sinusoidal J^c(t) at 7.945 GHz. A π/2 pulse of the Floquet drive field h_D(t) with a frequency of ω_D = δ_D + ∆ω is then applied to create excitations in a superposition of both the \({\hat{d}}_{+}\) and \({\hat{d}}_{-}\) modes, which is illustrated by State ‘I’ on the equator of the Bloch sphere (see Fig. 4a). Here, the duration of the π/2 pulse is 1/4^th of the Rabi-like oscillation period under the field amplitude |h_D | =5000 A/m, i.e.,τ₀ = 1/(4f_Rabi)=2.9 ns. After the first π/2 pulse is turned off (h_D(t) = 0), the excitation as superposition of \(\hat{{d}_{+}}\) and \({\hat{d}}_{-}\) modes would start to precess along the equator of the Bloch sphere (i.e., free evolution) with a precession frequency determined by the detuning amplitude ∆ω. Upon the completion of the free evolution period, a second π/2 pulse is applied to project the excitations (State ‘II’ in the Bloch sphere) into the \({\hat{d}}_{-}\) mode. The amplitude of the \({\hat{d}}_{-}\) mode after the completion of second π/2 pulse is determined by the relative phase between the two modes. The relative phase difference is proportional to both the frequency detuning ∆ω and the duration of free evolution τ. To computationally demonstrate this principle, we record the magnetization amplitude of the \({\hat{d}}_{-}\) mode \(({|\Delta {\bf{m}}|}_{{\hat{d}}_{-}})\) after the completion of the second π/2 pulse, as a function of the free evolution duration τ under ∆ω/2π = 10,15 and 20 MHz. As shown in Fig. 4b, the frequency for the variation of the \({|\Delta {\bf{m}}|}_{{\hat{d}}_{-}}\) with the duration τ is exactly equal to ∆ω. The result matches the analytical expression,\({|{\hat{d}}_{-}\left(t\right)|}^{2}=\frac{{|\alpha|}^{2}}{2}{\sin }^{2}\left(\frac{\Delta\omega}{2}\tau \right)\) (see Supplementary Note 6), where α is the initial mode amplitude in \({\hat{d}}_{+}\).

a (Top) the temporal waveform h_D(t)=|h_D|sin(ω∆t) when 0 ≤ ∆t ≤ τ₀ or τ₀ + τ ≤ ∆t ≤ 2τ + τ₀, and h_D(t)=0 otherwise; (Bottom) Schematic of operation sequences for Ramsey interference on a Bloch sphere. b The magnetization amplitude of the \({\hat{d}}_{-}\) mode of the magnon polariton obtained after the completion of the second π/2 pulse. Each data point in (b) was obtained from an independent simulation, where the free evolution duration τ and detuning amplitude ∆ω are different in each simulation.

Full size image

Coupled mode dynamics under the triple phonon-magnon-photon resonance condition

To simulate the coupled mode dynamics under the triple resonance condition, we design an electro-magno-mechanical system which contains a YIG/SiN bilayer membrane placed in a 3D photon cavity, as shown in Fig. 5a. The cavity hosts a nominal TE₁₀₀ mode photon (i.e., the EM wave has a half-wavelength profile along the x axis while is spatially uniform along y and z), which is a smaller portion of a larger-scale TE₁₀₁ mode cavity. Details of the system design and simulation set-up are shown in the Methods section. The resonant frequency of the cavity photon (~9.1 GHz) is the same as the frequency of the Kittel mode magnon to enable the formation of magnon polaritons. In the present set-up, the frequencies of the two magnon polariton modes are found to be 9.028 GHz \(({\hat{d}}_{-})\) and 9.14 GHz \(({\hat{d}}_{+})\), respectively. Based on the magnetoelastic backaction, the dynamically processing magnetization of the Kittel mode magnon in the YIG layer will generate chiral transverse acoustic (TA) phonons that has a wavevector along the thickness direction (y) of the bilayer, as has been demonstrated experimentally in a similar YIG/GGG bilayer^22,48. To obtain a large profile overlap between phonons and magnons, the layer thicknesses of the YIG (10 nm)/SiN (270 nm) bilayer membrane are designed to host a fundamental (n = 1) TA phonon of 9.13 GHz. Because this frequency is close enough to the frequency of the \({\hat{d}}_{+}\) mode magnon polariton, resonant interaction between the chiral TA phonon and the \({\hat{d}}_{+}\) mode magnon polariton can be enabled. Such triple resonance among the fundamental TA phonon, the Kittel mode magnon, and the k ≈ 0 mode photon is similar to the experiment by Zhang et al. ¹⁵ where the TA phonon of a much higher order (n = 2960) interacts with the Kittel mode magnon, resulting in a smaller mode profile overlap and hence a lower magnon-phonon coupling strength than the present design.

a A cavity electromagnonic system that contains a YIG/SiN bilayer membrane inside a 3D photon cavity. The arrow indicates the direction (mostly along z) and the vector length indicates the magnitude of the local cavity magnetic field H^EM. The schematic on the right illustrates the YIG/SiN bilayer membrane (not to scale), which occupies the space where x∈{14.85 mm, 15.4 mm}, y∈{40 mm, 280 nm}, and z∈{2.5 mm, 3 mm}. The bias magnetic field H^bias is applied along +y. The lower left corner of the cavity is defined as coordinate of the origin, i.e., (x, y, z) = (0,0,0). b Spatial distribution of the \({H}_{z}^{{\rm{EM}}}\) at t = 50 ns in the xz plane of the 3D system at y = 50 nm. c Profile of the strain component ε_xy along the thickness direction (y) of the YIG/SiN bilayer membrane at t = 50 ns. d Evolution of the mechanical displacement u_x and u_z from t = 50-51 ns at 50 nm above the bottom of the YIG/SiN bilayer membrane (this ___location is indicated by the filled circles in c). Under the triple phonon-magnon-photon resonance condition, evolution of (e), the TE₁₀₀ mode cavity photon, represented by the \({H}_{z}^{{\rm{EM}}}\) at the point (x, y, z) = (1.65 mm, 50 nm, 3 mm), (g), the Kittel mode magnon, represented by Δm_x= m_x(t)-m_x(t = 0), and (i), the standing chiral TA phonon mode at the detection point, represented by the ε_xy at the point indicated by the filled circle in (c). t = 0 is the moment the planar current pulse is injected to the cavity. f, h, j Frequency spectra of the cavity photon, the Kittel mode magnon, and the chiral TA phonon, respectively.

Full size image

The vectors in Fig. 5a show the direction and the magnitude of the local H^EM in the cavity. The hybridization of the Kittel mode magnons and cavity photons alters the local EM fields in the vicinity of the YIG resonator, as shown more clearly in Fig. 5b. Figure 5c shows the spatial profile of the fundamental TA phonon across the YIG/SiN bilayer. The right-handed phonon chirality is shown in Fig. 5d. Figure 5e, g, i shows the evolution of the \({H}_{z}^{{\rm{EM}}}\) in the cavity, the \(\Delta {m}_{x}\) in the YIG, and the local ε_xy in the SiN, and their frequency spectra are shown in Fig. 5f, h, j, respectively. As shown, the precessing magnetization of both the \({\hat{d}}_{-}\) and \({\hat{d}}_{+}\) mode magnon polaritons will generate chiral TA phonons of the same frequencies at 9.028 GHz and 9.14 GHz, respectively. The populations of the \({\hat{d}}_{-}\) and \({\hat{d}}_{+}\) modes are similar, as indicated by the similar spectral amplitudes of these two peaks in the frequency spectra of photons and magnons (see Fig. 5f, h). However, the population of the 9.14 GHz phonon, due to its proximity to the intrinsic phonon resonance frequency of 9.13 GHz, is significantly larger than that of the 9.028 GHz phonon (see Fig. 5j). Interestingly, there exists a strong peak at 9.13 GHz in the phonon frequency spectrum, even though the population of the driving \({\hat{d}}_{+}\) mode magnon polariton is low at 9.13 GHz. More interestingly, there exist small peaks at 9.13 GHz in the frequency spectra of both the \({H}_{z}^{{\rm{EM}}}\) and the ∆m_x, indicating an energy backflow from the magnetically excited 9.13 GHz phonon mode to both the magnon and photon systems. This energy backflow, which becomes pronounced in this case mainly because the damping terms for all three systems are turned off (i.e., no energy dissipation), is clear evidence of the triple phonon-magnon-photon resonance. As a control simulation (see Supplementary Fig. S3 in Supplementary Note 7), we found that turning on the elastic damping for the YIG and SiN leads to a significantly lower population for the magnetically existed 9.13 GHz phonon mode, and that there are no additional peaks at 9.13 GHz in the spectra of magnon and photons.

We have developed a 3D dynamical phase-field model that incorporates the coupled dynamics of strain, magnetization, and EM wave in a cavity electromagnonic system, which integrate magnon/phonon resonator(s) in a bulk 3D photon cavity. By solving the coupled equations of motion for these quantities under appropriate magnetic, mechanical, and EM boundary conditions, our computational model allows predicting the spatiotemporal evolution of strain, magnetization, and EM fields under various operating conditions directly from the fundamental material parameters. As examples, time-___domain dynamics of relevant modes in typical coherent gate operations (Rabi oscillation and Ramsey interference) are simulated. The physical validity and high numerical accuracy of the solvers for coupled magnon-photon dynamics in our dynamical phase-field model were demonstrated by understanding the simulation results with analytical fitting and rigorous Hamiltonian-based Floquet theory. We have also applied the dynamical phase-field model to design a cavity electro-magno-mechanical system that enables the triple phonon-magnon-photon resonance, and computationally demonstrate the resonant excitation of a chiral, fundamental (n = 1) TA phonon mode by magnon polaritons under such triple resonance condition.

In combination with the high throughput resulting from the GPU acceleration, the present 3D dynamical phase-field model can be used to guide the experimental design of the microwave photon cavity, magnon and phonon resonator(s) as well as the operating condition for the discovery of new physical phenomena as well as the optimization of key device features such as the coupling strength, mode swapping rate (e.g., Rabi flopping frequency), and cooperativity. The present dynamic phase-field model can also be extended to design and simulate cavity electromagnonic, magnomechanical, and electro-magno-mechanical systems with more complex structures, such as the on-chip systems integrating a coplanar microwave resonator and a magnon resonator^11,12,13,14, by implementing a more detailed treatment of the current dynamics in the normal/superconducting metal components (e.g., see a relevant recent modeling work⁴⁹).

Part 1: Description of the 3D dynamical phase-field model incorporating coupled dynamics of strain, magnetization, and EM waves in multiphase systems

A phase-field model leverages the symmetry-consistent use of continuum physical order parameters and their gradients to describe the total free energy of a spatially inhomogeneous system. The functional derivative of the total free energy \(({F}_{{\rm{tot}}})\) with respect to a specific order parameter yields the thermodynamic driving force that drives the evolution of the order parameter. For example, the effective magnetic field that drives the evolution of M is calculated as \({{\bf{H}}}^{{\rm{eff}}}=-\frac{1}{{\mu }_{0}}\frac{\delta {F}_{{\rm{tot}}}}{\delta {\bf{M}}}\), where μ₀ is vacuum permeability. In a dynamical phase-field model, equations of motion for all key order parameters are typically solved in their exact forms^36,37,50. This is different from conventional phase-field model, where faster-evolving order parameters are often assumed to reach steady or equilibrium state instantaneously⁵¹. A dynamical phase-field model is particularly necessary for hybrid systems featuring bidirectional dynamical energy exchange and conversion between different physical subsystems, as in cavity electromagnonic systems.

As an example, we consider a commonly used system that contains a bulk magnon resonator in a 3D microwave cavity. The evolution of the normalized magnetization m = M/M_s in the magnon resonator, where M_s is the saturation magnetization, is governed by the LLG equation, i.e.,

where \(\gamma\) is the gyromagnetic ratio; \({\alpha }_{{\rm{eff}}}\) is the effective magnetic damping coefficient; H^eff=H^app+H^anis+H^exch+H^mel + H^d + H^EM is the total effective magnetic field, where the externally applied magnetic field H^app(t) includes both the static bias magnetic field H^bias and a dynamic Floquet driving magnetic field h_D(t); H^anis and H^exch are the effective magnetic field resulting from the functional derivatives of the magnetocrystalline anisotropy energy and exchange coupling energy with respect to m, respectively, and their expressions are provided in our previous work³⁶; H^mel is the effective magnetoelastic field, H^d is the demagnetization field, and H^EM is the magnetic field of the cavity electromagnetic wave. H^mel is calculated as \({{\bf{H}}}^{{\rm{mel}}}=-\frac{1}{{\mu }_{0}}\frac{\partial {f}_{{\rm{elas}}}}{\partial {\bf{M}}}\). Here the elastic free energy density \({f}_{{\rm{elas}}}=\frac{1}{2}{c}_{{ijkl}}({\varepsilon }_{{kl}}-{\varepsilon }_{{kl}}^{0})({\varepsilon }_{{ij}}-{\varepsilon }_{{ij}}^{0})\), with i, j = x, y, z. Here the \({c}_{{ijkl}}\) is the elastic stiffness tensor; for magnets of cubic symmetry, the stress-free strain \({\varepsilon }_{{ii}}^{0}=\frac{3}{2}{\lambda }_{100}\left({m}_{i}^{2}-\frac{1}{3}\right)\) and \({\varepsilon }_{{ij}}^{0}=\frac{3}{2}{\lambda }_{111}{m}_{i}{m}_{j}\), where λ₁₀₀ and λ₁₁₁ are magnetostriction coefficients. The local total strain ε_ij can be written as \({\varepsilon }_{{ij}}\left(t\right)={\varepsilon }_{{ij}}^{{\rm{eq}}}+\Delta {\varepsilon }_{{ij}}\left(t\right)\), where the \({\varepsilon }_{{ij}}^{{\rm{eq}}}\) is the total strain at the initial equilibrium state and can be obtained by solving the mechanical equilibrium equation \(\nabla \cdot {\sigma }_{{ij}}^{{\rm{eq}}}=\nabla \cdot [{c}_{{ijkl}}({\varepsilon }_{{ij}}^{{\rm{eq}}}-{\varepsilon }_{{ij}}^{0,{\rm{eq}}})]=0\). For a stress-free, uniformly magnetized magnetic material, \({\varepsilon }_{{ij}}^{{\rm{eq}}}={\varepsilon }_{{ij}}^{0,{\rm{eq}}}\). The dynamical strain \(\Delta {\varepsilon }_{{ij}}\) is calculated via \(\Delta {\varepsilon}_{{ij}}=\frac{1}{2}(\frac{\partial \Delta {u}_{i}}{\partial j}+\frac{\partial \Delta {u}_{j}}{\partial i})\), and the time-varying local mechanical displacement ∆u = u-u^eq is obtained by solving the elastodynamic equation,

where ∆σ = σ-σ^eq is the dynamical stress; ρ is the mass density, and β is the stiffness damping coefficient. Since Eq. (4) is solved for the entire system, the material parameters \({c}_{{ijkl}}\), ρ, and β vary in different phases, where the continuity boundary condition for u and σ⁵² are applied at the interface between the magnon resonator and the microwave cavity. In this regard, by setting the \({c}_{{ijkl}}\) of the microwave cavity to be zero, the stress-free surface of the magnon resonator is automatically considered. The magnon resonator is cube-shaped, which is a computationally more tractable geometry because the entire simulation system is discretized by cube-shaped cells. Based on the coupling to the LLG equation solver, this numerical solver of the elastodynamic equation has previously been applied to simulate coupled magnon-phonon dynamics both the 1D^36,37 and 2D³⁹ magnetic multilayer system. The high numerical accuracy of this elastodynamic solver have been demonstrated through comparison to the analytical solutions in 1D system or the 2D simulations performed via the commercial COMSOL Multiphysics® software. Here, we demonstrate the numerical accuracy of this elastodynamic solver in a 3D elastically inhomogeneous multiphase system through the comparison to the simulation results obtained from the COMSOL Multiphysics® (Supplementary Note 8).

The demagnetization (stray) field H^d can be expressed as \({H}_{i}^{{\rm{d}}}\left(t\right)={H}_{i}^{{\rm{d}},{\rm{eq}}}+\Delta {H}_{i}^{{\rm{d}}}\left(t\right)\). The \({H}_{i}^{{\rm{d}},{\rm{eq}}}\) is produced by the magnetization m⁰ = m(t = 0) at the initial equilibrium state inside the magnon resonator, and can be obtained by solving the continuity equation for magnetic flux \(\nabla \cdot {{\bf{B}}}^{{\rm{eq}}}{\boldsymbol{=}}\nabla \cdot \left[{\mu }_{0}\left({{\bf{H}}}^{{\rm{d}},{\rm{eq}}}\,{\boldsymbol{+}}\,{{\bf{m}}}^{{{0}}}{M}_{\text{s}}\right)\right]=0\), which is part of the Maxwell’s equations. The magnetic boundary condition ∂m/∂n = 0 is applied on the surfaces of the magnon resonator, where n is normal vector to the surface. For a cubic or spheric magnet with spatially uniform magnetization, \({{\bf{H}}}^{{\rm{d}},{\rm{eq}}}=-\frac{1}{3}{M}_{s}({m}_{x}^{0},{m}_{y}^{0},{m}_{z}^{0})\). The dynamically changing \(\Delta {H}_{i}^{{\rm{d}}}\left(t\right)\), which emerges when m starts to evolve, does not need to be calculated separately. Rather, the magnetic-field component of the EM wave H^EM(t), which is obtained by solving the two dynamical equations in the Maxwell’s equations via the finite-difference time-___domain (FDTD) solver on Yee grid (to be discussed below), can automatically satisfy the magnetic flux continuity equation \(\nabla \cdot {\bf{B}}=0\) across the heterointerfaces⁵³. Specifically, the spatiotemporal evolution of the H^EM and the associated electric field component of the EM wave E^EM are simulated by solving the Maxwell’s equations,

Equations (5)–(6) indicate that the H^EM is produced by both the free charge current pulse J^c(t) via electric dipole radiation and the precessing m(t) via magnetic dipole radiation. The perfect electric conductor (PEC) boundary condition is applied on all surfaces of the microwave cavity for reflecting the EM wave without loss. Specifically, E_i = 0 and E_j = 0 on the ij surfaces of the cavity for PEC, with i = x, y, z, and j ≠ i. The elastic stiffness coefficients, the damping coefficient, and the mass density of YIG are listed in Supplementary Note 8. Other material parameters of YIG used in the simulations, including the magnetocrystalline anisotropy and the magnetoelastic coupling coefficients of YIG can be found in ref. ³⁷. Central finite difference is used for calculating spatial derivatives with a midpoint derivative approximation. Conventional Yee grid and the 3D FDTD method⁵³ are used to numerically discretize the EM wave and solve Eqs. (5)–(6). All dynamical equations are solved in a coupled fashion using the classical Runge-Kutta method with a time step ∆t = 5 × 10^–14s. The choice of ∆t is subjected to the Courant condition for numerical convergence in conventional FDTD algorithm, which requires \(\Delta t\le {l}_{0}/(\sqrt{3}v)\)³³, where l₀ is the simulation cell size and v is the EM wave velocity in the medium. Since the use of a larger ε_r leads to \(\sqrt{{\varepsilon }_{{\rm{r}}}}\) times smaller v compared to speed of light in vacuum, a larger ∆t can be used in this work, which significantly reduces the computational time in long-term dynamics simulation.

Part 2: detailed discussion of the size scaling method

Under the same excitation charge current, a larger ε_r would also lead to an E^EM that is \(\sqrt{{\varepsilon }_{{\rm{r}}}}\) times smaller, because the amplitude of E^EM is inversely proportional to the angular wavenumber k of the EM wave, with \(k=\omega \sqrt{{\varepsilon }_{0}{\varepsilon }_{{\rm{r}}}{\mu }_{0}}=\omega /v\). This relationship between E^EM and k can be quantitatively understood by rewriting J^c(t)=∂P/∂t and then analytically solving the wave equation of the E^EM under the plane-wave assumption (see details in⁵⁴).

Despite the smaller E^EM, the amplitude of H^EM would remain unchanged because the ratio of E^EM and the B field is also related by the EM wave velocity v. For example, let us we focus on the B_x \(({H}_{x}^{{\rm{EM}}})\) component, which is the main component that interacts with the magnon mode in the system shown in Fig. 1a. We also consider that only the dominant \({E}_{y}^{{\rm{EM}}}\) component is nonzero, which is consistent with Fig. 1c. In this case. Equation (6) can be rewritten as,

If we further write \({B}_{x}={B}_{x}^{0}{e}^{{\bf{i}}\left(\omega t-{kz}\right)}\) and \({E}_{y}^{{\rm{EM}}}={E}_{y}^{{\rm{EM}},0}{e}^{{\bf{i}}\left(\omega t-{kz}\right)}\) under the plane-wave assumption, Eq. (7) can be further rewritten into,

Therefore, although a larger ε_r reduces the \({E}_{y}^{{\rm{EM}}}\) by \(\sqrt{{\varepsilon }_{{\rm{r}}}}\) times, as shown in Eq. (8), the multiplication by \(\sqrt{{\varepsilon }_{{\rm{r}}}}\) leaves the B_x unchanged. The \({H}_{x}^{{\rm{EM}}}\) also remains unchanged, because (i) \({B}_{x}={\mu }_{0}{H}_{x}^{{\rm{EM}}}\) in the cavity and \({B}_{x}={\mu }_{0}\left(1+{\chi }_{{\rm{m}}}\right){H}_{x}^{{\rm{EM}}}\) in the YIG resonator \(({\bf{M}}={{\boldsymbol{\chi }}}_{{\rm{m}}}{{\bf{H}}}^{{\rm{EM}}})\); and (ii) the magnetic susceptibility tensor \({{\boldsymbol{\chi }}}_{{\rm{m}}}\) is independent of the ε_r.

To demonstrate the applicability of the conclusion above to cases that are more general than plane-wave assumption, we also scale down the size of the cavity electromagnonic system shown in Fig. 1a from the original size of 45 × 9 × 21 mm³ to a few different sizes, in addition to the 180 × 36 × 84 nm³ used in the main paper (for which ε_r was increased from 1 to 6.25 × 10¹⁰). These additional sizes include 225 × 45 × 105 nm³, 300 × 60 × 140 nm³, and 450 × 90 × 210 nm³, where the ε_r was increased from 1 to 4 × 10¹⁰, 2.25 × 10¹⁰, and 1 × 10¹⁰, respectively. We then perform dynamical phase-field simulations to model the excitation of magnon polaritons in these three systems in a similar manner to those in Fig. 1. As expected, the amplitude of \({E}_{y}^{{\rm{EM}}}\) is \(\sqrt{{\varepsilon }_{{\rm{r}}}}\) times smaller while the \({H}_{x}^{{\rm{EM}}}\) remains unchanged with the increasing ε_r, as shown by Supplementary Fig. S5 in Supplementary Note 9. By proportionally scaling down the size of the magnon resonator, the magnon-photon coupling strength would remain unchanged.

The main reason for scaling down a mm-scale system to the nm scale is to ensure that the magnon resonator only accommodates the Kittel mode magnon due to the dominant Heisenberg exchange coupling at the nm-scale. An alternative approach is to use one single cell to represent the magnon resonator (i.e., the macrospin approximation, as in ref. ⁵⁵). We adopt this approach in the design of cavity electro-magno-mechanical system with triple phonon-magnon-photon resonance, as will be discussed below. In this approach, the Heisenberg exchange coupling does not need to be considered, and the spatial discretization of the system is mainly determined by the need to discretize an EM wave. This approach allows for simulating a large-scale photon cavity with a relatively small number of cells and for designing the size and shape of the cavity, but does not permit studying the effect of the size (e.g., as in Fig. 1e) and the shape of the magnon resonator on the coupled magnon-photon dynamics. If using multiple simulation cells yet turning off the exchange coupling, then there would be no force to lock the discrete local magnetization vectors into one giant spin. As time goes, what we have found for the present cavity electromagnonic system (Fig. 1a) is that the numerical error will accumulate (i.e., the values of m_i will become more and more spatially nonuniform) and eventually lead to numerical divergence. Applying a stronger bias magnetic field would only alleviate this numerical issue by delaying, rather than preventing, the numerical divergence, and would also impose a constraint on the frequency range of the magnons that can be excited. This issue of numerical instability is even worse when the size of the magnon resonator is big enough such that the cavity magnetic field inside the resonator becomes highly inhomogeneous.

Part 3: Design and simulation set-up of a 3D cavity electro-magno-mechanical system enabling triple phonon-magnon-photon resonance

To excite a nominal TE₁₀₀ photon mode, we consider a cavity that is 16.5-mm-long (the half-wavelength of the 9.1 GHz EM wave at ε_r = 1) in the x axis and apply periodic boundary conditions to the xy surfaces as well as the PEC boundary condition to all the other surfaces. This set-up is equivalent to placing an array of YIG/SiN bilayer membranes in a large cavity and then probing the phonon-magnon-photon coupling in one of the repeating units (see Supplementary Fig. S6a in Supplementary Note 10 for the full-scale cavity design). To minimize the interaction between the neighboring YIG resonators, we chose a length of 5.5 mm along the z axis, which is sufficiently long to ensure that the EM field remains largely uniform in the xz plane except the region near the YIG, as shown by the spatial distribution of the \({H}_{z}^{{\rm{EM}}}\) in Fig. 5b. The simulated photon profile in Fig. 5a can be considered as a smaller portion of the profile of the TE₁₀₁ mode photon in a larger-scale cavity (see Supplementary Fig. S6b in Supplementary Note 10). Along the y axis of the cavity, which is parallel to the wavevector of the acoustic phonons and the thickness direction of the YIG/SiN bilayer, we consider a length of 400 nm to ensure the formation of the fundamental (n = 1) TA phonon mode at the GHz frequency. The EM field is always uniform along y because the EM wavelength (33 mm) is far larger than the 400 nm. Taken together, the dimensions of the 3D cavity are 16.5 mm (x) × 400 nm (y) × 5.5 mm (z). The size of the YIG/SiN bilayer membrane is 0.55 mm (x) × 280 nm (y) × 0.5 mm (z), where the thicknesses of the YIG and SiN layer are 10 nm and 270 nm, respectively.

Simulating the GHz phonon dynamics along the y axis requires that the cell size is nm-scale along the y axis. According to the Courant condition for the FDTD algorithm mentioned above, a time step ∆t on the order of 10⁻¹⁸s would be required to maintain numerical convergency. It would be computationally prohibitive to simulate the GHz mode dynamics over a time frame over hundreds of ns using such a small ∆t. To address this issue, we increase the ε_r from 1 to 10⁸, which allows for using a 10⁴ \((=\sqrt{{\varepsilon }_{{\rm{r}}}})\) times larger ∆t due to the slower EM wave velocity as discussed earlier. In the meantime, the EM wavelength and hence the size of the photon cavity are reduced by 10⁴ \((=\sqrt{{\varepsilon }_{{\rm{r}}}})\) times. Given that the original cavity is only 400-nm-long along the y axis, this would result in an impractical size of 0.04 nm. Alternatively, we consider (i) a 3D photon cavity of 1650 nm (x) × 400 nm (y) × 550 nm (z) with ε_r = 10⁸, where the dimensions of the cavity along the x and z axis were reduced by 10⁴ times yet the dimension along the y axis is kept the same as the original one; and (ii) a YIG/SiN bilayer of 55 nm (x) × 280 nm (y) × 50 nm (z) in the simulations, where the in-plane dimensions of the bilayer were reduced by 10⁴ times yet the thicknesses of the YIG (10 nm) and SiN (270 nm) layers remain unchanged. This treatment ensures that both the magnon-phonon and the magnon-photon coupling strength remain to be the same as those in the original-sized cavity for two reasons. First, the magnon-phonon coupling strength is determined only by the thicknesses of the YIG and SiN layer (along y) rather than their in-plane dimensions (x and z) because both the magnetization and strain vary only along the y axis. Second, the magnon-photon coupling strength is predominantly determined by the dimension ratio of the YIG-to-photon-cavity along the x axis because the cavity magnetic field H^EM varies largely along the x axis only (see Fig. 5a).

Because our focus is on the Kittel mode magnon and we are not interested in studying the size and shape effect of the YIG in this application example, we use one single cuboid-shaped cell with a size of (∆x, ∆y, ∆z) = (55 nm, 10 nm, 50 nm) to represent the YIG (i.e., the macrospin approximation) and a chain of 27 cells of the same dimension aligning along the y axis to discretize the juxtaposed SiN layer. Cells of the same dimension are also used to discretize the reduced-sized photon cavity, resulting in a total number of cells (N_x, N_y, N_z)=(30, 40, 11). A planar source current \({J}_{y}^{c}\left(t\right)={J}_{y}^{{\rm{c}},0}t{e}^{-{t}^{2}/2{\sigma }_{0}^{2}}\), which is spatially uniform along the yz plane, is injected at x = 275 nm to excite the cavity photon, with \({J}_{y}^{{\rm{c}},0}={10}^{12}{{\rm{A}}/{\rm{m}}}^{2}\) and σ₀ = 7 × 10^–11 s. The bias magnetic field \({H}_{y}^{\text{bias}}\) is applied along the +y direction, causing the magnetization to precess around the y axis, as shown in Fig. 5a. To set up the magnon-photon resonance, we identify the \({H}_{y}^{\text{bias}}\) that makes the FMR frequency ω_m equal the ω_c = 9.1 GHz by numerically simulating the mode splitting spectra of the magnon polaritons as a function of the \({H}_{y}^{\text{bias}}\) in the absence of coupling to phonons. A magnon-photon coupling strength g_cm of 2π×56 MHz is obtained from the mode splitting spectra (see Supplementary Fig. S6c in Supplementary Note 10). We have found that the g_cm remains unchanged when using other ε_r values to reduce the dimensions of the cavity along x and z to other values (while fixing the y-axis cavity dimension to 400 nm) and proportionally reducing the x-axis and z-axis dimension of the YIG/SiN bilayer (while fixing the y-axis dimension of the bilayer to 280 nm). The results are shown in Supplementary Fig. S6d in Supplementary Note 10.

Regarding the set-up of the phonon resonator, since the acoustic wave is spatially uniform in the xz plane, we use one single cell to represent the YIG/SiN bilayer along their x and z axes. Considering that the phonons are excited by a uniform magnetization precessing around the y axis and omitting the elastic damping (β = 0), Eq. (4) can be written as,

where B₁ and B₂ are the magnetoelastic coupling coefficients of the YIG resonator. Equation (8a-c) indicate that the Kittel mode magnons will excite chiral TA phonons of the same frequency which have a wavevector along the y axis (see Fig. 5d). Using procedures similarly to those described in ref. ³⁹, we analytically derive the frequencies of the standing TA phonon frequencies in the YIG/SiN bilayer as a function of the YIG and SiN layer thickness \(({d}^{{\text{YIG}}}{\space}{\rm{and}}\,{d}^{{\text{SiN}}})\), i.e.,

The first nonzero nontrivial solution of Eq. (12) yields the angular frequency of the fundamental (n = 1) acoustic phonon mode (ω_n=1), and so forth for the higher-order modes. Here \({v}^{{\text{YIG}}}=\sqrt{{c}_{44}^{{\text{YIG}}}/{\rho }^{{\text{YIG}}}}\) is the velocity of the TA phonon in YIG and likewise for the \({v}^{{\text{SiN}}}\). When \({d}^{{\text{YIG}}}=10{\space}{\rm{nm}}\) and \({d}^{{\text{SiN}}}=270{\space}{\rm{nm}}\), the resonant acoustic frequency ω_n=1/2π = 9.13 GHz, which is very close to the frequency (9.14 GHz) of the \({\hat{d}}_{+}\) mode magnon polariton. The elastic parameters of the YIG and SiN are provided in Supplementary Note 8.