Magnetic Susceptibility and Simulated Neutron Signal in the 2D Hubbard Model

TL;DR

This paper computes and compares dynamic spin susceptibilities in the 2D Hubbard model using advanced methods, relating findings to neutron experiments on cuprates, and explores how spin excitations evolve with interaction strength and doping.

Contribution

It introduces a comprehensive analysis of spin susceptibilities in the 2D Hubbard model using Dual Fermions, comparing results with other methods and experimental data, highlighting the evolution of spin excitations.

Findings

Identification of spin modes via peaks in Imχ(ω,q)

Observation of a 'Y'-shaped dispersion similar to experiments

Analysis of the evolution of spin excitations with doping and interaction strength

Abstract

We compute dynamic spin susceptibilities in the two-dimensional Hubbard model using the method of Dual Fermions and provide comparison to lattice Monte Carlo and cluster dynamical mean field theory. We examine the energy dispersion identified by peaks in ${\rm Im}\chi(\omega,q)$ which define spin modes and compare the exchange scale and magnon dispersion to neutron experiments on the parent La$_2$CuO$_4$ cuprate. We present the evolution of the spin excitations as a function of Hubbard interaction strengths and doping and explore the particle-hole asymmetry of the spin excitations. We also study the correlation lengths and the spin excitation dispersion peak structure and find a `Y'-shaped dispersion similar to neutron results on doped HgBa$_2$CuO$_{4+\delta}$.

Equations

$$H=\sum _{k,\sigma }\left({ϵ}_k-\mu \right)c_{k\sigma }^{†}{c}_{k\sigma }+U\sum _i{n}_{i\uparrow }{n}_{i\downarrow },$$

Full text

I Introduction

The interplay between charge, spin and superconducting fluctuations in the cuprate high-temperature superconductors has been experimentally documented with nuclear magnetic resonance (NMR)Fujita et al. (2012), inelastic neutron scattering (INS)Coldea et al. (2001); Lipscombe et al. (2009, 2007); Vignolle et al. (2007), resonant ultrasound spectroscopy (RUS)Shekhter et al. (2013); LeBoeuf et al. (2013), thermal probesGrissonnanche et al. (2015); Badoux et al. (2016), resonant inelastic X-ray scattering (RIXS)Le Tacon et al. (2011); Braicovich et al. (2010), as well as RamanDevereaux and Hackl (2007) and optical spectroscopies.Basov and Timusk (2005) All together, they paint a picture of a complex phase diagram comprised of many competing states where antiferromagnetic spin fluctuations play a crucial role.

Providing a theoretical description of even the normal high-temperature state of the cuprates has proven to be a formidable challenge. Much has been understood on the level of single particle properties including the pseudogap in cluster DMFT studies of the Hubbard modelHuscroft et al. (2001); Parcollet et al. (2004); Macridin et al. (2006); Kyung et al. (2006); Werner et al. (2009); Gull et al. (2009); Lin et al. (2009); Liebsch and Tong (2009); Sakai et al. (2009); Lin et al. (2010); Gull et al. (2010); Sordi et al. (2012); Lin et al. (2012); Merino and Gunnarsson (2014); Gunnarsson et al. (2015a); Chen et al. (2017) and its interplay with superconductivity.Lichtenstein and Katsnelson (2000); Maier et al. (2008); Civelli (2009); Gull and Millis (2012); Gull et al. (2013); Gull and Millis (2013); Sordi et al. (2013); Gull and Millis (2014, 2015); Chen et al. (2015); Maier et al. (2015); Maier and Scalapino (2018) However, making a direct connection to many experimental measurements requires the knowledge of two-particle susceptibilities in addition to single-particle quantities. Among these, magnetic excitation effects have been a central focus.

The theoretical understanding of magnetic excitations remains an ongoing challenge.Chen et al. (2015); Gunnarsson et al. (2015b); Chen et al. (2017) In the case of the insulating high-Tc parent compounds, experimental data has primarily been fit using a linear spin-wave theory in the Heisenberg limit at zero temperature.Coldea et al. (2001); Delannoy et al. (2009) This has allowed for a qualitative description of the excitations in terms of spin models and the determination of the strength of the exchange $J/t$. However, spin models do not describe itinerant fermion systems, such as the doped compounds. In addition, the Heisenberg parameters determined via spin model fits lie outside the regime where spin-models are a valid low-energy approximation of the Hubbard model. What is known to date about the magnetic excitations of the Hubbard model and how they connect to experimental measurements on the cuprates Jia et al. (2014); Zheng et al. (2017); Huang et al. (2017) comes from numerical studies of finite size systems. These methods accurately resolve short-range correlations but are limited in their temperature, doping, and momentum resolution and their ability to resolve long-ranged low-energy spin excitations.

In this paper we address this deficiency. We perform calculations on the Hubbard model in two dimensions using the technique of Dual Fermions (DF),Rubtsov et al. (2008) which is an approximate extension of the non-perturbative dynamical mean field theoryGeorges et al. (1996) and is believed to be accurate at high temperature.Iskakov et al. (2016) The method’s primary advantage is that it can recover continuous momentum dependence, and therefore that it does not suffer from the finite system size effects that limit determinantal lattice quantum Monte Carlo (DQMC) or cluster dynamical mean field theory results. In addition, the computational time is not substantially greater than that of single-site dynamical mean field theory, allowing us to probe the entire phase diagram with reasonable computational expense. We validate our results by comparing to cluster dynamical mean field theory in the dynamical cluster approximation (DCA) variant and DQMC in areas of parameter space that are accessible to those methods. We show that for a certain interaction strength, the dispersion of spin-waves closely resembles that observed experimentally at high temperatures (300K) in La2CuO4.Coldea et al. (2001) We further examine the effects of both electron and hole doping, providing a complete picture of the particle-hole asymmetry of spin excitations in the model.

The remainder of this paper proceeds as follows. In Section II we discuss the DF method and provide comparisons to other methods. Section III contains our main results from the DF method while section IV concludes.

II Methods

We study the single orbital Hubbard model in two dimensions on a square lattice.Scalapino (2007); LeBlanc et al. (2015) The Hamiltonian is

[TABLE]

where $\mu$ is the chemical potential, $k$ is momentum, $i$ labels sites in real-space, $U$ is the onsite Coulumb interaction, and the dispersion is given by $\epsilon_{k}=-2t\left[\cos(k_{x})+\cos(k_{y})\right]-4t^{\prime}\cos(k_{x})\cos(k_{y})$, in which $t$ and $t^{\prime}$ are the nearest and next-nearest neighbor hopping integrals.

II.1 Dual Fermions

We solve Eq. 1 in the Dual Fermion approximationRubtsov et al. (2008) using the open source Dual Fermion code of Ref. Antipov et al., 2015. This method treats all local correlations in a non-perturbative manner and perturbatively includes non-local correlations. In our implementation, self-consistent ladder diagrams for the non-local (‘dual’) self energy are summed in the charge and spin channels.Hafermann et al. (2009) This method is accurate at high temperatureLeBlanc et al. (2015) but uncontrolled, in the sense that contributions from higher order vertices and non-ladder diagrams are not included.Iskakov et al. (2016) Results must therefore be carefully benchmarked against other techniques, both for single- and for two-particle properties. To this end, in this section we present high temperature results from dynamical cluster approximation LeBlanc and Gull (2013); LeBlanc et al. (2015); Iskakov et al. (2016) and from determinantal quantum Monte Carlo.Blankenbecler et al. (1981)

The quantity of interest is the spin susceptibility on the real frequency axis. First, we numerically compute spin susceptibilities on the Matsubara axis as a function of transfer momentum $q$ and bosonic Matsubara frequency $i\nu_{n}=2\pi inT$. Representative simulation data obtained for a half-filled system at intermediate interaction $U/t=8$, temperature $T/t=0.2$, and next-nearest neighbor hopping strength $t^{\prime}/t=-0.3$ are shown in the inset of Fig. 1. Different curves denote susceptibilities at the momentum points indicated along a path from $X(\pi,0)$ to $M(\pi,\pi)$ in the Brillouin zone. Symmetry properties imply that only the real part of the susceptibility, $\text{Re}\chi(i\nu_{n},q)$, is non-zero. Conversion of $i\nu_{n}$ to the real frequency, $\omega$, is needed to obtain the frequency- and momentum-dependent susceptibility $\text{Im}\chi(\omega,q)$. This process relies on numerical methods of analytical continuation. Those results can be seen in the main panel of Fig. 1. $\text{Im}\chi(\omega,q)$ shows a single peak indicated by arrows. Its amplitude is largest at $q=(\pi,\pi)$ and quickly decays away from this point, while its peak position moves to higher frequencies. In the spirit of other works that examine spin structure factors and susceptibilities we define the frequency for this peak at each $q$-vector to be the spin-wave dispersion, $\omega_{s}(q)$.Jia et al. (2014); Zheng et al. (2017)

Analytic continuation of bosonic spectral functions, here based on the ALPSGaenko et al. (2017) open source maximum entropy codeLevy et al. (2017) with a Gaussian default model, is required to obtain the spin-wave dispersion. The analytic continuation is numerically ill posed, in the sense that many bosonic spectral functions will yield the same bosonic Matsubara response within error bars. We believe continuation uncertainties to be larger than any finite size error or stochastic uncertainty. These uncertainties are independent of the approximation error of the self-consistent ladder DF method, which we analyze in detail below.

II.2 Comparison to DCA

We first consider the accuracy of the DF method for single-particle quantities. Fig. 2 compares results for the single-particle self-energy to DCA calculations on an $8\times 8$ cluster for the nodal ($k=(\pi/2,\pi/2)$) and antinodal ($k=(\pi,0)$) points. The Dual Fermion result provides the majority of the momentum dependent contribution and is, at these temperatures and interaction strengths, comparable to DCA at a small fraction of the computational expense.

In Fig. 3 we compare the local spin susceptibility on the Matsubara axis for the lowest and first Matsubara frequencies obtained from Dual Fermions to those obtained in 8-site DCA clusters.Chen et al. (2017) Since DCA is limited in k-space resolution and computationally limited to small clusters when computing two-particle observables we compare only the local susceptibility $\chi_{loc}$. We see agreement at high temperature where the spin susceptibility is expected to be unstructured and the DCA result is close to the dynamical mean field result upon which the Dual Fermion treatment is built.

As temperature is lowered, the lowest Matsubara frequency for DF is consistently larger than DCA until around $T/t=0.2$, below which it declines at $T/t=0.18$ and $0.19$. At even lower $T/t$ the self-consistent ladder summation fails to converge. For this reason we consider $T=0.2t$ to be the lowest accessible temperature at $n=1$, $U/t=6$, and present results at this $T$ in section III.

II.3 Comparison to DQMC

At high temperature, we can compare the dynamical spin susceptibility $\chi(q,\omega)$ on the real frequency axis obtained from DF to those obtained from DQMC calculations. In Fig. 4 we plot the peak energy $\omega_{s}(q)$ of ${\rm Im}\chi(q,\omega)$ determined as in Fig. 1 as a function of momentum for various charge densities. The momenta, along high symmetry cuts in the Brillouin zone, start from $(\pi/2,\pi/2)$ and move diagonally to the $M$-point at $(\pi,\pi)$, via $X=(\pi,0)$ and $(\pi/2,\pi/2)$ to $\Gamma$ and back to $X$. Here, the parameters are set as $U/t=8$, $T/t=1/3$, and $t^{\prime}=-0.3$. Dual Fermion calculations are evaluated on a $16\times 16$ grid (solid symbols), and DQMC calculations are performed on the $8\times 8$ cluster (open symbols). Both of these calculations are converged in their $k$-space discretization. At the lowest density ($n=0.6$), Dual Fermion results are consistent with DQMC results. As we increase the density towards half-filling the spin excitation energy from Dual Fermion calculations are systematically lower than those from DQMC, but the momentum dependence is the same. To illustrate this, we show in magenta on each panel the DQMC curve rescaled by the factor $r$ needed to match $\omega_{s}^{DF}$ at $q=(\pi/2,\pi/2)$ and quote the value of $r$ in each case. After rescaling, we see that the overall structure of $\omega_{s}(q)$ is remarkably similar and that the approximate DF method is qualitatively correct within a single prefactor for each of the parameters examined, except near the $\Gamma$ point, as we comment on in section II.4.

Continuing to increase charge density to $n=1.15$, the spin dispersion around the $\Gamma$ point becomes identical, while the spin excitation at other momenta in Dual Fermion calculations has lower energy. We perform a similar analysis in Fig. 5 at fixed density, $n=1$, for variation in $U/t$. We find that for small values of $U/t$ again the DQMC and DF are quantitatively in agreement while for other values they have a distinct energy scale with the DF result again systematically lower than the DQMC result but with the same momentum dependence. The agreement in both highly doped and weak coupling cases but disagreement for larger interactions near half filling suggests that the method works best where the physics is close to that of the underlying DMFT approximation but may not be able to correctly represent the energy scales in the pseudogap and Mott insulating states.

To further analyze these conclusions, we plot the normalized magnetic susceptibility on the imaginary time axis, $\text{Re}\chi(q,\tau)$, for both calculations in Fig. 6. These data do not suffer from continuation errors. Panel (a) shows $\text{Re}\chi(q,\tau)$ at $n=0.6$ and $q=X$. It is found that $\text{Re}\chi(q,\tau)$ for both calculations are the same, resulting in the same $\omega_{s}(q)/t$ in Fig. 4 after analytic continuation. Panel (b) shows $\text{Re}\chi(q,\tau)$ at $n=1.0$ and $q=X$ where the two methods deviate. To assess the finite lattice size effects in DQMC, we show results for two cluster sizes ($4\times 4$ and $8\times 8$).

Figs. 4 - 5 together indicate that, at the high temperature tested here, the DF and DQMC are qualitatively the same. They differ quantitatively by a momentum independent prefactor that depends on both doping and interaction strength, $r\equiv r(U,\mu)$. As doped results from DCA or DQMC are difficult to obtain at lower temperature because of the fermion sign problem, we leave the temperature dependence of $r$ for future study. We continue under the assumption that the $\omega_{s}(q)$ from DF is correct up to an overall doping-, temperature-, and interaction dependent prefactor. This holds for all momenta except for the area near the $\Gamma$ point, which we address next.

II.4 Discrepancy at $\Gamma$ point

Figure 4 shows that $\omega_{s}(q=\Gamma)$ at $n=1$ does not touch to zero in DF calculations. However, since the spin operator commutes with the Hamiltonian, the spin susceptibility on the real frequency at the $\Gamma$ point should be uniformly zero, i.e. $\chi(q=\Gamma,\omega)=0$ for all values of $\omega$. On the Matsubara axis, this condition enforcing a spin excitation at zero frequency is that $\chi(q=\Gamma,i\Omega_{0})$ is a constant while $\chi(q=\Gamma,i\Omega_{n})=0$ as $n\neq 0$. To understand why DF results violate this condition, we plot spin susceptilibilty $\chi(q,\Omega_{n})$ and its dressed bubble (GG) and vertex terms in Fig. 7. The spin susceptibility equals the sum of the GG and the vertex terms. The right panel of Fig. 7 shows that vertex contribution does not exactly cancel the bubble contribution at $i\Omega_{n}$ for $n\neq 0$. Other methods of restoration exist.Rubtsov et al. (2012); van Loon et al. (2014) The magnitude of this violation is temperature dependent. As temperature increases, the value of the first Matsubara frequency decreases, see Fig. 8, and at high temperature ($T/t=0.5$) becomes essentially uniformly zero, satisfying the expectation that $\chi(q=\Gamma,\omega)=0$ for $\omega=0$.

The violation of this cancellation is induced by the approximate nature of the Dual Fermion procedure, which is based on a vertex function expansion of an auxiliary Anderson impurity coupled to a bath,Rubtsov et al. (2008) a problem for which spin is not conserved. While the exact summation of all contributions to all ordersIskakov et al. (2016) would restore this symmetry, the perturbative ladder Dual Fermion vertex only partially cancels the $GG$ bubble contribution. We expect $\omega_{s}(0,0)$ to decrease in energy as the system becomes metallic. This is reflected in both our high temperature and doped data.

III Results

We present DF results for the dynamical spin susceptibility of the single band Hubbard model obtained in the intermediate interaction strength regime where both weak coupling perturbative methods and strong coupling expansions fail. We have shown that the DF results are qualitatively correct at high temperatures and note that our calculations have the advantage of a fine-grained momentum resolution of up to $64\times 64=4096$ points in the Brillouin zone.

III.1 Spin-wave dispersion at half filling

Fig. 9 shows the momentum dependence of the spin wave dispersion $\omega_{s}(q)$ from the DF method of the half-filled system along high symmetry paths in the first Brillouin zone (see inset), for interaction strengths ranging from $U/t=3$ to $U/t=8$. Results are obtained at temperature $T/t=0.2$ and next-nearest neighbor hopping $t^{\prime}/t=-0.3$. At weak interaction $U/t=3$ (orange circles), the spin mode, $\omega_{s}(q)$, at $M$ and $X$ points is at higher energy than other interaction strengths. Increasing the interaction towards $U/t=8$ results in a decrease in $\omega_{s}(q)$ at these two points. This reduction in energy is accompanied by an increase in amplitude of the susceptibility, indicating strong low-energy spin fluctuations.

We compare our numerical data to the experimental spin susceptibilities obtained on LaCuO4Coldea et al. (2001) (black points, right axis, see also Ref. Headings et al., 2010) in the upper panel of Fig. 10. Also shown is a linear spin-wave fit with parameters determined by Coldea et al. (2001) Our simulations were obtained for $U/t=7.6$ and $t^{\prime}/t=-0.3$. Ref. Coldea et al., 2001 fitted the linear spin-wave results by a set of nearest and further Heisenberg exchange constants, $J,J^{\prime},J^{\prime\prime}$ as well as a substantial ring exchange $J_{c}$, and found $J=138$meV, $J^{\prime}=J^{\prime\prime}=2$meV and $J_{c}=40$meV.

From a comparison between the overall shape of the spin excitation dispersion in experiment, Fig. 10, and calculation, and in particular from the behavior between $M$, $X$, and $(\pi/2,\pi/2)$, we can conclude that an appropriate Hubbard interaction strength for modeling the data should be between $U=7t$ and $U=8t$. Specifically, the weak momentum dependence observed along the $M$-$X$-$(\pi/2,\pi/2)$ lines is inconsistent with an interaction strength much less than $U=7t$.

The lower panel of Fig. 10 shows calculated and experimental data for the imaginary part of the susceptibility at its maximum point $\omega_{s}(q)$. Spin excitations in both calculation and experiment are dominated by the strong low-energy excitation near $(\pi,\pi)$ and proportional to each other.

An obvious point of disagreement between our data for $\omega_{s}(q)$ and the experimental results is the overall magitude of the spin-wave dispersion. For generally accepted values of $U/t\sim 7-8$, and using a value of $t\sim 0.35eV$,Andersen et al. (1995) we find that $\omega_{s}(q=(\pi/2,\pi/2))=0.4t=142meV$. This value is a factor of $2.36$ below the experimental value. Since we expect the DF method to be correct only up to a prefactor (see Fig. 4 where $1/r\approx 2$ at the higher $T/t=1/3$) this scaling factor is consistent.

III.2 Doping dependence of the spin-wave dispersion

The spin-wave dispersion is highly doping dependent (Fig. 11), and data shows a pronounced difference between electron and hole doping, with the spin wave energy increasing much faster for electron than for hole doping, except very near the $M$ point. The interpretation of DF data is complicated by the fact that the rescaling prefactor is strongly doping and interaction dependent (Fig. 4). Also, the temperature dependence of the rescaling factors is unknown (data is only available down to $T/t=1/3$). We show our results without rescaling $r$, but we expect the general conclusions we draw here to be robust based on $r$ at higher temperature. At both intermediate ($U/t=5$) and strong ($U/t=8$) interaction strength, $\omega_{s}(q)$ near $(\pi,\pi)$ quickly increases with doping. At the $X$ point, hole doping reduces $\omega_{s}$ and electron doping increases $\omega_{s}$ at intermediate interaction, whereas for strong coupling electron doping leads to a sharp increase and only a moderate increase and eventual saturation for hole doping. No change is visible upon hole doping at $(\pi/2,\pi/2)$, whereas electron doping increases $\omega_{s}(q)$ by a factor of at least two.

Fig. 12 shows a detailed dependence of the spin-wave dispersion on doping, plotted as a function of occupation, $n$, for $U/t=8$ (left panel) and $U/t=5$ (right panel). In particular a strong difference between electron ($n>1$) and hole ($n<1$) doping is visible in the data at $(3\pi/4,0)$ and $(\pi,0)$, which shows no doping dependence on the hole-doped side but a rather strong momentum dependence on the electron-doped side. The data shown here are not rescaled to DQMC data. This result is qualitatively consistent with determinant quantum Monte Carlo calculations at higher temperature,Jia et al. (2014) and with LSCO data from Dean et al. (2013) (see also Ref. Jia et al., 2014). At $q=(\pi,\pi)$ both electron and hole doping increase $\omega_{s}$ shown in blue squares, which rises to higher energy for doping in each direction.

Finally, as interaction strength decreases from $U/t=8$ to 5 we see that $q=(\pi,\pi)$ maintains a minimal value near $n=1$ but the mode energy increases; and that $q=(\pi,0)$ and $(3\pi/4,0)$ continues to show little change in hole doping and increases upon electron doping.

Fig. 13 shows the amplitude of $\rm{Im}\chi(\omega,q)$ for fixed $q$ as a function of density for $U/t=8$, 5, and 3. The amplitudes of $\rm{Im}\chi(\omega,q)$ at $(3\pi/4,0)$, $(\pi,0)$, and $(\pi,\pi)$ are suppressed by both hole and electron dopings. For $q=(\pi/2,\pi/2)$, the amplitude of $\rm{Im}\chi(\omega,q)$ is enhanced by initial hole doping and then suppressed by doping more holes; on the electron doping side, the amplitude of $\rm{Im}\chi(\omega,q)$ is suppressed monotonically by electron doping at $q=(\pi/2,\pi/2)$.

III.3 Spin excitation dispersion near M point

The left panel of Fig. 14 shows a simulated spin excitation dispersion near $(\pi,\pi)$ for momenta $q=(\pi-\delta,\pi)$ and interaction strengths $U/t=8$ at half filling (black) and $n=0.95$ (red). At low frequencies, the dispersion peaks are located at $\delta=0$. At higher frequencies, the dispersions become incommensurate and begin to disperse at a characteristic frequency $\omega_{c}$, creating a ‘Y’-shape. As the system is hole doped from $n=1$ to $n=0.95$, $\omega_{c}$ increases substantially.

Data are consistent with a resolution effect, where a single peak at the $(\pi,\pi)$ point (for $\omega=0$) with a constant width disperses linearly away from $(\pi,\pi)$ for higher frequencies. Uncertainties of the analytical continuation procedure and the dual fermion approximation make it difficult to rule out alternative explanations.

For comparison, the right panel of Fig. 14 shows experimental data reproduced from Ref. Chan et al., 2016. A $t=460$ meVDas (2012); Vishik et al. (2014) allows us to establish the right hand axis of the left hand frame in Fig. 14. We can see that $\omega_{c}=40$ meV for $n=1$ is smaller than the experimental observation of $\omega_{c}=60$ meV while $n=0.95$ is much larger with $\omega_{c}=150$ meV.

To further examine the spin excitation spectrum around the $(\pi,\pi)$ region, we extract the correlation length $\xi$ by fitting the peak structure along the $(\pi,\delta)$ direction, see Fig. 15, at $T/t=0.2$ as a function of doping for $U/t=8$, $5$, and $3$. We observe that a reduction in interaction strength or doping in either direction causes a sharp decrease of the correlation length, which in turn results in a broadening of the spin structures at low energies. $\xi$ is marginally larger on the electron doped side (a result of the broken particle-hole symmetry) and this difference increases with the interaction strength $U/t$.

IV Conclusions

In summary, we have presented results from the Dual Fermion approximation for the spin excitation spectrum in the Hubbard model. The results show striking qualitative similarities with experiments in cuprate materials. The Dual Fermion results have been crosschecked both against DCA and against DQMC results at high temperature. Away from the $\Gamma$ point, the spectrum is believed to be accurate up to a momentum-independent prefactor and has a precise momentum resolution. As the prefactor is doping, interaction, and presumably temperature dependent, and as control calculations from DCA and DQMC are restricted to small systems or high temperature, a precise quantitative comparison between numerics and experimental spin excitations is not possible at this time. In addition, a quantitative comparison near the $\Gamma$ point requires a solution that respects spin conservation.

Nevertheless, from the overall shape of the spin excitation dispersion and in particular the behavior between $M$, $X$, and $(\pi/2,\pi/2)$, we can conclude that an appropriate Hubbard interaction strength for modeling experimental data should be between $U/t=7$ and $U/t=8$.

The existing literature shows that, on the single-particle level, agreement between model systems calculations and experiment is remarkable and is often not just qualitative but quantitative. On the two-particle level, similar progress has been made primarily for quantities that are either local (such as the NMR probe) or contain a slowly varying matrix element. One important exception is the work by Jia and collaborators Jia et al. (2014) on simulations of the resonant $X$-ray scattering cross section. These authors presented $\chi^{\prime\prime}(q,\omega)$ obtained from DQMC calculations for the Hubbard model with $U=8t$ and $t^{\prime}/t=-0.3$ at the temperature $T=t/3\sim 1200K$, and reported quantitative agreement with magnetic neutron scattering measurements. Our DQMC results at these parameters are in quantitative agreement, as are our DF results (up to a rescaling); our studies of other $U$ values reinforce the conclusion that the Hubbard model with $U\sim 7-8t$ is a good description of cuprate physics. The continuous momentum resolution of the DF technique now shows that, in addition, many of the momentum dependent features seen in experiment are also observed in calculations.

V Acknowledgments

We thank T. Devereaux and A.-M. Tremblay for fruitful discussions. Dual Fermion methods development has been supported by the Simons collaboration on the many-electron problem, application to susceptibilities and neutron spectroscopy by NSF-DMR-1606348. JPFL acknowledges the support of the Natural Sciences and Engineering Research Council of Canada (NSERC), RGPIN-2017-04253. Computational resources were provided by TG-DMR130036, ACENET and Compute Canada.