Extreme Terahertz Magnon Multiplication Induced by Resonant Magnetic Pulse Pairs

Nonlinear interactions of spin-waves and their quanta, magnons, have emerged as prominent candidates for interference-based technology, ranging from quantum transduction to antiferromagnetic spintronics. Yet magnon multiplication in the terahertz (THz) spectral region represents a major challenge. Intense, resonant magnetic fields from THz pulse-pairs with controllable phases and amplitudes enable high order THz magnon multiplication, distinct from non-resonant nonlinearities such as the high harmonic generation by below-band gap electric fields. Here, we demonstrate exceptionally high-order THz nonlinear magnonics. It manifests as 7${}^{\text{th}}$-order spin-wave-mixing and 6${}^{\text{th}}$ harmonic magnon generation in an antiferromagnetic orthoferrite. We use THz multi-dimensional coherent spectroscopy to achieve high-sensitivity detection of nonlinear magnon interactions up to six-magnon quanta in strongly-driven many-magnon correlated states. The high-order magnon multiplication, supported by classical and quantum spin simulations, elucidates the significance of four-fold magnetic anisotropy and Dzyaloshinskii-Moriya symmetry breaking. Moreover, our results shed light on the potential quantum fluctuation properties inherent in nonlinear magnons.

Recently, there is growing evidence that non-thermal control of magnetism by light is feasible [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17]. The THz frequency magnons [18, 19, 20, 21, 22, 23, 24] is central for emergent spin-wave interference-based devices [25] at least 1,000 times faster than the magneto-optical recording and device technologies that use thermal-magnetic excitations. Furthermore, THz control of magnetism is cross-cutting for coherent magnonics [26, 10, 11, 12, 27, 28], quantum magnetism [29, 1, 9], antiferromagnetic (AFM) spintronics [26, 2, 3], and quantum computing [30].

Microwave spin-wave excitation has been shown to induce megahertz frequency magnon multiplication [31]. However, THz frequency magnon multiplication remains an as yet unobserved nonlinear process, in which multiple THz magnons can fuse together to form a high-order magnon coherence. Recent experiments have applied THz frequency, multi-dimensional coherent nonlinear spectroscopy (THz-MDCS) to quantum systems, such as low-dimensional semiconductors [33, 32] and superconductors [35, 34, 36]. THz-MDCS study of magnetic materials have also yielded interesting results such as revealing magnon second harmonic generation [37] , upconversion [38] and mixing [39]. Thus far the results obtained are explained by calculating low-order spin susceptibilities. We place emphasis on the higher-order nonlinearity and collectivity of THz magnons, supported by classical and quantum spin simulations.

To explore magnon multiplication and control, an AFM orthoferrite is driven and probed by THz-MDCS. We use the magnetic fields $\mathbf{B}(t)$ of intense THz pulse-pair to resonantly induce collective spin dynamics and correlation, as illustrated in Fig. 1a. The excitation by such pulse-pair enables us to measure the phase of the AFM magnonic coherent nonlinear emission, in addition to the amplitude as in regular pump-probe (PP) experiments.This approach opens the door to the possibility of driving high-order nonlinear wave-mixing, which is controlled by the interactions of multiple, long-lasting magnon excitations persisting well beyond the duration of the initial laser pulses. We emphasize three additional points here. First, by leveraging both the real-time and relative phase dependencies of coherent emission, we can effectively distinguish, in frequency space after Fourier Transform (FT), discrete spectral peaks resulting from magnon multiplications of varying orders. Second, it is important to distinguish the magnon multiplication from the conventional high harmonic generation of solids [40]. Magnon multiplication utilizes resonant excitations of collective modes leading to high harmonics that emerge after the pulses. In contrast, the conventional electric field-induced high harmonics utilize photoexcitation occurring in the transparency region which exists only during the period of laser excitation. Third, resonant THz photoexcitations have proven to be highly effective in sustaining long-lasting coherence of low-energy collective modes in quantum materials. This is evidenced by recent studies involving topological materials [41, 42], semiconductors [43] and superconductors [36].

In this Letter, we demonstrate extremely high-order nonlinear magnonics by performing THz-MDCS measurements of Sm_0.4Er_0.6FeO₃ using an intense, resonant THz laser pulse-pair. The sensitive detection of coherent spectra enables us to observe up to sixth harmonic generation (6HG) and seventh-order spin wave-mixing (7WM) signals, which significantly exceed what has been previously reported. These distinct magnon multiplication peaks demonstrate a “super” resolution tomography of long-lived interacting magnon quantum coherences and distinguish between different many-body spin correlation functions that describe strongly-driven magnon states. By simulating the coherent dynamics generated by large pulse-pair-driven spin deviations from the equilibrium magnetization orientation, we reproduce the experimental THz-MDCS spectra and identify multi-magnon interaction originating from the dynamical interplay of four-fold magnetic anisotropy and Dzyaloshinskii-Moriya (DM) symmetry breaking.

Figure 1a illustrates our THz-MDCS measurement of quasi-AFM magnons in $a$-cut rare-earth orthoferrite sample Sm_0.4Er_0.6FeO₃ (Methods) at room temperature. This sample exhibits canted AFM order in equilibrium (red arrows), with a small net magnetization pointing along the $c$-axis (blue arrow). The measured linear THz response clearly shows a quasi-AFM magnon mode with energy at $\omega_{\text{AF}}\sim$ 2.5 meV (Fig. 1b). In our THz-MDCS experiment, two collinear THz pulses $E_{1}$ and $E_{2}$ (red lines) with a time delay $\tau$ (red arrow) are used, whose magnetic fields are polarized parallel to the sample $c$-axis. The transmitted nonlinear emission excited from such intense THz pulse-pair at a fixed $\tau$ is detected as the function of the gate time $t$ through electro-optical sampling.

Figure 1c presents representative time-domain traces of the two THz pulses, $E_{1}$ and $E_{2}$ (red and blue solid lines), with electric field strength up to 47.5 MV/m each. We clearly observe the THz light-driven single-order, long-lasting magnon quantum coherence. We expect that magnon-phonon interactions mainly drive the oscillation decay. Fig. 1d shows the coherent emission trace when both pulses are present, $E_{12}$ (light blue line), and the extracted nonlinear signal, $E_{\mathrm{NL}}(t,\tau)=E_{12}(t,\tau)-E_{1}(t,\tau)-E_{2}(t)$ (black line), recorded as the function of gate time $t$ for various inter-pulse delay times $\tau$. The exemplary $E_{\mathrm{NL}}(t,\tau)$ at the fixed $\tau=4.8$ ps in Fig. 1d is generated by the constructive interference of THz-excited magnon with energy $\sim$ 2.5 meV measured in the static THz spectrum (Fig. 1b). This result illustrates the simultaneous amplitude- and phase-resolved detection of the coherent nonlinear response. The subtraction of the single pulse signals $E_{1}$ and $E_{2}$ from the total signal $E_{12}$ allows an extraction of pure coherent nonlinear emissions arising from magnonic nonlinear interactions. The measured long-time dynamics of $E_{\mathrm{NL}}(t,\tau)$ exhibits pronounced, anharmonic oscillations lasting for a long time, as seen within the temporal range marked by a double-arrowed solid line in Fig. 1d. This observation indicates that nonlinear effects of strongly-driven AFM magnons last well beyond the pulse duration. Fig. 1e shows a spectrum of $E_{\mathrm{NL}}(t,\tau)$ obtained over the long oscillation period, which clearly demonstrates the presence of high-order magnonic nonlinearities. This spectrum shows sharp magnon multiplication peaks at frequencies up to six times the magnon energy (Fig. 1b), as marked by the dashed lines. They have been absent in prior THz-MDCS studies of ErFeO₃ [37, 38]. The harmonic peaks are robust from interactions among multiple long-lasting magnon quantum excitations, rather than laser fields. They appear for magnetic fields parallel to the $c$-axis (black line) and are absent for perpendicular to the $c$-axis (blue line).

Figures 2a and 2c show two-dimensional (2D) temporal and spectral profiles of the experimentally measured coherent nonlinear signals $E_{\mathrm{NL}}(t,\tau)$. For comparison, Figs. 2b and 2d present the corresponding results of our numerical simulations, to be discussed below. The oscillations along the horizontal gate time $t$-axis represent the time evolution of magnon interaction for different pulse-pair time delays $\tau$. These gate time $t$-oscillations are pronounced even when two incident pulses do not overlap in time, e.g., for a long inter-pulse delay $\tau=12$ ps. The observation of strong $t$-oscillations at such large $\tau$ indicates that light-induced coherence is stored in the magnon system well after the laser pulse. The oscillatory time-dependent nonlinear signal along the $\tau$-axis for a fixed $t$ displays coherent enhancement due to constructive interference when $\tau$ corresponds to in-phase magnon excitation, while the signal diminishes due to destructive interference when $\tau$ corresponds to out-of-phase magnon excitation. Fig. 2c shows the THz-MDCS spectra obtained by 2D-FT of $E_{\mathrm{NL}}(t,\tau)$ with respect to gate time $t$ (frequency $\omega_{t}$) and inter-pulse delay $\tau$ (frequency $\omega_{\tau}$), where a series of discrete spectral peaks can be observed.

Next, we will first discuss the peaks originating from conventional magnon rectification (MR), pump-probe (PP), and phase-reversal (R) or non-reversal (NR), followed by the newly discovered nonlinear magnonic wave-mixing and harmonic generation processes produced by high-order magnetic interactions as summarized in Table 3 (Methods). Particularly, Fig. 2c clearly shows that strong THz field driving leads to the formation of high-order nonlinear peaks never seen before, including third-harmonic (3HG), fourth-harmonic (4HG), fifth-harmonic (5HG), and sixth-harmonic generation (6HG) as well as five-wave (5WM), six-wave (6WM), and seven-wave mixing (7WM) peaks, as discussed below.

The different THz-MDCS peaks correspond to different nonlinear processes, which can be observed simultaneously and identified according to their locations in $(\omega_{t},\omega_{\tau})$ space. To analyze the measured 2D spectrum, we introduce frequency vectors for incident pulse 1 and 2 excitations, $\omega_{1}=(\omega_{\mathrm{AF}},\omega_{\mathrm{AF}})$ and $\omega_{2}=(\omega_{\mathrm{AF}},0)$ (red and black arrows in Fig. 2d), where $\omega_{\mathrm{AF}}\sim$ 2.5 meV denotes to the AFM magnon mode energy. The several discrete peaks separated along the vertical axis (relative phase), as shown in Fig. 2c, demonstrate that the pulse-pair time delay $\tau$ controls phase coherence between magnons excited by two incident pulses. Figs. 2e and 2f further show a zoom-in view of some high harmonic and wave mixing peaks observed in the experiment, while Figs. 2g and 2h show the corresponding peaks reproduced by the theory.

In order to compare the measured and simulated nonlinear processes, we first identify the physical origins of the series of discrete peaks seen in Fig. 2c along the vertical direction $\omega_{t}=\omega_{\mathrm{AF}}$, located at $\omega_{\tau}=-\omega_{\mathrm{AF}}$, 0, $\omega_{\mathrm{AF}}$, and $2\omega_{\mathrm{AF}}$. The $(\omega_{\mathrm{AF}},0)$ peak arises from the PP process (Table 3) that reflects the dynamics of magnon coherent populations. In particular, two field interactions during pulse 1 create a magnon population via the difference frequency process $\omega_{1}-\omega_{1}$. This population evolves during time $\tau$ and then interacts with a magnon excitation during pulse 2, resulting in a single magnon 1-quantum coherence (1QC). The peak at $(\omega_{\mathrm{AF}},2\omega_{\mathrm{AF}})$ shows the dynamics of a magnon 2-quantum coherence (2QC), which is generated by two magnon excitations during pulse 1 via the sum–frequency process $\omega_{1}+\omega_{1}$. This 2-magnon coherence is probed through its interaction with a magnon excitation during pulse 2.

The two peaks located at $(\omega_{\mathrm{AF}},\omega_{\mathrm{AF}})$ and $(\omega_{\mathrm{AF}},-\omega_{\mathrm{AF}})$ arise from the interaction of a 1QC generated by pulse 1 and a 1QC generated by pulse 2. Their interaction creates a magnon population via the difference frequency processes $\pm(\omega_{1}-\omega_{2})$, whose interaction with a magnon excitation during pulse 2 leads to a 1QC which is either phase reversed (R) or not reversed (NR) with respect to the 1QC of pulse 1. We also observe peaks along $\omega_{t}=3\omega_{\mathrm{AF}}$ that correspond to third harmonic generation (THG) signals. Here, the interaction of a 2QC created by pulse 1(2) with a 1QC generated by pulse 2(1) produces a THG signal at $(3\omega_{\mathrm{AF}},2\omega_{\mathrm{AF}})$ and $(3\omega_{\mathrm{AF}},\omega_{\mathrm{AF}})$.

The measured THz-MDCS 2D spectra allow us to identify the role of the antisymmetric exchange DM interaction in the studied AFM system. This is possible through peaks generated by second-order magnon nonlinear processes facilitated by the symmetry-breaking DM interaction. The latter interaction is already known to induce a canted AFM order in the ground state. Such spin canting with respect to the AFM orientation of a spin-up and a spin-down sub-lattice results in a net magnetization along the $c$-axis (blue arrow, Fig. 1a). The coherent spin dynamics along this axis then becomes anharmonic, which results in the second-order nonlinear peaks observed in the 2D spectra at $\omega_{t}=2\omega_{\mathrm{AF}}$ and at $\omega_{t}=0$ (Fig. 2c) that involve two magnon excitations. The difference frequency processes $\omega_{1}-\omega_{2}$ involving two magnon 1QC generate magnon rectification (MR) signals at $(0,\pm\omega_{\mathrm{AF}})$ via the DM interaction. The latter antisymmetric exchange also generates second harmonic generation (SHG) peaks at $(2\omega_{\mathrm{AF}},\omega_{\mathrm{AF}})$ and $(2\omega_{\mathrm{AF}},2\omega_{\mathrm{AF}})$. These SHG peaks arise from the sum-frequency generation processes $\omega_{1}+\omega_{2}$ and $\omega_{1}+\omega_{1}$, respectively.

Finally, the THz-MDCS spectra at Figs. 2c and 2e reveal high-order nonlinear peaks involving up to six magnon quanta. In particular, we observe strong 4HG, 5HG, and 6HG peaks, and also strong 5WM, 6WM, and 7WM nonlinear signals, shown in Fig 2f. It is important to note that the 2D frequency resolution enabled by the THz-MDCS allows us to not only excite coherent magnons but also achieve a “super” resolution tomography different from the conventional high harmonic generation spectroscopies that only detect one dimensional spectra. In the 2D frequency plane, the same, high-order magnon harmonics will appear as different correlated THz peaks along multiple axes, e.g., 6HG and 5WM peaks separated in 2D space. Such strong correlation enables the deterministic assignment of higher–order nonlinear peaks, even in the presence of measurement noise or FT artifacts, which were very difficult to discern before. As discussed below, our numerical simulations shown in Figs. 2d, 2g and 2h further corroborate our discovery of the higher-order magnon multiplication peaks. We attribute these high-order nonlinear signals to the existence of both four-fold magnetic anisotropy and DM interaction in our AFM system. The origin of all observed THz-MDCS peaks is summarized in Table 3 (Methods). Note that additional experimental signals at frequencies less than $\omega_{\mathrm{AF}}$ are not captured by our mean field model. While we do not discuss such signals here, they may arise from parametric processes involving two magnons with finite momenta $\mathbf{q}$ and $-\mathbf{q}$ such that $\omega_{\mathrm{AF}}=\omega_{\mathbf{q}}+\omega_{-\mathbf{q}}$[44].

To corroborate the above interpretation of experimental results, we use the following two sub-lattice Hamiltonian [45, 46]:

$\displaystyle H=J\,\mathbf{S}_{1}\cdot\mathbf{S}_{2}-\mathbf{D}\cdot\left(\mathbf{S}_{1}\times\mathbf{S}_{2}\right)-\sum_{i=1,2}\left(K_{a}\,S_{i,a}^{2}+K_{c}\,S_{i,c}^{2}\right)-K_{4}\sum_{\begin{subarray}{c}i=1,2\\ j=a,b,c\end{subarray}}S^{4}_{i,j}-\gamma\,\mathbf{B}(t)\cdot\sum_{i=1,2}\mathbf{S}_{i}\,.$

(1)

This classical spin model effectively describes the AFM resonance mode and THz magnon multiplication behavior in orthoferrite systems where the DM interaction is the predominant canting mechanism [47]. The first term of Eq. (1) accounts for the AFM coupling between the neighboring spins $\mathbf{S}_{1}$ and $\mathbf{S}_{2}$ with exchange constant $J>0$. The second term describes the Dzyaloshinskii-Moriya (DM) interaction, with anti-symmetric exchange parameter $\mathbf{D}$ aligned along the $b$-axis. The orthorhombic crystalline anisotropy is described by the third and fourth terms of the Hamiltonian, with anisotropy constants $K_{a}$ and $K_{c}$ along the $a$ and $c$ axes, respectively. To explain the observed high-order nonlinearities, we include the four-fold anisotropy $K_{4}$. Finally, the driving of the spin system by the laser magnetic field $\mathbf{B}(t)$ is described by the Zeeman interaction (the last term of the above Hamiltonian), where $\gamma=g\mu_{\mathbf{B}}/\hbar$ is the gyromagnetic ratio, $g=2$ is the Landé $g$-factor, and $\mu_{\mathrm{B}}$ denotes the Bohr magneton.

We have simulated the nonlinear spin dynamics by solving the full Landau-Lifshitz-Gilbert equations [2, 37] (LLG, Eq. (2) in Methods) derived from the Hamiltonian Eq. (1). Non-interacting magnon excitations can be described by linearizing this LLG equation in the case of small amplitude oscillations around the equilibrium magnetization. This approximation is valid for weak $\mathbf{B}(t)$. However, the full nonlinearities described by the LLG equation, central in many spin-wave device concepts, become critical for strong fields, which drive large amplitude of magnon oscillations. Here we characterize such magnon nonlinearities by solving the full nonlinear LLG equation (Eq. (2), Methods) for two magnetic fields polarized along sample $c$-axis. We obtain THz-MDCS spectra by doing 2D-FT of the nonlinear magnetization $\mathbf{M}_{\mathrm{NL}}$, which is obtained by subtracting the magnetizations induced by two the individual pulses from the nonlinear response of the net magnetization $|\mathbf{M}|$ under intense pulse-pair excitation (Methods). Figs. 2b and 2d present the 2D temporal profile $\mathbf{M}_{\mathrm{NL}}(t,\tau)$ and spectrum $\mathbf{M}_{\mathrm{NL}}(\omega_{t},\omega_{\tau})$, respectively. The simulated $\mathbf{M}_{\mathrm{NL}}(\omega_{t},\omega_{\tau})$ is fully consistent with the measured spectrum in Fig. 2c. This overall agreement of the two-time/two-frequency dependencies between theory and experiment allows us to assign the different discrete peaks observed in the experiment to specific nonlinear processes and magnetic interactions. The agreement also corroborates the magnonic origin of the observed high-order nonlinearities and symmetry-breaking peaks. Note that the remaining discrepancy with experimental results underscores the significance of considering the quantum properties of the spin in describing THz-MDCS experiments discussed in more detail below.

Figures 3a - 3f show the field dependence of the nonlinear spectral signals $E_{\mathrm{NL}}(\omega_{t},\tau)$ at magnonic multiplication frequencies $\omega_{t}$= 0, 1, 2, 3, 4, 5$\omega_{\mathrm{AF}}$. The $E_{\mathrm{NL}}(\omega_{t},\tau)$ traces are taken at $\tau$ = 4.8 ps to maximize the sensitivity similar to Fig. 1e. These high-order magnon multiplication peaks clearly grow nonlinearly with increasing field strength and follow power-law dependencies. By fitting (solid lines) the 5$\omega_{\mathrm{AF}}$ peak amplitude in Fig. 3a, we see that is clearly scales as $E^{5}$. Note that 0$\omega_{\mathrm{AF}}$ and 1$\omega_{\mathrm{AF}}$ peaks scale as $E^{2}$ and $E^{3}$, respectively, confirming that they originate from second-order and third-order nonlinear processes as discussed above.

Figures 3g and 3h compare a pulse-pair and single pulse pumping for generating high-order magnon multiplication peaks at 2, 3, 4$\omega_{\mathrm{AF}}$ (dashed vertical lines). The extracted $E_{\mathrm{NL}}(\omega_{t},\tau)$ spectra at fixed $\tau=$ 4.8 ps in Fig. 3g is presented for 100% (purple line), 70% (green line), and 50% (orange line) of the maximum THz field strength and grow nonlinearly with increasing field. Compared to the $E_{\mathrm{NL}}(\omega_{t},\tau)$ spectra, the single–pulse signal, $E_{2}(\omega_{t})$ (Fig. 3h), does not show any significant high-order magnon multiplication peaks. Magnonic nonlinearity and interaction effects are hidden in the response to a single strong pulse, but become prominent in the THz-MDCS. In the latter case, the single-pulse responses are subtracted out, so the pure nonlinear signals are determined by correlations between two pulse excitations, rather than by single pulse nonlinearities. As a result, THz-MDCS allows us to resolve the nonlinear interactions between coherent magnons excited by different pulses and associate them with different peaks in 2D frequency space.

Figure 4 presents our numerical calculations in more detail. Particularly, it provides direct insight into the physical origin of the different THz-MDCS spectral peaks by comparing the contributions from different magnetic interactions and anisotropies. We compare the spectra $\mathbf{M}_{\mathrm{NL}}(\omega_{t},\omega_{\tau})$ calculated in an electric field strength of $\sim 45$ MV/m between (i) the full-term calculation (Fig. 4a), (ii) a calculation without DM interaction, $D=0$ , where the ground state has no spin canting (Fig. 4b), and (iii) a calculation without four-fold magnetic anisotropy, i. e., by setting $K_{4}=0$ in the Hamiltonian (1) (Fig. 4c). All even-order wave-mixing signals vanish when the DM interaction is switched off, consistent with earlier works [37]. In addition, the third and higher odd-order nonlinear magnon signals are strongly suppressed (Fig. 4b). The third- and higher-order nonlinear magnon signals are also strongly suppressed if $K_{4}=0$ (Fig. 4c). The suppression of high order signals in the simulations (Figs. 4b and 4c) indicates that the high-order nonlinear peaks observed in THz-MDCS experiment are generated by the four-fold magnetic anisotropy as well as the DM interaction. This interpretation can be seen more clearly by plotting a spectral slice along $\omega_{\tau}=\omega_{\mathrm{AF}}$ in Fig. 4d, where the full-term calculation (black line) is compared with simulations where (i) the DM interaction is switched off (blue line), or (ii) the four-fold magnetic anisotropy is switched off (red line). This observation reveals the existence of a four-fold magnetic anisotropy and also demonstrates that the high order magnonic nonlinearities arise from the dynamical interplay between DM interaction and four-fold magnetic anisotropy.

Figures 4e and 4f underscore the significance of considering the quantum fluctuation properties of spins in describing the THz-MDSC spectra. While the simulations based on the classical spin LLG model in Fig. 2d can replicate the positions of the peaks in the THz-MDCS spectra, the strengths of the high-order nonlinear signals using classical spins are much weaker compared to the experimental data in Fig. 1e and Fig. 2c. For instance, the experimental data (gray squares) at $\omega_{\tau}=2\omega_{\mathrm{AF}}$ and $\omega_{\tau}=3\omega_{\mathrm{AF}}$ exhibits substantially larger high-harmonic generation signals in Fig. 4f compared to the result of the classical simulations (blue diamonds). To understand this we have also calculated $\mathbf{M}_{\mathrm{NL}}(\omega_{t},\omega_{\tau})$ for the quantum spin version of Hamiltonian (1). This is done by replacing the classical spins $\mathbf{S}_{i}$ in Eq. (1) by quantum spin operators $\hat{\mathbf{S}}_{i}$ at site $i$ of spin magnitude $S=5/2$. The resulting two-site quantum spin Hamiltonian has $(2S+1)^{2}=36$ eigenstates $|\Psi_{n}\rangle$ and eigenenergies $E_{n}$. The inset of Fig. 4e shows the eigenenergy spectrum of the quantum spin model. We have used the same parameters for $K_{a}$, $K_{c}$, and $K_{4}$ as in the classical simulations and adjusted $J$ and $D$ such that the magnon frequency, defined by $\omega_{\mathrm{AF}}=E_{1}-E_{0}$ in the quantum spin model, matches the one observed in the experiment (Table 2, Methods). The energy spectrum is anharmonic, however, the difference between the eigenenergies is given by multiples of the magnon energys such that one expects a harmonic energy spectrum.

To model the transmitted magnetic field, we calculate the dynamics of the magnetization along $c$-direction, $M^{c}(t)\equiv\langle\Psi(t)|(\hat{S}_{1,c}+\hat{S}_{2,c})|\Psi(t)\rangle$. Here $|\Psi(t)\rangle$ is the quantum state at time $t$ which is obtained via exact diagonalization. The calculated nonlinear magnetization $M_{\mathrm{NL}}^{c}(t,\tau)=M_{12}^{c}(t,\tau)-M_{1}^{c}(t,\tau)-M_{2}^{c}(t)$ is presented in Fig. 4e (red line) at a fixed inter-pulse delay of $\tau=$ 4.8 ps as in Fig. 1e. We observe high-harmonic generation peaks (vertical dashed lines) up to seventh-order. Compared to that, the result of the corresponding classical simulation in Fig. 4e (blue line) exhibits about one-order of magnitude smaller second to fifth-order harmonic generation signals. To identify the origin of the strong magnonic high-harmonic generation signals of the quantum spin simulation, we express $|\Psi(t)\rangle$ after the magnetic field pulse-pair excitation in terms of the eigenstates $|\Psi_{n}\rangle$, yielding $|\Psi(t)\rangle=\sum_{n}a_{n}\,e^{-iE_{n}t}|\Psi_{n}\rangle$ with $a_{n}\equiv\langle\Psi_{n}|\Psi(t=0)\rangle$. The dynamics of magnetization can then be written as $M^{c}(t)=\sum_{j,k}M^{c}_{j,k}\,a_{j}^{*}a_{k}e^{-i(E_{k}-E_{j})t}$ with magnetic dipole matrix elements $M^{c}_{j,k}\equiv\langle\Psi_{j}|(\hat{S}_{1,c}+\hat{S}_{2,c})|\Psi_{k}\rangle$. As a result, the peaks in the spectrum arise from transitions between different eigenstates. The contribution of the transitions between the different eigentstates to $M_{c}(t)$ is determined by strength of the magnetic dipole matrix elements $M^{c}_{j,k}$. Since the difference between eigenenergies (inset of Fig. 4e) corresponds to multiples of the magnon frequency, we obtain a harmonic $M^{c}$ spectrum as observed in Fig. 4e.

Figure 4f directly compares the magnonic high harmonic generation observed in the simulations and experiment by plotting the ratio between the peak strengths of the $n$th harmonic and the fundamental harmonic. The result of the experiment (gray squares) is compared with the ratios deduced from the quantum-spin simulation (red circles) and classical simulation (blue diamonds). The quantum spin simulation agrees well with the ratios from experiment while the classical simulation produces substantially smaller high-harmonic generation peaks. Our results from the quantum spin calculation depicted in Fig. 4f emphasize the importance of further investigation into the pronounced higher-order magnon nonlinearities observed in our THz-MDCS experiment, which notably exceed the predictions of the classical spin model. Additional factors also need to be considered in the future quantitative analysis of high-order magnon nonlinear peaks, including the THz electro-optic sampling response functions [51] and THz propagation effects [52, 53]. But we expect neither of these effects to substantially modify the peak strengths of the 2$\omega_{\mathrm{AF}}$ and 3$\omega_{\mathrm{AF}}$ nonlinear magnon peaks, which constitute the focus of the quantum spin calculations.

At last, we elaborate three key distinctions of the THz magnon multiplication observed in this work. First, magnon multiplication are generated by utilizing resonant excitation of collective spin wave modes with long-lasting quantum coherence that persists well after the laser pulses. High harmonic generation (HHG) in solids often occurs during the period of laser excitation. Second, THz-MDCS achieves super-resolution tomography of magnon interaction and nonlinearity. Our results reveal four-fold magnetic anisotropy and the Dzyaloshinskii-Moriya interaction in creating interactions among multiple long-lasting magnon quantum excitations. Such interaction leads to the high-order magnon multiplication. Third, the magnon multiplication peaks in the 2D spectrum can arise from resonant transitions between different quantum spin eigenstates [54]. Consequently, magnon multiplication nonlinearity can potentially provide direct information about quantum spin systems. Our results provide compelling implications of considering the quantum fluctuation properties of spins.

In summary, we demonstrate the extremely high-order magnonic multiplication and wave-mixing at THz frequency by driving a canted antiferromagnetic state. Our experimental results are reproduced by both classical and quantum spin simulations, signifying the importance of four-fold magnetic anisotropy and DM symmetry breaking, also the fluctuation properties of quantum spins in the higher order magnon multiplication. Our demonstration of magnonic nonlinearity and collectivity through THz coherent spectroscopy can be further extended to pursue nonlinear quantum magnonics in parallel with nonlinear quantum optics.

All the data generated in this study are provided in the Source Data file.

This work was supported by the Ames National Laboratory, the US Department of Energy, Office of Science, Basic Energy Sciences, Materials Science and Engineering Division under contract No. DEAC0207CH11358 (project design and THz-MDCS spectroscopy). Synthesis and basic characterizations of crystals was supported by National Natural Science Foundation of China (NSFC, No.52272014), Science and Technology Commission of Shanghai Municipality (No. 21520711300) and CAS Project for Young Scientists in Basic Research (YSBR-024). The quantum spin model simulations were supported by the U.S. Department of Energy, Office of Science, National Quantum Information Science Research Centers, Co-design Center for Quantum Advantage (C2QA) under contract number DE-SC0012704.

C.H. and L.L. performed the THz spectroscopy measurements with J.W.’s supervision. M.M., and Y.Y simulated the THz magnon nonlinearity using the quantum spin model. M.M., and I.E.P analyzed the THz spectra using the classical spin model. J.S., P.M., L.S. and A.W. grew the samples and performed crystalline quality and other basic characterizations. The paper is written by J.W., C.H., and M.M. with discussions from all authors. J.W. conceived and coordinated the project.

The authors declare no competing interests.