High-harmonic generation in quantum spin systems

TL;DR

This paper explores high-harmonic generation in one-dimensional quantum spin systems driven by magnetic fields, revealing how spin excitations influence the radiation spectrum and indicating potential for probing magnetic materials.

Contribution

It introduces a theoretical framework for understanding HHG in spin chains and links spectral features to spin excitations and phase transitions.

Findings

Plateau structures in HHG spectra explained by magnon processes.

Magnetic field induces phase transition affecting HHG cutoff energies.

High-harmonic signals can probe spin excitation spectra.

Abstract

We theoretically study the high-harmonic generation (HHG) in one-dimensional spin systems. While in electronic systems the driving by AC electric fields produces radiation from the dynamics of excited charges, we consider here the situation where spin systems excited by a magnetic field pulse generate radiation via a time-dependent magnetization. Specifically, we study the magnetic dipole radiation in two types of ferromagnetic spin chain models, the Ising model with static longitudinal field and the XXZ model, and reveal the structure of the spin HHG and its relation to spin excitations. For weak laser amplitude, a peak structure appears which can be explained by time-dependent perturbation theory. With increasing amplitude, plateaus with well-defined cutoff energies emerge. In the Ising model with longitudinal field, the thresholds of the multiple plateaus in the radiation spectra can be explained by the annihilation of multiple magnons. In the XXZ model, which retains the $\mathbf{Z}_{2}$ symmetry, the laser magnetic field can induce a phase transition of the ground state when it exceeds a critical value, which results in a drastic change of the spin excitation character. As a consequence, the first cutoff energy in the HHG spectrum changes from a single-magnon to a two-magnon energy at this transition. Our results demonstrate the possibility of generating high-harmonic radiation from magnetically ordered materials and the usefulness of high-harmonic signals for extracting information on the spin excitation spectrum.

Equations

$$\mathcal{H}(t)={\mathcal{H}}_{\mathrm{spin}}-\mathbf{B}(t)\cdot {\mathbf{S}}_{\mathrm{tot}},$$

$$I\propto |{\omega }^2\mathbf{M}(\omega ){|}^2.$$

$$\{\begin{array}{l}\mathcal{H}(t)\to \mathcal{H}(t+T_{\mathrm{per}}/2),\\ (S^x,S^y,S^z)\to (-S^x,-S^y,S^z),\end{array}$$

$$M^x(t+T_{\mathrm{per}}/2)=-M^x(t),$$

$$M^y(t+T_{\mathrm{per}}/2)=-M^y(t),$$

$$M^z(t+T_{\mathrm{per}}/2)=M^z(t),$$

$$M^x(2n\Omega )\propto {\int }_0^{T_{\mathrm{per}}}dt{e}^{i2n\Omega t}{M}^x(t)$$

$$={\int }_0^{\frac{T_{\mathrm{per}}}{2}}dt{e}^{i2n\Omega t}(M^x(t)+M^x(t+T_{\mathrm{per}}/2))=0.$$

$$M^z((2n+1)\Omega )\propto {\int }_0^{T_{\mathrm{per}}}dt{e}^{i(2n+1)\Omega t}{M}^z(t)$$

$$={\int }_0^{\frac{T_{\mathrm{per}}}{2}}dt{e}^{i(2n+1)\Omega t}(M^z(t)-M^z(t+T_{\mathrm{per}}/2))=0.$$

$${\mathcal{H}}_{\mathrm{Ising}}=-J\sum _j{S}_j^z{S}_{j+1}^z-H{S}_{\mathrm{tot}}^z,$$

$$\displaystyle B^{x}(t)=\left\{\begin{array}[]{ll}\displaystyle B\sin^{2}\Big{(}\frac{\Omega t}{2N_{\mathrm{cyc}}}\Big{)}\cos(\Omega t)&(0<t<T_{\mathrm{f}})\\ 0&(\mathrm{otherwise})\end{array}\right.,$$

$${\chi }^{\alpha \beta }(r,t)=-i\vartheta (t)⟨[S_r^{\alpha }(t),S_0^{\beta }(0)]⟩,$$

$${\chi }^{\alpha \beta }(q,\omega )={\int }_{-\infty }^{\infty }dt\sum _r{e}^{i(\omega t-qr)}{\chi }^{\alpha \beta }(r,t).$$

$$\mathcal{H}={\mathcal{H}}_{\mathrm{Ising}}-B{S}_{\mathrm{tot}}^x$$

$$E_2(q)=$$

$$=$$

$$2c_1-2c_2|\cos (q/2)|\le E_2(q)\le 2c_1+2c_2|\cos (q/2)|.$$

$$m_{\mathrm{sub}}^{x(z)}(\omega ;t_{\ast })\equiv \int dt{e}^{i\omega t}{m}^{x(z)}(t)W_{\mathrm{G}}(t;t_{\ast }),$$

$$M^{x(z)}(t)=\sum _{m,n}{\alpha }_m^{\ast }(t){\alpha }_n(t)⟨{\Phi }_m(B(t))|S^{x(z)}|{\Phi }_n(B(t))⟩.$$

$$i{\partial }_t{\alpha }_n(t)=$$

$$-i\sum _{m\ne n}({\partial }_{t}B(t))F_{nm}(B(t)){\alpha }_m(t),$$

$${\tilde{\mathcal{H}}}_{\mathrm{Ising}}(t)=-2J{m}^z(t)S^z-H{S}^z-B(t)S^x,$$

$${\chi }_{\mathrm{F}}^{\alpha \beta }(r,t;{\alpha }_0)=-i\vartheta (t)⟨\Psi (0;{\alpha }_0)|[S_r^{\alpha }(t;{\alpha }_0),S_0^{\beta }]|\Psi (0;{\alpha }_0)⟩$$

$${\chi }_{\mathrm{F}}^{\alpha \beta }(q,\omega )=⟨{\chi }_{\mathrm{F}}^{\alpha \beta }(q,\omega ;{\alpha }_0){⟩}_{{\alpha }_0}.$$

$${\mathcal{H}}_{\mathrm{XXZ}}=\sum _j[J_{xy}(S_j^x{S}_{j+1}^x+S_j^y{S}_{j+1}^y)-J_z{S}_j^z{S}_{j+1}^z],$$

$$\mathcal{H}={\mathcal{H}}_{\mathrm{XXZ}}-B{S}_{\mathrm{tot}}^x.$$

$$A_n=|{\alpha }_n^{\ast }{\alpha }_0⟨{\Phi }_n|S_{\mathrm{tot}}^x|{\Phi }_0⟩|,$$

$$\mathcal{H}$$

$$+\frac{J_{xy}+J_z}{4}\sum _j({\tilde{S}}_j^{+}{\tilde{S}}_{j+1}^{+}+{\tilde{S}}_j^{-}{\tilde{S}}_{j+1}^{-})-B(t)S_{\mathrm{tot}}^x.$$

$${\otimes }_j|\uparrow {⟩}_j={\otimes }_j\frac{|S^x=1/2{⟩}_j+|S^x=-1/2{⟩}_j}{\sqrt{2}}.$$

$${\tilde{\mathcal{H}}}_{\mathrm{XXZ}}(t)=$$

$$-B(t)S^x.$$

$$\mathcal{H}(t)={\mathcal{H}}_{\mathrm{spin}}+V(t),$$

Full text

I Introduction

The dynamics induced by light-matter coupling is an important problem in optical physics as well as nonequilibrium condensed matter and statistical physics. The application of strong laser pulses to a broad range of materials, including metals, semiconductors, and superconductors, results in rich physics and new phenomena, such as collective excitations Matsunaga et al. (2014); Lu et al. (2017), the control of order parameters Kirilyuk et al. (2010); Kampfrath et al. (2013), and fundamental changes in material properties Nasu (2004); Oka and Aoki (2009); Lindner et al. (2011). In particular, the high-harmonic generation (HHG), which is a nonlinear optical phenomenon observed in periodically driven systems, is attracting interest because of the underlying nontrivial charge dynamics and its technological relevance for attosecond laser science and the spectroscopy of charge dynamics Krausz and Ivanov (2009); Cavalieri et al. (2007).

HHG has originally been observed and studied in atoms and molecular gases McPherson et al. (1987); Ferray et al. (1988). Its mechanism can be understood by the so-called three step model, where tunnel-ionization occurs in the presence of a strong electric field, the released electrons are accelerated by the periodic field and eventually recombine with the ionized atoms by emitting the high-harmonic light Corkum (1993); Lewenstein et al. (1994). Recently the interest in this field has been renewed by the observation of HHG in various solids, in particular band insulators Ghimire et al. (2011); Schubert et al. (2014); Luu et al. (2015); Vampa et al. (2015a); Langer et al. (2016); Hohenleutner et al. (2015); Ndabashimiye et al. (2016); Liu et al. (2017); You et al. (2017); Kaneshima et al. (2018); Luu and Wörner (2018). Although the HHG in this case also originates from the dynamics of excited charges, the spatially periodic arrangement of atoms in solids leads to qualitative differences compared to atomic gases. Theoretical studies assuming weak correlations or employing an effective single particle picture have been performed to discuss the origin of the HHG in these band insulators Golde et al. (2008); Ghimire et al. (2011); Kemper et al. (2013); Higuchi et al. (2014); Vampa et al. (2014); Wu et al. (2015); Tamaya et al. (2016); Vampa et al. (2015b); Luu and Wörner (2016); Otobe (2016); Ikemachi et al. (2017); Osika et al. (2017); Hansen et al. (2017); Tancogne-Dejean et al. (2017a, b); Ikemachi et al. (2018); Ikeda et al. (2018). (For recent reviews, see Refs. Huttner et al. (2017); Kruchinin et al. (2018); Ghimire and Reis (2019).) It has been revealed that HHG originates from the intraband charge dynamics reflecting the non-parabolic shape of the bands Luu et al. (2015); Golde et al. (2008); Ghimire et al. (2011) and the interband dynamics corresponding to the recombination of excited charges Vampa et al. (2014, 2015b); Otobe (2016); Ikemachi et al. (2017); Osika et al. (2017). Furthermore, the existence of multiple bands and the interference between different excitation paths can play an important role Wu et al. (2015); Luu and Wörner (2016); Hansen et al. (2017); Ikemachi et al. (2017). Even though the details of its mechanism are still actively discussed, HHG in solids can be used to obtain important information about these solids, such as band and lattice structures Ghimire et al. (2011); Ndabashimiye et al. (2016); Liu et al. (2017); Hohenleutner et al. (2015); You et al. (2017). In addition, potential applications in new high-frequency laser sources are expected due to the high concentration of atoms compared to atomic gases Ndabashimiye et al. (2016). Stimulated by these developments and prospects, both experimentalists and theorists are making intensive efforts to understand the mechanism of HHG in greater detail and to explore new classes of materials, e.g., liquids Heissler et al. (2014); Luu and Wörner (2018), graphene Yoshikawa et al. (2017); Hafez et al. (2018), topological systems Chacón et al. (2018), strongly correlated systems Silva et al. (2018); Murakami et al. (2018); Murakami and Werner (2018); Tancogne-Dejean et al. (2018); Zhu et al. (2018), impurity-doped systems Yu et al. (2019) and magnetic metals Zhang et al. (2018).

In this paper, we explore a new avenue for HHG in solids, by considering the dynamics of the spin degrees of freedom in magnetic insulators, i.e. quantum spin systems. We theoretically study the excitation of these systems by time-periodic external magnetic fields, and evaluate the HHG signal resulting from the change of the magnetic moments [Fig. 1(a)]. This setup is relevant for materials with a large charge excitation gap whose low energy excitations are governed by the spin degrees of freedom. Recent developments in the field of metamaterials Mukai et al. (2014) and plasmonics Ciappina et al. (2017) enable the generation of strong magnetic field pulses with small electric field, which can be used to realize the setup considered in our study. The nonequilibrium dynamics of quantum spin systems, especially the dynamical control of the magnetization by laser fields, has been intensively studied both in the experimental Hohlfeld et al. (1997); Kimel et al. (2005); Kirilyuk et al. (2010) and theoretical communities Takayoshi et al. (2014a, b). In Refs. Takayoshi et al. (2014a, b), the magnetization dynamics in antiferromagnets has been calculated for laser fields with a frequency comparable to the exchange coupling. On the other hand, for the study of HHG, a lower photon energy is advantageous since it results in spectra with higher energy resolution and thus allows to elucidate the excitation structure. In this paper, we reveal that the HHG signal from spin systems can be associated with elementary spin excitations like magnons, just as the HHG in electronic systems reflects the dynamics of excited charges [see Fig. 1(b)]. These results suggest that the spin HHG can be potentially used as a probe of spin dynamics as well as for new laser sources in the THz regime. In Fig. 1(c), we summarize the similarities and differences to HHG in electronic systems, which are useful to keep in mind in the following discussion.

The present study focuses on one-dimensional ferromagnetic quantum spin systems described by the Ising model with longitudinal field and the XXZ model. These models are simple but fundamental, and can be realized in materials such as $\mathrm{Dy(C_{2}H_{5}SO_{4})_{3}\cdot 9H_{2}O}$, $\mathrm{LiTbF_{4}}$, $\mathrm{LiHoF_{4}}$  Wolf (2000) and $\mathrm{CoNb_{2}O_{6}}$ Coldea et al. (2010). We numerically investigate the nonequilibrium dynamics and the radiation spectrum resulting from the time-dependent magnetization by means of the infinite time-evolving block decimation (iTEBD) Vidal (2007), exact diagonalization (ED) calculations, and time-dependent mean-field theory (tdMF). To understand the relation between the HHG signal and elementary spin excitations, we also calculate the low-energy excitation structure of these systems by combining the density matrix renormalization group (DMRG) White (1992) and the time-evolving block decimation (TEBD) Vidal (2003). When the laser field is weak, a peak appears around the energy of the single-magnon excitation in both models, which can be explained by time-dependent perturbation theory. With increasing strength of the laser field, this peak structure changes to a plateau. We also find indications for multiple plateaus, whose thresholds are associated with the annihilation of (multiple) elementary spin excitations (magnons).

This paper is organized as follows. In Sec. II, we discuss general properties of the HHG in quantum spin systems. Section III presents the HHG signals resulting from the application of a linearly polarized laser to Ising models with longitudinal static field. Section IV is devoted to an analysis of the HHG signal from the laser application to the XXZ models. We summarize our results and discuss future extensions in Sec. V.

II HHG in quantum spin systems

In this section, we present the theory of HHG in quantum spin systems. In usual HHG, the electric field of a laser pulse induces a change of the electric polarization, which in turn produces electromagnetic waves. The total instantaneous radiated power is proportional to $|d\boldsymbol{j}(t)/dt|^{2}$, where $\boldsymbol{j}(t)$ is the electric current. If $\boldsymbol{j}(t)$ is a polarization current $d\boldsymbol{P}(t)/dt$ (with $\boldsymbol{P}$ the electric interband polarization), the power is proportional to $|d^{2}\boldsymbol{P}(t)/dt^{2}|^{2}$. In a similar way, we can consider the radiation of electromagnetic waves from a time-dependent magnetic dipole. The total instantaneous radiated power from the change of a localized magnetic dipole $\boldsymbol{M}(t)$ is proportional to $|d^{2}\boldsymbol{M}(t)/dt^{2}|^{2}$ Jackson (1998).

To study quantum spin systems in the presence of a time-dependent magnetic field $\boldsymbol{B}(t)$ we consider the Hamiltonian

[TABLE]

where $\mathcal{H}_{\mathrm{spin}}$ is the spin Hamiltonian and the last term represents the Zeeman coupling of the spins in the material with the magnetic field produced by the laser. We calculate the time evolution of the magnetization $\boldsymbol{M}(t)\equiv\langle\boldsymbol{S}_{\mathrm{tot}}(t)\rangle$, where $\boldsymbol{S}_{\mathrm{tot}}=\sum_{j}\boldsymbol{S}_{j}$ represents the summation over all spins and $\boldsymbol{S}_{\mathrm{tot}}(t)\equiv U^{-1}(t)\boldsymbol{S}_{\mathrm{tot}}U(t)$ ($U(t)=\mathcal{T}\int_{0}^{t}dt^{\prime}\exp[-i\mathcal{H}(t^{\prime})t^{\prime}]$ is the time evolution operator with $\mathcal{T}$ the time ordering). From this we obtain the Fourier transform of the magnetization $\boldsymbol{M}(\omega)=\int dte^{i\omega t}\boldsymbol{M}(t)$ and the radiation power

[TABLE]

The symmetry of the system may impose constraints on the structure of the HHG signal. For example, the inversion symmetry limits the HHG signal in electronic systems to odd harmonics. Now, let us consider the case when the time dependent Hamiltonian has a symmetry which can be represented as the combination of time translation and $\pi$ rotation around the $S^{z}$ axis

[TABLE]

where $T_{\mathrm{per}}\equiv 2\pi/\Omega$ is the period of the laser. Then, the magnetization satisfies

[TABLE]

if we assume a time-periodic steady state having the same symmetry as the Hamiltonian. In this case, the temporal Fourier transform of $M^{x}$ becomes 0 for $\omega=2n\Omega$ ($n$ is an integer) since

[TABLE]

In the same way, the temporal Fourier transform of $M^{z}$ becomes 0 for $\omega=(2n+1)\Omega$ ($n$ is an integer) since

[TABLE]

In the case of a finite pulse width, these arguments are strictly speaking not valid. Still we will see that in practice, these rules are satisfied except around the HHG peak in the case of weak laser fields.

In the following two sections, we will use numerical calculations to study the HHG in specific one-dimensional quantum spin systems.

III HHG in Ising models

Let us start by investigating the HHG in Ising models, which are among the simplest and most important models of magnets. In this case, the spin Hamiltonian $\mathcal{H}_{\mathrm{spin}}$ in Eq. (1) explicitly reads

[TABLE]

where $J>0$ is the ferromagnetic exchange coupling and $H>0$ is a static external magnetic field. $S^{x}$, $S^{y}$, and $S^{z}$ are spin-1/2 operators. The ground state of $\mathcal{H}_{\mathrm{Ising}}$ is a ferromagnetic state ($\langle S_{j}^{z}(0)\rangle=1/2$ for all $j$ ) and this state is perturbed by the application of a linearly polarized pulse laser $\boldsymbol{B}(t)=(B^{x}(t),0,0)$ in the $x$ direction. Because of the longitudinal field $H>0$, the $\mathbf{Z}_{2}$ symmetry of the system is broken. Hence, there is no quantum phase transition as a function of the transverse magnetic field, i.e., the ground state of the snapshot Hamiltonian $\mathcal{H}(t)=\mathcal{H}_{\mathrm{Ising}}-B^{x}(t)S_{\mathrm{tot}}^{x}$ remains gapped at any time. Though the main objective of this section is the theoretical analysis of the magnetization dynamics and HHG mechanism, the obtained results are relevant for materials having ferromagnetic dipole order such as $\mathrm{Dy(C_{2}H_{5}SO_{4})_{3}\cdot 9H_{2}O}$, $\mathrm{LiTbF_{4}}$ and $\mathrm{LiHoF_{4}}$ Wolf (2000).

We consider a magnetic field pulse of the form

[TABLE]

where $T_{\mathrm{f}}=2\pi N_{\mathrm{cyc}}/\Omega$, $\Omega$ is the laser frequency, $N_{\mathrm{cyc}}$ the number of laser cycles, and $B\sin^{2}(\frac{\Omega t}{2N_{\mathrm{cyc}}})$ the envelope of the pulse. In this paper, the parameters are fixed as $N_{\mathrm{cyc}}=9$ and $\Omega=1$ ($\Omega$ is also used as the energy scale by employing the units $\hbar=c=1$). The magnetic field pulse with $B=4$ is shown in the inset of Fig. 1(b).

In this section, the other parameters are set to $J=2$ and $H=6$, so that the gap is much larger than $\Omega=1$ and heating effects are suppressed. In addition, since we anticipate that the width of the plateau in the HHG signal is of the order of the characteristic energy scales of the spin system, we expect to observe several harmonics if $J$ and $H$ are chosen large compared to $\Omega$. The Ising model with smaller longitudinal field $H$, where the lifting of the two-fold degeneracy and hence the gap is smaller, is discussed in Appendix C.1. We numerically calculate the magnetization dynamics and report hereafter the normalized magnetizations $m^{x,y,z}\equiv M^{x,y,z}/N$, where $N$ is the number of spins. As explained in Sec. II, the radiation power of a magnetic dipole is proportional to $|\omega^{2}m^{\alpha}(\omega)|^{2}$.

Before studying the dynamics induced by the laser field, we investigate the excitation structure of the equilibrium system. To study excitations, we numerically calculate the dynamical structure factor (DSF), which is the imaginary part of the dynamical susceptibility. The method is as follows. We first obtain the ground state of the system by the DMRG White (1992), and then calculate the retarded correlation function

[TABLE]

where $\vartheta(t)$ is the step function, by the TEBD method Vidal (2003) for finite size systems. The dynamical susceptibility is the Fourier transform of the retarded correlation function,

[TABLE]

In this paper, we consider systems with size $N=120$, which are large enough that finite size effects can be neglected.

The DSFs $|\mathrm{Im}\chi^{xx}(q,\omega)|$ and $|\mathrm{Im}\chi^{zz}(q,\omega)|$ for the ground state of the Ising model in both longitudinal and transverse fields

[TABLE]

with $H=6$ are shown in Fig. 2. Equation (10) represents the snapshot Hamiltonian of $\mathcal{H}(t)=\mathcal{H}_{\mathrm{Ising}}-B^{x}(t)S_{\mathrm{tot}}^{x}$ at some fixed time $t$ corresponding to $B^{x}(t)=B$. If the transverse field is not present ($B=0$), the spins are completely localized and there is no dispersion since the Hamiltonian only contains $S^{z}$. The elementary excitation corresponds to a single spin flip, which has a gap $J+H$. In the presence of a nonzero transverse field, this flipped spin can propagate and transform into a magnon. The DSF shown in Fig. 2 represents the magnon dispersion. In Figs. 2(a)-2(c), a weak intensity is seen at twice of the energy of the lowest band (single-magnon dispersion). This corresponds to the two-magnon band. Since the single-magnon band has a cosine structure $E_{1}(q)=c_{1}+c_{2}\cos(q)$, the two-magnon band can be represented as

[TABLE]

by considering the momentum conservation, and we obtain

[TABLE]

This feature of the two-magnon band is observed more evidently when the longitudinal field $H$ is weak as mentioned in Appendix C.1.

In Fig. 3, we show the time evolution of $m^{x}$ and $m^{z}$ for the Hamiltonian $\mathcal{H}(t)=\mathcal{H}_{\mathrm{Ising}}-B^{x}(t)S_{\mathrm{tot}}^{x}$ with different values of the laser amplitude $B=2,4,6,8$. As the numerical method, we use the iTEBD Vidal (2007), which utilizes a matrix product state (MPS) representation. This method enables the simulation of infinite size systems, i.e., without finite-size effects, by assuming the translational invariance of the system. In this paper, we take the matrix dimension of the MPS as 100 and the time evolution is performed by the fourth-order Trotter decomposition with the time step $\Delta t=0.05$. The shape of the time evolving $m^{x}$ is similar to that of the applied laser magnetic field (Eq. (8)) for all values of the laser amplitude. The value of $m^{z}$ drops when $|m^{x}|$ grows, but otherwise the magnetization in the $z$ direction recovers to $m^{z}=1/2$. This demonstrates that the state of the system closely follows the ground state of the instantaneous Hamiltonian at each time. In other words, the laser frequency is slow enough for an adiabatic time evolution of the magnetization. In the present case of $H\gg J$, the gap is large ($\gg\Omega$) even for $B=0$, and it increases monotonically with increasing $B$ (Fig. 2). Thus, transitions to excited states through the Landau-Zener process are suppressed and the state remains in the snapshot ground state. However, if $B$ is further increased, the chain will eventually be disordered after the laser application, similarly to what is shown in Fig. 15(b) in Appendix C.1.

To investigate the HHG, we plot $|\omega^{2}m^{x}(\omega)|^{2}$ and $|\omega^{2}m^{z}(\omega)|^{2}$ on a logarithmic scale in Fig. 4. These spectra were obtained by first differentiating $m^{x(z)}(t)$ numerically as ${m^{\prime\prime}}^{x(z)}(t)=[m^{x(z)}(t+\Delta t)+m^{x(z)}(t-\Delta t)-2m^{x(z)}(t)]/(\Delta t)^{2}$, where $\Delta t$ is the time step, and then performing the Fourier transform. In the Fourier transform, we apply the Blackman window $W_{\mathrm{B}}(t)=0.42-0.5\cos(2\pi t/T_{\mathrm{f}})+0.08\cos(4\pi t/T_{\mathrm{f}})$ $(0<t<T_{\mathrm{f}})$ and $W_{\mathrm{B}}(t)=0$ (otherwise). The result in Fig. 4 clearly demonstrates the HHG for all values of $B$ in both magnetization components $m^{x}$ and $m^{z}$. Since the system satisfies the symmetry Eq. (2), $m^{x}(\omega)$ and $m^{z}(\omega)$ become 0 at $\omega=2n\Omega$ and $\omega=(2n+1)\Omega$ ($n$ is an integer), respectively, for steady states. Although the presented results are for the transient case, the magnitudes of $|\omega^{2}m^{x}(\omega)|^{2}$ and $|\omega^{2}m^{z}(\omega)|^{2}$ drop at $\omega=2n\Omega$ and $\omega=(2n+1)\Omega$, respectively. An exception occurs when $B$ is small and $\omega$ around the value corresponding to the excitation gap, as can be seen in Figs. 4(a), 4(b), 4(e), and 4(f) where the spectra exhibit peaks at $\omega=6,8,10$ in $m^{x}(\omega)$ and at $\omega=7,9$ in $m^{z}(\omega)$ for $B=2,4$. Since we consider the application of a laser pulse, the system is in a transient regime and does not reach a nonequilibrium steady state. Hence the conditions Eqs. (3) and (4) are not necessarily satisfied. The result for the peak position of the HHG spectra is supported by time-dependent perturbation theory (Appendix A). The peaks resulting from the perturbation theory are located at $\omega=H+J$ for $m^{x}$ and at $\omega=H+J\pm\Omega$ for $m^{z}$, i.e., they can appear at an arbitrary frequency (not necessarily an integer multiple of $\Omega$) depending on the values of $H$ and $J$. The validity of the time-dependent perturbation theory is also confirmed by the scaling of the radiation intensity with the laser amplitude $B$. In Fig. 5, we plot $|m^{x}(\omega)|$ and $|m^{z}(\omega)|$ in the region of small $B$. $|m^{x}(\omega)|$ and $|m^{z}(\omega)|$ at $\omega=n\Omega$ scale as $B^{n}$, while at $\omega=H+J\pm n\Omega$ they scale as $B^{n+1}$, which agrees with the prediction from the perturbation theory presented in Appendix A. This result indicates that $B\leq 2$ is in the perturbative regime.

In Fig. 4, when the field strength is sufficiently large, we can identify a frequency above which the intensity drops rapidly as well as multiple plateau structures. We can connect these cut-off energies with the excitation structures of the snapshot Hamiltonians, in particular those with the maximum value of $B$. In the dispersion relation obtained from the data in Fig. 2, the energy has a minimum (maximum) at $q=0$ ($q=\pi$), and the excitation gap corresponds to the mass of a magnon at $q=0$. We see that the intensity of $|\omega^{2}m^{x}(\omega)|^{2}$ and $|\omega^{2}m^{z}(\omega)|^{2}$ drops above the energies corresponding to integer multiples of the magnon mass at $q=0$, as indicated by the dashed lines in Fig. 4. This result suggests that for sufficiently large laser field amplitude, there occurs a spontaneous annihilation of $n(=1,2,3,\ldots)$ magnons, which leads to the emission of light with the frequency $n\Delta(B(t))$ at time $t$, where $\Delta$ represents the single-magnon energy gap. This situation is analogous to electron-hole or doublon-holon recombination in electron systems such as Mott insulators Vampa et al. (2014); Murakami et al. (2018), where the radiation originates primarily from the interband transitions.

Further insights can be obtained from a subcycle analysis. The subcycle Fourier transform of the magnetic moment is defined as

[TABLE]

where W_{\mathrm{G}}(t;t_{*})=\exp\big{[}-\frac{(t-t_{*})^{2}}{2\sigma^{2}}\big{]} ($\sigma=T_{\mathrm{per}}/8$) is a Gaussian window function. (An alternative way to compute time-dependent spectra is the wavelet analysis. We discuss the result of the wavelet analysis and the difference to the window Fourier transform in Appendix E.) In Fig. 6, we show the subcycle radiation spectrum $\log_{10}|\omega^{2}m_{\mathrm{sub}}^{x(z)}(\omega;t_{*})|^{2}$ for $B=8$ as a colormap and the multiple magnon excitation energies of the snapshot Hamiltonian at $t_{*}$ by the solid lines. In the low-energy region ($\omega<10$), although $|\omega^{2}m_{\mathrm{sub}}^{x(z)}(\omega;t_{*})|^{2}$ does not much depend on $t_{*}$, one can roughly identify an enhanced HHG signal following the one-magnon energy. This is due to the fact that the single-magnon band changes only little as a function of the transverse field (see Fig. 2). On the other hand, in the high-energy region, a high intensity signal is produced when the magnetic field is strong. In particular, we can clearly identify an enhanced HHG signal tracking the two-magnon and three-magnon lines, both in the radiation produced by the $x$ and $z$ magnetization components. These observations support the interpretation that the plateaus and their thresholds in the spin HHG originate from the annihilation of magnons.

We note that our discussion of the spin HHG so far has been based on the eigenstates or the energy structure of the snapshot Hamiltonians, as has been done for electronic systems using the Houston basis Wu et al. (2015) or assuming a slowly changing field Higuchi et al. (2014); Murakami et al. (2018). To be more specific, let us expand the wave function as $|\Psi(t)\rangle=\sum_{n}\alpha_{n}(t)|\Phi_{n}(B(t))\rangle$, where $|\Phi_{n}(B(t))\rangle$ is an eigenstate of the snapshot Hamiltonian with the eigenenergy $E_{n}(B(t))$, and express the magnetization as

[TABLE]

We can then classify the contributions to the magnetization dynamics according to the character of $|\Phi_{m}(B(t))\rangle$ and $|\Phi_{n}(B(t))\rangle$. The time dependence of the coefficients $\alpha_{n}$ follows from

[TABLE]

where $F_{nm}(B)=\langle\Phi_{n}(B)|\partial_{B}|\Phi_{m}(B)\rangle$. If the variation of $B(t)$ (with excitation frequency $\Omega$) is slow enough, $\partial_{t}B(t)$ is small and $E_{n}(B(t))$ can be approximated as a constant for a certain time interval. Hence the second term on the right hand side of Eq. (14) can be neglected and we can write $\alpha_{n}(t)\propto e^{-iE_{n}(B(t_{*}))t}$ for $t$ around $t_{*}$. If these approximations hold and the time-dependence of $|\Phi_{n}(B(t))\rangle$ (and hence that of $\langle\Phi_{m}(B(t))|S^{x(z)}|\Phi_{n}(B(t))\rangle$) is also small enough, the main contribution to $M^{x(z)}(t)$ [Eq. (13)] is proportional to $e^{-i[E_{n}(B(t_{*}))-E_{m}(B(t_{*}))]t}$ for $t$ around $t_{*}$, which oscillates with (multiple) magnon energies. If $|\Phi_{n}(B)\rangle$ and $|\Phi_{m}(B)\rangle$ differ by $l$ magnons, the radiation can be interpreted as originating from an $l$-magnon annihilation. However, in practice, there may be contributions from the second term on the right hand side of Eq. (14) and the time-dependence of $|\Phi_{n}(B(t))\rangle$, which leads to deviations from the simple magnon picture. Furthermore, the magnetization curve of the ground states for the Hamiltonian Eq. (10) is a nonlinear function of $B$. Since $m^{x}$ for the ground state with the field $B$ is an odd function, we see, by replacing $B$ in this equation by $B\cos(\Omega t)$, that the Fourier component of $n\Omega$ (with $n$ an odd integer) appears in $m^{x}(\omega)$ and its leading order is $B^{n}$. This partially explains the appearance of well-defined frequency components even in the energy region lower than the excitation gap seen in Fig. 4.

The above results suggest that for the parameters chosen in this study, the magnon picture is essentially valid and the dynamics is described in terms of well-ordered magnetic moments, i.e. the effect of quantum fluctuations is small. To confirm this point, we perform a tdMF analysis. The approximation $\sum_{j}S_{j}^{z}S_{j+1}^{z}\simeq 2m^{z}\sum_{j}S_{j}^{z}-N{m^{z}}^{2}$ leads to the tdMF Hamiltonian

[TABLE]

where $m^{z}(t)\equiv\langle S^{z}(t)\rangle$. We solve the Schrödinger equation with the Hamiltonian (15) by the fifth order Runge-Kutta method with the Cash-Karp parameters, and calculate the dynamics of $m^{x}(t)\equiv\langle S^{x}(t)\rangle$ and $m^{z}(t)$. The discretized time step is $\Delta t=0.05$. The result is also shown in Fig. 4. The curves of $|\omega^{2}m^{x}(\omega)|^{2}$ and $|\omega^{2}m^{z}(\omega)|^{2}$ calculated by the single spin dynamics agree well with those calculated by iTEBD up to the first HHG threshold. Note that there is no rescaling of the results and the agreement is quantitative. The deviations become larger above the first threshold. This indicates that correlations between magnons beyond mean-field theory are essential for the spontaneous recombination of multiple magnons.

Another useful perspective on HHG can be obtained from the Floquet picture Ikeda et al. (2018); Murakami and Werner (2018). The spectrum in the Floquet theory is derived from the Floquet DSF $|\chi_{\mathrm{F}}^{\alpha\beta}(q,\omega)|$, which is calculated in a similar way as $|\chi^{\alpha\beta}(q,\omega)|$. Let us consider the time-dependent Hamiltonian $\mathcal{H}(t;\alpha_{0})=\mathcal{H}_{0}-B\sin(\Omega t+\alpha_{0})S_{\mathrm{tot}}^{x}$ and represent the ground state of $\mathcal{H}(0;\alpha_{0})$ by $|\Psi(0;\alpha_{0})\rangle$, where $\alpha_{0}$ is the phase shift. We calculate the Floquet retarded correlation function

[TABLE]

[cf. Eq. (9)], where $S_{r}^{\alpha}(t;\alpha_{0})=U^{-1}(t;\alpha_{0})S_{0}^{\alpha}U(t;\alpha_{0})$ and $U(t;\alpha_{0})=\mathcal{T}\int_{0}^{t}dt^{\prime}e^{-i\mathcal{H}(t^{\prime};\alpha_{0})t^{\prime}}$. The DSF $\chi_{\mathrm{F}}^{\alpha\beta}(q,\omega;\alpha_{0})$ is defined as the Fourier transform of this correlation function, and we take the average relative to the phase shift $\alpha_{0}$ over a single cycle as

[TABLE]

Here we take $\alpha_{0}=n\pi/8$ ($n=0,1,\ldots,15$). In Fig. 7, we show the Floquet DSF $|\mathrm{Im}\chi_{\mathrm{F}}^{xx(zz)}(q,\omega)|$ for $B=8$. We can see the appearance of Floquet subbands with an energy splitting of $2\Omega$ rather than $\Omega$. The subbands of $|\mathrm{Im}\chi_{\mathrm{F}}^{xx}(q,\omega)|$ are located at $(\mathrm{one\;magnon\;band})\pm(\mathrm{odd\;integer})\Omega$ while those of $|\mathrm{Im}\chi_{\mathrm{F}}^{zz}(q,\omega)|$ are located at $(\mathrm{one\;magnon\;band})\pm(\mathrm{even\;integer})\Omega$. In $|\mathrm{Im}\chi_{\mathrm{F}}^{xx}(q,\omega)|$, the Floquet subbands of the negative energy magnon dispersion appear around $\omega\simeq 1$. These Floquet DSFs suggest that we can also interpret the high-harmonic peaks with energy below the magnon mass in terms of transitions between Floquet sidebands of the magnon spectrum.

In this section, we have focused on the model with strong longitudinal field ($H\gg J$) and large gap. With decreasing $H$, the gap decreases and the HHG behavior changes. We discuss the results of a weak $H$ model ($J=4$, $H=2$) in Appendix C.1. The emergence of HHG plateaus with a close relation to magnon energies can also be observed there.

IV HHG in XXZ models

In this section we consider another fundamental model of quantum magnets, the ferromagnetic XXZ model. The spin Hamiltonian is

[TABLE]

where $J_{z}>J_{xy}>0$. The difference from the Ising model is the term $S_{j}^{x}S_{j+1}^{x}+S_{j}^{y}S_{j+1}^{y}=\frac{1}{2}(S_{j}^{+}S_{j+1}^{-}+S_{j}^{-}S_{j+1}^{+})$, which acts as a kinetic term for the magnons. Note that the spins are completely frozen in the Ising model without laser. The ground state of Eq. (17) is a ferromagnetic state for $J_{z}>J_{xy}>0$ while it is a gapless Luttinger liquid for $|J_{z}|<|J_{xy}|$ Giamarchi (2004). The low energy excitations of Eq. (17) are magnons with dispersion $E(q)=J_{xy}\cos(q)+J_{z}$. Since the Hamiltonian Eq. (17) does not include the longitudinal static field $HS_{\mathrm{tot}}^{z}$, the system has a $\mathbf{Z}_{2}$ symmetry, and thus a quantum phase transition can be induced by applying a transverse field.

Here we consider the case where $J_{xy}$ is weak, $J_{xy}\ll J_{z}$, which is relevant for the modeling of quasi-one-dimensional magnetic insulators such as $\mathrm{CoNb_{2}O_{6}}$ Coldea et al. (2010). The parameters are set to $J_{xy}=2$, $J_{z}=10$, and $\Omega=1$. We take both $J_{xy}$ and $J_{z}$ to be larger than $\Omega$ so that the HHG plateau contains several harmonics. For the analysis of the model with strong $J_{xy}$ ($J_{xy}\lesssim J_{z}$), see Appendix C.2. In Fig. 8, we show the DSF in the ground state of the XXZ model with a transverse field $B$, which corresponds to the snapshot Hamiltonian of the system under laser irradiation,

[TABLE]

In contrast to the case of the Ising model, the low-energy excitation spectrum is continuous due to the existence of the kinetic term. The lower bound of the dispersion at $q=0$ decreases with increasing $B$. The gap closes and a phase transition happens at $B_{\mathrm{c}}\simeq 6$. Before the transition ($B<B_{\mathrm{c}}$), $\chi^{xx}$ shows a stronger intensity than $\chi^{zz}$, because the spins are primarily aligned in the $z$ direction in the ground state. When $B$ is small enough, the DSF has a strong intensity near the one magnon dispersion for $B=0$ (i.e., $E(q)=J_{xy}\cos(q)+J_{z}$), and in particular the strongest intensity is found at $q=\pi$. On the other hand, after the transition ($B>B_{\mathrm{c}}$), the intensity of $\chi^{zz}$ becomes much stronger than $\chi^{xx}$, because the spins are mainly aligned in the $x$ direction in the ground state, and the strongest intensity is observed at $q=0$. The dispersion captured by $\chi^{zz}$ is sharp, and it can be interpreted as a single-magnon band in terms of the spin wave theory (see Appendix B). We also note that the upper bound of the continuous dispersion at $q=0$ captured by $\chi^{xx}$ corresponds to a two-magnon state since its energy is twice the excitation energy at $q=\pi$ captured by $\chi^{zz}$ for $B>B_{\mathrm{c}}$ .

The time evolution of $m^{x}$ and $m^{z}$ calculated by iTEBD is shown in Fig. 9. The time evolution of $m^{x}$ essentially tracks the laser magnetic field Eq. (8) for small $B$, but the shape changes especially near the peaks of the intensity as $B$ is increased. Higher frequency components than $\Omega$ appear near the peaks, and these contribute to the HHG (see the sub-cycle analysis below). The time evolution of $m^{z}$ drastically changes its behavior depending on whether $B$ is smaller or larger than $B_{\mathrm{c}}$. For $B<B_{\mathrm{c}}$, the magnitude of $m^{z}$ decreases when the laser intensity is strong, otherwise $m^{z}\simeq 1/2$, which demonstrates that the state follows the ground state of the snapshot Hamiltonian, i.e., the time evolution is almost adiabatic. However, for $B>B_{\mathrm{c}}$, $m^{z}$ suddenly decreases from $1/2$, which shows that the system makes transitions to excited states of the snapshot Hamiltonian.

The HHG spectra $|\omega^{2}m^{x}(\omega)|^{2}$ and $|\omega^{2}m^{z}(\omega)|^{2}$ are shown in Fig. 10. Here the same Blackman window is used as in the Ising case. The HHG structure is clear for the weak field $B$ while it is noisier after the transition. Since the system satisfies the symmetry Eq. (2), the magnitudes of $|\omega^{2}m^{x}(\omega)|^{2}$ and $|\omega^{2}m^{z}(\omega)|^{2}$ drop at $\omega=2n\Omega$ and $\omega=(2n+1)\Omega$, respectively, except that $m^{x}(\omega)$ has a peak and $m^{z}(\omega)$ has a dip around $\omega=12$ for $B=2$. This energy corresponds to the upper bound of the single-magnon band $J_{xy}+J_{z}$, and we can explain the peaks at $\omega=J_{xy}+J_{z}$ for $m^{x}(\omega)$ and at $\omega=J_{xy}+J_{z}\pm\Omega$ for $m^{z}(\omega)$ in the small $B$ region in terms of the time-dependent perturbation theory as shown in Appendix A.

As we increase $B$ and leave the perturbative regime, plateau structures develop in the low-energy region. Again we can connect these cut-off energies (threshold energies) with the spin excitation structure. As depicted in Fig. 10, for $B<B_{\mathrm{c}}$, the threshold of the HHG plateau corresponds to $\Delta_{q=0}$, which is the upper bound of the dispersion obtained from $\chi^{xx}$ at $q=0$. For $B>B_{\mathrm{c}}$, the threshold of the first HHG plateau is determined by $2\tilde{\Delta}_{q=\pi}$, where $\tilde{\Delta}_{q=\pi}$ is the excitation gap corresponding to $\chi^{zz}$ at $q=\pi$. This energy scale is not very apparent in $|\omega^{2}m^{z}(\omega)|^{2}$ but we can see that $|\omega^{2}m^{x}(\omega)|^{2}$ is larger than $|\omega^{2}m^{z}(\omega)|^{2}$ by several orders near the threshold energy [dashed-dotted lines in Figs. 10(d), 10(e), 10(i), and 10(j)] and dominates the HHG. Note that for $B>B_{\mathrm{c}}$, the spins are mostly aligned in the $x$ direction in the ground state and $S^{z}$ works as a spin-flip (magnon generation) operator. Even though $\Delta_{q=0}=2\tilde{\Delta}_{q=\pi}$ and this mode can also be excited by the $S^{x}$ operator, the intensity of $\chi^{zz}$ is much larger than that of $\chi^{xx}$ as seen in Fig. 8. Thus it is more natural to regard it as a two-magnon process.

In the same way as we have done for the Ising model, we can obtain further insight into the origin of the HHG by performing a subcycle analysis for the XXZ model. In Fig. 11, we show the subcycle radiation spectrum Eq. (12) for $B=4$ and $B=10$ (below and above the critical field, respectively). In the case of $B=4$ [Fig. 11(a)], the strong intensity in the HHG signal follows the single-magnon and two-magnon excitation energy ($\Delta_{q=0}$ and $2\Delta_{q=0}$) of the snapshot Hamiltonian at each time, which suggests that again the threshold can be associated with the annihilation of multiple magnons at $q=0$ for $B<B_{\mathrm{c}}$. In the case of $B=10$ [Fig. 11(b)], the strong intensity in the HHG signal follows the two-magnon excitation energy ($2\tilde{\Delta}_{q=\pi}$). There is also some additional intensity in the energy range $\omega=25$ – 40 in Fig. 11(b), which may correspond to higher order excitations such as four magnon processes.

To confirm that the threshold of the HHG plateaus corresponds to the magnetic excitation structure, especially magnon modes, we perform an ED calculation for a system of $N=8$ sites. The system size is small, but the ED calculations reproduce quantitatively the behavior of the HHG spectra for small $B$ as can be seen in Fig. 10. Although there is a quantitative deviation from the iTEBD results in the case of strong $B$, the HHG signals show a qualitative agreement. In particular, the threshold energy of the first plateau is the same for ED and iTEBD. We denote the eigenstates of the snapshot Hamiltonian at the time when the laser intensity takes the maximum ($t_{\mathrm{peak}}=\pi N_{\mathrm{cyc}}/\Omega$) by $|\Phi_{n}\rangle$ and their eigenenergies by $E_{n}$. In Fig. 12(a), we show $\langle\Phi_{n}|S_{\mathrm{tot}}^{x}|\Phi_{n}\rangle$ calculated by ED for large $B$. In the present model, even though $S_{\mathrm{tot}}^{x}$ is not a conserved quantity, the spins basically align in the $S^{x}$ direction in the ground state for large $B$ and the expectation values $\langle\Phi_{n}|S_{\mathrm{tot}}^{x}|\Phi_{n}\rangle$ are almost discretized and distributed around integer values. The expectation values near $-2$ are highlighted with cross markers in Fig. 12. From Fig. 12(a), the energy threshold of the first HHG plateau corresponds to the upper bound of the $\langle\Phi_{n}|S_{\mathrm{tot}}^{x}|\Phi_{n}\rangle\simeq-2$ sector ($E_{n}-E_{0}=24.5$). Since the ground state is in the $\langle\Phi_{n}|S_{\mathrm{tot}}^{x}|\Phi_{n}\rangle\simeq-4$ sector, two spins are flipped, i.e., two magnons are generated. In Fig. 12(b), we plot the quantity

[TABLE]

where $\alpha_{n}\equiv\langle\Phi_{n}|\Psi(t_{\mathrm{peak}})\rangle$ represents the overlap between the state at $t=t_{\mathrm{peak}}$ and the $n$-th excited state $|\Phi_{n}\rangle$ of the snapshot Hamiltonian ($|\Phi_{0}\rangle$ is the ground state). This quantity is directly related to $m^{x}$ through Eq. (13). We see that there is a strong intensity at the energy $E_{n}-E_{0}=24.5$, which agrees with the threshold energy in Fig. 10(e). Hence we can conclude that the threshold of the first HHG plateau is dictated by the two-magnon mode $2\tilde{\Delta}_{q=\pi}$. In addition, Fig. 12(b) suggests that the contribution to the HHG signal mainly comes from the two-magnon sector ($\langle\Phi_{n}|S_{\mathrm{tot}}^{x}|\Phi_{n}\rangle\simeq-2$).

Further insight into the HHG signal with large $B$ can be obtained by rewriting the Hamiltonian. Since the spin alignment axis is $S^{x}$ in the case of very strong laser field $B$, the magnon creation and annihilation operators correspond to $\tilde{S}^{\pm}=S^{y}\pm iS^{z}$. Using these operators, the Hamiltonian (Eq. (1) with Eq. (17)) becomes

[TABLE]

The $\tilde{S}_{j}^{+}\tilde{S}_{j+1}^{+}+\tilde{S}_{j}^{-}\tilde{S}_{j+1}^{-}$ term creates and annihilates magnons (at large $B$) in pairs. In the Hamiltonian Eq. (20), the Hilbert space is separated into the sectors with $S_{\mathrm{tot}}^{x}=(\mathrm{even\;integer})$ and $S_{\mathrm{tot}}^{x}=(\mathrm{odd\;integer})$ since the parity of the magnon number is a conserved quantity. Hence the state remains in the same sector during the time evolution. The initial state is the ferromagnetic state, which corresponds to a Schrödinger cat state in the $S^{x}$ basis,

[TABLE]

Thus, this state has weight in both $S_{\mathrm{tot}}^{x}=(\mathrm{even\;integer})$ and $S_{\mathrm{tot}}^{x}=(\mathrm{odd\;integer})$ sectors. $M_{x}$ has nonzero expectation values for states within the same sector, $\langle\Psi_{\mathrm{even}}|S_{\mathrm{tot}}^{x}|\Psi_{\mathrm{even}}\rangle+\langle\Psi_{\mathrm{odd}}|S_{\mathrm{tot}}^{x}|\Psi_{\mathrm{odd}}\rangle$, while $M_{z}$ has nonzero expectation values for the states between different sectors, $\langle\Psi_{\mathrm{odd}}|S_{\mathrm{tot}}^{z}|\Psi_{\mathrm{even}}\rangle+\langle\Psi_{\mathrm{even}}|S_{\mathrm{tot}}^{z}|\Psi_{\mathrm{odd}}\rangle$. This expression explains why the two-magnon mode is evident in $m^{x}(\omega)$ [Figs. 10(d) and 10(e)] while it is not apparent in $m^{z}(\omega)$ [Figs. 10(i) and 10(j)].

We also analyze the dynamics of this system by means of the tdMF theory. The mean field Hamiltonian is

[TABLE]

The radiation power spectrum calculated by the Hamiltonian Eq. (21) is shown in Fig. 10. The tdMF result shows a peak or plateau structure in the HHG spectrum, but quantitatively it deviates strongly from the iTEBD and ED results, in contrast to the case of the Ising model. This is due to the strong quantum fluctuations induced by the $S_{j}^{+}S_{j+1}^{-}+S_{j}^{-}S_{j+1}^{+}$ term, and implies that the tdMF theory does not provide a good description of the XXZ model.

V Summary and discussions

In this paper, we studied HHG in quantum spin systems driven by a laser magnetic field. When the laser is applied to magnetic insulators, it drives the magnetic dipole which generates electromagnetic radiation with power proportional to $|\omega^{2}M(\omega)|^{2}$. We considered two specific but fundamental quantum spin chain models, the Ising model with static longitudinal field and the XXZ model. In both cases, when the magnetic field is strong enough, the spin HHG shows a (multiple-)plateau structure, which is associated with the annihilation of (multiple) magnons.

To be more specific, in the Ising model case, the excitation gap does not close in the presence of a transverse field since the $\mathbf{Z}_{2}$ symmetry is explicitly broken. When the laser amplitude is weak enough, the time-dependent perturbation theory is valid, which explains the appearance of a peak around the frequency $J+H$. With increasing laser amplitude, the shape of the HHG spectrum changes from a peak structure to a plateau structure. The subcycle analysis suggests that the HHG originates from the annihilation of magnons. The cutoff energies, above which the radiation intensity drops, correspond to integer multiples of the single-magnon excitation energy at $q=0$. Since the magnetic field is stronger than the interaction, the tdMF theory provides a quantitative description.

In the XXZ model without longitudinal field, the system has a $\mathbf{Z}_{2}$ symmetry and a phase transition happens at a critical value of the transverse field. The structure of the HHG spectrum changes depending on whether the peak amplitude of the laser magnetic field is below or above the critical field. Similarly to the Ising case, when the laser amplitude is small, the time-dependent perturbation theory is valid and explains the appearance of a peak around the frequency $J_{xy}+J_{z}$. As the laser amplitude increases, the peak structure transforms into a plateau structure. The cutoff energy of this plateau corresponds to the single-magnon mass at $q=0$ below the critical field. When the laser amplitude is larger than the critical field, the threshold is determined by the two-magnon excitation at $q=\pi$. The subcycle analysis and the ED analysis suggest that also in the XXZ model case, the annihilation of magnons leads to the HHG signal. The tdMF approach is not effective in this model due to the quantum fluctuation caused by the $J_{xy}$ term.

Now let us discuss the similarities and differences between the HHG from spin systems and that from insulating electron systems such as semiconductors and Mott insulators Vampa et al. (2014, 2015b); Murakami et al. (2018). In the latter case, a periodic electric field creates charge carriers (electrons and holes in semiconductors, and doublons and holons in Mott insulators) and these carriers move around in response to the applied electric field. The HHG originates from the dynamics of these charge carriers, which can be separated into the interband and intraband current. The interband current corresponds to the creation and recombination of charge carriers, while the intraband current represents the contribution from hopping processes which do not change the number of charge carriers, i.e. where the carriers remain in the same conduction/valence or Hubbard band. In contrast, in the spin systems, the magnetic field can excite magnetic excitations (magnons) but there is no preferable direction to move since the homogeneous magnetic field, unlike the electric field, does not produce a spatially dependent potential. Hence, the HHG signal originating from the dynamics of the magnetization is analogous to the interband current, while there is no counterpart to the intraband current. Our finding that the spin HHG is associated with the annihilation of magnons is reminiscent of the electron HHG which is dominated by the recombination of charge carriers Vampa et al. (2014); Murakami et al. (2018).

Experimentally, the HHG from spins excited by time-periodic magnetic fields can be realized by choosing large gap insulating materials, and by taking advantage of metamaterials to selectively enhance the magnetic field Mukai et al. (2014). For example, $\mathrm{CoNb_{2}O_{6}}$ Coldea et al. (2010) can be represented as a ferromagnetic XXZ chain with $J_{xy}\ll J_{z}$, and therefore the discussion in Sec. IV is relevant for this material, while examples of Ising magnets (Sec. III) such as $\mathrm{Dy(C_{2}H_{5}SO_{4})_{3}\cdot 9H_{2}O}$, $\mathrm{LiTbF_{4}}$ and $\mathrm{LiHoF_{4}}$ are discussed in Ref. Wolf (2000). For $\mathrm{CoNb_{2}O_{6}}$, since the value of $J_{z}(=10)$ is 1.94 meV Coldea et al. (2010), the energy unit is $0.194\;\mathrm{meV}=1.67\;\mathrm{T}=2\pi\times 0.0469\;\mathrm{THz}$ by noting that $g\mu_{\mathrm{B}}B$ and $\hbar\Omega$ have the dimension of energy, where $g\simeq 2$ is the Landé $g$ factor for electron spins, $\mu_{\mathrm{B}}=0.0579\;\mathrm{meV/T}$ is the Bohr magneton, and $\hbar=6.58\times 10^{-13}\;\mathrm{meV}\cdot\mathrm{s}$. $\Omega=1$ and $B=2$ thus correspond to $\Omega=2\pi\times 0.0469\;\mathrm{THz}$ and $B=3.34\;\mathrm{T}$, respectively.

Our results demonstrate the possibility of generating high-harmonic signals in spin systems, which may be utilized for new laser sources in the THz regime or to obtain information about the magnetic excitations of these spin systems under strong fields. In the present work, we focused on one-dimensional ferromagnets but the fact that the tdMF results show a similar HHG spectrum strongly suggests that the HHG signal can also be produced in higher dimensional magnets. Radiation from the magnetic dipole should be possible also in ferrimagnets and antiferromagnets. Although the total magnetization is zero in antiferromagnets, the laser magnetic field produces a net magnetization and a HHG signal can be expected. Since there exist various kinds of quantum spin systems, studying these other types of magnetic insulators is an interesting direction for future research.

Acknowledgements.

We acknowledge fruitful discussions with Takashi Oka, Tran Trung Luu and Gabriel Aeppli. ST is supported by the Swiss National Science Foundation under Division II and ImPact project (No. 2015-PM12-05-01) from the Japan Science and Technology Agency. YM and PW are supported by ERC Consolidator Grant No. 724103.

Appendix A Time-dependent perturbation theory

In this appendix we analyze the spin system in the presence of a laser field using the time-dependent perturbation theory. The Hamiltonian is

[TABLE]

where $V(t)$ represents the laser-matter interaction which is assumed here for simplicity to have the form

[TABLE]

We switch to the interaction picture. The state and operator are represented as $|\Psi(t)\rangle_{\mathrm{I}}=e^{+i\mathcal{H}_{\mathrm{spin}}t}|\Psi(t)\rangle$ and $\mathcal{O}_{\mathrm{I}}=e^{i\mathcal{H}_{\mathrm{spin}}t}\mathcal{O}e^{-i\mathcal{H}_{\mathrm{spin}}t}$, respectively, where $|\Psi(t)\rangle$ and $\mathcal{O}$ are the state and operator in the Schrödinger picture. The equation of motion becomes

[TABLE]

where $V_{\mathrm{I}}(t)=e^{i\mathcal{H}_{\mathrm{spin}}t}V(t)e^{-i\mathcal{H}_{\mathrm{spin}}t}$. From Eq. (25), we derive

[TABLE]

We denote the eigenenergy and eigenstate of $\mathcal{H}_{\mathrm{spin}}$ by $E_{n}$ and $|\varphi_{n}\rangle$, respectively. Let us expand $|\Psi(t)\rangle_{\mathrm{I}}$ in the basis of $|\varphi_{n}\rangle$,

[TABLE]

We substitute (27) into (26) and take the inner product with $\langle\varphi_{n}|$, to obtain

[TABLE]

where $V_{nm}(t)=\langle\varphi_{n}|V_{\mathrm{I}}(t)|\varphi_{m}\rangle=e^{-i(E_{m}-E_{n})t}\langle\varphi_{n}|V(t)|\varphi_{m}\rangle.$ In the present case,

[TABLE]

for $t\geq 0$. At $t=0$, the system is in the ground state $c_{0}(0)=1$ and $c_{n}(0)=0$ ($n\geq 1$), thus Eq. (28) becomes

[TABLE]

Physical observables are calculated as

[TABLE]

and specifically for the magnetization as

[TABLE]

First we consider the Ising model $\mathcal{H}_{\mathrm{spin}}=\mathcal{H}_{\mathrm{Ising}}$ [Eq. (5) in the main text]. The ground state is the configuration with all spins up $|\varphi_{0}\rangle=|\!\!\!\uparrow\uparrow\cdots\uparrow\rangle$ and the first excited states $|\varphi_{n}\rangle$ ($n=1,\ldots,N$) are single spin flipped states $|\varphi_{n}\rangle=S_{n}^{-}|\varphi_{0}\rangle$. Since the excitation gap is $E_{n}-E_{0}=H+J$ ($n=1,\ldots,N$), we can calculate

[TABLE]

Hence the magnetization (29) becomes

[TABLE]

For $M_{x}$, the order $B$ term contains components with frequency $\Omega$ and $H+J$, while for $M_{z}$, the order $B^{2}$ term contains components with frequency $2\Omega$ and $H+J\pm\Omega$. The full calculation of the $\mathcal{O}(B^{3})$ terms is difficult, but we can see that the $c_{n}^{*}(t)c_{0}(t)e^{i(H+J)t}$ term contains $e^{3i\Omega t}$ and $e^{i(H+J\pm 2\Omega)t}$. Thus, we can surmise that for $M_{x}$, the leading order of frequency $n\Omega$ is $B^{n}$ ($n$: odd) and that of frequency $H+J\pm n\Omega$ is $B^{n+1}$ ($n$: even) while for $M_{z}$, the leading order of frequency $n\Omega$ is $B^{n}$ ($n$: even) and that of frequency $H+J\pm n\Omega$ is $B^{n+1}$ ($n$: odd).

Next we consider the XXZ model $\mathcal{H}_{\mathrm{spin}}=\mathcal{H}_{\mathrm{XXZ}}$ [Eq. (17) in the main text], where the laser field is again assumed to be Eq. (24). The ground state of $\mathcal{H}_{\mathrm{XXZ}}$ is the fully polarized ferromagnetic state $|\varphi_{0}\rangle=|\!\uparrow\uparrow\cdots\uparrow\rangle$ and its eigenenergy is $E_{0}=-\frac{NJ_{z}}{4}$. Due to the symmetry breaking, $|\!\downarrow\downarrow\cdots\downarrow\rangle$ is also a ground state, but we assume that the initial state is $|\varphi_{0}\rangle$. The low-energy excited states are single-magnon states $|\varphi_{n}\rangle=\frac{1}{\sqrt{N}}\sum_{j=1}^{N}e^{i\frac{2\pi n}{N}j}S_{j}^{-}|\varphi_{0}\rangle$ ($n=1,\ldots,N$) and their eigenenergy is $E_{n}=-\frac{(N-4)J_{z}}{4}+J_{xy}\cos(\frac{2\pi n}{N})$ ($n=1,\ldots,N$). Noting that

[TABLE]

we obtain

[TABLE]

for $n=1,\ldots,N$, where $\delta_{nN}$ is the Kronecker delta. Thus we derive

[TABLE]

Therefore the magnetization (29) becomes

[TABLE]

Similarly to the case of the Ising model, we can surmise that for $M_{x}$, the leading order of frequency $n\Omega$ is $B^{n}$ ($n$: odd) and that of frequency $J_{xy}+J_{z}\pm n\Omega$ is $B^{n+1}$ ($n$: even) while for $M_{z}$, the leading order of frequency $n\Omega$ is $B^{n}$ ($n$: even) and that of frequency $J_{xy}+J_{z}\pm n\Omega$ is $B^{n+1}$ ($n$: odd).

Appendix B Spin wave theory

We consider the system

[TABLE]

where $J_{z}>J_{xy}>0$ with general spin-$S$. The number of sites is $N$ and we consider periodic boundary conditions.

First, let us determine the classical ground state. For $B=0$ (large $B$), the spin is polarized along the $S^{z}$ ($S^{x}$) axis, thus we can assume that the direction of the spins is in the $xz$ plane, $S_{j}=S(\sin\phi,0,\cos\phi)$. The energy is

[TABLE]

where the constant term is neglected. The configuration minimizing $E$ is

[TABLE]

We introduce new spin axes $\tilde{S}_{j}^{x,y,z}$ as $S_{j}^{x}=\cos\phi\tilde{S}_{j}^{x}+\sin\phi\tilde{S}_{j}^{z}$, $S_{j}^{y}=\tilde{S}_{j}^{y}$ and $S_{j}^{z}=-\sin\phi\tilde{S}_{j}^{x}+\cos\phi\tilde{S}_{j}^{z}$, so that the spin is polarized along the $\tilde{S}_{j}^{z}$ axis. Next we perform the Holstein-Primakoff transformation,

[TABLE]

where $a_{j}$ and $a_{j}^{\dagger}$ are annihilation and creation operators for bosons (magnons), and $n_{j}\equiv a_{j}^{\dagger}a_{j}$ is the number operator. Expanding in powers of $1/S$ and retaining terms up to second order in $a_{j}$ and $a_{j}^{\dagger}$ yields an expression of the Hamiltonian in terms of magnon operators,

[TABLE]

The first order term of $a_{j}$, $a_{j}^{\dagger}$ vanishes if one imposes the condition (33). After the Fourier transform $a_{k}=\frac{1}{\sqrt{N}}\sum_{j}e^{-ikj}a_{j}$, $a_{k}^{\dagger}=\frac{1}{\sqrt{N}}\sum_{j}e^{ikj}a_{j}^{\dagger}$, we obtain

[TABLE]

where

[TABLE]

and $E_{\mathrm{CL}}=NS^{2}(J_{xy}\sin^{2}\phi-J_{z}\cos^{2}\phi)-NSB\sin\phi$ is the classical ground state energy (see Eq. (32)). Note that $g(k,\phi)=g(-k,\phi)$. We then perform the Bogoliubov transformation $b_{k}=a_{k}\cosh\theta_{k}+a_{-k}^{\dagger}\sinh\theta_{k}$, $a_{k}=b_{k}\cosh\theta_{k}-b_{-k}^{\dagger}\sinh\theta_{k}$, where $\tanh 2\theta_{k}=\frac{f(k,\phi)}{g(k,\phi)}$ ($\theta_{k}=\theta_{-k}$). Finally the Hamiltonian becomes

[TABLE]

where

[TABLE]

is the quantum correction to the classical ground state energy which is a constant.

The magnon band structure from the spin wave theory $-f(k,\phi)\sinh 2\theta_{k}+g(k,\phi)\cosh 2\theta_{k}$ is shown in Fig. 13. The excitation gap closes and the transition happens at $B=J_{xy}+J_{z}$.

Appendix C Additional analysis of models with different parameters

C.1 Ising model with weak static field

In the main text, we considered the Ising model with a strong longitudinal field $H\gg J$. In this subsection, we study how the radiation spectrum and the excitation structure are changed if the static field is weak $H\ll J$. The parameters are set to $J=4$, $H=2$, and $\Omega=1$.

The DSFs $|\mathrm{Im}\chi^{xx}(q,\omega)|$ and $|\mathrm{Im}\chi^{zz}(q,\omega)|$ for the ground states of Eq. (10) are shown in Fig. 14. The low energy excitation is again a magnon and the shape of the dispersion is similar to the high field case, but the size of the excitation gap decreases at first with the introduction of $B$ and then increases. This behavior is caused by the weak $\mathbf{Z}_{2}$ symmetry breaking due to the small longitudinal field $H$. For $H=0$, the $\mathbf{Z}_{2}$ symmetry is recovered and a gap closing (i.e., a quantum phase transition) happens at the critical field $B=B_{\mathrm{c}}$. For $H>0$, the $\mathbf{Z}_{2}$ symmetry is explicitly broken and the excitation gap opens at $B_{\mathrm{c}}$. However the gap size is small for $H\ll J$, and thus the gap size becomes a nonmonotonous function of $B$. The continuous spectrum corresponding to the two-magnon mode [Eq. (11) in the main text] appears more evidently in Figs. 14(b)-14(d) compared with the strong field case. We can see an additional excitation between the single-magnon band and the two-magnon continuum, which is a two-magnon bound state. When $B$ is small, the energy of this state (two-spin flips on nearest neighbor sites) is $\simeq J+2H(=8)$. With increasing $B$, this bound state is strongly hybridized with the two-magnon continuum and is finally merged into it.

In Fig. 15, we show the time evolution of $m^{x}$ and $m^{z}$. For $B=6$, the shape of $m^{x}(t)$ is clearly different from the sinusoidal curve of the laser field [Eq. (8) in the main text] especially near the peaks, which gives rise to a strong HHG signal. In contrast to the case of strong static fields, the final value of $m^{z}$ deviates from the value $1/2$ for large $B$. (This deviation will also happen in the strong longitudinal field case for large $B/H$.) As discussed above, the gap of the system first decreases and then increases as a function of $B$, while the Landau-Zener tunneling happens mainly near the minimum of the gap. Thus, transitions to excited states of the snapshot Hamiltonian occur, which results in the drop of the final value of $m^{z}$ from 1/2.

In order to investigate the HHG, we show $|\omega^{2}m^{x}(\omega)|^{2}$ and $|\omega^{2}m^{z}(\omega)|^{2}$ in Fig. 16. When $B$ is small, the behavior of the radiation spectrum is similar to the case of high static field. The intensity of $|\omega^{2}m^{x}(\omega)|^{2}$ and $|\omega^{2}m^{z}(\omega)|^{2}$ generically peaks at $\omega=(2n+1)\Omega$ and $\omega=2n\Omega$ ($n$: integer), respectively, but at $\omega=J+H=6$, $|\omega^{2}m^{x}(\omega)|^{2}$ exhibits a local maximum and $|\omega^{2}m^{z}(\omega)|^{2}$ shows a dip. This is consistent with the time-dependent perturbation theory. When $B$ becomes larger, a plateau structure appears in the HHG signal and its threshold corresponds to the single-magnon excitation energy (the dashed lines in Fig. 16) or the energy of the two-magnon bound state (the dotted lines in Fig. 16). As $B$ is further increased, the threshold of the HHG plateau changes from the single-magnon energy to twice of the magnon energy at $q=\pi$ (dashed-dotted lines in Fig. 16). This behavior is similar to the XXZ model with small $J_{xy}$ (see Sec. IV), where the threshold corresponds to the single-magnon energy before the transition and the two-magnon energy after the transition. In the present case, due to the existence of the longitudinal field $H$, the change of the threshold energy scale is not a transition but a crossover.

We also show the analysis by the tdMF theory with the Hamiltonian Eq. (15) in Fig. 16. For the weak laser amplitude $B$, the agreement between the iTEBD and tdMF theories is quantitatively good. When $B$ becomes large, the spectra start to deviate above the single-magnon energy but the threshold of the HHG plateau is almost the same ($B=2$). For $B=6$ [Fig. 16(d) and 16(h)], the HHG signal calculated by the tdMF theory becomes less prominent above the single-magnon energy ($\omega\simeq 5$) while the threshold is the two-magnon energy for iTEBD. This result implies that the tdMF theory can reproduce the single-magnon dynamics but fails to capture multiple-magnon processes.

In Fig. 17, we show the subcycle radiation spectrum $\log_{10}|\omega^{2}m_{\mathrm{sub}}^{x}(\omega;t_{*})|^{2}$ for $B=2$ and $B=6$. Green and purple solid lines show the energy of the single-magnon at $q=0$ and of two magnons at $q=\pi$ for the snapshot Hamiltonian, respectively. In Fig. 17(a), some intensity exists between the two lines, which may be associated with the two-magnon bound state seen in Fig. 14(a) and 14(b) (around $\omega=8$ – 9). For large laser field amplitude ($B=6$), the intensity around the two-magnon energy becomes prominent. This subcycle analysis supports the interpretation that the crossover of the HHG signal is caused by a change of the dynamics from a single-magnon to a two-magnon dominated process.

C.2 XXZ model with strong $J_{xy}$

We next consider the XXZ model with $J_{xy}$ stronger than that in the main text, i.e. $J_{xy}\lesssim J_{z}$. The parameters are set to $J_{xy}=8$, $J_{z}=10$, and $\Omega=1$. In Fig. 18, we show the DSF of the Hamiltonian Eq. (18) (in the main text). The single-magnon dispersion ($E(q)=J_{xy}\cos(q)+J_{z}$ for $B=0$) splits by the introduction of the $B$ field. The lower bound of the spectrum at $q=0$ decreases with increasing $B$, and the gap closes at $B_{\mathrm{c}}\simeq 6$, where a phase transition happens. After the transition, the intensity of $\chi^{zz}$ is stronger than $\chi^{xx}$, but both are still comparable. The dispersion captured by $\chi^{zz}$ has a dip around $q=\pi$, a property which is reproduced by the spin wave theory (see Appendix B). However, the relation $\tilde{\Delta}_{q=0}=2\tilde{\Delta}_{q=\pi}$ (for $B>B_{\mathrm{c}}$) does not hold in contrast to the weak $J_{xy}$ case.

The time evolution of $m^{x}$ and $m^{z}$ calculated by iTEBD is shown in Fig. 19. The behavior of $m^{x}(t)$ and $m^{z}(t)$ is similar to that in the weak $J_{xy}$ case. The time evolution of $m^{z}$ is different depending on whether $B$ is smaller or larger than $B_{\mathrm{c}}$. In particular, $m^{z}$ suddenly decreases from $1/2$, when $B$ exceeds $B_{\mathrm{c}}$.

In Fig. 20, we show the HHG spectra $|\omega^{2}m^{x}(\omega)|^{2}$ and $|\omega^{2}m^{z}(\omega)|^{2}$. When $B$ is small, the result is again described by the time-dependent perturbation theory (Appendix A), and there is a peak at $\omega=J_{xy}+J_{z}=18$ for $m^{x}(\omega)$ and $\omega=J_{xy}+J_{z}\pm\Omega$ for $m^{z}(\omega)$ in the case of $B=2$. In contrast to the weak $J_{xy}$ case, the threshold of the plateau corresponds to $m_{q=0}$ for both $B<B_{\mathrm{c}}$ and $B>B_{\mathrm{c}}$. Since there is a dip around $q=\pi$ for $B>B_{\mathrm{c}}$ as is seen from the dispersions in Figs. 18(i) and 18(j), the relation $\Delta_{q=0}=2\tilde{\Delta}_{q=\pi}$ does not hold. Hence the energy scale of the threshold of the plateau corresponds to the mode excited by the operator $S_{\mathrm{tot}}^{x}$.

Figure 21 shows the subcycle radiation spectrum Eq. (12) (in the main text) for $B=4$ and $B=10$ (below and above the critical field, respectively). In the case of $B=4$ [Fig. 21(a)], the strong intensity in the HHG signal follows the single-magnon excitation energy $\Delta_{q=0}$ of the snapshot Hamiltonian, which indicates that the threshold is related to the annihilation of single magnons at $q=0$ for $B<B_{\mathrm{c}}$. In the case of $B=10$ [Fig. 21(b)], the strong intensity in the HHG signal still roughly follows the energy of $\Delta_{q=0}$ and $2\Delta_{q=0}$.

To obtain more information on the relation between the HHG spectra and the excitation structure, we perform ED calculations for a system with $N=8$ sites. Although there is a quantitative deviation from the iTEBD results, the ED calculations qualitatively reproduce the behavior of the HHG spectra, especially the peaks and plateaus as shown in Fig. 20. In Fig. 22(a), we show $\langle\Phi_{n}|S_{\mathrm{tot}}^{x}|\Phi_{n}\rangle$ calculated for the eigenstates of the snapshot Hamiltonian at $t=t_{\mathrm{peak}}$. Although the discretization of $\langle\Phi_{n}|S_{\mathrm{tot}}^{x}|\Phi_{n}\rangle$ is not as clear as in the weak $J_{xy}$ case and the values are not necessarily close to integers, the eigenstates can be roughly classified into sectors. In Fig. 22(b), we plot the quantity $A_{n}=|\alpha_{n}^{*}\alpha_{0}\langle\Phi_{n}|S_{\mathrm{tot}}^{x}|\Phi_{0}\rangle|$ ($\alpha_{n}\equiv\langle\Phi_{n}|\Psi(t_{\mathrm{peak}})\rangle$) [Eq. (19) in the main text]. We see that there is a strong intensity at the energy $E_{n}-E_{0}=20.9$, which agrees with the threshold energy in Figs. 20(e) and 20(j). The eigenstate at $E_{n}-E_{0}=20.9$ belongs to the $\langle\Phi_{n}|S_{\mathrm{tot}}^{x}|\Phi_{n}\rangle\simeq-1.5$ sector (depicted by the cross marks in Fig. 22), and this sector is connected to the $\langle\Phi_{n}|S_{\mathrm{tot}}^{x}|\Phi_{n}\rangle\simeq-2$ sector in the weak $J_{xy}$ case. As is seen from Eq. (20) in the main text, the hybridization between two sectors characterized by different eigenvalues of $S_{\mathrm{tot}}^{x}$ is caused by the term $\frac{J_{xy}+J_{z}}{4}\sum_{j}(\tilde{S}_{j}^{+}\tilde{S}_{j+1}^{+}+\tilde{S}_{j}^{-}\tilde{S}_{j+1}^{-})$, which becomes stronger as $J_{xy}$ is increased. This strong hybridization explains the results that the values of $\langle\Phi_{n}|S_{\mathrm{tot}}^{x}|\Phi_{n}\rangle$ deviate from integer and that the state with the energy $E_{n}-E_{0}=20.9$ is strongly excited by the $S_{\mathrm{tot}}^{x}$ operator. By recalling that $\Delta_{q=0}$ is not equal to $2\tilde{\Delta}_{q=\pi}$, this excitation of $E_{n}-E_{0}=20.9$ cannot be regarded as two free magnons created by the $S^{z}$ operator, which implies that magnon-magnon interaction effects are important.

The radiation power spectrum calculated by the tdMF Hamiltonian Eq. (21) (in the main text) is also shown in Fig. 20. The plateau structure of the radiation spectrum does not appear in the tdMF analysis and the high harmonic signals decay exponentially as the frequency becomes larger. Due to the strong quantum fluctuations induced by $J_{xy}=8$, the tdMF theory does not give a good description in this case.

Appendix D Subtraction of linear response

In this section, we show the time evolution of $m^{x}$ after the subtraction of the linear response component in order to illustrate the origin of high harmonic generation. First we calculate the magnetization dynamics $m_{0}(t)\equiv\langle S^{x}(t)\rangle$ under the weak laser field (8) with $B=0.2$ (i.e., in the linear response regime). Then we subtract this linear response component from $m^{x}(t)$. In Fig. 23, $m^{x}(t)-(B/B_{0})m_{0}(t)$ is shown for both the Ising model with $J=2$ and $H=6$ and the XXZ model with $J_{xy}=2$ and $J_{z}=10$. In both models, the discrepancy from the linear response becomes large near the maxima of the laser field amplitude, and the nonlinear component increases with increasing laser intensity. In particular, $m^{x}(t)-(B/B_{0})m_{0}(t)$ has a non-sinusoidal shape, which is a manifestation of strong nonlinearity.

Appendix E Wavelet analysis

In this section, we show the time-resolved radiation spectra obtained by a wavelet analysis. This approach is similar to the subcycle analysis, but in contrast to the latter, the time and energy resolution depends on $\omega$. In the low energy regime, the time ($t_{*}$) resolution is low and the energy ($\omega$) resolution is high, while it is the opposite in the high energy regime. The wavelet transform for the second derivative of the magnetic moment is defined as

[TABLE]

where

[TABLE]

is a mother function for the Gabor wavelet. In Fig. 24, we show the wavelet spectrum $|{m^{\prime\prime}}_{\mathrm{W}}^{x}(\omega;t_{*})|^{2}$ for the Ising model with $J=4$, $H=2$, and $B=8$ and the XXZ model with $J_{xy}=2$, $J_{z}=10$, and $B=10$. In the low-energy region, the signal is smeared out in the time direction while there are clearly resolved peaks along the $\omega$ direction. In the high energy region, on the other hand, the structures are smeared out along the $\omega$ axis, while the time evolution of the spectral features can be well captured.