Persistent spin excitations in doped antiferromagnets revealed by resonant inelastic light scattering
Abstract
How coherent quasiparticles emerge by doping quantum antiferromagnets is a key question in correlated electron systems underlying an understanding of the phase diagram of copper oxides. Recent resonant inelastic xray scattering (RIXS) experiments in holedoped cuprates have purported to measure high energy collective spin excitations that persist well into the overdoped regime and bear a striking resemblance to those found in the parent compound, challenging the perception that spin excitations should weaken with doping and have a diminishing effect on superconductivity. Here we show that RIXS at the Cu edge indeed provides access to the spin dynamical structure factor once one considers the full influence of light polarization. Further we demonstrate that highenergy spin excitations do not correlate with the doping dependence of T, while lowenergy excitations depend sensitively on doping and show ferromagnetic correlations. This suggests that highenergy spin excitations are marginal to pairing in cuprate superconductors.

Stanford Institute for Materials and Energy Sciences, SLAC National Accelerator Laboratory and Stanford University, Menlo Park, California 94025, USA

Department of Applied Physics, Stanford University, Stanford, California 94305, USA

Department of Physics, Stanford University, Stanford, California 94305, USA

Advanced Photon Source, Argonne National Laboratory, Lemont, Illinois 60439, USA

Department of Physics and Astronomy, University of British Columbia, Vancouver, British Columbia, Canada V6T 1Z1

Quantum Matter Institute, University of British Columbia, Vancouver, British Columbia, Canada V6T 1Z4

Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto, 6068502, Japan

Department of Physics and Astrophysics, University of North Dakota, Grand Forks, North Dakota 58202, USA
Introduction
Initial Cu edge RIXS measurements on undoped and weakly underdoped cuprates [Lucio2010, Guarise2010] complemented earlier neutron and Raman scattering experiments, seeming to favor a spinfluctuation scenario as a viable explanation of superconductivity [Tacon2011]. However, more recent RIXS measurements on overdoped cuprates [Tacon2011, Tacon2013, Dean2013a, Dean2013b] have shown persistent highenergy spin excitations to very high doping levels where superconductivity disappears. In contrast neutron and Raman measurements display an absence of robust spin excitations in the overdoped regime [Wakimoto2007, Yuan2012]; this conflict undermines an understanding of unconventional superconductivity in the cuprates, making an investigation of how spin excitations manifest in the RIXS cross section a crucial component to its resolution.
In this letter, we reconcile these seemingly incompatible experimental results by computing Cu edge RIXS spectra using exact diagonalization (ED), capable of reproducing major experimental features. We demonstrate with light polarization analysis that the RIXS crosssection in a crossedpolarization geometry can be interpreted simply in terms of the spin dynamical structure factor which enables a comparison between different scattering experiments. Utilizing determinant quantum Monte Carlo (DQMC), we study in detail the momentum and doping dependence of , finding strong changes near both and with relatively insensitive antiferromagnetic zone boundary (AFZB) paramagnons upon hole doping. Moreover, with electron doping these same AFZB paramagnons harden significantly, providing a testable experimental prediction from this work. Underlying this observed behavior is a framework of local spin exchange which remains robust even with significant doping away from the parent antiferromagnet. In contrast, our calculations show a sensitive evolution of lowenergy paramagnons near and which give evidence for the predominance of ferromagnetic correlations. These results highlight the importance of spectral weight and dispersion at low energies in establishing a relevant energy scale and strength of spin fluctuations for pairing rather than higherenergy AFZB paramagnons.
Results
0.1 The relationship between RIXS and
RIXS is a resonant technique and its sensitivity to magnetic excitations arises as a result of corelevel spinorbit interactions in the intermediate state [see Fig. 1(a)]. Although in Mott insulators it has been shown that the RIXS crosssection can be approximated by when the charge excitations are gapped, it is not clear whether the same approximation carries over to doped systems where the ground state is no longer that of a Mott insulator with commensurate filling [Ament2009, MauritsPRL].
To answer this question, we numerically evaluate the RIXS crosssection as a function of momentum [Fig. 1(c)] directly from the KramersHeisenberg formula [RevModPhys.83.705, Jia2012] using small cluster ED of an effective singleband Hubbard Hamiltonian (including both nearest and nextnearest neighbor hopping and onsite Coulomb repulsion ) at various electron concentrations ; details are given in the Methods section. Figure 1(c) displays the RIXS spectra calculated for the experimental geometry discussed in Ref. [Tacon2011]. The RIXS spectra, even without outgoing polarization discrimination, agree well with for different electron concentrations at the chosen momentum space points accessible on the same finite size clusters for each calculation. This is particularly true at halffilling where the charge gap ensures that only spinexcitations can be visible in the given energy range. The main differences occur in the doped systems at higher energy (close to ) which are of less interest for our spin analysis. The important result shown in Fig. 1(c) concerns the suppression of these higher energy peaks in the crosspolarized geometry [see Fig. 1(c) RIXS] which results in a significant improvement in the comparison between the RIXS crosssection and . This indicates that Cu edge RIXS with crossed polarizations (a fourparticle correlator) provides access to the spin excitation spectrum encoded in (a twoparticle correlator) for doped as well as undoped cuprates. [Ro further confirm this agreement between RIXS and for more momentum points, we can manually adjusted the Cu edge energy to be so that momentum points up to can be reached. For more details see Supplementary Note 1]
0.2 The momentum and doping dependence of
Having established a relationship between RIXS and , we focus now on the momentum and doping dependence of for the singleband Hubbard model. Here we employ the numerically exact DQMC method [see Methods] with maximum entropy analytic continuation on larger lattices with fine control of the electron concentration through the chemical potential. As shown in Fig. 2 the DQMC calculations qualitatively reproduce both the momentum and doping evolution of the RIXS measurements found in Refs. [Lucio2010, Guarise2010, Tacon2011, Tacon2013, Dean2013a, Dean2013b]. A comparison between the intensity and dispersion of lowenergy magnetic excitations near and shows a transition from antiferromagnetic to ferromagnetic spin correlations with increasing doping. In an antiferromagnetic system, the dynamical spin structure factors show gapless excitations at both and , with strong intensity around . In a ferromagnetic system, the dynamical spin structure factors show strong intensity with gapless excitations at and much weaker intensity with gapped excitations when approaching . However, the spectra of higherenergy AFZB paramagnons show relatively little change with hole doping other than a general decrease of intensity, suggesting that spin excitations do not soften even in the heavily overdoped regime. For electron doping, AFZB paramagnons surprisingly harden by 50% at 15% doping which has been observed recently in the prototypical electrondoped cuprate NdCeCuO. [WeiSheng] For additional analysis and discussion, see Supplementary Discussion.
0.3 Theory to understand AFZB paramagnons
The behavior of these AFZB paramagnons stands in stark contrast to naive expectations of spin softening with either hole or electron doping: (i) longrange AF order collapses quickly upon doping with a small (intermediate) concentration of holes (electrons), and one would expect spin excitations to soften accordingly [Vladimirov2009, Kar2011]; (ii) shortrange AF correlations are further weakened due to a dilution of AF bonds [see Fig. 3(a)] in the locally static spin picture. Indeed, the nearest neighbor spinspin correlations from the DQMC calculations with electron doping decrease in a manner surprisingly well described by a locally static spin picture where the doped electrons are immobile, as shown in Fig. 3(b).
We can address these points by considering the role of threesite exchange [Daghofer2008] which lowers the system energy when both doped carriers and AF correlations are present [see Fig. 3(a)]. A local spin flip in an otherwise AF background produces a ferromagnetic alignment of nearest neighbor spins which costs additional energy (of the same order as the spin exchange ) by suppressing hole (or doubleoccupancy) delocalization represented by the threesite terms. In fact, if we consider only the spin exchange contributions, the combined energy of a single spin flip in the doped system (breaking both spin exchange and threesite bonds) is larger than that of the undoped system by [see Fig. 3(a)]. This hardening has been observed in ED calculations of for the , , and Hamiltonians [see Methods] as shown in Fig. 3(c) upon electron doping.
The situation is more subtle with hole doping, because this “locally static model” no longer completely applies as seen in Fig. 3(b). The negative nextnearest neighbor hopping (positive for electron doping) promotes magnetic sublattice mixing and a much larger destruction of the AF correlations [Bala1995]. With hole doping the trend observed in RIXS is fully recovered only in the Hubbard model, as shown in Fig. 2, implying that higher order processes absent in type models become crucial in quantitatively reproducing the spin wave dispersion [Dellanoy2009].
Discussion
How do these results reconcile the seemingly contradictory observations between RIXS, neutron and Raman scattering? First, inelastic neutron scattering probes spin excitations particularly well around momentum transfer, showing a vanishing spectral weight in the regime of large hole doping [Wakimoto2007, JPSJNeutronReview]. This behavior is also visible in the numerical results presented in Fig. 2(a) which suggest that the impact of doping on the intensity and dispersion of excitations near and is not symmetric. The decreasing correlation length with doping, evidenced by the spin gap at and the weak dispersion towards , thus impacts these momentum points more strongly than the AFZB paramagnons, in accordance with a locally static picture. Second, Raman scattering [Yuan2012, Muschler2010, Onose2004] shows a softening of the socalled bimagnon (double spinflip or twomagnon) response upon both hole and electron doping. This trend has been reproduced by our ED calculations of the Raman response shown in Fig. 4 [see Methods for the calculation details and the verification of bimagnon peaks]. Strong magnonmagnon interactions reduce the bimagnon Raman peak energy from twice that of the single magnon bandwidth as determined by AFZB magnons and quickly reduce the overall intensity. Taken as a whole, our results provide a qualitative, and in some cases quantitative, agreement with the salient experimental features of neutron scattering, Raman, and RIXS measurements, suggesting that coherent propagating spin waves quickly disappear with the destruction of longrange AF order upon doping, while shortrange, single spinflip processes can survive to high doping levels as reflected in the evolution of .
Full polarization control will allow RIXS to become an effective tool for directly observing spin dynamics along the AFZB, particularly noting the electron/hole doping differences. Together with the dome shaped superconducting phase diagram, these results imply that AFZB spin fluctuations might play a relatively minor role in the pairing mechanism, consistent with established experimental and numerical observations [Scalapino2012, Tallon, White_Scalapino_1999]. This calls into question a simple view of pairing which emphasizes only the spin exchange energy scale . However, we suggest that a definitive resolution to this issue would come from future RIXS experiments along the BZ diagonal (out to ) to illuminate the evolution from antiferro to ferromagnetic correlations, compare with neutron scattering results, and ultimately shed additional light on the intriguing mystery of cuprate hightemperature superconductivity.
Methods
0.4 Numerical techniques
We use exact diagonalization (ED) to evaluate the RIXS crosssection from the KramersHeisenberg formula [RevModPhys.83.705], spin dynamical structure factor , and Raman scattering crosssection [RevModPhys.79.175] on small clusters with periodic boundary conditions. We employ a 12site Betts cluster in evaluating the RIXS crosssection and shown in Fig. 1. The Raman scattering response shown in Fig. 4 has been evaluated on 16 and 18site square (or diamondshaped) clusters and the 18site cluster was employed to evaluate for , and shown in Fig. 3. The ED calculations for are performed with the Parallel ARnoldi PACKage (PARPACK) and the crosssections obtained by use of the biconjugate gradient stabilized method and continued fraction expansion [Jia2012]. The ED calculations on and models are performed using the Lanczos algorithm. Finite temperature determinant quantum Monte Carlo (DQMC) simulations [Blank1981, White1989] were performed on to obtain the imaginary time spinspin correlation function from which the real frequency response function was obtained by analytic continuation using the maximum entropy method (MEM). [Jarrell1996, MEM2] These simulations were performed on square lattice clusters with periodic boundary conditions at an inverse temperature for the same Hubbard Hamiltonian parameter values utilized in the ED studies. For this set of parameters, the DQMC method exhibits a significant fermion sign problem [SignProblemPRB] over the entire doping range which we address in the MEM [MEM2] (see Supplementary Methods and Supplementary Fig. 6). MEM requires the use of a model function for determining an entropic prior in the analytic continuation routine. We utilize a Lorentzian model whose peak as a function of is determined from a simple spin wave dispersion at small out to the AFZB; however, beyond the AFZB the model assumes no softening as expected for longrange antiferromagnetism with the top of the magnon band set by approximations for the spin exchange and an assumed reduction of the spin moment by quantum fluctuations. While some quantitative changes occur with significant changes to these default models, we have checked that the qualitative behavior remains robust. The MEM routine returns the real frequency spin susceptibility from which is obtained from the fluctuationdissipation theorem. More details about the models and numerical algorithms can be found in the following Methods and Supplementary Methods.
0.5 Rixs
The Cu edge RIXS crosssection is calculated using the KramersHeisenberg formula [RevModPhys.83.705] for the singleband Hubbard model
(1)  
in which
(2)  
where is the momentum transfer; and are the incident photon energy (in our study the Cu edge) and photon energy transfer, respectively; is the ground state energy of the system in the absence of a corehole; is the ground state wave function; (and h.c.) dictates the dipole transition process from Cu to the level (or from Cu to ), with the xray polarization either or (the polarization vector parallel or perpendicular to the scattering plane); and is the inverse corehole lifetime (see Supplementary Note 2). In , and represent a sum over the nearest and next nearest neighbor sites respectively. The Hamiltonian for the intermediate state also involves the onsite energy for creating a core hole, Coulomb interaction induced by the corehole and spinorbit coupling , all denoted as in . represents the spinorbital coupling coefficients. The angle between the incident and the scattered photon propagation vectors is set to be . The parameters used in the RIXS calculation are eV, eV, eV, eV, eV, eV and eV [Tsutsui2000, Kourtis2012]. RIXS spectra at halffilling are taken only at the Cu resonance, and upon doping at the resonance closest to the halffilling Cu edge resonant energy. The RIXS results were obtained for a Lorentzian broadening with half width at half maximum (HWHM) = () and a Gaussian broadening with HWHM = () on the energy transfer. The spin dynamical structure factor , discussed in the next section, for was calculated using the same parameters to make comparison to our RIXS results.
0.6 Spin Dynamical Structure Factor
The spin dynamical structure factor is defined as
(3) 
where we have studied the , and Hamiltonians:
(4) 
(5) 
is the corresponding ground state energy of the model Hamiltonian; , for ; , for and ; ; and is restricted in the subspace without double occupancy and , in which the operator annihilates a dressed electron whose hopping conserves the number of effective doubly occupied sites [Harris1967]. To explore the similarities and differences between and (with and without ), we calculate on the three model Hamiltonians with the parameters , and (corresponding to by the relation ).
0.7 Locally Static Model
In the “locally static model” [see Fig. 3(b)] it is assumed that the holes destroy the shortrange spinspin correlations solely by the effect of ‘static’ doping, i.e. by removing the spins and thus cutting spin bonds. In this case the nearest neighbor spinspin correlation can be calculated in the following way:
(6) 
where is the abbreviation of for two neighboring sites and , and is the concentration of either doped holes () or doped electrons ().
0.8 Raman Scattering
We calculate the Raman scattering crosssection in the channel using the nonresonant response function for [RevModPhys.79.175]:
(7)  
in which is the bare band dispersion, with parameters and . is taken as eV to make comparison with experimental data.
This twoparticle response also has been studied recently in cluster dynamical meanfield theory [Millis2012] showing that, if calculated fully gauge invariantly, the nonresonant Raman response shows the presence of a strong bimagnon peak at half filling. Nevertheless, the Raman spectrum calculated using this method for doped systems is sensitive to both charge and spin excitations in the low energy regime. Our identification of the bimagnon excitations in the Raman spectra relies primarily on the qualitative evolution of the peaks in agreement with experimental observations [Muschler2010]. We note that all of the excitations visible in our Raman spectra correspond to transitions which have symmetry. At halffilling, the energy of the excitation to which we assign bimagnon character lies within the charge gap which makes the bimagnon assignment clear. Upon either hole or electron doping, we expect to develop charge excitations in the Raman response at low energy and for the twomagnon response to soften and decrease in intensity. Our assignment corresponds to an upper bound for the bimagnon energy scale with doping where the additional structure at low energies signals either charge excitations or a substantial broadening of the bimagnon excitations now represented by multiple features in the exact diagonalization (ED) result (as discussed in connection with comparisons between determinant quantum monte carlo (DQMC) and ED results). However, the energy scale clearly softens and, more importantly, the intensity drops (significantly) in agreement with the experimental observations where the bimagnon “peak” becomes nearly impossible to distinguish from the charge background almost immediately upon crossing the AFM phase boundary.

We would like to thank G. Ghiringhelli, R. Hackl, J. P. Hill, B. J. Kim, M. Le Tacon, W.S. Lee and J. Tranquada for discussions. This work was supported at SLAC and Stanford University by the U.S. Department of Energy, Office of Basic Energy Sciences, Division of Materials Science and Engineering, under Contract No. DEAC0276SF00515 and by the Computational Materials and Chemical Sciences Network (CMCSN) under Contract No. DESC0007091. CJJ is also supported by the Stanford Graduate Fellows in Science and Engineering. CCC is supported by the Aneesur Rahman Postdoctoral Fellowship at Argonne National Laboratory, operated under the U.S. Department of Energy Contract No. DEAC0206CH11357. TT is supported by the GrantinAid for Scientific Research (Grant No. 22340097) from MEXT. TT and TPD acknowledge the YIPQS program of YITP, Kyoto University. A portion of the computational work was performed using the resources of the National Energy Research Scientific Computing Center (NERSC) supported by the U.S. Department of Energy, Office of Science, under Contract No. DEAC0205CH11231.

The authors declare no competing financial interests. Correspondence and requests should be addressed TPD (tpd@stanford.edu).
Supplementary Material: Persistent spin excitations in doped antiferromagnets revealed by resonant inelastic light scattering
Supplementary Figures
Supplementary Notes
0.9 Supplementary Note 1: The comparison of with RIXS at other momentum points
Showing more momentum points for comparison between RIXS and would better demonstrate the agreement between the two spectra. In the main text other momentum points provided by the largest available 12site cluster has not been shown as these other points are not “accessible” in Cu edge RIXS due to the limited momentum transfer from photons at the Cu edge energy (). Nevertheless, to answer solely questions about differences between and RIXS, one can “artificially” adjust the photon momentum by changing the Cu edge photon energy to cover the full Brillouin zone. Supplementary Fig. S1 shows the RIXS crosssection in the incoming/outgoing polarization (red) and the spin dynamical structure factor (blue), with the Cu edge energy tuned to for the remaining three momenta accessible on the 12site cluster: , , and (the latter is accessible in RIXS, but the intensity is nearly zero, so it was not shown in the main text). Each spectrum has been calculated using exact diagonalization (ED) for the Hubbard model on a 12site cluster for electron concentration . One can see that for each momentum, RIXS in the geometry approximates well the spin dynamical structure factor . These results reconfirm the agreement between and the RIXS crosssection in the crosspolarized scattering geometry.
0.10 Supplementary Note 2: The corehole lifetime dependence
Our choice of in the main text is very realistic for cuprates – in particular, it agrees with the most recent estimates in Ref. [Bisogni2013] and is also in agreement with the chosen value of in Ref. [Igarashi2012]. While the RIXS results might sensitively depend on , choosing it in agreement with experiments properly captures all the effects related to finite corehole lifetime. Nevertheless, decreasing or increasing the assumed value of by a factor of two still does not lead to significant differences in the spectral weight of excitations observed in the calculated RIXS spectrum. More precisely, even for unphysical values of , RIXS in the crosspolarized channel and in the doped case approximates fairly well the spin dynamical structure factor , as shown in Supplementary Fig. S2. Naturally, further reduction of the value of (which definitely requires unrealistic values of at the Cu edge) leads to significant changes in the RIXS spectra – but the study of this phenomenon is beyond of the scope of this paper.
Supplementary Discussion
0.11 The behavior of at and
In Supplementary Fig. S3 we show the behavior of the spin dynamical structure factor at as a function of doping, in analogy with Fig. 2c in the main text. While the momentum point is hard to access in INS [JPSJNeutronReview] and has so far not been accessed in RIXS, the spectrum at this point could be accessed and verified against these findings in future RIXS experiments.
Let us note, that the behavior of at as a function of doping is very similar to the behavior at as a function of doping (cf. left panel of Supplementary Fig. S3 and Fig. 2c in the main text). However, unlike the point, it is clear that the origin of the hole doping dependence can be explained less accurately using the model with 3site terms (although the 3site terms contribute to the hardening of the excitations with electron doping), cf. right panel of Supplementary Fig. S3. One must also invoke higher order exchange processes (present in the Hubbard model) to fully explain the observed doping dependence of these spin excitations. We note that the “local picture” (see Fig. 3 of the main text) better represents the behavior of along the  direction rather than along the diagonal direction in the Brillouin zone due to a more significant destruction of the antiferromagnetic correlations close to with doping.
For completeness, we show the doping dependence of the integrated intensity and the half width at half maximum (HWHM) at the and points in Supplementary Fig. S4. Unlike the point (cf. Supplementary Fig. S5 and discussion in the next section) the variation of intensity with doping at these momenta is rather small. This agrees with the assumptions made when interpreting the RIXS results in terms of the “locally static model” and stating that one is primarily sensitive to local spin excitations at these momentum space points: the energy of these local excitations seems to be rather weakly affected by doping (see main text) and one also could expect that their intensity should not change very much with doping.
0.12 The behavior of near and the connection to INS
Our calculated is consistent with inelastic neutron scattering (INS), capturing important doping dependent characteristics of the measurements. We find qualitative agreement between the DQMC results, as shown in Supplementary Fig. S5, and the INS measurements on the following points:
On the hole doped side: (i) The momentumintegrated magnetic scattering intensity in a region around in the INS experiments strongly decreases with hole doping [cf. Fig. 7 (b) and Fig. 8 in Ref. [JPSJNeutronReview]]; similarly, in our calculations the intensity of the magnetic response at and nearby points decreases significantly with hole doping. (ii) The energy scale of magnetic excitations near increases with significant hole doping [JPSJNeutronReview], similar to the behavior we observe.
On the electron doped side: (i) The momentumintegrated spectral weight of the socalled “low energy peaks” in the INS experiments seems to decrease more slowly with electron doping than with hole doping [cf. inset of Fig. 18 in Ref. [JPSJNeutronReview]]; in our calculations we observe a similar trend – the intensity of the peak decreases far slower with electron doping than with hole doping. (ii) The slope of low energy magnetic excitations near (which provides an estimate for the antiferromagnetic spin stiffness) becomes less steep with doping [cf. Fig. 18 in Ref. [JPSJNeutronReview]]; we qualitatively observe a similar trend.
Nevertheless, the relatively high temperature and lack of a dense momentumspace mesh which are limiting factors of the DQMC method preclude any quantitative, substantive comparison between our results and some of the aspects of the INS spectra. In particular, our momentumspace mesh does not provide sufficient resolution to address the magnetic incommensurability which is visible as part of the INS “hourglass” structure around .
Supplementary Methods
0.13 Fermion sign problem in DQMC
The DQMC method for simulating the manybody Hubbard Hamiltonian exhibits a sign problem which restricts the minimum accessible temperatures (maximum ) and, to a lesser extent, the maximum accessible cluster sizes [SignProblemPRB]. The sign problem arises from a nonpositive definite measure of the probability distribution for different auxilliary field configurations in the DQMC method. The sign problem can be partially overcome by performing measurements during DQMC process using a weighted distribution which accounts for this nonpositive measure or “fermion sign”. We have accounted for this source of statistical variation between different estimates of the correlation function using the maximum entropy method (MEM) analytic continuation techniques based on Bayesian inference [MEM2], and the twoparticle response can be reliably extracted for the values of the sign encountered in the simulations. (see Supplementary Fig. 6)
0.14 ED vs DQMC
Both ED and DQMC are numerically exact methods. However, there are significant differences in the types of questions which can be better answered by each technique. Both methods are based on small, realspace cluster models of the full lattice Hamiltonian. More specifically:
ED is a wavefunctionbased method restricted by the exponential scaling of the Hilbert space dimension with cluster size which leads to two effects: (i) ED can be utilized only on small clusters. The intermediate states of the Cu edge RIXS process for a 12site cluster discussed in the main text, considering the Cu core orbitals, has a Hilbert space dimension of = 14,976,864 at halffilling (while the size of the Hamiltonian matrix is ). The Hamiltonian matrix is usually very sparse (usually a few hundreds of nonzero elements within each column or row), and can be handled by implementing massively parallel algorithms. By extension, the size of the realspace cluster limits the number and spacing of accessible momentum space points for any comparison between RIXS and in our study. (ii) Any spectral function, e.g. singleparticle addition or removal spectrum, spin and charge dynamical structure factors, RIXS crosssection, or Raman response function, only has finite spectral support over the discrete spectrum of eigenstates for the Hamiltonian represented in the given Hilbert space. The spectral features thus naturally evolve, often significantly, with changes in the cluster size. The advantages of ED over other techniques rest with the fact that it is a zero temperature formalism capable of addressing ground state properties and that it provides direct access to the wavefunction from which one can evaluate any correlation function or crosssection, i.e. RIXS and on equal footing (albeit with severe limitations on the Hilbert space size). In particular, RIXS at the transition metal edge has been calculated by several groups using the ED technique [RevModPhys.83.705].
DQMC is a finite temperature, imaginary time Green’s function or propagator based method limited mainly by the fermion sign problem. One can perform simulations on significantly larger clusters with a more dense momentum space mesh than that of ED; however, the fermion sign problem limits the lowest temperatures that can be accessed with this method and one needs to employ an additional numerical procedure for obtaining real frequency response functions from the imaginary time data. While we can obtain both single and twoparticle response functions using this method, there is to this pointintime no method for estimating the multiparticle RIXS crosssection with this technique which requires an explicit handling of the corehole intermediate state. So, the DQMC method can provide full realfrequency information about single and two particle response functions on a relatively dense mesh in momentum space, but limited to fairly high temperatures.
Despite their significant differences, these two methods give qualitatively similar results as shown in Supplementary Fig. S7. This figure displays a comparison between obtained from DQMC simulations on square clusters and that obtained from ED on both 12 and 18site small clusters at approximately the same doping level and momentum space position. One can see that the lowest energy magnon peak aligns in energy for both DQMC and ED. The presence of several discrete higher energy states in ED makes up part of the tail emanating on the high energy side of the main magnon peak in DQMC. In fact, increasing the ED cluster size shows that these states form a continuous band of excitations captured by DQMC. This behavior is a manifestation of the discrete energy domain which supports spectral functions obtained using the ED technique.
The origin of the “high energy” excitations, visible both in ED (as separate peaks at ca. energy transfer) and DQMC (as the high energy “tail” extending beyond ca. energy transfer), may be explained as follows: (i) Due to the strong interaction between holes/electrons and spin fluctuations, there are many different “spinflip” states in . (ii) These states can be understood as spin fluctuations dressed with holes/electrons (cf. Ref. [Martinez1991] for a somewhat related problem of the hole dressed by magnons which also leads to spectral functions with “many peaks”). It is remarkable that these dressed spin fluctuations have a very similar dispersion relation to that of the “pure” spin fluctuations in the undoped case. The findings in this manuscript point to the origin of this effect in the peculiar interplay between the 3site terms, softening of spin excitations due to spin removal with hole/electron doping, and, most importantly, rather small effects on dispersion, at least on the electron doped side, due to the dressing of spin fluctuations with charge.