Using RIXS to uncover elementary charge and spin excitations in correlated materials

Despite significant progress in resonant inelastic x-ray scattering (RIXS) experiments on cuprates at the Cu L-edge, a theoretical understanding of the cross-section remains incomplete in terms of elementary excitations and the connection to both charge and spin structure factors. Here we use state-of-the-art, unbiased numerical calculations to study the low energy excitations probed by RIXS in undoped and doped Hubbard model relevant to the cuprates. The results highlight the importance of scattering geometry, in particular both the incident and scattered x-ray photon polarization, and demonstrate that on a qualitative level the RIXS spectral shape in the cross-polarized channel approximates that of the spin dynamical structure factor. However, in the parallel-polarized channel the complexity of the RIXS process beyond a simple two-particle response complicates the analysis, and demonstrates that approximations and expansions which attempt to relate RIXS to less complex correlation functions can not reproduce the full diversity of RIXS spectral features.

Resonant inelastic x-ray scattering (RIXS) is an experimental technique in which the transferred energy, momentum and polarization associated with incident and scattered x-ray photons can be measured and analyzed to reveal information about the elementary excitations of a system [1]. In recent years RIXS has attracted considerable attention and positioned itself as a primary experimental technique to probe the excitations in correlated materials, especially transition-metal oxides [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17]. RIXS possesses atomic sensitivity with incoming photons resonantly tuned to a specific atomic absorption edge, making it a particularly unique and powerful tool for characterizing excitations across the Brillouin zone in these materials.
The direct RIXS process consists of two dipole transitions, as shown schematically in Fig. 1 for the Cu L-edge. In the first step an incoming photon excites the ground state by promoting an electron from a filled core shell (Cu 2p) into the valence shell (Cu 3d). An intermediate state manifold forms following the charge and spin shakeup which accompany the introduction of a local core-hole potential and this new carrier in the valence shell. In the second step an electron from the valence shell, possessing appropriate atomic character, fills the core-hole, accompanied by an outgoing photon, which leaves the system in an excited final state.
Due to the complexity of this process and the presence of the intermediate state, an interpretation of RIXS spectra has been hindered by an incomplete understanding of how it may be related to other, more fundamental, response functions governing spin, charge, lattice and orbital excitations. In addition, it is well known that polarization plays a key role in understanding selective excitations in Raman spectroscopy [18], but a full polarization analysis involving both incident and scattered x-rays has become possible only now in experiment and only shown theoretically to be important for spin-flip excitations [19]. One would expect that a more complete understanding, at a fundamental level, could be obtained by analyzing the full theoretical RIXS cross-section, accounting for the influence of both the incident and scattered light polarization for a given experimental scattering geometry.
To better understand the RIXS cross-section from a theoretical perspective, here we will explore a rather general question: what low energy excitations are measured by RIXS, particularly in cuprates at the Cu L-edge? There certainly are a number of ways to address this question analytically and numerically. Here, our interest is confined to the low energy excitations which involve spin and/or charge degrees of freedom, neglecting both lattice and higher energy orbital or charge transfer excitations in this work. Our aim will be to qualitatively understand which excitations are encoded in the RIXS cross-section, and any connections to both spin and charge dynamical structure factors, under different scattering conditions; and we will quantitatively compare numerical estimates of the full RIXS cross-section with various approximations to tease out this information.
Naively, one may expect the Cu L-edge RIXS crosssection to be the resonant analog of the hard x-ray nonresonant IXS, [20,21] whose response is the charge dynamical structure factor. It then may seem counterintuitive that RIXS at the Cu L-edge rose to promi- nence for the successful empirical measurement of the spin (magnon) dispersion in undoped cuprates, making it a complementary experimental probe to the wellestablished inelastic neutron scattering [4,5]. With theoretical support [19,22], experiments on doped cuprates proved sensitive to magnetic excitations, with a similar cross-section to the spin dynamical structure factor, regardless of doping level [7][8][9][10][11][12][13][14][15][16][17]. These observations were interpreted in terms of persistent magnetic excitations up to an unexpectedly high doping level. However, this semi-empirical connection to the spin dynamical structure factor has been based primarily on approximate theoretical and numerical treatments for the full RIXS crosssection. The fast collision approximation [23][24][25] and the effective operator approach [26] suggested that only single magnon excitations or S(q, ω) should be measured by RIXS at the Cu L-edge [22,26,27]. More sophisticated treatments -the ultrashort core-hole lifetime (UCL) expansion [28][29][30] and the UCL-inspired ansatz [31] -highlighted that RIXS should be sensitive also to bimagnon excitations. A further extension of UCL also pointed out the importance of three-magnon excitations [32]. However, these studies made no explicit comparison between the approximations and the full RIXS cross-section, nor have they concentrated on the sensitivity of RIXS to the charge dynamical structure factor. Perhaps more importantly, these approximations suffer from several severe limitations. First, the fast collision approximation (or the effective operator approach) [22,26,27] relies on an estimation of the dynamics only at the site where the core-hole has been created in the intermediate state. However, the intermediate state is drawn from a manifold of states, which differ from the ground state not only locally, at the site where a core-hole is created, but also on neighboring sites due to hopping or spin exchange associated with the dynamical screening process. As a result, and especially upon doping, this approximation fails to capture key elements of the full RIXS process. Second, the UCL approximation relies on the assumption that the energies of the intermediate state manifold are much smaller than the inverse corehole lifetime Γ, which allows an approximation based on only the first few (two) terms in the UCL Taylor-serieslike expansion. However, this assumption need not hold, especially in "itinerant", doped systems where a number of intermediate states may have energies ∝ t ∼ Γ, making the UCL a non-convergent approximation.
In this paper we study the low energy excitations of RIXS at the Cu L-edge in an unambiguous way by numerically evaluating the cross-section for a twodimensional Hubbard model comprising "effective" Cu 3d orbitals supplemented by local Cu 2p core levels. Using exact diagonalization, which previously was applied to the study of paramagnons in cuprates [19], we compare the exact RIXS cross-sections to approximate crosssections for the same model using the same method. The intrinsically correlated nature of our Hamiltonian distinguishes this study from the analytically exact RIXS calculations for the case of completely uncorrelated electrons [33]. In the next section, we present and compare numerical results for the exact and approximate RIXS cross-sections, and discuss content of the excitations and consequences for experiments. The paper ends with conclusions, and appendices which contain details and longer derivations.

A. Exact RIXS cross-section
The RIXS cross-section at the Cu L-edge is [1,27] (1) where |i (|f ) is the initial (final) state of the system in the RIXS process with energy E i (E f ), transfered mo- where q i and q f (ω i and ω f ) are the incoming and outgoing photon momentum (energy), and e = e i · (e f ) † is the tensor that describes the incoming (i) and outgoing (f ) photon polarizations. Here the operator O q,e = 1/ √ N j e iq·j O j,e and which describes the evolution of the system in the RIXS experiment from the initial state to the final state via the intermediate states accessible via the core-hole to valence band dipole transitions. N is the number of lattice sites in the system. The dipole transition operator D j,e = σ,α,β (A e α p † jασ d jσ + h.c.) with p jασ (d jσ ) annihilates a hole in the 2p (3d) shell with spin σ.A e α is the matrix element of the dipole transition between the 2p α orbitals and 3d x 2 −y 2 orbital written as A e α = d x 2 −y 2 ,σ |ˆ ·r|p ασ for polarizationˆ . Γ is the inverse core-hole lifetime. Note that with the exception of the schematic in Fig. 1, "hole notation" has been used throughout the paper, such that the dipole transitions at the Cu L-edge correspond to transitions from the initial 3d * 1 2p 0 state to a 3d * 0 2p 1 intermediate state, and finally from the intermediate state to the final state configuration. Here 3d * corresponds to holes nominally (Zhang-Rice singlets [34]) in the single band notation.
Here we use the single-band Hubbard model as it carries the key features of correlated materials in the low energy regime, which are the relevant energy regime for the study of charge and spin excitations for cuprates in RIXS measurements. The single-band results can be applied directly to cuprates and can be generalized to other multi-orbital correlated materials. The Hamiltonian H defined on a 2D square lattice describes the relevant interactions of the "effective" 3d and 2p orbitals, consisting of two parts: The first part, Eq. (4), is the well-known single-band Hubbard model with the nearest (next-nearest) neighbor hopping t (t ) and on-site Hubbard repulsion U . Here, the operator d † in the single-band Hubbard model creates a "3d * " hole, which should be distinguished from the actual Cu 3d x 2 −y 2 hole in the multiband model. The second part, Eq. (5), describes (i) the energy splitting between the 2p and the 3d shells through the difference in site energy d − p , (ii) the repulsion between the 2p and the 3d holes (the "core-hole potential") U c , and (iii) the spin-orbit coupling in the 2p shell λ with the matrix elements χ σσ αα ≡ p ασ |l · s|p α σ , where the l · s term represents the spin-orbit coupling operator.
The RIXS cross-section is calculated using exact diagonalization on a 12-site cluster. In this cluster, momentum points (2π/3,0) and (π/2, π/2) are accessible, providing information along both the Brillouin zone axis and diagonal, relevant for RIXS experiments. As the Cu 2p core levels and the spin-orbit coupling are included (which also means that the total spin is not a good quantum number in the intermediate state), the Hilbert space size is ∼ 10 7 . Ground state eigenvectors and eigenvalues were obtained using the implicitly restarted Arnoldi method encoded in the Parallel ARPACK [35] libraries. The cross-section itself was obtained using the biconjugate gradient stabilized method [36] and continued fraction expansion [37]. The numerical technique has been applied previously to calculate RIXS at the Cu K-edge and L-edges [19,20]. Numerical results were obtained for parameters which can relatively well reproduce the low energy physics for cuprates: d − p = 2325t (which gives the typical splitting between the Cu 2p and the Cu 3d shell of 930 For doped systems, we investigate at the resonance where the character of the intermediate state is similar to that in the undoped case on the core-hole site. The angle between the incident and the scattered photons is set to 50 • , with the scattering plane parallel to xz, i.e. perpendicular to the xy plane on which we define the 2D Hubbard Hamiltonian, and the incoming polarization is chosen to be π. This scattering geometry is consistent with that used most commonly in RIXS measurments for cuprates [8][9][10][12][13][14][15]. The relation between polarization and the transferred momenta [22,25,26,38] follows from this scattering geometry: e i = (sin θ, 0, cos θ), , and the angle θ is related to the transferred momentum via k x = 1.07π sin(θ − 65 • ). The calculated RIXS spectra will be presented at Figs. 2 -5. The comparison with the approximated spectra will be presented in Section II C below.

B. Approximate cross-section
The approximations follow from integrating out the core hole degrees of freedom using two approaches. We firstly assume that the energy of the incoming photon ω i is tuned to the main resonance at the Cu L 3 edge, as dictated by H c . [39] This is equivalent to a projection into a subspace with only one 3d * hole on the core-hole site in either the initial state |i or the final state |f . We then perform the UCL approximation by expanding the RIXS operator, Eq. (2), in a power series in H/Γ and only keeping the first two terms. See Appendix A for FIG. 2. Exact RIXS cross-sections and approximations calculated using exact diagonalization. Top (bottom) panels show spectra for the π − σ (π − π) polarization geometries. Left, middle, and right panels show RIXS spectra calculated for n = 0.83 ("hole-doping"), n = 1 ("half-filling"), and n = 1.17 ("electron doping") electron fillings, respectively. The spectra are summed over all momenta along the nodal direction and along the ky = 0 direction that are accessible in RIXS on a 12-site cluster, weighted with the RIXS form factors, i.e.: S(ω) = q |Wπ−σ| 2 S(q, ω),Ñ (ω) = q |Wπ−π| 2Ñ (q, ω) and similarly for S (ω) and N (ω). The elastic response has been removed in each panel.
details. These approximations can be written separately for the two polarization conditions as follows: π−σ RIXS in the UCL approximation: where the local RIXS form factor W π−σ = −ı2 sin θ/(3Γ), the spin operator is the 2D coordination number, and the 2D structure factor is ε k = γ k + t η k /t with γ k = (cos k x + cos k y )/2 and η k = cos k x cos k y . The structure factor ε k has A 1g symmetry. Note that the first term of the expansion S(q, ω) has the form of the spin dynamical structure factor, while the second term S (q, ω) is a rather complicated four particle response which probes spin and charge excitations. The latter corresponds toD  Fig. 2 RIXS spectra. The π − σ (π − π) relative polarizations are shown on the left (right) panels, respectively. Note that the spectral weightÑ = 5×10 −5 (eV ) −2 in the half-filled case and does not appear in the figure.All momenta along the nodal direction and along the qy = 0 direction accessible on a 12-site cluster are taken into account (elastic response excluded). the 3-spin Greens function (see Appendix B) and to joint spin and charge excitations.
π−π RIXS in UCL approximation: where the constrained density operator Before evaluating the correlation functions, we want to highlight that none of these approximations give the standard charge dynamical structure factors for the parallel-polarization channel. The first term of the expansion,Ñ (q, ω), is not the standard charge dynamical structure factor. It represents a more complicated fourparticle response function in the original space of unprojected fermions. In the following, it carries information about projected charge excitations, which should in no way be confused with information about the full charge response. That RIXS does not probe the standard charge response stems from the fact that initial states with double occupancy on the core-hole site can not be excited in the RIXS process due to the Pauli principle. The second termÑ (q, ω), a complicated four-particle response function as well, corresponds to 2-spin or bimagnon excitations (see Appendix B) and other symmetry projected charge excitations, represented again by a correlator beyond the familiar two-particle charge response.

C. Comparing exact and approximate results
In the following, we will present a systematic comparison of the full RIXS spectra with the approximations using the UCL expansion for the two polarization conditions. As a momentum resolved technique, RIXS has the power to measure the dispersion of elementary excitations, which is one of its main advantages compared to traditional optical or Raman scattering, where the momentum transfer is limited to q ∼ 0. [18] However, there are limitations in cluster size, as well as a limited number of poles using a finite-size cluster. Thus, summing over all the accessible momentum points (results as shown in Fig. 2) gives us a complete picture of the energies and the distribution of intensities for the excitations, so that the comparisons to the approximations can be made in a single shot. To better quantify this comparison, we also calculate and compare the total spectral weight carried by the excitations, cf. Fig. 3. Nevertheless, in section II D, we compare exact RIXS cross-sections and approxima-tions at the momentum points (2π/3,0) and (π/2, π/2), and discuss connections to experiments.
The full cross-sections are shown in Fig. 2 for two polarization channels: (i) the cross-polarized channel [π − σ with π (σ) incoming (outgoing) polarization] and (ii) the parallel-polarized channel [π − π with π (π) incoming (outgoing) polarization]. For each the spectra are calculated for three different doping levels (n = 0.83, "hole-doping", n = 1, "half-filling", and n = 1.17, "electron doping"). The π − σ RIXS spectra presented here agree with those presented in Ref. [19]. The approximate cross-sections S(ω) and S (ω) for cross-polarization and N (ω) andÑ (ω) for parallel-polarization are also shown in Fig. 2. The results for both full and approximated RIXS are shown for a Lorentzian broadening with half width at half maximum (HWHM) = 0.025t for the energy transfer. Note that the spectra in Fig. 2 correspond to a momentum summation over all points accessible in the 12-site cluster, to provide a holistic picture of the character of excitations probed by RIXS and the utility of various approximations (see Fig. 2

caption for details).
First note the results in the cross-polarized channel (the π −σ channel). On a qualitative level, the line shape of the full RIXS cross-section can be reproduced well by the spin dynamical structure factor S(q, ω) (the first term of the UCL approximation). This is true at halffilling, where all charge excitations have been gappedout, while in either the electron-or hole-doped cases one can observe some relatively small discrepancies between the two spectra. When adding higher order terms from the effective expansion, to a large extent this observation remains unchanged since these terms encode similar excitations to S(q, ω) together with excitations of mixed charge and spin character, as can be seen readily from the form of the operator in Eq. (8). However, on a quantitative level, this comparison breaks down, with discrepancies in the overall intensity and integrated spectral weight which can become very large (see Fig. 3). While one may have expected that higher order terms in the effective expansion should provide a more satisfactory qualitative and quantitative agreement, these do not help in reducing the differences; on the contrary, these additional terms actually enhance the quantitative mismatch. Note that the first two terms of the UCL expansion suggest larger spectral weight for the hole-doped case compared to that with electron-doping, in contrast to the behavior for the RIXS cross-section. This suggests that the differences cannot be attributed to a simple rescaling factor.
The quantitative mismatch between RIXS at the Cu L-edge and the approximations highlights the role that the intermediate state wavefunction plays in the RIXS spectra, just as in the case for RIXS at the Cu Kedge. [20] RIXS is an intrinsic four-particle process, where the wavefunction overlap between the ground state and Cross-sections at q = (2π/3, 0) with the elastic response removed.Top (bottom) panels show spectra for the π−σ (π−π) polarization. Left, middle, and right panels show RIXS spectra for n = 0.83, n = 1, and n = 1.17 electron fillings, respectively. The approximate spectra are weighted with the RIXS form factors, i.e.S(q, ω) → |Wπ−σ| 2 S(q, ω),Ñ (q, ω) → |Wπ−π| 2Ñ (q, ω) and similarly for S (q, ω) andÑ (q, ω). mation on the fundamental excitation energies, but unfortunately it cannot provide reasonable spectral weights on a quantitative level. This mismatch also illustrates the failure of these approximations. That being said, both numerically and empirically, in the cross-polarized scattering geometry the RIXS cross-section qualitatively corresponds to the spin dynamical structure factor which encodes information about spin excitations at the two particle level, which underscores RIXS utility as a complementary probe to inelastic neutron scattering.
An altogether different situation arises in the parallelpolarized channel (the π − π channel). While a comparison between the RIXS spectrum and the "projected" charge excitations produces a modest qualitative agreement between the line shapes, both missing peaks as well as significant differences in the spectral weights undermine any quantitative value in this comparison. The addition of the higher order terms seems to be needed for a better qualitative comparison of the line shapes although both terms support similar spectral excitations based on the form of the operators in Eqs. (10) and (11) and a spectral weight analysis precludes any quantitative agreement. In both cases, while the RIXS spectral lineshape may be approximated by the two expansion terms, neither provides a faithful representation for the proper two-particle charge response encoded in the simple dynamical structure factor, placing statements about the true charge excitation character of the RIXS crosssection on less solid footing.

D. Consequences for RIXS experiments
The preceding section presented a comparison between cross-sections integrated in momentum, as well as energy for a total spectral weight analysis. In this section we show spectra at two particular momenta: q = (2π/3, 0) and q = (π/2, π/2) (see Fig. 4 and Fig. 5) to underscore those results, shown in a context amenable to experiment. When doped, the Hubbard model, and by extension cuprates, will possess spin and charge excitations in a similar low energy regime which will appear, either directly or in a more complicated way reflecting the complexity of the cross-section, in the RIXS spectrum for the crossed-and parallel-polarization channels, respectively. Thus, to satisfactorily distinguish between the magnetic and charge channel, or two-spin excitations, one must perform measurements which can discriminate the outgoing polarizations. Unfortunately, to this point RIXS experiments have been unable to fully distinguish between the cross-polarized and parallel-polarized channels, making some statements with the help of a careful analysis of experimental RIXS scattering geometry. The newly constructed, state-of-the-art RIXS end-station at ESRF now FIG. 5. Cross-sections at q = (π/2, π/2) point with the elastic response removed.Top (bottom) panels show spectra for the π − σ (π − π) polarization. Left, middle, and right panels show RIXS spectra calculated for n = 0.83, n = 1, and n = 1.17 electron fillings, respectively. The approximate spectra are weighted with the RIXS form factors, i.e.S(q, ω) → |Wπ−σ| 2 S(q, ω), N (q, ω) → |Wπ−π| 2Ñ (q, ω) and similarly for S (q, ω) andÑ (q, ω).
provides an opportunity to perform such measurements (cf. Ref. [16]); and other end-stations (currently operational or to be commissioned in the coming years) also would allow for the differentiation between the crossedand parallel-polarized channels.
Even with outgoing polarization discrimination, in either the crossed-or parallel-polarized channels one needs to carefully invoke either S(q, ω) orÑ (q, ω) as approximations for the full RIXS cross-section. This especially may be true when analyzing RIXS spectral weights as a function of doping. None of these approximations address intermediate state effects, or equivalently differences in the cross-section with changes in the incoming photon energy. [20] Thus, to address the full resonant profile one always needs to calculate the full RIXS crosssection.

III. CONCLUSIONS
Overall, one must carefully apply approximations for calculating the RIXS cross-section. The nonlocal character of the intermediate state can become particularly important in correlated ground states with longer range entanglement. Therefore, we suggest the full RIXS simulations will be needed to verify the character of excitations.
In the cross-polarized channel, we have shown that on a qualitative level Cu L-edge RIXS line shapes correspond to the spin dynamical structure factor S(q, ω), consistent with lowest order approximations as postulated by the fast collision approximation (or the effective operator approach) [19,22,26,27], see Appendix C. As a consequence, we expect that the line shapes reported from cross-polarized RIXS experiments can be reproduced to some extent by theoretical modeling of the spin dynamical structure factors (or empirically through inelastic neutron scattering experiments when also considering differences in the effective matrix elements between the two techniques). However, the detailed analysis in this paper suggests that a quantitative comparison between RIXS and the two-particle spin and charge dynamical structure factors would be impractical. One should not expect a meaningful comparison between different spectral weights obtained from these different techniques, either experimentally or from simulation.
In the parallel-polarized channel, the situation is further complicated by the operator form taken by the approximations themselves. On a qualitative level the primary contributions seem to follow from higher order termsÑ (q, ω), with a notable exception at half-filling, see Appendix C. However, any precise quantitative comparison to experiment would require calculating the full RIXS cross-section. At the same time none of the terms in the approximations correspond to proper two-particle charge excitations, but rather inherently reflect the complexity of the RIXS process. Thus, while line shapes in RIXS should closely resemble the line shapes of the (projected) charge excitations, the spectral weights may be quite different, with a more complicated analysis required to tease out the character of various spectral peaks. The UCL approximation should be valid when all the relevant eigenenergies of the intermediate state Hamiltonian H are much smaller than the inverse core-hole lifetime Γ. We use the following spectral decomposition where {|N } are eigenstates of H with energy {E N }. We are interested in RIXS at the resonant edge between the 2p 0 3d * 1 initial configuration and the 2p 1 3d * 0 intermediate state configuration. (All expressions are presented in hole language.) This means that we need to exclude all intermediate states |N which contain one hole in the d orbital on site j, i.e. on the core-hole site. We add a projection operatorP j and rewrite the above expression as We define the following HamiltoniansH = H + U c iασσ n d iσ n p iασ andH c = H c − U c iασσ n d iσ n p iασ . Note thatH andH c commute (and H =H +H c ), since the intermediate states of RIXS are such that they always contain either a hole in the 2p shell or in the 3d shell (guaranteed by the projection operatorsP j ). Then we obtain: where |n are eigenstates ofH with energy ε n , |n c are eigenstates ofH c with energy ε nc . Consequently we can write Note that with a single hole in the p shell, there are just two eigenstates ofH c : n c ∈ |L 2 , |L 3 with energies ε L2 , ε L3 . They correspond to the two 'j' eigenstates j = 1/2 and j = 3/2, respectively, split by the spin-orbit coupling ∝ λ. (We note that j represents the angular momentum and j represents site j on the cluster.) This implies Resonance approximation We assume that the incoming x-ray photons are tuned to the L 3 resonance, i.e. ω i ε L3 . Since λ Γ (and all eigenenergies ε n λ), we can neglect the contribution from the L 2 resonance and obtain . (17) UCL expansion We adopt the UCL expansion [32,40,41] for Cu L-edge RIXS to obtain This means that the RIXS operator can be rewritten as Second order UCL approximation Keeping only terms with l = 0 and l = 1 gives The validity of this approximation has been discussed in detail in the main text of the paper.
Change of projection operators It is convenient to use another operator -P j -which projects to the sector with no doubleoccupancy in the d level on site j, i.e. where the core-hole is created. This gives A similar expression holds also for the D † j,e f dipole operator.
The l = 1 term Next we evaluate Since wherē due to the fact that all terms in the Hamiltonian which contain site j will vanish when "sandwiched"' between the dipole operators D and evaluated on the initial state |i . Here m(j) and m (j) are nearest and next nearest neighbors of site j. Thus, the following expression holds: Commuting the HamiltonianH −H mj (which does not contain operators on site j) with the operator D j,e i and P j to gives SinceH|i = H|i = E i |i = 0 to set the energy zero, we are left with the following expression P j D † j,e f |L 3 L 3 |D j,e i P jHmj |i . However, due to P j , the U and the U c terms will never contribute and P j D † j,e f |L 3 L 3 |D j,e i P j d † lσ d jσ )|i = 0 for any l = m, m . Thus, we obtain Note the asymmetry in the above expression, i.e. the lack of the hermitian conjugate terms ∝ d † jσ d mσ -this asymmetry expresses the fact that in RIXS we are only sensitive to sites on which the 3d holes reside.
Introducing so-called local matrix elements of RIXS, we obtain where the local RIXS form factors follow from e.g. Ref.
[27] (cf. Eq. (2) and Fig. 1 in Ref. [27]): Let P j n j P j =ñ j [whereñ j = σñ jσ = σd † jσd jσ and d † jσ = d † jσ (1 − n j,−σ )] and P j S z j P j = S z j . Combining the above equations, we finally arrive at the expression for the RIXS operators in the UCL approximation Substituting the above expressions into Eq. (1) and performing Fourier transformations, we obtain Eqs. (6) and (9). As the first (second) order UCL terms have real (imaginary) contributions, the interference terms vanish and the full RIXS cross-section consists of separate first and second order UCL terms. Note that if only these first two terms are considered in the UCL approximation, then the RIXS spectrum does not depend on the size of the core-hole potential U c (the latter will appear only in higher order corrections in the UCL approximation). double occupancies), it is expected that at relatively low doping levels, the first terms should be dominant (even though their amplitude scales with J/Γ and not with t/Γ). Hence, we compare the spectra of these first two terms, S 2 (q, ω) and S 3 (q, ω), in the second order of the UCL expansion with the spectra of S (q, ω) andÑ (q, ω) (cf. Fig. S1) from the Hubbard model. We see that at half-filling S (q, ω) can be approximated relatively well by the '3-spin' excitations probed by S 3 (q, ω) and that N (q, ω) can be approximated relatively well by the "2spin" or bi-magnon excitations probed by S 2 (q, ω). Most of the discrepancies between these two spectra can be found in the high energy regime and are therefore attributed to the failure of the t-J model expansion at higher energies. The electron-and hole-doped cases show much less pronounced agreement, where the projected spin and charge excitations S (q, ω) andÑ (q, ω) also probe the low energy charge excitations below the gap which can give a relatively large contribution in the spectrum.
The prediction that RIXS can probe the "2-spin" and the "3-spin" excitations at half-filling already has been put forward in Refs. [29-32, 42, 43] in the case of the Heisenberg model, consistent with our UCL expansions for the Hubbard and t-J models. Note that usually the Greens functions containing the '2-spin' and the '3-spin' operators are referred to as probing the 'two-magnon' and the 'three-magnon' spectrum, though this terminology may be used loosely in this context. Finally, the fact that the charge dynamical structure factor for the half-filled Hubbard model probes the 'two-magnon' spectrum also has been discussed in the context of nonresonant Raman scattering (cf. Ref. [18, 44, and 45]).  Fig. S2 shows the dominant Cu L-edge RIXS process in the cross-polarized channel in hole language. The hole in the 3d * orbital on a particular site j in the initial state of RIXS is transferred via the dipole operator D into the 2p orbital on the same site j with a non-uniquely defined spin in the intermediate state of RIXS due to spin-orbit coupling in the core. This is transferred back via the dipole operator D † to the d orbital on the same site j with a spin flip in the final state compared to the initial state of the RIXS process. While we demonstrate this process on a single site j in real space, in reality a coherent superposition of such excitations are created with phase factors e ıqj ; this leads to a well-defined, single spin flip with momentum q in the final state of RIXS, i.e. RIXS is sensitive to the spin dynamical structure factor S(q, ω). Fig. S3 shows the dominant Cu L-edge RIXS process in the parallel-polarized channel in hole language. The hole in the 3d * orbital on a particular site j in the initial state of RIXS is transferred via the dipole operator D into the 2p orbital on the same site j with a well defined spin in the intermediate state. In the intermediate state a "shakeup" happens, which creates a "2-spin" and/or charge excitation in the final state. While this process is shown on a single site j (and its neighbors) in real space, in reality a coherent superposition of such excitations is created with a phase factor e ıqj ; leading to a "2-spin" or charge excitation created with momentum q in the final state, i.e. RIXS has some sensitivity toÑ (q, ω), which unfortunately have no simple analog in the standard twoparticle charge response function.