Torque equilibrium spin wave theory study of anisotropy and Dzyaloshinskii-Moriya interaction effects on the indirect K$-$ edge RIXS spectrum of a triangular lattice antiferromagnet

TL;DR

This paper applies torque equilibrium spin wave theory to analyze how anisotropy and Dzyaloshinskii-Moriya interactions influence the RIXS spectra of a triangular lattice antiferromagnet, providing theoretical predictions for experimental testing.

Contribution

The study extends TESWT to include DM interactions and zero-point fluctuations, enabling detailed analysis of RIXS spectra in noncollinear quantum magnets with experimental relevance.

Findings

Anisotropy causes a downshift of the roton RIXS peak at the M point.

DM interaction stabilizes the spiral phase and causes an upshift of RIXS peaks.

Distinct spectral behaviors at M and M' points reveal effects of anisotropy and DM interaction.

Abstract

We apply the recently formulated torque equilibrium spin wave theory (TESWT) to compute the $1/S$-order interacting $K$ -edge bimagnon resonant inelastic x-ray scattering (RIXS) spectra of an anisotropic triangular lattice antiferromagnet with Dzyaloshinskii-Moriya (DM) interaction. We extend the interacting torque equilibrium formalism, incorporating the effects of DM interaction, to appropriately account for the zero-point quantum fluctuation that manifests as the emergence of spin Casimir effect in a noncollinear spin spiral state. Using inelastic neutron scattering data from Cs$_2$CuCl$_4$ we fit the 1/S corrected TESWT dispersion to extract exchange and DM interaction parameters. We use these new fit coefficients alongside other relevant model parameters to investigate, compare, and contrast the effects of spatial anisotropy and DM interaction on the RIXS spectra at various points across the magnetic Brillouin zone. We highlight the key features of the bi- and trimagnon RIXS spectrum at the two inequivalent rotonlike points, $M(0,2 \pi/\sqrt{3})$ and $M^{\prime}(\pi,\pi/\sqrt{3})$, whose behavior is quite different from an isotropic triangular lattice system. While the roton RIXS spectrum at the $M$ point undergoes a spectral downshift with increasing anisotropy, the peak at the $M^\prime$ location loses its spectral strength without any shift. With the inclusion of DM interaction the spiral phase is more stable and the peak at both $M$ and $M^\prime$ point exhibits a spectral upshift. Our calculation offers a practical example of how to calculate interacting RIXS spectra in a non-collinear quantum magnet using TESWT. Our findings provide an opportunity to experimentally test the predictions of interacting TESWT formalism using RIXS, a spectroscopic method currently in vogue.

Tables

Table 1. Table 1: Parameter values of Cs 2 CuCl 4 using different methods. The first line is our TESWT fitting results. The second line is our 1/S-SWT fitting parameters. The third line gives the fitting parameters of series expansion (SE) method  Fjærestad  et al. ( 2007 ) . The last line gives the parameters measured by Electron-Spin-Resonance (ESR)  Zvyagin  et al. ( 2014 ) .

Method	$J$(meV)	$J^{\prime }$(meV)	$D$(meV)
TESWT	$0.480\pm 0.009$	$0.152\pm 0.015$	$0.012\pm 0.002$
$1/S$-LSWT	$0.573\pm 0.009$	$0.239\pm 0.014$	$0.012\pm 0.001$
SE	$0.374\pm 0.005$	$0.128\pm 0.005$	$0.020\pm 0.002$
ESR	$0.41\pm 0.02$	$0.122\pm 0.006$	$-$

Equations

$$\mathcal{H}=\sum _{⟨ij⟩}{J}_{ij}{\mathbf{S}}_i\cdot {\mathbf{S}}_j+H_{\mathrm{DM}},$$

$$H_{\mathrm{DM}}=-\sum _i\mathbf{D}\cdot [{\mathbf{S}}_i\times ({\mathbf{S}}_{i+{\mathbf{\delta }}_1}+{\mathbf{S}}_{i+{\mathbf{\delta }}_2})],$$

$${\mathbf{S}}_i/S=\cos (\mathbf{Q}\cdot {\mathbf{r}}_i){\hat{z}}_0+\sin (\mathbf{Q}\cdot {\mathbf{r}}_i){\hat{x}}_0,$$

$$E_0(\mathbf{Q})=3NJ{S}^2({\lambda }_{\mathbf{Q}}-{\eta }_{\mathbf{Q}})=3NJ{S}^2{\gamma }_{\mathbf{Q}},$$

$${\lambda }_{\mathbf{k}}=\frac{1}{3}(\cos k_x+2\alpha \cos \frac{k_x}{2}\cos \frac{\sqrt{3}}{2}{k}_y),$$

$${\eta }_{\mathbf{k}}=\frac{2}{3}\eta \sin \frac{k_x}{2}\cos \frac{\sqrt{3}}{2}{k}_y,$$

$${\nabla }_{\mathbf{Q}}{E}_0(\mathbf{Q})=0,$$

$$\left\{\begin{array}[]{cc}\sin Q_{x}+\alpha\sin\frac{Q_{x}}{2}\cos\frac{\sqrt{3}}{2}Q_{y}+\eta\cos\frac{Q_{x}}{2}\cos\frac{\sqrt{3}}{2}Q_{y}=0,\\ \alpha\cos\frac{Q_{x}}{2}\sin\frac{\sqrt{3}}{2}Q_{y}-\eta\sin\frac{Q_{x}}{2}\sin\frac{\sqrt{3}}{2}Q_{y}=0.\end{array}\right.$$

$$Q_{0}=\left\{\begin{array}[]{cc}2\arccos(-\frac{\alpha}{2}),&\alpha<2,\\ 2\pi,&\alpha\geq 2.\end{array}\right.$$

$$S_i^{x_0}=S_i^z\sin {\theta }_i+S_i^x\cos {\theta }_i,$$

$$S_i^{z_0}=S_i^z\cos {\theta }_i-S_i^x\sin {\theta }_i,$$

$$\displaystyle\mathcal{H}=\sum_{\langle ij\rangle}\Big{[}J_{ij}S_{i}^{y}S_{j}^{y}+J_{ij}^{+}(S_{i}^{z}S_{j}^{z}+S_{i}^{x}S_{j}^{x})$$

$$\displaystyle+J_{ij}^{-}(S_{i}^{z}S_{j}^{x}-S_{i}^{x}S_{j}^{z})\Big{]},$$

$$J_{ij}^{+}=$$

$$J_{ij}^{-}=$$

$$S_i^z=S-n_i,S_i^{-}=a^{†}\sqrt{2S-n_i},S_i^{+}=(S_i^{-}{)}^{†},$$

$$\mathcal{H}=E_0(\mathbf{Q})+H_2+H_3+H_4,$$

$$H_{2}=\sum_{\mathbf{k}}\Big{[}A_{\mathbf{k}}a_{\mathbf{k}}^{\dagger}a_{\mathbf{k}}+\frac{B_{\mathbf{k}}}{2}(a_{\mathbf{k}}^{\dagger}a_{-\mathbf{k}}^{\dagger}+a_{-\mathbf{k}}a_{\mathbf{k}})\Big{]},$$

$$A_{\mathbf{k}}=3JS[{\lambda }_{\mathbf{k}}+{\xi }_{\mathbf{k}}-2{\gamma }_{\mathbf{Q}}],$$

$$B_{\mathbf{k}}=3JS[{\xi }_{\mathbf{k}}-{\lambda }_{\mathbf{k}}],$$

$${\xi }_{\mathbf{k}}=$$

$$a_{\mathbf{k}}=u_{\mathbf{k}}{b}_{\mathbf{k}}+v_{\mathbf{k}}{b}_{-\mathbf{k}}^{†},$$

$$u_{\mathbf{k}}=\sqrt{\frac{A_{\mathbf{k}}+{\epsilon }_{\mathbf{k}}}{2{\epsilon }_{\mathbf{k}}}},v_{\mathbf{k}}=-\frac{B_{\mathbf{k}}}{|B_{\mathbf{k}}|}\sqrt{\frac{A_{\mathbf{k}}-{\epsilon }_{\mathbf{k}}}{2{\epsilon }_{\mathbf{k}}}}.$$

$${\epsilon }_{\mathbf{k}}=\sqrt{A_{\mathbf{k}}^2-B_{\mathbf{k}}^2}.$$

$$H_2=E_2(\mathbf{Q})+\sum _{\mathbf{k}}{\epsilon }_{\mathbf{k}}{b}_{\mathbf{k}}^{†}{b}_{\mathbf{k}},$$

$$E_2(\mathbf{Q})=3NSJ{\gamma }_{\mathbf{Q}}+\frac{1}{2}\sum _{\mathbf{k}}{\epsilon }_{\mathbf{k}},$$

$${\nabla }_{\mathbf{Q}}[E_0(\mathbf{Q})+E_2(\mathbf{Q})]=0.$$

$${\nabla }_{\mathbf{Q}}[E_0({\mathbf{Q}}_1)+E_2({\mathbf{Q}}_1)]={\nabla }_{\mathbf{Q}}{E}_2({\mathbf{Q}}_0)+\Delta \mathbf{Q}\cdot K=0,$$

$$K_{\alpha ,\beta }=\frac{{\partial }^2{E}_0({\mathbf{Q}}_0)}{\partial Q_{\beta }{Q}_{\alpha }}.$$

$$\Delta \mathbf{Q}=-\mathbf{w}\cdot K^{-1},$$

$$w_{\alpha }=\frac{\partial E_2({\mathbf{Q}}_0)}{\partial Q_{\alpha }}.$$

$$H_3=-\sqrt{\frac{S}{2}}\sum _{⟨ij⟩}$$

$$H_3=\frac{3JSi}{2}\sqrt{\frac{3}{2SN}}\sum _{1+2=3}({\bar{\gamma }}_1+{\bar{\gamma }}_2)(a_1^{†}{a}_2^{†}{a}_3-a_3^{†}{a}_1{a}_2),$$

Full text

PACS number(s)

78.70.Ck, 75.25.-J, 75.10.Jm

I Introduction

In a recent publication Cheng et. al., Ref. Luo et al., 2015, highlighted the features of the indirect $K$ -edge resonant inelastic x$-$ray scattering (RIXS) bi- and trimagnon spectrum of an isotropic triangular lattice antiferromagnet (TLAF). The TLAF is known to possess a $120^{\circ}$ long range ordered state even after quantum fluctuations are considered Singh and Huse (1992); Huse and Elser (1988); Bernu et al. (1994); Li et al. (2015); Ono et al. (2003); Kadowaki et al. (1987); Ishii et al. (2011); Poienar et al. (2010); Toth et al. (2011, 2012); Shirata et al. (2012); Koutroulakis et al. (2015); Susuki et al. (2013). The authors considered the self-energy corrections to the spin-wave spectrum to pinpoint the nontrivial effects of magnon damping and very weak spatial anisotropy on RIXS. It was shown that for a purely isotropic TLAF model, a multipeak RIXS spectrum appears which is primarily guided by the damping of the magnon modes. Interestingly enough it was demonstrated that the roton momentum point is immune to magnon damping (for the isotropic case) with the appearance of a single-peak RIXS spectrum. It was suggested that this feature could be utilized as an experimental signature to search for or detect the presence of roton like excitations in the lattice. However, including XXZ anisotropy leads to additional peak splitting, including at the roton wave vector.

At present no theoretical guidance exists for experimentalists on how to interpret the RIXS spectrum of the ordered phase in a geometrically frustrated triangular lattice quantum magnet, though a proposal has been put forward to detect spin-chirality terms in triangular-lattice Mott insulators via RIXS Ko and Lee (2011). Furthermore, as discussed in this article the existing spin wave theory formulation used for the isotropic case fails beyond the isotropic point and with Dzyaloshinskii-Moriya (DM) interaction included in the model.

Lately, the nature of the ground and excited states of the TLAF has garnered some attention Chen et al. (2013); Schmidt and Thalmeier (2014); Hauke et al. (2011); Kohno et al. (2007); Swanson et al. (2009); Fishman and Okamoto (2010); Ghioldi et al. (2015); Suzuki et al. (2014); Weichselbaum and White (2011); Hauke (2013). A high magnetic field phase diagram study of the TLAF has also been performed Starykh et al. (2014). An appropriate theoretical treatment of interactions in a TLAF must consider spin wave quantum fluctuation effects Fjærestad et al. (2007). Zero-point quantum fluctuations of a noncollinear ordered quantum magnet gives rise to spin Casimir effect Du et al. (2015, 2016). As a spin analog of the Casimir effect in vacuum, the spin Casimir effect describes the various macroscopic Casimir forces and torques that can potentially emerge from the quantum spin system. The physical consequence of the Casimir torque, generated due to the underlying lattice anisotropy, is the modification of the ordering wave vector, which is much smaller than the classical value. The modification in the ordering wave vector can cause the spin spiral state to become unstable, in turn rendering the standard spin wave theory expansion (1/S-SWT) approach inapplicable. Thus, the generic interacting spin wave theory is not appropriate.

To remedy the effect of singular behavior (which is not a precursor to the onset of quantum disordered phases) that naturally arises in noncollinear systems due to the presence of spin Casimir torque, Du et. al. Du et al. (2015, 2016), proposed the torque equilibrium spin-wave theory (TESWT). The regularization scheme of TESWT formalism removes the naturally occuring divergences within the interacting 1/S-SWT formalism of the anisotropic quantum lattice model. It was shown that TESWT gives a much closer final ordering vector to the results of series expansion (SE) and modified spin wave theory (MSWT) method Hauke et al. (2011); Weihong et al. (1999). Furthermore, its prediction of the phase diagram is consistent with the previous numerical studies Hauke et al. (2011); Weihong et al. (1999).

Historically, the concept of a roton minimum and a rotonlike point in the TLAF was introduced by Zheng et. al. Zheng et al. (2006a, b). Using SE method the authors identified a local minimum in the magnon dispersion at the high symmetry $M^{\prime}$ point, ($\pi$, $\pi/\sqrt{3}$). Drawing analogy with the appearance of a similar dip (local minimum) that is observed in the excitation spectra of superfluid 4He Feynman (1998) and the fractional quantum Hall effect Girvin et al. (1986), the authors proposed the “roton” nomenclature to describe the minimum in the magnon dispersion. The dip in the spectrum is also present at the other high symmetry $M$ point, $(0,2\pi/\sqrt{3})$, in the middle of the Brillouin zone (BZ) face edge. Zheng et. al. noted that a roton minimum is absent in the linear spin wave theory (LSWT) spectrum. Thus, the occurence of the rotonlike point is a consequence of quantum fluctuations arising in a frustated magnetic material Kubo and Kurihara (2014); Powalski et al. (2015). In a subsequent publication the concept of the rotonlike point was extended to the case of an anisotropic lattice by Fjaerestad et.al. Fjærestad et al. (2007). Additionally, a square lattice system with $J^{\prime}/J>2$ has also been predicted to support the roton minima Zheng et al. (2006a); Kubo and Kurihara (2014).

Further support of the roton feature was provided by the 1/S-SWT study of Starykh et.al. Starykh et al. (2006a). Based on their work it was proposed that rotons are part of a global renormalization (weak local minimum), with large regions of (almost) flat dispersion. The appearance of rotonlike minima and what was dubbed as a roton excitation has also been studied in an anisotropic spin-1/2 TLAF from the perspective of an algebraic vortex liquid theory Alicea et al. (2006); Alicea and Fisher (2007). Several anomalous roton minima were predicted in the excitation spectrum in the regime of lattice anisotropy where the canted Neel state appears. From the perspective of the algebraic vortex liquid theory formulated in terms of fermionic vortices in a dual field theory, it was proposed that the roton is a vortex anti-vortex excitation, thereby, lending credence to use of the word roton as an apt description. Rotons have also been predicted to exist in field induced TLAF magnetic systems Maksimov et al. (2016). The field-induced transformations in the dynamical response of the XXZ model create the appearance of rotonlike minima at the $K$ point. Experimental evidence of the rotonlike point can be found in recent inelastic neuron scattering (INS) spectrum of the $\alpha$-CaCr2O2 system Toth et al. (2011, 2012). Examples of TLAF where anisotropy and DM interaction are present are plethora Perkins and Brenig (2008); Perkins et al. (2013); Fjærestad et al. (2007); Shirata et al. (2012); Ono et al. (2003); Kadowaki et al. (1987); Ishii et al. (2011); Poienar et al. (2010); Toth et al. (2011, 2012); Coldea et al. (2001).

With advancements in instrumental resolution of the next-generation synchrotron radiation sources, RIXS spectroscopy presents itself as a novel experimental tool to investigate the nature of the bimagnon RIXS spectrum and the influence of the roton Dean (2015). As a spectroscopic technique RIXS has the ability to probe both single-magnon and multimagnon excitations across the entire BZ Ament et al. (2011); Nagao and Igarashi (2007); Haverkort (2010). Using RIXS it is possible to probe high energy excitations in cuprates Chen and Sushkov (2013); Pakhira et al. (2012). Considering the physical behavior that has been studied within the context of RIXS TLAF and the fact that departures from the isotropic triangular lattice geometry is a norm in a frustrated TLAF, this begs the question $-$ “What is the influence of spatial anisotropy and DM interaction on the bi- and trimagnon $K$ -edge indirect RIXS bimagnon spectrum at the rotonlike points and the other BZ points of an anisotropic triangular lattice ?”

In this article, we utilize material parameters relevant to Cs2CuCl4 to elucidate the $K$ -edge RIXS behavior of the rotonlike points and also the bimagnon behavior at the $Y$ point. We apply TESWT to our quantum Heisenberg model with spatial anisotropy and DM interaction on a triangular lattice. Using a TESWT upto first order in $1/S$, we compute the final ordering vector, the spin-wave energy, and phase diagram with different anisotropy parameters. We find the phase diagram has a physically consistent behavior in the ordering wave vector $\mathbf{Q}$. We find that the presence of a relatively small DM interaction can make the spiral state more stable. We calculate the interplay of x-ray scattering and bi- and trimagnon excitation. We find that the evolution of the RIXS spectra at rotonlike points is non-trivial. In the isotropic case all the rotonlike points are identical due to the 60*∘* rotation symmetry of the underlying isotropic triangular lattice. However, in the presence of symmetry breaking DM interaction terms the equivalence breaks down to give rise to two distinct points $-$ $M$ and $M^{\prime}$, see Fig. 9. Thus we investigate and track the evolution of the spectra at these two points separately. With increasing anisotropy the spectral weight at these points are subdued, even though the rotonlike points lie outside the region of magnon damping. Additionally, we find that the RIXS spectrum at the rotonlike $M$ point undergoes a spectral downshift. However, for the $M^{\prime}$ point the location of the peak is stable, albeit suppressed as the strength of the perturbation is increased. We also track the bimagnon RIXS evolution at the $Y$ point in the Brillouin zone to compare and contrast with the behavior at the rotonlike points. The spectrum at Y shows more peaks than at $M$ or $M^{\prime}$. Thus, the roton excitation spectrum is more stable Luo et al. (2015).

This paper is organized as follows. In Sec. II we introduce the model spin-1/2 anisotropic TLAF with DM interaction. In Sec. III.1 we state the spin wave formalism required to compute the wave vector renormalization (Sec. III.2) and renormalized dispersion (Sec. III.3). In Sec. IV, we extend the applicability of the TESWT formalism to include the effects of DM interaction. In Sec. IV.1, we elaborate on the TESWT method, compute the ordering vector and dispersion, and perform a TESWT INS fitting (Sec. IV.2). We then calculate the phase diagram in Sec. IV.3. In Sec. V we compute the indirect RIXS spectra. In Sec. V.1 we compute the non-interacting bi- and trimagnon spectrum. In Sec. V.2 we outline the formalism to compute the interacting bimagnon RIXS spectrum by including the quartic interactions. In Sec. V.3 we track the evolution of the roton energy to provide a physical explanation of the trend exhibited by the RIXS spectrum with anisotropy and DM interaction. In Sec. V.4 we state the results for the total indirect $K$ -edge RIXS intensity. Finally, in Sec. VI we provide our conclusions.

II Model Hamiltonian

The antiferromagnetic Heisenberg model on the anisotropic triangular-lattice is widely believed to be well described by Cs2CuBr4 and Cs2CuCl4 Ghorbani et al. (2016). While Cs2CuCl4 exhibits spin-liquid behavior over a broad temperature range Coldea et al. (2001, 2003), the Cs2CuBr4 compound exhibits a magnetically ordered ground state with spiral order in zero magnetic field Ono et al. (2003). For \alpha$$-CaCr2O4, though it is reported to have two inequivalent Cr3+ ions and four different exchange interactions, the nature of the distortion is such that the average of the exchange interactions along any direction is approximtely equal.

We consider the spin-1/2 antiferromagnetic Heisenberg model on the anisotropic triangular-lattice perturbed by a DM interaction, described by the Hamiltonian

[TABLE]

where $\langle ij\rangle$ refers to nearest-neighbor bonds on the triangular lattice, $J_{ij}=J$ denotes the exchange constants along the horizontal bonds and $J_{ij}=J^{\prime}$ the diagonal bonds, see Fig. 1. The asymmetric DM interaction between neighboring spins is given by

[TABLE]

where $\mathbf{D}=(0,D,0)$ with $(D>0)$ and $\boldsymbol{\delta}_{1,2}$ are the nearest neighbor vectors along the diagonal bonds as shown in Fig. 1. In the classical limit, the spin operators are replaced by the three-component vectors

[TABLE]

where the spin forms a spiral with the ordering vector $\mathbf{Q}$. The classical ground state energy is given by

[TABLE]

with

[TABLE]

where the dimensionless ratios $\alpha=J^{\prime}/J$ and $\eta=D/J$ denote the relative interaction strengths. For the determination of the ordering vector $\mathbf{Q}$ we have to minimize the classical ground state energy

[TABLE]

which amounts to finding the roots of the equations

[TABLE]

Anticipating that this condition leads to a spiral along the $x$ axis $\mathbf{Q}=(Q_{0},0)$, we obtain the solution in the absence of DM interaction as

[TABLE]

Apriori, it is not clear whether the classical ordering vector correctly describes the long-ranger order in the quantum frustrated system. In fact, the classical wave vector will be renormalized by quantum fluctuations as will be discussed in Sec. IV.1.

III Linear Spin-wave theory

III.1 $1/S$ expansion

Before we set up the spin-wave expansion, it is convenient to transform the spin components from the laboratory frame $(x_{0},z_{0})$ to the rotating frame $(x,z)$ through

[TABLE]

where $\theta_{i}=\mathbf{Q}\cdot\mathbf{r}_{i}$. The rotating Hamiltonian takes the form

[TABLE]

where we have defined

[TABLE]

SWT amounts to applying the Holstein-Primakoff (HP) transformation to bosonize the rotating Hamiltonian (12)

[TABLE]

where $n_{i}=a_{i}^{\dagger}a_{i}$ and $a_{i}^{\dagger}$ ($a_{i}$) is the magnon creation (annihilation) operator for a given site $i$. Under the assumption of diluteness of the HP boson gas, $n_{i}/(2S)\ll 1$, one arrives at the interacting spin-wave Hamiltonian to the first order expansion of the square root

[TABLE]

where the first term is the classical energy and $H_{n}$ denotes terms of the $n^{th}$ power in the HP boson operators $a^{\dagger}(a)$.

III.2 Quadratic terms: first-order corrected LSWT

After Fourier transformation we obtain the quadratic Hamiltonian in momentum space as

[TABLE]

with

[TABLE]

where

[TABLE]

Diagonalization of $H_{2}$ is performed with the canonical Bogoliubov transformation

[TABLE]

with the parameters $u_{\mathbf{k}}$ and $v_{\mathbf{k}}$ defined as

[TABLE]

As a result we obtain the linear spin-wave dispersion

[TABLE]

It is noted that the magnon spectrum has zeros at $\mathbf{k}=0$ while a gap is opened at $\mathbf{k}=\mathbf{Q}$ in the presence of DM interaction. The diagonalized Hamiltonian $H_{2}$ is given by

[TABLE]

where the zero-point energy

[TABLE]

is the $1/S$ correction to the classical ground-state energy.

Generally, the first-order correction of LSWT $\mathbf{Q}_{1}=\mathbf{Q}_{0}+\Delta\mathbf{Q}$ is determined by minimizing the sum $E_{0}(\mathbf{Q})+E_{2}(\mathbf{Q})$

[TABLE]

Neglecting higher order terms, we obtain

[TABLE]

with

[TABLE]

A straightforward calculation gives $1/S$ correction to the classical wave vector

[TABLE]

where

[TABLE]

In Fig. 2 we display the variation of the ordering wave vector renormalization against lattice anisotropy computed using LSWT, 1/S corrected LSWT, and TESWT. It is clear that while the LSWT formulation extends the spiral phase region, the first-order correction from 1/S-LSWT gives an unphysical result as $\alpha\rightarrow 2$ while $\eta=0$. Inclusion of DM interaction rounds the singularity with an angle that is greater than $2\pi$. The root cause of this divergence originates from spin Casimir torque Du et al. (2015, 2016). In a frustrated spiral system, the strong quantum fluctuation effect leads to failure in the first-order correction. In Sec. IV we will discuss and implement the TESWT approach which offers a solution to this issue. The equations to generate the TESWT results are reported in that section.

III.3 Cubic and quartic terms: renormalized dispersion

The 1/S correction to the spin wave dispersion has to be accounted for in a non-collinear structure. The interplay of magnon decay as it arises from the non-collinear structure is also considered Zhitomirsky and Chernyshev (2013); Chernyshev and Zhitomirsky (2006); Mourigal et al. (2013). The three-boson term that arises from the coupling between transverse and longitudinal fluctuations in the noncollinear spin structure takes the form Chernyshev and Zhitomirsky (2009),

[TABLE]

In momentum space, we obtain

[TABLE]

where we have defined

[TABLE]

In the above we have adopted the convention that $1=\mathbf{k}_{1}$, $2=\mathbf{k}_{2}$, etc. For example, $a_{1}\equiv a_{\mathbf{k}_{1}}$. Performing the Bogoliubov transformation in $H_{3}$ we obtain the interaction terms expressed via the magnon operators as

[TABLE]

The three-boson vertices are given by

[TABLE]

with $\bar{V}_{a,b}$ given by

[TABLE]

We notice that the three-magnon vertices are of order $1/\sqrt{S}$ relative to the linear spin-wave Hamiltonian and they must occur in pairs in any self-energy or polarization diagram. The quartic term $H_{4}$ in the interacting spin-wave Hamiltonian (16) reads

[TABLE]

To obtain the explicit forms of the quasiparticle representation of $H_{4}$, we introduce the following mean-field averages

[TABLE]

The Hartree-Fock decoupling of the $H_{4}$ yields the quadratic Hamiltonian

[TABLE]

where

[TABLE]

We then obtain the Hartree-Fock corrected $H_{2}$ term as

[TABLE]

where

[TABLE]

Finally, the normal-ordered quartic term $\tilde{H}_{4}$ in the quasiparticle representation describes the multi-magnon interactions. In the hierarchy of $1/S$ expansion, terms relevant for our calculations are the lowest order irreducible two-magnon scattering amplitude

[TABLE]

with the vertex function given by

[TABLE]

where we have defined

[TABLE]

The effective $1/S$ interacting spin$-$wave Hamiltonian in terms of the magnon operators reads

[TABLE]

At zero temperature the bare magnon propagator is defined as

[TABLE]

The first order $1/S$ correction to the magnon energy is determined by the Dyson equation

[TABLE]

with the one-loop self-energy $\Sigma(\mathbf{k},\omega)=\Sigma_{a}(\mathbf{k},\omega)+\Sigma_{b}(\mathbf{k},\omega)+\Sigma_{c}(\mathbf{k})$, where $\Sigma_{c}(\mathbf{k})=\delta\varepsilon_{\mathbf{k}}$ is a frequency-independent Hartree-Fock correction, while $\Sigma_{a,b}(\mathbf{k},\omega)$ are calculated as

[TABLE]

The on-shell solution consists of setting $\omega=\varepsilon_{\mathbf{k}}$ in the self-energy Eqs. (51) and (52) leads to the following expression for the $1/S$ renormalized spectrum

[TABLE]

where $\bar{\omega}_{\mathbf{k}}=\mathrm{Re}[\omega_{\mathbf{k}}]$ is the renormalized spin$-$wave energy and $\Gamma_{\mathbf{k}}=-\mathrm{Im}[\omega_{\mathbf{k}}]$ represents the magnon decay rate. In Fig. 3, we plot the $1/S$ LSWT dispersion of Cs2CuCl4 Fjærestad et al. (2007).

IV Torque Equilibrium Spin Wave Theory

Zero-point quantum fluctuation in a non-collinear ordered spin structure can lead to deviations in the measured ordering wave vector compared to the classical one. The correction emerging from the spin Casimir effect is usually neglected, but it was recently shown that this is not a bonafide assumption. In Du et. al. Du et al. (2015, 2016) it was clearly established that in certain situations a standard spin wave theory is no longer applicable due to the spin Casimir quantum effect, even when the system is long-range ordered. An important consequence of these quantum fluctuations is on the spiral state which can become unstable, which is different from the case of long-range-order melting. As mentioned earlier the classical signatures of these instabilities are the divergences of the ordering wave vector at the quantum critical point and the strongly singular one-loop expansions of the energy spectrum and the sublattice magnetization. In this section, we extend the applicability of the TESWT formalism to include the effects of DM interaction in an anisotropic TLAF. Using INS experimental data from Cs2CuCl4 Coldea et al. (2003), we obtain fitting parameters for the exchange constants and DM interactions utilized in subsequent indirect $K$ -edge RIXS calculations.

IV.1 TESWT formalism

Spin Casimir effect will change the classical ground state to a new saddle point. This new ground state can be unambiguously determined once we compute the value of $\mathbf{Q}$. An ordinary approach is considering the $1/S$ correction $\Delta Q$, as we show in Sec. III.2. However, such a method gives an unphysical result, see Fig. 2. As $\alpha\rightarrow 2$, the $1/S$ correction $\Delta Q$ becomes infinites.

The basic idea of TESWT is to minimize the ground state energy. The spin Casimir torque is defined as

[TABLE]

where $|\Psi_{vac}\rangle$ represents the quasiparticle vacuum state. Then the torque equilibrium condition is

[TABLE]

where $\mathbf{Q}$ is the final ordering vector, $H_{cl}=E_{0}(\mathbf{Q})$ is the classical energy. Using the fact that the spin-wave spectrum function $\varepsilon_{\mathbf{k}}$ is only well defined at $\mathbf{Q}_{cl}$, we try to find a system whose classical ordering vector is $\mathbf{Q}$ for convenience of calculation. Thus we shift the function depending on classical ordering vector $\mathbf{Q}_{cl}$ to $\mathbf{Q}$ by

[TABLE]

where $\widetilde{H}_{2},\widetilde{A}_{\mathbf{k}}$ and $\widetilde{B}_{\mathbf{k}}$ are functions of another spin system whose classical ordering vector $\widetilde{\mathbf{Q}}_{cl}$ equals $\mathbf{Q}$. The counterterm is given by $H^{c}_{2}$ whose effects are considered in the $A_{\mathbf{k}}(B_{\mathbf{k}})$ coefficients through $A^{c}_{\mathbf{k}}(B^{c}_{\mathbf{k}})$. In principle, we have many combinations of $(\widetilde{\alpha},\widetilde{\eta})$ that satisfy this condition. As $\eta/\alpha$ is small, within perturbation theory, we believe $\widetilde{\eta}=\eta$ is a reasonable choice. Thus the new parameters can be deduced by solving the following self-consistent equations

[TABLE]

The spin Casimir torque is then expressed approximately as $\mathbf{T}_{sc}(\mathbf{Q})=\widetilde{\mathbf{T}}_{sc}(\mathbf{Q})$. Thus the torque equilibrium equation in Eq.( 55) can be written as

[TABLE]

Note, the exchange parameters on the left-hand side of the equation are exact as $\alpha,\eta$. While the parameters on the right-hand side approximate as $\widetilde{\alpha}=-2\cos{\frac{Q}{2}}-\eta\cot{\frac{Q}{2}}$. We solve this equation numerically and give the results in Fig. 2. If there is no DM interaction, TESWT gives $Q=2\pi$ for $\alpha\geq 1.2$, which are similar to the results of numerical methods Hauke et al. (2011); Weihong et al. (1999). The LSWT, however, gives a wider region for spiral order phase, can’t describe the region for $1.2\leq\alpha\leq 2$. As anticipated, even a small DM interaction, $\eta=0.05$, changes our final ordering vector. The DM interaction improves the spiral order stabilization and enlarges it’s region of validity.

We diagonalize $\widetilde{H}_{2}(\widetilde{\alpha},\widetilde{\eta},\mathbf{Q})$ and treat $H_{2}^{c}$ as a counterterm. Since we are considering a $1/S$ theory, we neglect the counterterm contributions from $H_{3}^{c}$ and $H_{4}^{c}$ Du et al. (2015, 2016). Thus, we can write the Hamiltonian as

[TABLE]

Following the procedure outlined in Sec. III, the effective TESWT Hamiltonian now reads

[TABLE]

where $\widetilde{F}$ means $F(\widetilde{\alpha},\widetilde{\eta},\mathbf{Q})$ ($F$ is an arbitrary operator) and

[TABLE]

Thus, we shifted the classical ordering vector $\mathbf{Q}_{cl}$ to the final ordering vector $\mathbf{Q}$ using TESWT. Therefore, the first order $1/S$ corrected magnon dispersion can now be changed to

[TABLE]

IV.2 INS fitting

As discussed above, with anisotropy the application of 1/S-LSWT formalism is tricky. But, application of TESWT requires magnetic interaction parameters computed within that formalism. The most direct way to achieve this goal is to compare the theoretical dispersion with the experimental data. We fit the INS data of Cs2CuCl4 Coldea et al. (2003) to Eq. (64) using iterative least squares estimation both by TESWT and $1/S$-LSWT. Our fitting parameters along with results from other sources are reported in Table. 1. Our dispersion line fits are reported in Fig. 3. The absence of higher order terms within our TESWT could be a source of disagreement with the series expansion results Fjærestad et al. (2007), which is an all numerical method that considers higher order terms Zheng et al. (2006b). As the fitted dispersion by TESWT gives a reasonable comparison with the experimentally fitted SE method parameters, we believe that our TESWT can capture the essential physical behavior. While it maybe fruitful to investigate the above mentioned discrepancy, within the context of our RIXS calculation we do not expect the improved interaction constants to bring about much qualitative or quantitative differences.

IV.3 Sublattice magnetization

Next, we study the phase diagram of the anisotropic triangular-lattice. In a spin system, the sublattice magnetization can describe the phase transition behavior. The second-order correction of the sublattice magnetization contributes little to the result. Thus, we only consider the first order correction to the sublattice magnetization as

[TABLE]

where

[TABLE]

In Fig. 4 we plot the sublattice magnetization $\langle S\rangle$ variation with spatial anisotropy. Our result without DM interaction is consistent with previous numerical studies Hauke et al. (2011); Weihong et al. (1999). Consistent with our previous analysis of Fig. 2, the spiral order is destroyed at $\alpha\geq 1.2$. In addition, the spiral order is unsafe at $\alpha\leq 0.5$, consistent with modified spin wave results Hauke et al. (2011). The DM interaction, which originates from spin-orbit coupling, helps to generate a non-collinear spin ground state. It is evident from Fig. 4, as $\eta$ gets bigger, the phase transformation point in the region $\alpha\leq 0.5$ diminshes until it disappears. On the opposite end, the sublattice magnetization recovers thereby making the $\alpha\geq 1.2$ zone less susceptible to drastic effects of quantum fluctuation. These findings suggest that the DM interaction enlarges the region of the spiral state. Our focus in this article is on the multimagnon RIXS spectrum in the spiral phase. Thus, we can use the computed phase diagram to extract the appropriate choice of parameters. We find that TESWT not only gives a consistent physical estimate of the final ordering vector, but also correctly predicts the phase diagram of an anisotropic TLAF, helping to better understand the behavior of the spiral ground state of such a geometrically frustrated system.

V Indirect RIXS spectra

V.1 Noninteracting bi- and trimagnon RIXS

In this section we calculate the bi- and trimagnon RIXS spectrum. The results in this part use TESWT while the LSWT approach is shown in Appendix A. The indirect RIXS scattering operator, is given by van den Brink (2007); Forte et al. (2008)

[TABLE]

where $\mathbf{q}$ is the scattering momentum. In quasiparticle representation, the magnon creation parts of the RIXS scattering operator can be given by

[TABLE]

where the bimagnon scattering matrix element is

[TABLE]

and the trimagnon scattering matrix element is

[TABLE]

We neglect the corrections from magnon interactions for the trimagnon intensity, which appear at $1/S^{2}$ order. Next, using Eqs. (92) and  (93) stated in Appendix A we obtain the following expressions for $I_{2}(\mathbf{q},\omega)$ (noninteracting bimagnon) and $I_{3}(\mathbf{q},\omega)$ (trimagnon) scattering intensity

[TABLE]

where $\omega_{\mathbf{k}}^{(0)}=\widetilde{\varepsilon}_{\mathbf{k}}+\varepsilon_{\mathbf{k}}^{c}$.

In Fig. 5 we display our results of the noninteracting bi- and trimagnon RIXS spectra at various points across the BZ. Overall the agreement between the LSWT and the TESWT formalism is reasonable. Our TESWT result generates more peaks for the bimagnon intensity. We note that in the isotropic regime $\alpha=1$, our TESWT results are identical with the LSWT formalism since the final ordering vector Q equals the classical vector Qcl, see Fig. 11. As discussed earlier, the TESWT is the physically correct formalism in the presence of anisotropy.

V.2 Interacting bimagnon RIXS spectra

We now proceed with the analysis of $1/S$ correction to the two-magnon Green’s function by taking into account both the self-energy correction to the single magnon propagator $G$ according to the Dyson equation and the vertex insertions to the two-particle propagator $\Pi$ which satisfies the Bethe-Salpeter (BS) equation Canali and Girvin (1992); Luo et al. (2015)

Using the procedure outlined in our prior work Luo et al. (2015) and Feynman rules in momentum space, we obtain the following equations for the two$-$particle propagator and the associated vertex function as

[TABLE]

where the basic one-magnon propagator up to $1/S$ order is now given by

[TABLE]

The lowest order two-particle irreducible interaction vertex in Fig. 6(c) reads

[TABLE]

in which the frequency-independent four-point vertex $\mathcal{V}_{4}$ coming from the quartic Hamiltonian can be written as

[TABLE]

and the other four vertices $\mathcal{V}_{3}^{(a-d)}$ in the same $1/S$ order which are assembled from two three-point vertices and one frequency-dependent propagator can be written as

[TABLE]

In the above we have retained only the bare propagator $\mathrm{G}_{0}$ for each intermediate line in $\mathcal{V}_{3}^{(a-d)}$ in the spirit of $1/S$ expansion. Note, the vertex expressions here are different from those stated within the traditional 1/S-SWT approach Luo et al. (2015). The vertex expressions here are shifted by the correct TESWT wave vector as represented by the tilde notation. Based on the above generalization, we now derive the final solution of the interacting RIXS intensity from the ladder approximation BS equation.

We adopt a numerical approach to compute the interacting bimagnon RIXS intensity. We assume that two on-shell magnons are created and annihilated in the repeated ladder scattering process with $\omega^{\prime}\approx-\omega_{\mathbf{k}}^{(0)}=-\widetilde{\varepsilon}_{\mathbf{k}}-\varepsilon_{\mathbf{k}}^{c}$ and $\omega_{1}\approx-\omega_{\mathbf{k}_{1}}^{(0)}=-\widetilde{\varepsilon}_{\mathbf{k}_{1}}-\varepsilon_{\mathbf{k}_{1}}^{c}$. We substitute (73) and (V.2) into (94) to obtain

[TABLE]

where $\Pi_{\mathbf{k}}=2[\omega-\omega_{\mathbf{k}+\mathbf{q}}-\omega_{\mathbf{k}}+i0^{+}]^{-1}$ is the renormalizated two-magnon propagator in the absence of vertex correction. To proceed further we divide the BZ into $N$ points and replace the continuous momenta $(\mathbf{k},\mathbf{k}^{\prime},\mathbf{k}_{1})$ with discrete variables $(\mathrm{m,n,l})$. Thus, we can write

[TABLE]

where

[TABLE]

Adopting the matrix notation $\Gamma=(\hat{1}-V\Pi)^{-1}$ we obtain the final form of the $\mathrm{\hat{\chi}}$ matrix as

[TABLE]

where we have defined the following $N\times N$ matrices,

[TABLE]

The interacting bimagnon RIXS susceptibility is computed as

[TABLE]

We use Eqs. (73) - (88) and Eq. (92) stated in Appendix A to numercially compute our interacting bimagnon RIXS intensity at $M$, $M^{\prime}$, and $Y$ BZ points.

In Fig. 7 we show the spectra at the $Y$ point. The first panel is a reproduction of our previous result reported in Ref. Luo et al., 2015. In Fig. 7(b) we display the result of TESWT Cs2CuCl4 RIXS. Compared to the isotropic case or to the other anisotropic situations, panels (c) and (d), this spectrum is substantially broadened. With enhanced anisotropy the lattice can be envisioned as disintegrating into a set of loosely coupled chains. Thus, instead of bimagnons one can expect the emergence of spinons as is expected in 1d systems. 1d RIXS has been able to capture multi-spinon excitations Schlappa et al. (2009, 2018). Thus, the predicted RIXS spectrum feature could be used to confirm quasi-1d to 2d dimensional crossover features of Cs2CuCl4 Starykh et al. (2006b). In Figs. 7(c) or 7(d) we can compare the effects of including a tiny DM interaction. We find that there is a prominent low energy peak with a relatively muted higher energy response. This tiny DM interaction does not bring about any spectral down- or up- shift. The spectral weight is simply redistributed.

V.3 RIXS signatures at roton points

In Fig. 8 we display the interacting RIXS intensity variation at the two anisotropic roton points $\mathbf{q}=M$ and $\mathbf{q}=M^{\prime}$ with varying lattice anisotropy and DM interaction. The anisotropy parameter choices ensure that the TLAF does not decouple into a set of loosely coupled 1d chains, where the bosonization description has been shown to apply  Starykh et al. (2006b). The upper panel Figs. 8(a) and 8(b) are results for zero DM interaction. Note, the two spectrum coincide in the isotropic limit since the two roton points are equivalent due to $C_{3v}$ symmetry of the isotropic triangular lattice Luo et al. (2015), while they evolve differently in the presence of spatial anisotropy. In particular, we find that the roton spectra at $\mathbf{q}=M$ point (the roton point along the $k_{y}$ direction in BZ) is very sensitive to anisotropy. Though the single-peak structure is stable against $J^{\prime}/J$, the peak position undergoes a spectral downshift with increased anisotropy. On the other hand, for the $q=M^{\prime}$ point (along the diagonal BZ direction), the peak location of the spectra does not change much, in comparison to the $M$ point, in the presence of anisotropy. In the lower panel, Figs. 8(c) and 8(d), we display the behavior of the RIXS spectra with DM interaction. Contrary to the isotropic case, the presence of DM interaction introduces a spectral upshift at both $\mathbf{q}=M$ and $\mathbf{q}=M^{\prime}$. The dashed line in the lower panel is the result of using realistic parameters generated from the Cs2CuCl4 INS data fit based on TESWT.

To gain insight into the roton behavior of the RIXS spectra we track the evolution of the roton minimum in the single magnon dispersion along $\Gamma\rightarrow M$ and $M^{\prime}$, both parallel and perpendicular to the BZ path, see Fig. 9. A bimagnon excitation requires $\omega_{\mathbf{k}+\mathbf{q}}+\omega_{\mathbf{k}}$ amount energy. We notice that the one magnon dispersion along $M$ displays more sensitivity compared to that along $M^{\prime}$. The asymmetrical sensitivity to the dispersion stiffness explains the origins of the differing roton RIXS spectra behavior. Increasing anisotropy reduces the one magnon energy (softening) near the $M$ point (the first column in Fig. 9), thus leading to a spectral downshift in Fig. 8. Whereas for the $M^{\prime}$ point, the overall energy scale of the dispersion is not affected (the third column in Fig. 9). We observe neither a drastic hardening nor softening. Thus, the RIXS spectrum holds steady without any shift. The softening and subsequent flattening of the dispersion at the $M$ point suggests that for the anisotropic TLAF, the roton feature is retained more at the $M$ point compared to the $M^{\prime}$. However, inclusions of the DM interaction increases the one magnon energy both near $M$ and $M^{\prime}$ points (the second and fourth column in Fig. 9), introducing a spectral upshift. This could be understood by the fact that DM interaction introduces a gap, thus it requires more energy to create a single magnon and in turn a bimagnon excitation.

The evolution of the spectral height in Fig. 8 can also be explained. As anisotropy weakens the coupling between the TLAF spins to transform the material to a quasi-1d spin chain, it is more difficult to create a bimagnon excitation. In RIXS, this will cause a decrease in the value of the bimagnon scattering matrix element $|\widetilde{M}(\mathbf{k}+\mathbf{q},-\mathbf{k})|$ in turn leading to a reduction in the spectral weight, see Fig. 8(a) and 8(b). On the contrary, the presence of the DM interaction encourages interactions beyond the traditional Heisenberg type. Thus, it assists with the creation of bimagnons, see Fig.8(c), where the spectral weight increases. But for the ${\bf q}=M^{\prime}$ point, the actual nature of the magnon bands is not affected by the DM interaction, see Fig. 9 fourth column. Thus, the height of the RIXS spectrum does not change with DM interaction. Note, in all the above discussion we have assumed that the triangular lattice does not break down to a set of coupled 1d spin chains. The $\alpha=0.5$ RIXS spectra could well describe the Cs2CuBr4 compound.

V.4 Total RIXS

In Fig. 10 we report the total RIXS spectrum for Cs2CuCl4 with TESWT fitting parameters. The total RIXS spectrum comprises of the bi- and trimangon response. We use Eqs. (92) and  (93) to compute the spectrum. The interacting bimagnon (Eq. (88)) and noninteracting trimagnon intensity (Eq.(72)) are summed over to get the total RIXS spectrum. As expected, the trimagnon peak is located at a higher energy than the bimagnon response. In the response for the $M$ and $Y$ points, the main peaks are separated, see Figs. 10(a) and 10(c). At the $M^{\prime}$ point in Fig. 10(b), a small bigmagnon peak is obvious while the main peaks of bi- and trimagnon are mixed.

We note that the spectrum height of the bimagnon undergoes a special evolution. Bimagnon has a height near the boundary of BZ ($M$ and $M\prime$ points) but vanishes when it is close to the center of BZ ($Y$ point). A similar trend for the bigmanon can also be observed in Figs. 5 and 11. This is due to the behavior of the RIXS scattering element from the indirect $K$ -edge RIXS scattering operator in Eq. (67). For wave vector choice $\mathbf{q}$ close to the high symmetry $\Gamma$ point, the RIXS bimagnon matrix element occuring from $\mathcal{R}_{\mathbf{q}}$ gives a vanishingly small contribution. Thus the spectral weight of the bimagnon is substantially weakened near the $\Gamma$ point. Without DM interaction, the contribution is purely from the trimagnon excitations at the $\Gamma$ point in the isotropic TLAF, see Fig. 11(a). The above observations on the total RIXS spectrum should be helpful in distinguishing the contributions of the two different multimagnon excitations.

VI Conclusion

Due to the possible realization of various unusual ordered or disordered phases, frustrated magnetism is an active area of research in condensed matter physics Starykh (2015). Traditionally, information on the magnetic ground state and single magnon excitations is inferred from inelastic neutron scattering (INS) experiments Rønnow et al. (2001); Coldea et al. (2001). However, with the advent of RIXS spectroscopy experimentalists now have a probe that can comprehensively investigate a wide range of energy and momentum values in BZ.

In this article, we have demonstrated the application of a recently proposed spin-wave theory scheme called TESWT to the indirect $K$ -edge RIXS. As highlighted in this paper it is not a trivial matter to ensure that the sanctity of the spin spiral state is preserved. We performed a TESWT fitting of Cs2CuCl4 INS data, which gives $\alpha\approx 0.316$ and $\eta\approx 0.025$. Using these realistic parameters we computed the indirect $K$ -edge bi- and trimagnon RIXS spectra within TESWT formalism. Our results allow us to confirm that in contrast to the isotropic model, quantum fluctuations in the noncollinear anisotropic TLAF can generate divergent fluctuations with drastic effects on the magnetic phase diagram. We find that the behavior of the RIXS spectra is influenced with the occurence of two inequivalent rotonlike points, $M(0,2\pi/\sqrt{3})$ and $M^{\prime}(\pi,\pi/\sqrt{3})$. While the roton RIXS spectra at the $M$ point undergoes a spectral downshift with increasing anisotropy, the peak at the $M^{\prime}$ is not affected. However, the peak at $M^{\prime}$ does not exhibit any downshift. We believe in the anistorpic case the $M$ point retains more of the roton feature. Finally, we find that in the total RIXS spectra, the features of the bimagnon and the trimagnon are certainly different and thus can be easily distinguished within an experimental setting. While resolution and intensity issues may plague the $K$ -edge, we hope the calculation in this paper and our past publication Luo et al. (2015) will inspire experimentalists to improve resolution to test our predicted $K$ -edge RIXS behavior.

In conclusion, our theoretical investigation of the indirect RIXS intensity in the spiral antiferromagnets on the anisotropic triangular lattice demonstrates that RIXS has the potential to probe and provide a comprehensive characterization of the dispersive bimagnon and trimagnon excitations in the TLAF across the entire BZ, which is far beyond the capabilities of traditional low$-$energy optical techniques Devereaux and Hackl (2007); Vernay et al. (2007); Perkins and Brenig (2008); Perkins et al. (2013).

Acknowledgements.

We thank Radu Coldea for sharing with us the INS data for Cs2CuCl4. T D. acknowledges invitation, hospitality, and kind support from Sun Yat-Sen University Grant No. OEMT–2017–KF–06. T. D. acknowledges funding support from Augusta University Scholarly Activity Award. S. J., C. L and D. X. Y. are supported by NKRDPC Grants No. 2018YFA0306001, No. 2017YFA0206203, NSFC-11574404, NSFG-2015A030313176, National Supercomputer Center in Guangzhou, and Leading Talent Program of Guangdong Special Projects.

Appendix A Isotropic TLAF RIXS spectra

In this Appendix we compare the results of LSWT and TESWT for the isotropic lattice. We apply linear spin wave theory to the calculation of indirect $K$ -edge RIXS spectrum in this section. After the usual HP and Bogoliubov transformation application, the magnon creation parts of the RIXS scattering operator can be expressed as

[TABLE]

where the bimagnon and trimagnon scattering matrix element expression are given by

[TABLE]

Note that all the coefficients and functions are defined at the classical ordering vector $\mathbf{Q}_{cl}$ in LSWT. The frequency and momentum dependent magnetic scattering intensity is related to the multimagnon RIXS response function via the fluctuation-dissipation theorem

[TABLE]

where the total indirect $K$ -edge RIXS susceptibility is given by

[TABLE]

In the above $\chi_{2}(\mathbf{q},\omega)$ could be either a noninteracting or interacting two$-$magnon susceptibility, but $\chi_{3}(\mathbf{q},\omega)$ is the non$-$interacting three$-$magnon susceptibility. The susceptibilities can be expressed explicitly from the corresponding multi-magnon Green’s function defined as

[TABLE]

where $\Pi$ and $\Lambda$ denote the bi- and trimagnon propagator, respectively. The momentum-dependent two-magnon and three-magnon Green’s function in terms of Bogoliubov quasiparticles are defined as

[TABLE]

where $\mathcal{T}$ is the time-ordering operator and $\langle\cdot\rangle$ is the average of the ground state. Using Eq. (96) and Eq. (97), we can compute the noninteracting and the interacting RIXS spectra. The non-interacting spectrum can be calculated by applying Wick’s theorem to Eq. (96) and Eq. (97). The final expressions are stated in Eqs. (71) and  (72).