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 7th-order spin-wave-mixing and 6th harmonic magnon generation in an antiferromagnetic orthoferrite. We use THz two-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.

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 highorder magnon coherence.Recent experiments have applied THz frequency, multi-dimensional coherent nonlinear spectroscopy (THz-MDCS) to quantum systems, such as low-dimensional semiconductors [32,33] and superconductors [34][35][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 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.4 Er 0.6 FeO 3 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.

Results
Figure 1a illustrates our THz-MDCS measurement of quasi-AFM magnons in a-cut rare-earth orthoferrite sample Sm 0.4 Er 0.6 FeO 3 (Methods) at room temperature.This sample exhibits canted AFM order in equilibrium (red arrows), with a small net magnetization pointing along the caxis (blue arrow).The measured linear THz response clearly shows a quasi-AFM magnon mode with energy at ω AF ∼ 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 τ (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 τ 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 line), recorded as the function of gate time t for various inter-pulse delay times τ .The exemplary E NL (t, τ ) at the fixed τ = 4.8 ps in Fig. 1d is generated by the constructive interference of THz-excited magnon with energy ∼ 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 NL (t, τ ) 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 NL (t, τ ) 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 3 [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 NL (t, τ ).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 τ .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 τ = 12 ps.
The observation of strong t-oscillations at such large τ indicates that light-induced coherence is stored in the magnon system well after the laser pulse.The oscillatory time-dependent nonlinear signal along the τ -axis for a fixed t displays coherent enhancement due to constructive interference when τ corresponds to in-phase magnon excitation, while the signal diminishes due to destructive interference when τ corresponds to out-of-phase magnon excitation.Fig. 2c shows the THz-MDCS spectra obtained by 2D-FT of E NL (t, τ ) with respect to gate time t (frequency ω t ) and inter-pulse delay τ (frequency ω τ ), 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 (ω t , ω τ ) space.To analyze the measured 2D spectrum, we introduce frequency vectors for incident pulse 1 and 2 excitations, ω 1 = (ω AF , ω AF ) and ω 2 = (ω AF , 0) (red and black arrows in Fig. 2d), where ω AF ∼ 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 τ 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 ω t = ω AF , located at ω τ = −ω AF , 0, ω AF , and 2ω AF .The (ω 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 ω 1 − ω 1 .This population evolves during time τ and then interacts with a magnon excitation during pulse 2, resulting in a single magnon 1-quantum coherence (1QC).The peak at (ω AF , 2ω 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 ω 1 + ω 1 .This 2-magnon coherence is probed through its interaction with a magnon excitation during pulse 2.
The two peaks located at (ω AF , ω AF ) and (ω AF , −ω 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 ±(ω 1 −ω 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 ω t = 3ω 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ω AF , 2ω AF ) and (3ω AF , ω 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 ω t = 2ω AF and at ω t = 0 (Fig. 2c) that involve two magnon excitations.The difference frequency processes ω 1 − ω 2 involving two magnon 1QC generate magnon rectification (MR) signals at (0, ±ω AF ) via the DM interaction.The latter antisymmetric exchange also generates second harmonic generation (SHG) peaks at (2ω AF , ω AF ) and (2ω AF , 2ω AF ).These SHG peaks arise from the sum-frequency generation processes ω 1 + ω 2 and ω 1 + ω 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 ω 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 q and −q such that ω AF = ω q + ω −q [44].
To corroborate the above interpretation of experimental results, we use the following two sub-lattice Hamiltonian [45,46]: 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 S 1 and S 2 with exchange constant J > 0. The second term describes the Dzyaloshinskii-Moriya (DM) interaction, with anti-symmetric exchange parameter 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 Finally, the driving of the spin system by the laser magnetic field B(t) is described by the Zeeman interaction (the last term of the above Hamiltonian), where γ = gµ B /h is the gyromagnetic ratio, g = 2 is the Landé g-factor, and µ 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).Noninteracting magnon excitations can be described by linearizing this LLG equation in the case of small amplitude oscillations around the equilibrium magnetization.This approximation is dependencies.By fitting (solid lines) the 5ω AF peak amplitude in Fig. 3a, we see that is clearly scales as E 5 .Note that 0ω AF and 1ω 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ω AF (dashed vertical lines).The extracted E NL (ω t , τ ) spectra at fixed τ = 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 NL (ω t , τ ) spectra, the single-pulse signal, E 2 (ω 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.

Discussion
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 M NL (ω t , ω τ ) calculated in an electric field strength of ∼ 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 ω τ = ω 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.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 Here |Ψ(t)⟩ is the quantum state at time t which is obtained via exact diagonalization.The calculated nonlinear magnetization 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 nth 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ω AF and 3ω 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 wavemixing 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.(99.99%) as the starting materials.The mixture was pressed into pellets and sintered in air at 1400 K for 20 h.After confirmation by X-ray diffraction that the products had converted to the orthoferrite structure, the powder was then compacted into the rods.The rods were sintered in a vertical tube furnace at 1500 K for 15 h in air.Single crystal growth was performed in the optical floating-zone furnace (FZ-T-10000-H-VI-P-SH, Crystal System Corp).

Methods Sample
The sample measured in the experiment is a-cut Sm 0.4 Er 0.6 FeO 3 single crystal with size around 10 mm×5 mm×1 mm.At room temperature, the canted spin antiferromagnetic order leads to a small net magnetization along c-axis (Fig. 1a).The magnetic resonance mode (quasi-AFM mode) can be excited if B THz is parallel to this net magnetization.The in-plane axis is determined by rotating the sample to achieve the strongest magnon oscillation signal.In the experiment, we fixed sample c-axis perpendicular to THz electric field.

THz multi-dimensional coherent nonlinear spectroscopy
To generate and characterize THz-MDCS spectra, the Sm 0.4 Er 0.6 FeO 3 single crystal sample was excited by two broadband THz laser pulses.The measured coherent differential transmission To acquire the THz nonlinear signal, a double modulation scheme has been used.Two optical choppers modulate two intense THz beams at a frequency of 500 Hz (pulse 1) and 250 Hz (pulse 2), respectively.Such a scheme allows to isolate nonlinear THz signals.The FT is performed at 7 -28 ps range of the long lasting time-domain oscillatory traces as shown in Figs.1c and 1d.The range is chosen to ensure avoidance of the stimulating THz pulses and echo pulses.The high order THz magnon multiplication obtained were possible by (1) choosing the Sm-doped magnetic material with the bigger four-fold magnetic anisotropy compared to un-doped ones, (2) applying an experimental scheme of 2 pulses vs 1 pulse driving and (3) optimizing the detection sensitivity of 2D coherent nonlinear signals.

THz-MDCS signal simulations
To simulate the THz-MDCS signals, we solve Landau-Lifshitz-Gilbert (LLG) equations derived from the Hamiltonian (1) [37]: which describe the nonlinear coherent spin dynamics induced by the magnetic field B(t).Here, ∂S i corresponds to the time-dependent effective local magnetic field at site i and α is the Gilbert damping coefficient.We numerically solve the full nonlinear LLG equations (2) with input THz magnetic field pulses 1 and 2 polarized along c-axis for realistic orthoferrite material parameters summarized in Table 1.The parameters used match the magnon mode frequency, free-induction decay, and nonlinear magnon signals observed in the experiment.We then calculate the net magnetization M ∝ S 1 +S 2 of the two AFM sub-lattices.To directly simulate our THz-MDCS experiment, we calculate the nonlinear differential magnetization M NL , which corresponds to the measured E NL .M NL is obtained by computing the net magnetization induced by both pulses, M 12 (t, τ ), as a function of the gate time t and the inter-pulse time τ , as well as the net magnetizations resulting from pulse 1, M 1 (t, τ ), and pulse 2, M 2 (t).The nonlinear differential magnetization is then given by M NL (t, τ ) = M 12 (t, τ ) − M 1 (t, τ ) − M 2 (t) for the collinear geometry used in the experiment.The THz-MDCS spectra are obtained by Fourier transform of M NL (t, τ ) with respect to both t (frequency ω t ) and τ (frequency ω τ ).To analyze these 2D spectra, we introduce the frequency vector (ω t , ω τ ).Since both pulses resonantly drive magnon 1-quantum coherences (1QC) characterized by the magnon frequency ω AF ∼ 2.5 meV (Fig. 1b), the frequency vectors of the two pulses 1 and 2 can be written as ω 1 = (ω AF , ω AF ) and ω 2 = (ω AF , 0).The observed spectral peaks and corresponding nonlinear processes are summarized in Table 3 Table 3. Nonlinear magnon processes contributing to the THz-MDCS spectra 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 con-

Figures 4e and 4f
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 ω τ = 2ω AF and ω τ = 3ω 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 M NL (ω t , ω τ ) for the quantum spin version of Hamilto- is presented in Fig. 4e (red line) at a fixed interpulse delay of τ = 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 fifthorder harmonic generation signals.To identify the origin of the strong magnonic high-harmonic generation signals of the quantum spin simulation, we express |Ψ(t)⟩ after the magnetic field pulse-pair excitation in terms of the eigenstates |Ψ n ⟩, yielding |Ψ(t)⟩ = n a n e −iEnt |Ψ n ⟩ with a n ≡ ⟨Ψ n |Ψ(t = 0)⟩.The dynamics of magnetization can then be written as M c is plotted as a function of the gate (real) time t and the delay time between pulses 1 and 2, τ , which controls the relative phase accumulation between the two phase-locked THz fields.The time-resolved coherent nonlinear dynamics is then explored by varying the inter-pulse delay τ .By measuring the electric fields in the time-domain through electro-optic sampling (EOS) by a third pulse, we achieve phase-resolved detection of the sample response as a function of gate time t for each time delay τ .The nonlinear response signals, separated from the linear response, leads to enhanced resolution that enables us to observe high order magnonic nonlinearity.Experimental detailsOur THz-MDCS setup is driven by a 800 nm Ti: Sapphire amplifier with 4 mJ pulse energy, 40 fs pulse duration, and 1 kHz repetition rate.The majority part of laser output (2.7 W) is split into two beams with similar power (1.35 W each), and the inter-pulse delay between them is controlled by a mechanical delay stage.Two beams are recombined and focused into a MgO doped LiNbO 3 crystal to generate intense THz pulses through the tilted pulse front scheme.The intense THz pulse has a peak field strength ∼ 47.5 MV/m, central frequency ∼ 1 THz (4.1 meV), and bandwidth ∼ 1.5 THz [48-50].The focused THz spot size at the sample position is ∼ 1.3 mm.The rest part of the laser output (∼ 1 mW) is used to detect THz fields in time-domain by EOS using a 2 mm thick (110) oriented ZnTe crystal.

Fig. 2 .Fig. 3 .
Fig.2.Super-resolution tomography of magnonic multiplication and correlation.a, Twodimensional false-color plot of the measured coherent nonlinear transmission signals E NL (t, τ ) as a function of the gate time t and the delay τ induced with the electric field strength 47.5 MV/m.b, The simulated E NL (t, τ ).c and d, The E NL (ω t , ω τ ) spectra from 2D-FT of time domain images a and b.The THz-MDCS spectra show PP, R, NR, and 2QC peaks, which are generated by third-order nonlinear processes.Second-order nonlinear processes lead to MR and SHG peaks.Importantly, strong THz driving also leads here to formation of high-order nonlinear signals including THG, 4HG, 5HG, 6HG as well as 5WM, 6WM and 7WM peaks.The frequency vectors for two incident pulses 1 and 2, ω 1 = (ω AF , ω AF ) and ω 2 = (ω AF , 0), are marked as red and black arrows.e and f, Zoom-in view of new high harmonic and wave mixing peaks observed in our experiments which involve up to six magnon quanta (main text).g and h, Corresponding zoom-in peaks reproduced by theory.

Fig. 4 .
Fig.4.Simulations of high-order nonlinear magnonic correlations.ac, M NL (ω t , ω τ ) for (a) the full-term calculation, (b) a calculation without DM interaction, D = 0, where the ground state has no spin canting, i. e., M = 0 in the ground state, and (c) a calculation without fourfold anisotropy in the Hamiltonian, K 4 = 0. d, M NL (ω t , ω τ ) at fixed ω τ = ω AF .The full-term calculation (black line) is compared with simulations where the DM interaction is switched off (blue line) and the four-fold magnetic anisotropy is switched off (red line).(e) M c NL (ω t , τ ) at a fixed inter-pulse delay of τ = 4.8 ps.The result of the full quantum spin simulation (red line) is shown together with the result of the classical simulation (blue line).The quantum spin simulation yields stronger high-harmonic generation peaks (vertical dashed lines) compared to the classical simulation.Inset: Eigenenergies E n of the unperturbed quantum spin Hamiltonian which has (2s + 1) 2 = 36 eigenstates for the studied two-site spin-s = 5/2 system.(f) Ratio between the strength of nth harmonic and fundamental harmonic generation peaks.The ratios extracted from the experimental mHHG spectrum in Fig.1e(squares) are compared with the results of the quantum spin simulation (circles) and the classical simulation (rectangles).
Polycrystalline sample Sm 0.4 Er 0.6 FeO 3 was prepared by the conventional solid state reaction method, using the high purity oxide powder Sm 2 O 3 (99.99%),Er 2 O 3 (99.99%)and Fe 2 O 3

Table 1 .
. Parameters used in the classical simulations

Table 2 .
Parameters used in the quantum spin simulations signal nonlinear process frequency space MR ω 1