Introduction

Angle-resolved photoemission spectroscopy (ARPES) experiments show, quite universally in various cuprates1,2,3,4,5,6,7,8,9, a high energy anomaly in the form of a waterfall-like structure. The onset of these waterfalls is between 100 and 200 meV, at considerably higher energy than the distinctive low-energy kinks10,11,12, and they end at even much higher binding energies around  ~1 eV3. Also, their structure is qualitatively very different: an almost vertical and smeared-out waterfall and not a kink from one linear dispersion to another that is observed at lower binding energies. Akin waterfalls have been reported most recently in nickelate superconductors13,14, there starting at around 100 meV. This finding puts the research focus once again on this peculiar spectral anomaly. With the close analogy between cuprates and nickelates15,16 the observation of waterfalls in nickelates gives fresh hope to eventually understand the physical origin of the waterfalls.

Quite similar as for superconductivity, various theories have been suggested for waterfalls in cuprates, including: the coupling to hidden fermions17, the proximity to quantum critical points18, and multi-orbital physics19,20. The arguably most widespread theoretical understanding is the coupling to a bosonic mode, such as phonons21 or spin fluctuations (including spin polarons)22,23,24,25,26. Here, in contrast to the low energy kinks, the electron-phonon coupling appear a less viable origin for waterfalls, simply because the phonon energy is presumably too low. Also the spin coupling J in cuprates is below 200 meV, which however might concur with the onset of the waterfall. But, its ending at 1 eV is barely conceivable from a spin fluctuation mechanism, as it is almost an order of magnitude larger than J. Even the possibility that waterfalls are matrix element effects that are not present in the actual spectral function has been conjectured27.

The simplest model for both, superconducting cuprates and nickelates, is the one-band Hubbard model for the Cu(Ni) 3\({d}_{{{\rm{x}}}^{2}-{{\rm{y}}}^{2}}\) band. In the case of cuprates, the more fundamental model might be the Emery model which also includes the in-plane oxygen orbitals. However, with some caveats such as doping-depending hopping parameters, a description by the simpler Hubbard model is qualitatively similar28,29. In the case of nickelates, these oxygen orbitals are lower in energy, but instead rare earth 5d orbitals become relevant and cross the Fermi level30,31,32,33,34. Still, the simplest description is that of a one-band Hubbard model plus largely detached 5d pockets35,36. This simple description is confirmed by ARPES that shows no additional Fermi surfaces and only 5d A pockets for SrxLa(Ca)1−xNiO213,14.

In this paper, we show that waterfalls naturally emerge when a Hubbard band splits off from the central quasiparticle band. This splitting-off is sufficient for, and even necessitates a waterfall-like structure. Using dynamical mean-field theory (DMFT)37 we can exclude that spin fluctuations are at work, as the feedback of these on the spectrum would require extensions of DMFT38. For the doped model, the waterfall prevails in a large range of interactions, which explains its universal occurrence in cuprates and nickelates. A one-on-one comparison of experimental spectra to those of the Hubbard model with ab initio determined parameters for cuprates and nickelates also shows good agreement. Previous papers pointing toward a similar mechanism9,23,39,40,41,42,43,44,45 have, to the best of our knowledge, been quite general, without the more detailed analysis or understanding which the present paper provides. Among others, Macridin et al.23 noted a positive slope of the DMFT self-energy at intermediate frequencies, but eventually concluded that spin fluctuations lead to waterfalls; Moritz et al.9,41 emphasized that waterfalls simply connect Hubbard and quasiparticle bands; and Sakai et al.40 pointed out the importance of the quasiparticle renormalization and vicinity to a Mott transition, advocating the momentum dependence of the self-energy. All these publications use similar numerical quantum Monte Carlo simulations for the Hubbard model either directly for a finite  lattice or for lattice extensions of DMFT. In such calculations, it is difficult to track down whether spin fluctuations23 or other mechanisms9,40,41 are in charge.

Results

Waterfalls in the Hubbard model

Neglecting matrix elements effects, the ARPES spectrum at momentum k and frequency ω is given by the imaginary part of Green’s function, i.e., the spectral function

$$A({\bf{k}},\omega )=-\frac{1}{\pi }{\rm{Im}}\underbrace{\frac{1}{\omega -{\varepsilon }_{{\bf{k}}}-{{\Sigma }}(\omega )+i\delta }}_{\equiv G({\bf{k}},\omega )}.$$
(1)

Here, δ is an infinitesimally small broadening and εk the non-interacting energy-momentum dispersion. For convenience, we set the chemical potential μ ≡ 0. The non-interacting εk is modified by electronic correlations through the real part of the self-energy ReΣ(ω) while its imaginary part describes a Lorentzian broadening of the poles (excitations) of Eq. (1) at

$$\omega={\varepsilon }_{{\bf{k}}}+\,{\mathrm{Re}}{{\Sigma }}(\omega ).$$
(2)

Please note that we have here omitted the momentum dependence of the self-energy which holds for the DMFT approximation, while non-local correlations can lead to a k-dependent self-energy. This k-dependence can, e.g., arise from spin fluctuations and lead to a pseudogap. In Supplementary Figs. 1 and 2, we compare DMFT to an extension of DMFT, the dynamical vertex approximation (DΓA46) that includes such non-local correlations. We show that although non-local correlations can further corroborate the presence of waterfall-like structures, their underlying origin remains tied to local correlations.

Figure 1 shows our DMFT results for the Hubbard model on the two-dimensional square lattice at half-filling with only the nearest neighbor hopping t. We go from the weakly correlated regime (left) all the way to the Mott insulator (right). The spectrum then evolves from the weakly broadened and renormalized local density of states (LDOS) resembling the non-interacting system in panel (a) to the Mott insulator with two Hubbard bands at  ± U/2 in panel (d). In-between, in panel (c), we have the three-peak structure with both Hubbard bands and a central, strongly-renormalized quasiparticle peak in-between; the hallmark of a strongly correlated electron system that DMFT so successfully describes37. Panel (b) is similar to panel (c), with the difference being that the Hubbard bands are not yet so clearly separated. This is the situation where waterfalls emerge in the k-resolved spectrum shown in Fig. 1(j).

Fig. 1: LDOS, graphical solution of the pole equation, and spectral function.
figure 1

Top (ad): DMFT LDOS for the two-dimensional Hubbard model with nearest neighbor hopping t = 0.3894 eV at half-filling and—from left to right–increasing U. The shaded areas denote the filled states. Temperature is room temperature T = t/15 except for the last column (U = 15t) where T = t. Energies are in units of eV. Middle (eh): Graphical solution for the poles of the Green’s function in Eq. (2) as the crossing point (colored circles) between εk + ReΣ(ω) (solid lines in three colors for the three k points indicated by vertical lines in the bottom panel) and ω (black dashed line); the colored dashed lines denote εk for the same three momenta. Note that in the insulating case [panel (h)] the dashed lines are shifted by U/2 and the results scaled by a factor of 1/11. Also shown is the imaginary part of the self-energy (light gray; right y-axis). Bottom (il): k-resolved spectral function A(kω) along the nodal direction Γ = (0, 0) to M = (ππ), showing a waterfall for U = 5t in panel (j). Also plotted are energy distribution curve maxima (EDC MAX, gray circles) and momentum distribution curve maxima (MDC MAX, gray squares) defined as the maxima of A(kω) as a function of ω and k, respectively.

Waterfalls from \(\partial {\rm{Re}}{{\Sigma }}(\omega )/\partial \omega=1\)

To understand the emergence of this waterfall feature, we solve in Fig. 1(e–h) the pole equation (2) graphically. That is, we plot the right-hand side of Eq. (2), εk + ReΣ(ω), for three different momenta (colored solid lines), with each momentum indicated by a vertical line of the same color in panels (i-l). The left-hand side of Eq. (2), ω, is plotted as a black dashed line. Where they cross, indicated by circles in panels (i-l), we have a pole in the Green’s function and a large spectral contribution.

For U = 2t (leftmost column), the excitations are essentially the same as for the non-interacting system ω ≈ εk, with the self-energy only leading to a minor quasiparticle renormalization and broadening. In the Mott insulator at large U and zero temperature, on the other hand, Σ(ω) = U2/(4ω). Finite T and hopping t regularize this 1/ω pole seen developing in Fig. 1(h), but then turning into a steep positive slope of ReΣ(ω) around ω = 0. Instead of a delta-function, ImΣ(ω) becomes a Lorentzian (light gray curve; note the rescaling). Thus, while there is an additional pole-like solution around ω = 0, it is completely smeared out.

Now for U = 8t in Fig. 1(g) we have for large ω the same pole-like behavior as in the Mott insulator, though of course with a smaller U2 prefactor. On the other hand, at small frequencies ω we have the additional quasiparticle peak which corresponds to a negative slope ∂ReΣ(ω)/∂ωω=0 < 0 that directly translates to the quasiparticle renormalization or mass enhancement m*/m = 1 − ∂ReΣ(ω)/∂ωω=0 > 1. Altogether, εk + ReΣ(ω) must hence have the form seen in Fig. 1(g): we have one solution of Eq. (2) at small ω in the range of the negative, roughly linear ReΣ(ω), which corresponds to the quasiparticle excitations. We have a second solution at large ω, where we have the 1/ω self-energy as in the Mott insulator, which corresponds to the Hubbard bands. For a chosen k, there is a third crossing in-between, where the self-energy crosses from the Mott like 1/ω to the quasiparticle like—ω behavior. Here, the self-energy has a positive slope. This pole is however not visible in A(kω) [Fig. 1(k)], simply because the smearingImΣ(ω) is very large. It would not be possible to see it in ARPES.

However, numerically, one can trace it as the maximum in the momentum distribution curve (MDC), i.e, maxkA(kω) along Γ to M, shown as squares in Fig. 1k. This MDC shows an S-like shape since the positive slope of ReΣ(ω) in this intermediate ω range is larger than one (dashed black line). Consequently, for εk at the bottom of the band (orange and blue lines) this third pole in panel (g) is close to the quasiparticle pole, while for εk closer to the Fermi level μ ≡ 0 (green line) it is close to the pole corresponding to the Hubbard band.

For the smaller U of Fig. 1(e), on the other hand, Σ is small and thus also the positive slope in the intermediate ω range must be smaller than one (dashed black line). Together with the continuous evolution of the self-energy from (e) to (j), this necessitates that for some Coulomb interaction in-between, the slope close to the inflection point in-between Hubbard and quasiparticle band equals one: \(\partial {\rm{Re}}{{\Sigma }}(\omega )/\partial \omega=1\). That is the case for U ≈ 5t shown in Fig. 1(k).

Now there is only one pole for each momentum. For the momentum closest to the Fermi level μ (green line), it is in the quasiparticle band where ReΣ(ω) ~ − ω at small ω. When we reduce εk, i.e., shift the εk + ReΣ(ω) curve down, there is one momentum (blue curve) where the crossing is not in the quasiparticle band nor in the Hubbard band but in the crossover region between the two, with the positive slope of ReΣ(ω). As this slope is one, the blue and black dashed lines are close to each other in a large energy region. That is, we are close to a pole for many different energies ω. Given the finite imaginary part of the self-energy, we are thus within reach of an actual pole. Consequently, we get a waterfall in Fig. 1(j) with spectral weight in a large energy range for this blue momentum. Finally, for εk’s at the bottom of the band (orange line), the crossing point is in the lower Hubbard band. Altogether this leads to a waterfall as a crossover from the quasiparticle to the Hubbard band.

Doped Hubbard model

Next, we turn to the doped Hubbard model in Fig. 2. The main difference is that now for the U →  limit, we do not get a Mott insulator, but keep a strongly correlated metal. As a consequence the U-range where we have waterfall-like structures is much wider, which explains that they are quite universally observed in cuprates and nickelates. Strictly speaking, an ideal vertical waterfall again corresponds mathematically to a slope one close to the inflection point of ReΣ(ω). This is the case for U ≈ 8t in Fig. 2c. However, with the much slower evolution with U at finite doping, we have a large U range with waterfall-like structures, first at small U in the form of moderate slopes as in panel (b), and then for large U in form of an S-shape-like structure as in panel (d) that are akin to waterfalls. The survival of the S-shape structure even at U = 15t strongly suggests that it extends up to U → .

Fig. 2: Spectral functions and second derivatives of MDCs for 20% hole doping.
figure 2

Top (ad): Momentum-resolved spectral function for the two-dimensional Hubbard model with nearest neighbor hopping t = 0.3894 eV—from left to right—increasing U, room temperature T = t/15, and 20% hole doping, showing waterfall-like structures in a large interaction range. Bottom (eh): Second momentum derivative of the spectral function, which is usually employed in experiments to better visualize the waterfalls. Besides the MDC maxima (MDC MAX, gray squares), we also plot the minima of the second derivative of the MDCs (SD MDC MIN, gray diamonds).

Figure 2 (e-h) shows the second derivative (SD) of the MDC, i.e., ∂2A(kω)/∂k2 along the momentum line Γ to M. This SD MDC is usually used in an experiment to better visualize the waterfalls; and indeed we see in Fig. 2e–h that the waterfall becomes much more pronounced and better visible than in the spectral function itself.

Connection to nickelates and cuprates

Let us finally compare our theory for waterfalls to ARPES experiments for nickelates and cuprates. We here refrain from adjusting any parameters and use the hopping parameters of the Hubbard model that have been determined before ab initio by density functional theory (DFT) for a one-band Hubbard model description of Sr0.2La0.8NiO235, La2−xSrxCuO447, and Bi2Sr2CuO6 (Bi2201)48. Similarly, constrained random phase approximation (cRPA) results are taken for the interaction U34,48. CRPA U values of ref. 47 are rounded up similarly to nickelates as in refs. 34,35 to mimic the frequency dependence of U. The parameters are listed in the captions of Figs. 3 and 4.

Fig. 3: DFT+DMFT calculations of waterfalls in nickelates.
figure 3

Waterfalls in the one-band DMFT spectrum (a; top) and its second derivative (c; bottom) for Sr0.2La0.8NiO2, compared to experiment13 (exp, golden circles), together with the MDC maxima (MDC MAX, gray squares), and the minima of the second derivative of the MDCs (SD MDC MIN, gray diamonds). In the right column (b; d), we added a broadening Γ = 1 eV to the DMFT self-energy to mimic disorder effects. The ab initio determined parameters of the Hubbard model for nickelates are35: t = 0.3894 eV, \({t}^{{\prime} }=-0.25t,{t}^{{\prime\prime} }=0.12,U=8t\), 20% hole doping, and we take a sufficiently low temperature T = 100/t.

Figure 3 compares the waterfall structure in the one-band Hubbard model for nickelates to the ARPES experiment13 (for waterfalls in nickelates under pressure cf.45). The qualitative agreement is very good. Quantitatively, the quasiparticle renormalization is also well described without free parameters. The onset of the waterfall is at a similar binding energy as in ARPES, though a bit higher, and at a momentum closer to Γ.

This might be due to different factors. One is that nickelate films still have a high degree of disorder, especially stacking faults. We can emulate this disorder by adding a scattering rate Γ to the imaginary part of the self-energy. For Γ = 1 eV, we obtain Fig. 3b, d which is on top of experiment also for the waterfall-like part of spectrum, though with an adjusted Γ. Indeed, we think that this Γ is a bit too large, but certainly disorder is one factor that shifts the onset of the waterfall to lower binding energies. Other possible factors are (i) the ω-dependence of U(ω) in cRPA which we neglect, and (ii) surface effects on the experimental side to which ARPES is sensitive. Also (iii) a larger U would according to Fig. 2 result in an earlier onset of the waterfall. At the same time, it would however also increase the quasiparticle renormalization which is, for the predetermined U, in good agreement with the experiment.

Figure 4 compares the DMFT spectra of the Hubbard model to the energy-momentum dispersions extracted by ARPES for two cuprates. Panels (a-d;g-j) show the comparison for four different dopings x of La2−xSrxCuO4. Again, we have a good qualitative agreement including the change of the waterfall from a kink-like structure at large doping x = 0.3 in panels (d,j) to a more S-like shape at smaller doping x = 0.12 in panels (a,g). The same doping dependence is also observed for Bi2201 from panels (f,l) to (e,k). Note, lower doping effectively means stronger correlations, similar to increasing U in Fig. 2, where we observe the same qualitative change of the waterfall. Altogether this demonstrates that even changes in the form of the waterfall from kink-like to vertical waterfalls to S-like shape can be explained. Quantitatively, we obtain a very good agreement at larger dopings, while at lower dopings there are some quantitative differences. However, please keep in mind that we did not fit any parameters here.

Fig. 4: DFT+DMFT calculations of waterfalls in cuprates.
figure 4

Waterfalls in the DMFT spectrum (top) and its second derivative (bottom) for (ad;gj) La2−xSrxCuO4 at four different x and (e, f;k, l) Bi2Sr2CuO6 (Bi2201) for x = 0.16 and x = 0.26 hole doping compared to experiment (golden circles, ref. 2), (green circles, ref. 5). Also shown are the MDC maxima (MDC MAX, gray squares), and the minima of the second derivative of the MDCs (SD MDC MIN, gray diamonds). The ab initio determined parameters of the one-band Hubbard model for La2−xSrxCuO4 are47: t = 0.4437 eV, \({t}^{{\prime} }=-0.0915t,{t}^{{\prime\prime} }=0.0467t,U=7t\). CRPA U values of ref. 47 are rounded up similarly to nickelates as in refs. 34,35 to mimic the frequency dependence of U. Those for Bi2201 are48: t = 0.527 eV, \({t}^{{\prime} }=-0.27t,{t}^{{\prime\prime} }=0.08t,U=8t\) (rounded). Again we set T = 100/t.

Discussion

Umbilical cord metaphor

At small interactions U all excitations or poles of the Green’s function are within the quasiparticle band; for very large U and half-filling all are in the Hubbard bands; and for large, but somewhat smaller U we have separated quasiparticle and Hubbard bands. We have proven that there is a qualitatively distinct fourth “waterfall" parameter regime. Here, the Hubbard band is not yet fully split off from the quasiparticle band, and we have a crossover in the spectrum from the Hubbard to the quasiparticle band in the form of a waterfall. This waterfall must occur when turning on the interaction U and is, in the spirit of Ockham, a simple explanation of the waterfalls observed in cuprates, nickelates, and other transition metal oxides. Even the change from a kink-like to an actual vertical waterfall to an S-like shape with increasing correlations agrees with the experiment.

As Supplementary Movie 1, we provide a movie of the spectrum evolution with increasing U. Figuratively, we can call this evolution the “birth of the Hubbard band", with the quasiparticle band being the “mother band" and the Hubbard band the “child band". The waterfall is then the “umbilical cord" connecting the “mother band" and “child band" before the latter becomes fully disconnected from the former. As a matter of course, such metaphors are never perfect. Here, e.g., we rely on the time axis being identified with increasing U. However, one could also interpret it vice versa, that is, as the quasiparticle band disconnects from the Hubbard band as U decreases.

Methods

In this section, we outline the model and computational methods employed. The two-dimensional Hubbard model for the 3\({d}_{{{\rm{x}}}^{2}-{{\rm{y}}}^{2}}\) band reads

$${\mathcal{H}}=\sum _{ij\sigma }{t}_{ij}{\hat{c}}_{i\sigma }^{\dagger }{\hat{c}}_{j\sigma }+U\sum _{i}{\hat{n}}_{i\uparrow }{\hat{n}}_{i\downarrow }.$$
(3)

Here, tij denotes the hopping amplitude from site j to site i, which we restrict to nearest neighbor t, next-nearest neighbor \({t}^{{\prime} }\), and next-next-nearest neighbor hopping t; \({\hat{{c}_{i}}}^{\dagger }\) (\(\hat{{c}_{j}}\)) are fermionic creation (annihilation) operators, and σ marks the spin; \({\hat{n}}_{i\sigma }={\hat{c}}_{i\sigma }^{\dagger }{\hat{c}}_{i\sigma }\) are occupation number operators; U is the Coulomb interaction.

DMFT calculations were done using w2dynamics49 which uses quantum Monte Carlo simulations in the hybridization expansion50. For the analytical continuation, we employ maximum entropy with the chi2kink method as implemented in the ana_cont code51.