Two-magnon excitations in resonant inelastic x-ray scattering studied by spin-density-wave formalism

TL;DR

This paper models two-magnon excitations in RIXS at the transition-metal K-edge using a Hubbard-type model and SDW formalism, providing results that align well with experiments and revealing new spectral features.

Contribution

It introduces a Hubbard-type $d$-$p$ model combined with SDW formalism for RIXS, offering a more detailed approach than previous Heisenberg models.

Findings

Calculated RIXS spectra match experimental momentum dependence.

Spectra show large two-magnon weight near the zone center.

Differences from Heisenberg models suggest new experimental tests.

Abstract

We study two-magnon excitations in resonant inelastic x-ray scattering (RIXS) at the transition-metal $K$-edge. Instead of working with effective Heisenberg spin models, we work with a Hubbard-type model ($d$-$p$ model) for a typical insulating cuprate La$_2$CuO$_4$. For the antiferromagnetic ground state within the spin-density-wave (SDW) mean-field formalism, we calculate the dynamical correlation function within the random-phase approximation (RPA), and then obtain two-magnon excitation spectra by calculating the convolution of it. Coupling between the $K$-shell hole and the magnons in the intermediate state is calculated by means of diagrammatic perturbation expansion in the Coulomb interaction. Calculated momentum dependence of RIXS spectra agrees well with that of experiments. A notable difference from previous calculations based on the Heisenberg spin models is that RIXS spectra have a large two-magnon weight near the zone center, which may be confirmed by further careful high-resolution experiments.

Equations

$$H_{FL}=\alpha \sum _{⟨i,j⟩}[\mathbf{e}\cdot {\mathbf{r}}_{ij}][{\mathbf{e}}^{\prime }\cdot {\mathbf{r}}_{ij}]{\mathbf{s}}_i\cdot {\mathbf{s}}_j,$$

$$H=H_d+H_s+H_{sd}+H_p+H_x.$$

$$H_d$$

$$H_0$$

$$H^{\prime }$$

$$H_s=\sum _i\sum _{\sigma }{\epsilon }_{1s}{s}_{i\sigma }^{†}{s}_{i\sigma }=\sum _{\mathbf{k}}\sum _{\sigma }{\epsilon }_{1s}{s}_{\mathbf{k}\sigma }^{†}{s}_{\mathbf{k}\sigma },$$

$$H_{sd}$$

$$H_p=\sum _{\mathbf{k}}\sum _{\mu }\sum _{\sigma }{\epsilon }_{4p\mu }(\mathbf{k})p^{\prime }{}_{\mathbf{k}\mu \sigma }^{†}{p}^{\prime }{}_{\mathbf{k}\mu \sigma },$$

$$H_x=\sum _{\mathbf{k},\mathbf{q}}\sum _{\mu }\sum _{\sigma }{w}_{\mu }(\mathbf{q},\mathbf{e}){\alpha }_{\mathbf{q}\mathbf{e}}{p}^{\prime }{}_{\mathbf{k}+\mathbf{q}\mu \sigma }^{†}{s}_{\mathbf{k}\sigma }+h.c.,$$

$$w_{\mu }(\mathbf{q},\mathbf{e})=-\frac{e}{m}\sqrt{\frac{2\pi }{|\mathbf{q}|}}\mathbf{e}\cdot ⟨4p_{\mu }|\mathbf{p}|1s⟩\propto \mathbf{e}\cdot {\mathbf{e}}_{\mu },$$

$$H_{MF}^{\prime }$$

$$H_{d,MF}=\sum _{\mathbf{k}}\sum _j{E}_j(\mathbf{k})a_{j\mathbf{k}}^{†}{a}_{j\mathbf{k}},$$

$$d_{{\mathbf{k}}_{\sigma }\sigma }$$

$$d_{{\mathbf{k}}_{\sigma }\sigma }^{†}$$

$$\hat{G}({\bf k},\xi)=\left[\begin{array}[]{cc}G_{\uparrow\uparrow}({\bf k},\xi)&G_{\uparrow\downarrow}({\bf k},\xi)\\ G_{\downarrow\uparrow}({\bf k},\xi)&G_{\downarrow\downarrow}({\bf k},\xi)\end{array}\right].$$

$$G_{{\sigma }_1{\sigma }_2}(\mathbf{k},\xi )=\sum _j{u}_{{\sigma }_1,j}(\mathbf{k})u_{{\sigma }_2,j}^{\ast }(\mathbf{k})G_j(\mathbf{k},\xi ),$$

$$G_j(\mathbf{k},\xi )=\frac{1}{\xi -E_j(\mathbf{k})}.$$

$$G_j^R(\mathbf{k},\xi )$$

$$G_j^A(\mathbf{k},\xi )$$

$$G_j^C(\mathbf{k},\xi )$$

$$n_d$$

$$\mathbf{m}$$

$$\Delta ={\epsilon }_d+U-{\epsilon }_p\approx 4\mathrm{eV},$$

$$J\approx\frac{4t_{dp}^{4}}{\Delta^{2}}\biggl{[}\frac{1}{U}+\frac{1}{2\Delta}\biggr{]}\approx 0.15\,{\rm eV}.$$

$${\Pi }_{{\sigma }_1^{\prime }{\sigma }_2^{\prime },{\sigma }_2{\sigma }_1}(\mathbf{q},t^{\prime }-t)=⟨{\rho }_{{\sigma }_1^{\prime }{\sigma }_2^{\prime }}(\mathbf{q},t^{\prime }){\rho }_{{\sigma }_2{\sigma }_1}(-\mathbf{q},t)⟩,$$

$${\rho }_{{\sigma }_1{\sigma }_2}(\mathbf{q},t)$$

$${\rho }_{{\sigma }_1{\sigma }_2}(\mathbf{q})$$

$${\Pi }_{{\sigma }_1^{\prime }{\sigma }_2^{\prime },{\sigma }_2{\sigma }_1}(\mathbf{q},\omega )=\frac{{\chi }_{{\sigma }_1^{\prime }{\sigma }_2^{\prime },{\sigma }_2{\sigma }_1}^R(\mathbf{q},\omega )-{\chi }_{{\sigma }_1^{\prime }{\sigma }_2^{\prime },{\sigma }_2{\sigma }_1}^A(\mathbf{q},\omega )}{i(1-e^{-\omega /T})},$$

$${\chi }_{{\sigma }_1^{\prime }{\sigma }_2^{\prime },{\sigma }_2{\sigma }_1}^R(\mathbf{q},t^{\prime }-t)$$

$${\chi }_{{\sigma }_1^{\prime }{\sigma }_2^{\prime },{\sigma }_2{\sigma }_1}^A(\mathbf{q},t^{\prime }-t)$$

$${\chi }_{{\sigma }_1{\sigma }_2,{\sigma }_3{\sigma }_4}^{R,A}(\mathbf{q},\omega )={\chi }_{{\sigma }_1{\sigma }_2,{\sigma }_3{\sigma }_4}^{(0)R,A}(\mathbf{q},\omega )$$

$$-\sum _{{\sigma }_i^{\prime }}{\chi }_{{\sigma }_1{\sigma }_2,{\sigma }_3^{\prime }{\sigma }_4^{\prime }}^{(0)R,A}(\mathbf{q},\omega ){\Gamma }_{{\sigma }_3^{\prime }{\sigma }_1^{\prime },{\sigma }_2^{\prime }{\sigma }_4^{\prime }}^{(0)}{\chi }_{{\sigma }_1^{\prime }{\sigma }_2^{\prime },{\sigma }_3{\sigma }_4}^{R,A}(\mathbf{q},\omega ),$$

$${\Gamma }_{{\sigma }_3^{\prime }{\sigma }_1^{\prime },{\sigma }_2^{\prime }{\sigma }_4^{\prime }}^{(0)}=U({\delta }_{{\sigma }_3^{\prime }{\sigma }_4^{\prime }}{\delta }_{{\sigma }_1^{\prime }{\sigma }_2^{\prime }}-{\delta }_{{\sigma }_3^{\prime }{\sigma }_2^{\prime }}{\delta }_{{\sigma }_1^{\prime }{\sigma }_4^{\prime }})$$

$${\chi }_{{\sigma }_1{\sigma }_2,{\sigma }_3{\sigma }_4}^{(0)R}(\mathbf{q},\omega )=\frac{1}{N}\sum _{\mathbf{k}}{\int }_{-\infty }^{\infty }\frac{d\xi }{\pi }\sum _{j,j^{\prime }}{u}_{{\sigma }_4,j}(\mathbf{k})u_{{\sigma }_1,j}^{\ast }(\mathbf{k})u_{{\sigma }_2,j^{\prime }}(\mathbf{k}+\mathbf{q})u_{{\sigma }_3,j^{\prime }}^{\ast }(\mathbf{k}+\mathbf{q})$$

$$\displaystyle\times\Bigl{[}f(\xi){\rm Im}\{G_{j}^{R}({\bf k},\xi)\}G_{j^{\prime}}^{R}({\bf k}+{\bf q},\xi+\omega)+f(\xi+\omega)G_{j}^{A}({\bf k},\xi){\rm Im}\{G_{j^{\prime}}^{R}({\bf k}+{\bf q},\xi+\omega)\}\Bigr{]},$$

Full text

I Introduction

Resonant inelastic x-ray scattering (RIXS) at the transition-metal absorption edges is a promising powerful tool to detect various elementary excitations in strongly correlated electron systems [1, 2]. Particularly, the RIXS technique utilizing the transition-metal $K$- or $L$-edge x-rays can probe momentum dependence of electronic excitations such as charge [3, 4, 5], orbital [6, 7], and magnon excitations [8, 9]. To analyze theoretically those momentum-dependent excitations in RIXS, various effective theoretical methods have been adopted, e.g., exact diagonalization [10, 11, 12], perturbation expansion[13, 14, 15], ultra-short life-time expansion [16, 17], dynamical mean-field theory (DMFT) [18].

In x-ray scattering with linearly polarized x-rays, the total spin moment is conserved, if effects of the spin-orbit coupling are negligible. Therefore, in contrast to neutron scattering, only an even number of magnons can be excited, whereas an odd number of magnons are prohibited to be excited. In fact, excitations with 500 meV energy and characteristic momentum dependence were observed in the Cu $K$-edge RIXS for La2CuO4  [19, 20], and have been identified as two-magnon excitations, based on the agreement with theoretical calculations [21, 22, 23, 24, 25]. In those previous theoretical works they adopted Heisenberg spin Hamiltonians and the spin-wave (SW) approximation.

A large number of theoretical studies on two-magnon excitations have been done in the context of Raman light scattering over several decades [26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37]. The essential microscopic process in most of those studies is the inter-site spin exchange, which can be described effectively by the so-called Fleury-Loudon (FL) Hamiltonian [26]:

[TABLE]

where $\alpha$ is a constant, ${\bf s}_{i}$ is the spin operator at magnetic ion site $i$, ${\bf e}$ and ${\bf e}^{\prime}$ are the electric-field vectors of emitted and absorbed rays, and ${\bf r}_{ij}$ is the coordinate vector connecting sites $i$ and $j$. The FL Hamiltonian can be incorporated into Heisenberg spin Hamiltonians as perturbation with ease. In calculations with the FL Hamiltonian, one can reproduce the observed characteristic lineshapes of two-magnon Raman spectra by taking account of magnon-magnon interactions, which crucially evidenced the importance of magnon-magnon interactions [27, 28, 29, 30, 31]. The FL Hamiltonian has been derived microscopically from a Hubbard Hamiltonian in the large-$U$ limit, being intended for the high-$T_{c}$ cuprates [33, 35].

Two-magnon excitations in RIXS should be distinguished from those in Raman light scattering, since a $K$-shell (i.e., $1s$ core) hole created in the intermediate state can play an essential role. We need to consider excitation processes in the presence of the core hole at a transition-metal site, which are not involved in Raman light scattering. A microscopic mechanism of two-magnon excitations in $K$-edge RIXS was before proposed by van den Brink [22]. He took account of virtual inter-site hopping processes to screen the $1s$ hole in the intermediate state of RIXS, and thereby calculated the modification to the antiferromagnetic Heisenberg spin exchange $J$ around the $1s$ hole. The modified exchange integral between the excited site and neighboring sites induces inter-site spin exchange excitations between those sites, and consequently two magnons are created before the $1s$ hole is finally annihilated. In Ref. 21, Nagao and Igarashi incorporated this mechanism into the previous perturbative framework developed by the author and Igarashi. They replaced a ladder of electron and hole propagators in Ref. 14 with a ladder of magnon propagators, and calculated the dynamical correlation function by using the SW approximation and including the magnon-magnon interaction within $1/S$ expansion. These calculations seem to explain well the experimental observations so far in La2CuO4 [21, 22, 25].

A central aim of the present work is to describe two-magnon excitations in $K$-edge RIXS using an itinerant Hubbard-type model ($d$-$p$ model) and the spin-density-wave (SDW) formalism, instead of using a Heisenberg spin model and the SW approximation. To our knowledge, studies of two-magnon RIXS based on the SDW formalism are still lacking, although some studies of two-magnon Raman light scattering based on the SDW formalism have been done [36]. We deal with the coupling between the $1s$ hole and a pair of magnons within perturbation expansion in the Coulomb interaction, and calculate the magnon propagators within the random-phase approximation (RPA). This approach is another natural way of extending our previous perturbative formulation to explain well the two-magnon RIXS spectra, as we shall demonstrate below.

II Theoretical Framework

II.1 Schematic illustration of two-magnon RIXS and Hamiltonian

Before presenting theoretical details, we schematically illustrate typical $K$-edge RIXS processes involving two-magnon creation: An incident x-ray near the $K$-edge promotes a $1s$ electron resonantly to an empty $p$ band above the Fermi energy ($E_{F}$). Subsequently, correlated $d$ electrons near $E_{F}$ are excited to screen the created $1s$ hole, due to the Coulomb interaction between the $1s$ and $d$ orbitals. Such a screening process can be expressed by excitation of an electron-hole pair, as illustrated in Fig. 1 (i). This electron-hole excitation in the intermediate state can decay into a pair of magnons, due to the spin degree of freedom of the excited electron and hole, as shown in Fig. 1 (ii) and (iii). After two magnons are created, the initially-promoted $p$ electron goes back to the $1s$ state, emitting an x-ray, in the final state. The energy loss and momentum change between the incident and emitted x-rays are carried away by the two magnons.

To be specific, hereafter we restrict our discussion to a typical copper oxide La2CuO4, although our discussion below is applicable also to other transition-metal compounds. To describe the electron dynamics illustrated above, we use the following Hamiltonian:

[TABLE]

This is the same Hamiltonian that we used in Ref. 14. We present each term explicitly in the following.

$H_{d}$ describes the electronic states near $E_{F}\equiv 0$. We take a Hubbard-type Hamiltonian ($d$-$p$ model) for the Cu$3d_{x^{2}-y^{2}}$ and O$2p_{x,y}$ orbitals in a single CuO2 layer:

[TABLE]

Here $d_{i\sigma}$ and $p_{a\ell\sigma}$ ($d_{i\sigma}^{{\dagger}}$ and $p_{a\ell\sigma}^{{\dagger}}$) are the annihilation (creation) operators for the Cu$3d_{x^{2}-y^{2}}$ and O$2p_{\ell}$ electrons with spin $\sigma$, where the Cu$3d_{x^{2}-y^{2}}$ and O$2p_{\ell}$ orbitals form a $\sigma$-bond. $n_{di\sigma}$ is the number operator for the $3d_{x^{2}-y^{2}}$ electrons with spin $\sigma$ at Cu site $i$. Summation with $\langle i,a\rangle$ ($\langle a,b\rangle$) is over nearest-neighbor Cu-O (O-O) bonds. We take $t_{dp}=1.3$ eV, $t_{pp}=0.65$ eV [38], and $U=11$ eV as in our previous work [14].

For the $1s$ electrons, we assume a completely localized orbital at each Cu site:

[TABLE]

where $\varepsilon_{1s}$ is the one-particle energy, $s_{i\sigma}$ ($s_{i\sigma}^{{\dagger}}$) is the annihilation (creation) operator for the $1s$ electrons with spin $\sigma$ at Cu site $i$. $s_{{\bf k}\sigma}(s_{{\bf k}\sigma}^{{\dagger}})$ is the momentum representation of $s_{i\sigma}(s_{i\sigma}^{{\dagger}})$.

$H_{sd}$ is the core-hole potential, i.e., the Coulomb interaction between the $1s$ and $3d$ electrons:

[TABLE]

where $n_{di}$ and $n_{si}$ are the number operator for the $3d_{x^{2}-y^{2}}$ and $1s$ electrons with spin $\sigma$ at Cu site $i$, respectively. Since the $3d$ and inner-shell $1s$ orbitals are both strongly localized, the core-hole potential $V_{sd}$ is usually chosen to be comparable to the Cu$3d$ Coulomb interaction $U$ for La2CuO4.

$H_{p}$ describes the conduction $4p$ band electrons:

[TABLE]

where $p^{\prime}{}_{{\bf k}\mu\sigma}$ ($p^{\prime}{}_{{\bf k}\mu\sigma}^{{\dagger}}$) is the annihilation (creation) operator of the $4p_{\mu}$ electron ($\mu=x,y,z$) with momentum ${\bf k}$ and spin $\sigma$. $H_{x}$ describes the resonant $1s$-$4p$ electric-dipole transition induced by x-rays:

[TABLE]

where $\alpha_{\bf{q}\bf{e}}$ is the annihilation operator of a photon with momentum $\bf{q}$ and polarization $\bf{e}$. The electric-dipole transition matrix $w_{\mu}(\bf{q},\bf{e})$ is given by

[TABLE]

in natural units ($c=\hbar=1$), where $e$ and $m$ are the elementary charge and the electron mass, ${\bf e}_{\mu}$’s are the orthonormal basis vectors.

II.2 Two-magnon RIXS formula

Nozières and Abrahams (NA) developed a theoretical framework of electron Raman scattering by means of Keldysh perturbation theory, and discussed the threshold singularity in metals [39]. Extending the NA’s framework, we derived a formula for RIXS intensity at the transition-metal $K$ edge [13, 14, 15], and thereby we explained experimental observations on charge-transfer and orbital excitations in $K$-edge RIXS for several transition-metal compounds [40, 41, 42]. To calculate two-magnon RIXS spectra, we further extend our previous framework by a way different from Ref. 21, fully based on the itinerant picture.

We adopt the SDW mean-field approach to describe the antiferromagnetic (AF) ground state [43, 44]. Within the SDW mean-field formalism, the Coulomb interaction part $H^{\prime}$ is approximated by

[TABLE]

where ${\bf k}_{\sigma}$ is defined as ${\bf k}_{\uparrow}={\bf k}$ and ${\bf k}_{\downarrow}={\bf k}+{\bf q}_{AF}$ with ${\bf q}_{AF}$ the magnetic-ordering vector. $n_{d}$ and ${\bf m}$ are the mean-fields for the $d$-electron number and spin moment to be determined self-consistently, and ${\bf\sigma}$ is the Pauli matrix vector. In our calculation, we set the $z$-axis of spin along the crystallographic [001] (i.e., the $c$-axis) direction, and assume the commensurate AF ground state with ${\bf q}_{AF}=(\pi,\pi)$ and ${\bf m}\parallel[110]$ as observed in neutron scattering [46]. Introducing new fermion annihilation and creation operators, the mean-field Hamiltonian $H_{d,MF}\equiv H_{0}+H^{\prime}_{MF}$ is diagonalized as:

[TABLE]

where $j$ is band index. $E_{j}({\bf k})$ is the diagonalized band energy which the chemical potential is already subtracted from. The chemical potential is always determined so that the total electron number equals five per unit cell. The original $d$-electron annihilation and creation operators in the momentum representation are related to $a_{j{\bf k}}$ and $a_{j{\bf k}}^{{\dagger}}$ by

[TABLE]

where $u_{\sigma,j}({\bf k}_{\sigma})$’s are the $d$-electron elements of the diagonalization matrix. The $d$-electron Green’s function has a $2\times 2$ matrix form:

[TABLE]

Each element of $\hat{G}({\bf k},\xi)$ is expressed as

[TABLE]

with

[TABLE]

The advanced, retarded and causal branches of the Green’s functions are

[TABLE]

where $\gamma$ is a small positive, and $f(\xi)=[e^{\xi/T}+1]^{-1}$ is the Fermi distribution function at temperature $T$. Using the Green’s function $G^{R}_{j}({\bf k},\xi)$, the self-consistency equations for $n_{d}$ and ${\bf m}$ become

[TABLE]

Throughout our study, numerical integration in energy and momentum is carried out by discretizing the interval $-40$ eV $<\xi<40$ eV into 16000 energy points and the first Brillouin zone into $N=80\times 80$ k points. We take $\gamma=0.16$ eV. We choose $\varepsilon_{d}-\varepsilon_{p}\approx-7{\rm eV}$ so that $n_{d}=1.53$ and $\varepsilon_{d}-\varepsilon_{p}+Un_{d}/2=1.4$ eV for the non-magnetic state. The spin moment and insulating gap calculated for the AF ground state are $|{\bf m}|=0.61\mu_{B}$ and 2.55 eV, respectively. The charge-transfer energy between the Cu-$3d$ and O-$2p$ states is

[TABLE]

and thus the AF Heisenberg exchange between nearest-neighbor Cu sites is evaluated as [45]:

[TABLE]

This value is consistent with that used in most Heisenberg-model calculations for La2CuO4.

To describe the electron dynamics in RIXS, we define the dynamical correlation function by

[TABLE]

where $\rho_{\sigma_{1}\sigma_{2}}({\bf q},t)$ is the density operator in the Heisenberg representation:

[TABLE]

Fourier transform of the dynamical correlation function is related to the linear-response susceptibility in terms of the fluctuation-dissipation theorem:

[TABLE]

where $\chi_{\sigma^{\prime}_{1}\sigma^{\prime}_{2},\sigma_{2}\sigma_{1}}^{R,A}({\bf q},\omega)$ are the retarded and advanced branches of the susceptibility, and are the Fourier transforms of

[TABLE]

We calculate $\chi_{\sigma_{1}\sigma_{2},\sigma_{3}\sigma_{4}}^{R,A}({\bf q},\omega)$ within RPA:

[TABLE]

where

[TABLE]

is the Coulomb interaction at Cu site. Diagrammatic expression of RPA is presented in Fig. 2(a). $\chi_{\sigma_{1}\sigma_{2},\sigma_{3}\sigma_{4}}^{(0)}({\bf q},\omega)$ is the bare susceptibility calculated by

[TABLE]

For numerical calculations, we use more convenient expressions:

[TABLE]

where

[TABLE]

Energy-momentum integration in Eq. (37) is numerically performed using fast Fourier transformation (FFT).

Below we need also the causal susceptibility, $\chi_{\sigma^{\prime}_{1}\sigma^{\prime}_{2},\sigma_{2}\sigma_{1}}^{C}({\bf q},\omega)$, which is the Fourier transform of

[TABLE]

Here T$[...]$ means the ordinary time-ordered product. This causal component can be calculated from the advanced and retarded components:

[TABLE]

where $b(\omega)=[e^{\omega/T}-1]^{-1}$ is the Bose distribution function.

Single-magnon spectrum to be observed in neutron scattering is not $\Pi_{\sigma^{\prime}_{1}\sigma^{\prime}_{2},\sigma_{2}\sigma_{1}}({\bf q},\omega)$ itself, but

[TABLE]

where $\sigma_{\mu}$ and $\sigma_{\nu}$ are the Pauli matrices, and ${\bf q}_{\sigma_{1}\sigma_{2}}\equiv{\bf k}_{\sigma_{1}}-{\bf k}_{\sigma_{2}}$. Summation in spin is restricted to the cases when ${\bf q}_{\sigma^{\prime}_{1}\sigma^{\prime}_{2}}$ and ${\bf q}_{\sigma_{1}\sigma_{2}}$ are equivalent up to reciprocal-lattice translations. This restriction is necessary, because neutron scattering experiments observe the components whose absorbed and emitted momenta are equivalent [47]. We should note that $\Pi_{\sigma^{\prime}_{1}\sigma^{\prime}_{2},\sigma_{2}\sigma_{1}}({\bf q},\omega)$ describes the dynamical processes where the absorbed and emitted momenta are ${\bf q}-{\bf q}_{\sigma_{1}\sigma_{2}}$ and ${\bf q}-{\bf q}_{\sigma^{\prime}_{1}\sigma^{\prime}_{2}}$, respectively, as illustrated in Fig. 2(b), where we express diagrammatically $\Pi_{\sigma^{\prime}_{1}\sigma^{\prime}_{2},\sigma_{2}\sigma_{1}}({\bf q},t^{\prime}-t)$ as well as $\chi_{\sigma^{\prime}_{1}\sigma^{\prime}_{2},\sigma_{2}\sigma_{1}}^{R,A}({\bf q},\omega)$ by a solid wavy line.

RIXS intensity from two-magnon excitations can be calculated in a similar way to that in Ref. 14 by means of the Keldysh diagrammatic technique [39]. In the Keldysh diagrammatic representation, probability of excitations where a pair of magnons are left in the final state is expressed by diagrams where a pair of magnon propagators bridge the upper normally and lower reversely time-ordered branches. As in our previous works, we adopt the Born approximation, i.e., the lowest-order perturbation with respect to the core-hole potential $V_{sd}$. Typical diagrams, which we shall calculate, are displayed in Fig. 2(c), where a $3d$ triangle loop connects three RPA propagators, giving the lowest-order coupling between a charge mode (c) and two magnons (m1 and m2). These diagrams give the transition probability of the physical processes of Fig. 1, where the charge mode excited by the $1s$ hole decays into a pair of magnons in the final state. In Fig. 2(c), we express $\chi_{\sigma_{1}\sigma_{2},\sigma_{3}\sigma_{4}}^{C}({\bf q},\omega)$ as well as $\Pi_{\sigma_{1}\sigma_{2},\sigma_{3}\sigma_{4}}({\bf q},\omega)$ by a solid wavy line. As known from Fig. 2(c), RIXS intensity from two-magnon excitations consists of two contributions:

[TABLE]

where the left [right] diagram including a parallel [crossed] pair of magnon propagators expresses $W_{p}(q{\bf e};q^{\prime}{\bf e}^{\prime})$ [$W_{c}(q{\bf e};q^{\prime}{\bf e}^{\prime})$]. $q=({\bf q},\omega)$ and ${\bf e}$ [$q^{\prime}=({\bf q}^{\prime},\omega^{\prime})$ and ${\bf e}^{\prime}$] are the four-momentum and polarization of the absorbed [emitted] x-ray. Analytic expressions of these diagrams are given in the form of the convolution of the dynamical correlation function:

[TABLE]

where $Q\equiv({\bf Q},\Omega)\equiv q-q^{\prime}\equiv({\bf q}-{\bf q}^{\prime},\omega-\omega^{\prime})$ is the momentum transfer and energy loss of x-rays, $p=({\bf p},\zeta)$, $V_{\sigma_{1}\sigma_{2},\sigma_{3}\sigma_{4}}({\bf e},{\bf e}^{\prime};\omega;Q;p)$ is the scattering function expressed by the product of the $1s$-$4p$ triangle loop (the resonance factor), the core-hole potential screened by the charge mode c, and the $3d$ triangle loop:

[TABLE]

where

[TABLE]

with

[TABLE]

$V_{sd}\Lambda_{\tau\tau}({\bf Q},\Omega)$ corresponds to the screened core-hole potential, which is reduced from the bare $V_{sd}$ due to electron-hole excitations $\chi_{\tau_{1}\tau_{2},\sigma\sigma}^{C}({\bf Q},\Omega)$. Energy-momentum integration in Eq. (47) is numerically performed using FFT. In terms of the density of states (DOS) of the $4p_{\mu}$ band,

[TABLE]

the resonance factor is written as

[TABLE]

For the below numerical calculations for RIXS spectra, we use $\rho_{4p\mu}(\varepsilon)$ obtained from a first-principles band calculation and $\varepsilon_{1s}=-8980$ eV, which reproduce well experimental x-ray absorption spectra (XAS) as shown in the Appendix. Rigorously speaking, the expressions (42) and (43) include not only two-magnon excitations but also purely charge excitations. However, as far as we see low-energy excitations below 1 eV, we may regard the two-magnon excitations completely dominate the spectral weights, since charge modes below the insulating gap ($\sim 2.5$ eV) should be absent in the final state.

The above expression of RIXS intensity contains not only magnon propagators, but the resonance factor and the $3d$ triangle loop. To see only the the contribution from two-magnon excitations, we define the bare part of two-magnon excitations:

[TABLE]

which is diagrammatically represented in Fig. 2(d). $\tilde{W}({\bf Q};\Omega)$ is obtained from Eq. (42) by maintaining total-spin conservation but neglecting the momentum-frequency and polarization dependences in $V_{\sigma_{1}\sigma_{2},\sigma_{3}\sigma_{4}}({\bf e},{\bf e}^{\prime};\omega;Q;p)$.

Finally, we define the integrated weights:

[TABLE]

where we choose $E_{c}=1$ eV, and the scaling constants $C$ and $\tilde{C}$ to satisfy $I({\bf e},{\bf e}^{\prime};{\bf Q})=\tilde{I}({\bf Q})=1$ at ${\bf Q}=(\pi,0)$, for the following numerical calculations.

III Numerical Results

III.1 Single-magnon excitations

Consistency between calculated single-magnon spectra and neutron scattering data is prerequisite for calculation of two-magnon RIXS spectra. Numerical results for the dynamical spin correlation function $S_{zz}({\bf q},\omega)$ are displayed in Fig. 3(a). As shown in Ref. 44, the SDW approach yields a spin-wave (i.e., single-magnon) dispersion relation precisely agreeing with neutron scattering. Although we use the $d$-$p$ model which includes the O-$2p$ orbitals and differs from the simple Hubbard model in Ref. 44, the dispersion is well reproduced again to agree with neutron scattering data. It is known that SW calculation in the simple nearest-neighbor Heisenberg model leads to $\omega(\pi,0)=\omega(\pi/2,\pi/2)$, being inconsistent with neutron scattering. An advantage of the SDW formalism is that the magnon excitation energies are correctly reproduced at those q points. One of notable features is that the strong divergent behavior around ${\bf q}=(\pi,\pi)$, as observed in neutron scattering [48], which is clearly seen in terms of the integrated intensity in Fig. 3(b). This divergent behavior becomes important for the two-magnon RIXS intensity around the zone center, as we shall see below.

In SW approaches, magnon DOS can be calculated straightforwardly, since the dispersion relation of magnons is expressed explicitly in terms of trigonometric functions. On the other hand, in our SDW approach, magnon DOS is difficult to calculate precisely, since the magnon energies are obtained only as numerical values and the calculated magnon peaks in $S_{\mu\nu}({\bf q},\omega)$ are accompanied by some broadness (i.e., damping).

III.2 Dependence on x-ray momentum transfer

Calculated two-magnon RIXS spectra $W(q{\bf e};q^{\prime}{\bf e}^{\prime})$ are compared with experimental data at various momentum transfers in Fig. 4. In Fig. 4, also the bare two-magnon spectra $\tilde{W}({\bf Q},\Omega)$ are drawn. In Fig. 4(e), vanishing of the two-magnon RIXS spectrum at ${\bf Q}=(\pi,\pi)$ is reproduced, being consistent with experiments and also with other previous theoretical calculations. Considering uncertainty due to elastic-line subtraction and limited resolution in the experiment, we can regard the theoretical curves as agree well with the experimental plots.

To see more detailed momentum dependence of the two-magnon RIXS spectrum, we present the intensity map along symmetry lines in the first Brillouin zone in Fig. 5(a), where we compare with experimental peak positions. Although the calculated intensity map shows broadness of spectra, the experimental peak positions fall well within an energy-loss range where the calculated intensity is relatively strong, except for ${\bf Q}=(\pi,\pi)$.

A notable feature is that the RIXS intensity becomes strong around ${\bf Q}=0$ and $\Omega=0$, which has not been pointed out in previous experiments and theoretical calculations. This notable feature is straightforwardly understood in the following way: Intensity of single-magnon excitations becomes divergent toward ${\bf q}=(\pi,\pi)$ as observed in neutron scattering and also as calculated in III.1. Therefore, excitation of two magnons with ${\bf p}_{1},{\bf p}_{2}\approx(\pi,\pi)$ contributes strongly to the two-magnon spectrum around ${\bf Q}={\bf p}_{1}+{\bf p}_{2}\approx(2\pi,2\pi)\equiv(0,0)$. This is a striking difference from the previous two-magnon RIXS calculations. The difference arises from the difference in starting models, rather than from that in the underlying microscopic mechanisms of two-magnon creation. As shown in Ref. 22, the magnetic scattering operator involving two spins is proportional to or commutable with the starting Heisenberg spin Hamiltonian at ${\bf Q}=0$, and therefore can never excite the ground eigenstate to any other states. This peculiarity of Heisenberg spin models leads to the vanishing of RIXS intensity at ${\bf Q}=0$. On the other hand, such vanishing does not occur in Hubbard-type Hamiltonians.

In Fig. 5(b), we compare the integrated weights $I({\bf e},{\bf e}^{\prime};{\bf Q})$ and $\tilde{I}({\bf Q})$ with experimental data. Integration of the intensities yields a large weight near the zone center, which contradicts the vanishing weight concluded in previous theoretical and experimental works. In Heisenberg spin models, the intensity vanishes at the zone center and therefore also the integrated weight vanishes there, unless ring spin exchange is taken into account [21, 22, 25]. As we pointed out already, this vanishing arises from the peculiarity of Heisenberg spin models. In experiment, the RIXS intensity around the zone center is overlaid by the much stronger elastic line, and is difficult to distinguish from it. In Ref. 20, the two-magnon RIXS weight is set to zero at the zone center after background subtraction. However, supposing that two-magnon RIXS intensity may be subtracted together with the elastic line in analyzing experimental data, still one cannot exclude the possibility that the two-magnon RIXS weight takes a finite value at the zone center.

In Fig. 5(b), the calculated weight at ${\bf Q}=(\pi,\pi)$ seems finite, which will be inaccurate due to the limited precision of numerical integration. One may expect that the weight approaches zero more closely, if we take smaller $\gamma$ and finer energy-momentum discretization for numerical integration.

In Fig. 5(c), the intensity map of the bare part of two-magnon excitations, $\tilde{W}({\bf Q},\Omega)$, is presented along symmetry lines in the first Brillouin zone. From comparison between $W(q{\bf e};q^{\prime}{\bf e}^{\prime})$ and $\tilde{W}({\bf Q},\Omega)$ in Fig. 4, Figs. 5(a) and (c), the bare two-magnon part $\tilde{W}({\bf Q},\Omega)$ well captures the overall spectral properties of RIXS spectra $W(q{\bf e};q^{\prime}{\bf e}^{\prime})$, except around the zone center ${\bf Q}=0$.

III.3 Dependence on incident x-ray energy and polarization

In experiments [19, 20], the two-magnon RIXS weight was not observed for the polarization geometry $E\parallel ab$. To explain this, we calculate the dependence of two-magnon RIXS spectra on x-ray polarization direction, whose results are shown in Fig. 6. According to the calculated results, the two-magnon RIXS intensity for $E\parallel ab$ will be enhanced at higher incident x-ray energies (by about 4 eV) than that for $E\parallel c$. However, the intensity for $E\parallel ab$ is much weaker even at the resonance than for $E\parallel c$, which explains the missing of this feature for $E\parallel ab$ in experiments. This weakness of the intensity for $E\parallel ab$ arises from the smallness of the resonance factor, and therefore, roughly speaking, can be attributed to the smallness and broadness of the $4p$ DOS $\rho_{4px,y}(\varepsilon)$ (see Fig. 7). More careful investigations may confirm a small but finite two-magnon weight for $E\parallel ab$.

IV Discussion and Summary

Comparing the RIXS intensity $W(q{\bf e};q^{\prime}{\bf e}^{\prime})$ and bare two-magnon part $\tilde{W}({\bf Q},\Omega)$ is quite illuminating. From Fig. 4, calculated $W(q{\bf e};q^{\prime}{\bf e}^{\prime})$ and $\tilde{W}({\bf Q},\Omega)$ exhibit quite similar spectral properties except around the zone center ${\bf Q}=0$. This similarity suggests that, except around ${\bf Q}=0$, the spectral properties are determined almost only by the two-magnon states in the final state, and do not depend on the details of the momentum-frequency dependences of scattering function $V_{\sigma_{1}\sigma_{2},\sigma_{3}\sigma_{4}}({\bf e},{\bf e}^{\prime};\omega;Q;p)$. On the other hand, difference around ${\bf Q}=0$ is striking, which is more clearly seen in the integrated weights $I({\bf e},{\bf e}^{\prime};{\bf Q})$ and $\tilde{I}({\bf Q})$ in Fig. 5(b): $I({\bf e},{\bf e}^{\prime};{\bf Q})$ is much suppressed around ${\bf Q}=0$. This suppression around ${\bf Q}=0$ arises from the momentum dependence of $\Lambda_{\tau\tau^{\prime}}({\bf Q},\Omega)$. As we mentioned, $\Lambda_{\tau\tau^{\prime}}({\bf Q},\Omega)$ reflects the strength of core-hole screening due to the electron-hole excitation mode with momentum Q. We show $|\Lambda_{\tau\tau}({\bf Q},\Omega)|$ in Fig. 5(d), where the smallness of $|\Lambda_{\tau\tau}({\bf Q},\Omega)|$ around ${\bf Q}=0$ indicates that the core-hole potential is well screened to be effectively weak there. Such an effect arising from local core-hole screening is characteristic of RIXS, and absent in Raman light scattering.

As we explained already, the strong intensity around ${\bf Q}=0$ is not obtained from Heisenberg spin models. Furthermore, we point out that this strong intensity around ${\bf Q}=0$ cannot be captured by the FL framework on Raman light scattering. The FL Hamiltonian Eq. (1) leads to a momentum dependent scattering vertex as a result from the nature of inter-site interactions, and fails to pick up the strong contribution from the two magnons ${\bf p}_{1},{\bf p}_{2}\approx(\pi,\pi)$. In fact, the dominant B1g Raman spectrum is calculated by integrating a function including a factor $\sim\cos p_{x}-\cos p_{y}$ [35, 37], which makes the above two magnons ${\bf p}_{1},{\bf p}_{2}$ much less important.

In Raman light scattering, characteristic spectral lineshapes confirmed the importance of magnon-magnon interactions [27, 28, 29, 30, 31]. Raman light scattering spectrum can be related to the ${\bf Q}=0$ weight, since both the absorbed and emitted rays take a negligibly small wavenumber. Therefore, as in Raman light scattering, magnon-magnon interactions may affect the RIXS intensity at the zone center ${\bf Q}=0$. However, in our calculation on RIXS, magnon-magnon interactions are not included. In our framework, magnon-magnon interactions are expressed by diagrams in which the two wavy lines, m1 and m2, in Fig. 2(c) are connected by some diagram elements. Unfortunately, magnon-magnon interactions are difficult to include with some closed form and feasible numerical computation. On the other hand, it is still controversial whether magnon-magnon interactions are important for ${\bf Q}\neq 0$ or not. One reason why magnon-magnon interactions may not be effective for ${\bf Q}\neq 0$ is that two magnons with finite ${\bf Q}$ have relative non-zero velocity and rapidly go away from each other, and therefore the interaction between them is effectively diminished [20, 23].

In summary, we have calculated the two-magnon RIXS spectra at the $K$-edge, using an itinerant Hubbard-type Hamiltonian and the SDW mean-field formalism. Single-magnon excitation has been described within RPA. Coupling between the $K$-shell hole and the magnons in the intermediate state has been dealt with by means of diagrammatic perturbation expansion in the Coulomb interaction. The calculated momentum and polarization dependences of two-magnon RIXS spectra agree well with those of experiments. A sharp contrast to previous studies based on Heisenberg spin models and the SW approximation is that the two-magnon excitations can take a large weight near the zone center. Further high-resolution experiments around the zone center may provide us with insights on: which of the itinerant Hubbard-type models and the Heisenberg spin models are more appropriate for the description of two-magnon RIXS near the zone center, how strongly the core-hole potential is screened, as well as how effective the magnon-magnon interactions are.

Acknowledgements.

The author would like to thank Prof. J. Igarashi, Dr. K. Ishii and Prof. T. Nagao for invaluable communications.

Appendix A X-ray absorption spectra

X-ray absorption spectra (XAS) at the $1s$-$4p$ resonance can be approximately calculated from the conduction $4p$-band DOS $\rho_{4p\mu}(\varepsilon)$ with neglecting core-hole bound states:

[TABLE]

where $q=({\bf q},\omega)$, and $\Gamma_{1s}=1$ eV in our numerical calculation. We calculate $\rho_{4p\mu}(\varepsilon)$ using the WIEN2k code [49]. In Fig. 7, calculated XAS are compared with the experimental data read from Ref. 50 for two polarization directions $E\parallel ab$ and $E\parallel c$. For the both cases of polarization, the main peak positions agree well with the experimental ones by setting $\varepsilon_{1s}=-8980$ eV.