Introduction

A key question in quantum materials research is whether or not the single-band Hubbard model describes the properties of the high-temperature (high-Tc) superconducting cuprates1,2,3. On the one hand, several studies have demonstrated a direct mapping between multi-orbital CuO models and effective single-band descriptions4,5,6,7. At the same time, quantum cluster methods have found evidence for a d-wave superconducting state6,8 in the Hubbard model, with a nonmonotonic Tc as a function of doping that resembles the dome found in real materials. On the other hand, a growing number of state-of-the-art numerical studies on extended Hubbard and t-J clusters have found evidence for stripe-ordered ground states for model parameters relevant to the cuprates9,10,11,12,13. While density matrix renormalization group (DMRG) simulations of multi-leg hole (h)-doped Hubbard ladders do obtain a superconducting ground state for nonzero values of the next-nearest-neighbor hopping \({t}^{{\prime} }\)14, its order parameter does not have the correct \({d}_{{x}^{2}-{y}^{2}}\) symmetry15 found in the cuprates16. Conversely, DMRG calculations for six- and eight-leg t-J cylinders obtain the correct order parameter but only on the electron (e)-doped side of the phase diagram12. These results cast significant doubt on the long-held belief that the Hubbard model describes the h-doped cuprates. Nevertheless, hope remains that it may capture the e-doped materials.

From an experimental perspective, charge-density-wave (CDW) correlations have been established as a ubiquitous feature of the underdoped cuprates17,18. Initially observed by inelastic neutron scattering in the form of intertwined spin and charge stripes19, short-range CDW correlations have now been reported in nearly all families of cuprates using scanning tunneling microscopy20,21 and resonant inelastic x-ray scattering (RIXS)22,23,24,25,26,27,28,29,30,31,32,33. Importantly, these CDW correlations persist up to high temperatures, particularly on the e-doped side of the phase diagram25,28,29,30.

Given their ubiquity, these CDW correlations must be accounted for by any proposed effective model for the cuprates. Evidence for charge modulations, both in the form of unidirectional stripe correlations or short-range CDW correlations, has now been found in a variety of finite-temperature quantum Monte-Carlo (QMC) simulations of the h-doped Hubbard model13,34,35,36,37,38. These simulations are generally restricted to high temperatures by the Fermion sign problem10,34,35 (except for very recent constrained path QMC calculations38) and focus on the h-doped model. However, the observed cuprate CDWs exhibit a significant electron-hole asymmetry. On the h-doped side, they can intertwine with spin-density modulations to form stripes while they coexist with uniform antiferromagnetic (AFM) correlations on the e-doped side25,28,30. These differences have raised questions on whether the e- and h-doped CDWs share a common origin29,30. Addressing these questions requires detailed calculations for the e-doped Hubbard model, which analyze the spin and charge correlations. To date, such calculations have not been performed.

Here, we provide insights into this question by studying the spin and charge correlations of the e-doped two-dimensional Hubbard model and contrasting them with the h-doped case, using the dynamical cluster approximation (DCA)39 and a nonperturbative QMC cluster solver40 (see Methods). Working on large (16 × 4) rectangular clusters that are wide enough to support large-period stripe correlations if they are present34, we vary the electron density 〈n〉 across both sides of the phase diagram to contrast the correlations in each case. Our calculations uncover robust two-component CDW correlations on the e-doped side, which consists of superimposed checkerboard (0.5, 0.5) (in the unit of 2π/a implied for all vectors in k space) and unidirectional QCDW = ( ± δc, 0) components. These CDW correlations appear to be decoupled from any stripe-like modulations of the spins and instead coexist with short-range antiferromagnetic correlations. This behavior is in direct contrast to the h-doped case, where we find evidence for fluctuating stripe correlations in both the charge and spin degrees of freedom34. Our results agree with experimental observations on the e-doped cuprates, including the observed doping dependence of QCDW. This supports the notion that the single-band Hubbard model captures the e-doped side of the high-Tc phase diagram.

Results

Figure 1 compares the static (ω = 0) charge χc(Q, 0) and spin χs(Q, 0) susceptibilities for the h- (〈n〉 = 0.8) and e-doped (〈n〉 = 1.2) Hubbard model with \(U/t=6,\, {t}^{{\prime} }/t=-0.2\), and varying temperature (see Supplementary Note 1 for error information of the lowest temperature results). In the h-doped case (Fig. 1a, c, e), unidirectional charge and spin stripes form as the temperature is lowered, consistent with prior finite-temperature studies34,35,38. These correlations manifest as incommensurate peaks in the static susceptibility at wave vectors Qc = ( ± δc, 0) and Qs = (0.5 ± δs, 0.5) for the charge and spin channels, respectively. For the spin channel in Fig. 1e, the dashed line shows a fit of the lowest temperature data using two Lorentzian functions. Figure 1c plots the charge correlations for the h-doped case along the (Qx, 0.5) direction, where we observe a weak double-peak structure emerging at the lowest accessible temperature. This modulation is weaker than the Qc structure in Fig. 1a, such that the charge correlations are predominantly stripe-like.

Fig. 1: Static spin and charge correlations in the doped single-band Hubbard model.
figure 1

a, c show the static charge susceptibility χc(Q, 0) along the (Qx, 0) and (Qx, 0.5) (in the unit of 2π/a) directions, respectively, for the h-doped system with \({t}^{{\prime} }/t=-0.2\) and 〈n〉 = 0.8. e shows the corresponding static spin susceptibility χs(Q, 0) along the Q = (Qx, 0.5) direction. The yellow dashed lines show incommensurate peaks obtained from fitting multiple Lorentzian functions plus a constant background to the βt = 6 data. (The constant contribution is not shown.) b, d, f show the corresponding results for the e-doped case with \({t}^{{\prime} }/t=-0.2\) and 〈n〉 = 1.2. These results were obtained using a 16 × 4 cluster, and the inverse temperatures β = 1/T are reported in units of 1/t.

Full size image

We observe qualitatively different correlations in the e-doped case shown in Fig. 1b, d, f. At high temperature (β ≤ 6/t), χc(Q, 0) has a single broad peak centered at q = (0, 0), which can again be decomposed into two incommensurate Lorentzian peaks centered at ± δc, indicative of a unidirectional charge stripe. As the temperature is lowered, these peaks sharpen and become discernible without fits while the q-independent background remains constant. The charge correlations also have a relatively temperature-independent (0.5, 0.5) component (Fig. 1e) of similar strength as the stripe-like charge correlations. In contrast to the h-doped case, we find no indication of a spin-stripe at this doping; the spin susceptibility has a single peak centered at (0.5, 0.5) for all accessible temperatures.

Comparing panels a and b, we see that the charge stripe correlations in the h- and e-doped systems develop differently as the temperature decreases. In the h-doped case (Fig. 1a), the incommensurate peaks grow while δc shifts to smaller values as the temperature decreases34. In the e-doped case, the incommensurate peaks grow (Fig. 1b) while the value of χc(Q, 0) near zone center is suppressed, resulting in well-defined peaks centered at ≈(±0.25, 0). In addition, the (0.5, 0.5) component is significantly stronger in the e-doped case (Fig. 1d). However, since the height of the incommensurate peaks in Fig. 1b does not appear to level off, the stripe correlations could dominate at lower temperatures.

The corresponding correlation functions in real-space, obtained at the lowest accessible temperatures (βt = 6 and 10 for the h- and e-doped cases, respectively) are plotted in Fig. 2. (The data for the h-doped case are regenerated from Fig. 1D of ref. 34 with error bars.) The charge correlations in the h-doped case (Fig. 2a) only show a stripe pattern. In contrast, the e-doped case (Fig. 2b) has a clear short-range (0.5, 0.5) checkerboard-like pattern near r = 0, superimposed over a stripe-like component. A similar (but weaker) pattern is also observed at a higher temperature in the determinantal quantum Monte-Carlo simulation for the e-doped Hubbard model (see Supplementary Note 2).

Fig. 2: Static spin and charge correlations in real-space.
figure 2

a χc(r, 0) for the e-doped system (\({t}^{{\prime} }=-0.2t,\langle n\rangle=1.2\)) at the lowest accessible inverse temperature βt = 10. b χc(r, 0) for the h-doped system (\({t}^{{\prime} }=-0.2t,\langle n\rangle=0.8\)) at the lowest accessible inverse temperature βt = 6. c, d show the staggered spin-spin correlations χs,stag(r, 0) for the h- and e-doped cases, respectively. + and − signs indicate the sign of the correlations whose absolute mean is larger than two standard errors. The dashed lines indicate the approximate nodes in the spin and charge stripe modulations.

Full size image

Figure 2c, d show the staggered spin correlation function for the h- and e-doped cases, respectively. Here, the blue region in the middle represents one AFM domain, while the red region on both sides indicates neighboring AFM domains with the opposite phase. (Note that the staggered spin correlation function that is plotted contains an additional negative sign on the B-sublattice as explained in the Methods section.) Consistent with Fig. 1 and as observed before34, the h-doped case has a clear stripe pattern whose period is roughly twice that of the charge modulations. In contrast, the e-doped model is dominated by short-range AFM correlations with only a single phase inversion appearing at longer distances. To summarize, the h-doped system has intertwined spin and charge stripe correlations, while the e-doped system manifests CDW correlation with stripe-like Q = (δc, 0) and checkerboard (0.5, 0.5) components and nearly uniform AFM spin correlations.

Figure 3 examines how the stripe component of the charge correlations develops for the e-doped case with density and \({t}^{{\prime} }\) for a fixed inverse temperature β = 8/t. Figure 3a shows χc(Q, 0) along the (Qx, 0) direction for various densities, while holding \({t}^{{\prime} }/t=-0.2\) fixed. At 〈n〉 = 1.125, the charge susceptibility has a single peak centered at Q = 0. With further electron doping, the peak splits into two well-defined peaks that become discernible without fitting for 〈n〉 ≥ 1.2. At the same time, the uniform background grows with doping due to the increased metallicity in the system35. Figure 3b presents the same quantity for various values of \({t}^{{\prime} }\), while fixing 〈n〉 = 1.2. At \({t}^{{\prime} }/t=-0.1\), the charge susceptibility again has a single broad peak at Q = 0. As \(|{t}^{{\prime} }|\) increases, the middle peak is suppressed, leading to the appearance of a pair of incommensurate peaks. At the same time, the uniform background remains almost unchanged.

Fig. 3: Evolution of the charge-density-wave correlations with model parameters.
figure 3

a, b show χc(Q, 0) for the e-doped case along the (Qx, 0) direction at β = 8/t. a shows the effect of different electron fillings 〈n〉 at \({t}^{{\prime} }=-0.2t\) while b shows the effect of varying \({t}^{{\prime} }\) for fixed 〈n〉 = 1.2. c The evolution of the incommensurate CDW peak Qcdw = (δc, 0) with doping, obtained from fitting the spectra in a with three Lorentzian functions and a constant background. For comparison, we also plot the measured values of Qcdw extracted from RIXS experiments28. d The evolution of the incommensurate CDW peak as a function of \({t}^{{\prime} }\), obtained from fitting the data shown in b. All results were obtained on a 16 × 4 cluster with U = 6t and βt = 8 except for the βt = 10 results in c. The error bars in c, d are standard deviations of errors from fitting the data to three Lorentzian functions plus a constant background.

Full size image

To assess these trends quantitatively, we fit the curves in Fig. 3a with three Lorentzian functions centered at Qx = 0 and ± Qc plus a constant background to extract the wave vector δc = Qc. As shown in Fig. 3c, δc increases approximately linearly as a function of the doping ρ = 〈n〉 − 1. This trend agrees with experimental observations for the e-doped cuprates28; however, the observed Qcdw(ρ) curve is shifted to lower doping levels relative to our results. We also extracted δc as a function of \({t}^{{\prime} }\), as shown in Fig. 3d, where we find that δc linearly shifts to larger values with \(|{t}^{{\prime} }|\). Extrapolating the experimental data in Fig. 3c to ρ = 0.2 yields Qcdw ≈ 0.32 r.l.u., which corresponds to \(|{t}^{{\prime} }/t|\, \approx \, 0.4\) in our model. This value is much larger than the typical values used to model the e-doped cuprates using the single-band model.

RIXS experiments have found that the CDW in the e-doped cuprates has a significant dynamical component25,28,30. To compare with these measurements, we show in Fig. 4 the dynamical spin and charge structure functions S(Q, ω) and N(Q, ω) obtained using a parameter-free differential evolution analytic continuation (DEAC) algorithm41. (A comparison of the DEAC results to those obtained with more conventional techniques is provided in Supplementary Note 3.) The dynamic spin structure factor S(Q, ω) displays the typical persistant antiferromagnetic paramagnon spectrum obtained with other QMC methods42, with large spectral weight at Q = (0.5, 0) at an energy near ω = t and a downward dispersion towards Q = 0. The dynamic charge structure factor exhibits similar behavior with a large spectral weight near the zone edge and a downward dispersion towards the zone center. Still, the spectral weight is concentrated at higher energies.

Fig. 4: Dynamical spin and charge structure factors for the e-doped model.
figure 4

a, b show the dynamical spin S(Q, ω) and charge N(Q, ω) structure factors, respectively, along the Q = (Qx, 0) direction. Results are shown for an e-doped model with \({t}^{{\prime} }/t=-0.2\) and 〈n〉 = 1.2, obtained on a 16 × 4 cluster with U = 6t and β/t = 10. c shows the sum of the two as a crude approximation for Cu L-edge RIXS spectra. The black circles and diamonds indicate the locations of the maxima in S(Q, ω) and N(Q, ω), respectively. All panels are plotted on the same color scale, as indicated on the right.

Full size image

As an approximation of the predicted RIXS intensity, Fig. 4c plots a superposition of S(Q, ω) and N(Q, ω). Due to the vastly different energy scales in the spin and charge correlations, one sees two distinct upward dispersing branches, i.e., a low-energy branch due to the spin correlations and a high-energy branch due to the charge correlations. Since the spin correlations have a much larger amplitude, the low-energy branch has a stronger overall intensity, while the higher energy (charge) branch is muted. The behavior for the dynamical correlations near zone center agrees well with the data reported in ref. 30. Our results for the higher energy charge excitations call for future RIXS measurements in this regime; however, the high-energy portion of the charge excitations may overlap with the intra-atomic orbital excitations on the Cu atom (the so-called dd-excitations), which are commonly observed at the Cu L-edge43.

Discussion

The correlations in the e-doped case cannot be attributed solely to the weak-coupling Lindhard physics44. For example, as discussed in Supplementary Note 4, if we adjust the effective U in a random-phase approximation (RPA) to match the charge correlations χc(Qc, 0) then the predicted correlations at (0.5, 0.5) are much stronger than those reported here. Similarly, we do not resolve any peak structure in χs(Q, 0) along the (Qx, 0) direction as one would expect in such a weak-coupling framework. While the peak positions predicted by RPA are close to the values reported here for 〈n〉 = 1.2, the resulting temperature, doping, and \({t}^{{\prime} }\) evolution is inconsistent with the observed behavior (see Supplementary Note 4).

We have demonstrated that the CDW correlations observed in DCA simulations of the single-band Hubbard model have a pronounced particle-hole asymmetry. Our results for the e-doped system are also in qualitative agreement with RIXS experiments on Nd2−xCexCuO425,30; however, there are some notable quantitative differences. For example, our predicted δc(ρ) dependence is shifted to higher electron doping in comparison to experiment. This discrepancy may be related to challenges in determining the carrier concentration in the CuO2 planes. ARPES measurements of Pr1.3−xLa0.7CexCuO4−δ (PLCCO), for example, have suggested that the doped electron concentration of CuO2 plane can be larger than the Ce concentration x by up to 0.08 e/Cu, depending on the annealing method45. This discrepancy is comparable to the shift observed in Fig. 3c. Adjusting the value of \({t}^{{\prime} }\) could also partially account for this difference; single-band fits to the measured Fermi surface of PLCCO estimate \(|{t}^{{\prime} }/t|\approx 0.34 - 0.4\). Finally, the periodicity of the charge modulations may be affected by the DCA mean-field34. In the h-doped case, DCA and DQMC predict different stripe periods for the same model parameters. Nevertheless, our calculations reproduce the qualitative doping dependence observed in the real materials and support the notion that the single-band Hubbard model captures the physics of the e-doped cuprates.

Methods

The model

We consider the single-band Hubbard model on a two-dimensional square lattice. The Hamiltonian is

$$H=-\mathop{\sum}\limits_{{{{{{{{\bf{i}}}}}}}},{{{{{{{\bf{j}}}}}}}},\sigma }{t}_{{{{{{{{\bf{i}}}}}}}},{{{{{{{\bf{j}}}}}}}}}^{}{c}_{{{{{{{{\bf{i}}}}}}}},\sigma }^{{{{\dagger}}} }{c}_{{{{{{{{\bf{j}}}}}}}},\sigma }^{}-\mu \mathop{\sum}\limits_{{{{{{{{\bf{i}}}}}}}},\sigma }{n}_{{{{{{{{\bf{i}}}}}}}},\sigma }+U\mathop{\sum}\limits_{{{{{{{{\bf{i}}}}}}}}}{n}_{{{{{{{{\bf{i}}}}}}}},\uparrow }{n}_{{{{{{{{\bf{j}}}}}}}},\downarrow }.$$
(1)

Here, \({c}_{{{{{{{{\bf{i}}}}}}}},\sigma }^{{{{\dagger}}} }\) (\({c}_{{{{{{{{\bf{i}}}}}}}},\sigma }^{}\)) creates (annihilates) a spin-σ ( = ↑, ↓) electron at site i = a(ix, iy), where a = 1 is the lattice constant, ti,j is the hopping integral between sites i and j, μ is the chemical potential, and U is the Hubbard repulsion. To model the cuprates, we set ti,j = t and \({t}^{{\prime} }\) for nearest and next-nearest-neighbor hopping, respectively, and zero otherwise, and take U/t = 6 unless stated otherwise.

The dynamical cluster approximation

We study the model in Eq. (1) using DCA8,39,46 as implemented in the DCA++ code47. The DCA represents the infinite bulk lattice in the thermodynamic limit by a finite-size cluster embedded in a self-consistent dynamical mean-field. The intra-cluster correlations are treated exactly, while the mean-field approximates the inter-cluster degrees of freedom. We use rectangular Nc = 16 × 4 clusters that are large enough to support spin and charge stripe correlations if they are present in the model34.

Assuming short-ranged correlations, the self-energy Σ(k, iωn) can be approximated by the cluster self-energy Σ(K, iωn), where K are the cluster momenta. We obtain the coarse-grained single-particle Green function

$$\bar{G}({{{{{{{\bf{K}}}}}}}},\, {{{{{{{\rm{i}}}}}}}}{\omega }_{n}) =\frac{{N}_{{{{{{{{\rm{c}}}}}}}}}}{N}\mathop{\sum}\limits_{{{{{{{{{\bf{k}}}}}}}}}^{{\prime} }}G({{{{{{{\bf{K}}}}}}}}+{{{{{{{{\bf{k}}}}}}}}}^{{\prime} },\, {{{{{{{\rm{i}}}}}}}}{\omega }_{n})\\ =\frac{{N}_{{{{{{{{\rm{c}}}}}}}}}}{N}\mathop{\sum}\limits_{{{{{{{{{\bf{k}}}}}}}}}^{{\prime} }}\frac{1}{{{{{{{{\rm{i}}}}}}}}{\omega }_{n}+\mu -\varepsilon ({{{{{{{\bf{K}}}}}}}}+{{{{{{{{\bf{k}}}}}}}}}^{{\prime} })-\Sigma ({{{{{{{\bf{K}}}}}}}},\, {{{{{{{\rm{i}}}}}}}}{\omega }_{n})},$$
(2)

by averaging the lattice Green function G(k, iωn) over the N/Nc momenta \({{{{{{{{\bf{k}}}}}}}}}^{{\prime} }\) in a patch about the cluster momentum K. (N and Nc are the number of site in the lattice and cluster, respectively.) In this way, the bulk problem is reduced to a finite-size-cluster problem.

We solve the cluster problem with the continuous-time, auxiliary-field quantum Monte-Carlo algorithm (CT-AUX)40,48. The expansion order of the CT-AUX QMC algorithm is typically 500–1500, depending on the temperature and value of \({t}^{{\prime} }\). Depending on the value of the average fermion sign for a given parameter set, we measure 1−5 × 108 samples for the correlation functions. Six to eight iterations of the DCA loop are typically needed to obtain good convergence for the DCA cluster self-energy.

To study the spin and charge correlations, we measure the static (ω = 0) real-space spin-spin correlation function

$${\chi }_{{{{{{{{\rm{s}}}}}}}}}({{{{{{{\bf{r}}}}}}}},\, 0)=\frac{1}{N}\mathop{\sum}\limits_{{{{{{{{\bf{i}}}}}}}}}\int\nolimits_{0}^{\beta }\left\langle {\hat{S}}_{{{{{{{{\bf{r}}}}}}}}+{{{{{{{\bf{i}}}}}}}}}^{z}(\tau ){\hat{S}}_{{{{{{{{\bf{i}}}}}}}}}^{z}(0)\right\rangle d\tau$$
(3)

and density-density correlation function

$${\chi }_{{{{{{{{\rm{c}}}}}}}}}({{{{{{{\bf{r}}}}}}}},\, 0)=\frac{1}{N}\mathop{\sum}\limits_{{{{{{{{\bf{i}}}}}}}}}\int\nolimits_{0}^{\beta }\left[\langle {n}_{{{{{{{{\bf{r}}}}}}}}+{{{{{{{\bf{i}}}}}}}}}(\tau ){n}_{{{{{{{{\bf{i}}}}}}}}}(0)\rangle -\langle {n}_{{{{{{{{\bf{r}}}}}}}}+{{{{{{{\bf{i}}}}}}}}}(\tau )\rangle \langle {n}_{{{{{{{{\bf{i}}}}}}}}}(0)\rangle \right]d\tau .$$
(4)

Here, r = (rx, ry)a and i = (ix, iy)a are lattice vectors and \({\hat{S}}_{{{{{{{{\bf{i}}}}}}}}}^{z}=\frac{1}{2}({n}_{{{{{{{{\bf{i}}}}}}}},\uparrow }-{n}_{{{{{{{{\bf{i}}}}}}}},\downarrow })\) and \({n}_{{{{{{{{\bf{i}}}}}}}}}={\sum }_{\sigma }{c}_{{{{{{{{\bf{i}}}}}}}},\sigma }^{{{{\dagger}}} }{c}_{{{{{{{{\bf{i}}}}}}}},\sigma }^{}\) are the local spin-z and density operators of the effective cluster problem. The correlation functions χc,s(r, 0) are measured directly by the cluster solver. The static momentum-space susceptibilities χc(Q, 0) and χs(Q, 0) are obtained by a Fourier transform. We also plot the staggered spin correlation function \({\chi }_{{{{{{{{\rm{s}}}}}}}},{{{{{{{\rm{stag}}}}}}}}}({{{{{{{\bf{r}}}}}}}})={(-1)}^{({r}_{x}+{r}_{y})}S({{{{{{{\bf{r}}}}}}}},\, 0)\) to highlight the relative phases of the antiferromagnetic domains.

Analytic continuation

The dynamical charge N(Q, ω) and spin S(Q, ω) structure factors are obtained from the imaginary part of the charge and spin susceptibilities using

$$N({{{{{{{\bf{Q}}}}}}}},\, \omega )=\frac{{{{{{{{\rm{Im}}}}}}}}{\chi }_{{{{{{{{\rm{c}}}}}}}}}({{{{{{{\bf{Q}}}}}}}},\, \omega )}{1-{{{{{{{{\rm{e}}}}}}}}}^{-\beta \omega }},$$
(5)

and

$$S({{{{{{{\bf{Q}}}}}}}},\, \omega )=\frac{{{{{{{{\rm{Im}}}}}}}}{\chi }_{{{{{{{{\rm{s}}}}}}}}}({{{{{{{\bf{Q}}}}}}}},\, \omega )}{1-{{{{{{{{\rm{e}}}}}}}}}^{-\beta \omega }}.$$
(6)

The real frequency susceptibility is related to the imaginary time susceptibility by the integral equation

$${\chi }_{{{{{{{{\rm{s,c}}}}}}}}}({{{{{{{\bf{Q}}}}}}}},\, \tau )=\int\nolimits_{0}^{\infty }\frac{d\omega }{\pi }\frac{{{{{{{{{\rm{e}}}}}}}}}^{-\tau \omega }+{{{{{{{{\rm{e}}}}}}}}}^{-(\beta -\tau )\omega }}{1-{{{{{{{{\rm{e}}}}}}}}}^{-\beta \omega }}{{{{{{\mathrm{Im}}}}}}}\,{\chi }_{{{{{{{{\rm{s,c}}}}}}}}}({{{{{{{\bf{Q}}}}}}}},\, \omega ).$$
(7)

Inverting this relationship to obtain \({{{{{{{\rm{Im}}}}}}}}\chi ({{{{{{{\bf{Q}}}}}}}},\omega )\) is an ill-posed problem. We used three independent methods to perform the analytic continuation to gauge our relative confidence in the various features. These include the method of Maximum Entropy49, a parameter-free differential evolution algorithm41, and stochastic optimization50. The results obtained using the differential evolution algorithm are shown in the main text, while the remaining results are provided in Supplementary Note 3.