Collective modes in the paramagnetic phase of the Hubbard model

TL;DR

This paper investigates the charge excitation spectra and collective modes of the Hubbard model in the paramagnetic phase, revealing a continuum, a gapless zero-sound mode, and a high-frequency mode, with detailed analysis of their dependence on interaction strength and temperature.

Contribution

It provides a detailed analytical and numerical study of charge excitations and collective modes in the Hubbard model using slave-boson methods, including an approximate expression for the high-frequency mode.

Findings

Identification of a gapless zero-sound mode with anisotropic velocity

Derivation of an analytical expression for the high-frequency mode

Finite temperature phase diagrams showing suppression of magnetic instabilities

Abstract

The charge dynamical response function of the Hubbard model is investigated on the square lattice in the thermodynamical limit. The obtained charge excitation spectra consist of a continuum, a gapless collective mode with anisotropic zero-sound velocity, and a correlation induced high-frequency mode at $\omega\approx U$. The correlation function is calculated from Gaussian fluctuations around the paramagnetic saddle-point within the Kotliar and Ruckenstein slave-boson representation. Its dependence on the on-site Coulomb repulsion $U$ and density is studied in detail. An approximate analytical expression of the high frequency mode, that holds for any lattice with one atom in the unit cell, is derived. Comparison with numerical simulations, perturbation theory and the polarization potential theory is carried out. We also show that magnetic instabilities tend to vanish for $T\gtrsim t/6$, and finite temperature phase diagrams are established.

Equations

$$H=\sum _{i,j}{t}_{ij}\sum _{\sigma ,{\sigma }^{\prime },{\sigma }^{\prime \prime }}{z}_{i{\sigma }^{\prime \prime }\sigma }^{†}{f}_{i\sigma }^{†}{f}_{j{\sigma }^{\prime }}{z}_{j{\sigma }^{\prime }{\sigma }^{\prime \prime }}+U\sum _i{d}_i^{†}{d}_i.$$

$${\underline{z}}_i=e_i^{†}{\underline{L}}_i{M}_i{\underline{R}}_i{\underline{p}}_i+{\underline{\tilde{p}}}_i^{†}{\underline{R}}_i{M}_i{\underline{L}}_i{d}_i$$

$$M_i$$

$${\text{cali}A}_i$$

$$\displaystyle{\cal L}^{\rm B}=\sum_{i}\bigg{[}e_{i}^{\dagger}\partial_{\tau}e_{i}$$

$${\mathcal{L}}^{\mathrm{F}}=$$

$$\displaystyle\psi(k)=\big{(}\delta e(k),\delta d^{\prime}(k),\delta d^{\prime\prime}(k),\delta p_{0}(k),\delta\beta_{0}(k),\delta\alpha(k),$$

$$\displaystyle\delta p_{1}(k),\delta\beta_{1}(k),\delta p_{2}(k),\delta\beta_{2}(k),\delta p_{3}(k),\delta\beta_{3}(k)\big{)}$$

$${\psi }_{\mathrm{MF}}=(e,d,0,p_0,{\beta }_0,\alpha ,0,0,0,0,0,0)$$

$$\int d\tau \mathcal{L}(\tau )={\mathcal{S}}_{\mathrm{MF}}+\sum _{k,\mu ,\nu }{\psi }_{\mu }(-k)S_{\mu \nu }(k){\psi }_{\nu }(k)$$

$${\chi }_s(k)=⟨\delta {\text{cali}S}_z(-k)\delta {\text{cali}S}_z(k)⟩=2p_0^2{S}_{11,11}^{-1}(k).$$

$${\chi }_c(k)$$

$$z_0=p_0(e+d)\sqrt{\frac{2}{1-{\delta }^2}}$$

$$E_{\mathbf{k}}=z_0^2{t}_{\mathbf{k}}-(\mu -{\beta }_0)$$

$$e=\frac{x^2+\delta }{2x},d=\frac{x^2-\delta }{2x},p_0^2=1-\frac{x^4+{\delta }^2}{2x^2},$$

$$\alpha =\frac{p_0^2{x}^2{U}_0}{2}\left(\frac{1}{x^2+\delta }+\frac{1}{1-\delta }\right),{\beta }_0=\alpha -\frac{x^2{U}_0}{4}\left(1+\frac{2p_0^2}{1-{\delta }^2}\right)$$

$$U_0=-8{\epsilon }_0/(1-{\delta }^2)$$

$${\epsilon }_0=\frac{2}{L}\sum _{\mathbf{k}}{t}_{\mathbf{k}}{n}_F(E_{\mathbf{k}})$$

$$\frac{2}{L}\sum _{\mathbf{k}}{n}_F(E_{\mathbf{k}})=1-\delta$$

$$\frac{(1-x^2)x^4}{x^4-{\delta }^2}=\frac{U}{U_0}.$$

$$F$$

$$\chi_{c}(k)=\frac{e^{2}S_{55}(k)\big{[}\tilde{S}_{33}\big{(}2p_{0}^{2}\Gamma_{1}(k)-8dp_{0}\Gamma_{2}(k)+8d^{2}\Gamma_{3}(k)\big{)}+2e^{2}p_{0}^{2}S_{55}(k)(\omega+{\rm i}0^{+})^{2}\big{]}}{\tilde{S}_{33}\big{[}\Gamma_{2}^{2}(k)-\Gamma_{1}(k)\Gamma_{3}(k)\big{]}-\frac{e^{2}}{(e+d)^{2}}S_{55}(k)\big{[}p_{0}^{2}\Gamma_{1}(k)+2(e-d)p_{0}\Gamma_{2}(k)+(e-d)^{2}\Gamma_{3}(k)\big{]}(\omega+{\rm i}0^{+})^{2}}$$

$${\tilde{S}}_{33}$$

$${\Gamma }_1(k)$$

$${\Gamma }_2(k)$$

$${\Gamma }_3(k)$$

$${\chi }_c(k)=\frac{{\chi }_0(k)}{1+f^s(k){\chi }_0(k)}$$

$$\displaystyle f^{s}(k)=\frac{U}{2}\bigg{[}1$$

$${\chi }_m(k)=\frac{2}{L}\sum _{\mathbf{q}}(t_{\mathbf{q}}+t_{\mathbf{q}+\mathbf{k}}{)}^m\frac{n_F(E_{\mathbf{q}+\mathbf{k}})-n_F(E_{\mathbf{q}})}{(\omega +\mathrm{i}{0}^{+})-(E_{\mathbf{q}+\mathbf{k}}-E_{\mathbf{q}})}.$$

$$\chi_{{\rm pp}}(k)=\frac{\chi^{\rm sc}(k)}{1+\big{[}f^{s}_{\rm pp}({\bf k})+(\omega^{2}/{\bf k}^{2})g^{s}_{\rm pp}({\bf k})\big{]}\chi^{\rm sc}(k)}$$

$$c_s({\theta }_{\mathbf{k}})=\frac{{\omega }_{\mathrm{ZS}}(\mathbf{k})}{|\mathbf{k}|}.$$

$${\epsilon }_{\mathbf{k}}=\frac{2}{L}\sum _{\mathbf{q}}{t}_{\mathbf{q}+\mathbf{k}}{n}_F(E_{\mathbf{q}}).$$

$${\omega }_{\mathrm{UHB}}(\mathbf{k})\approx \frac{U_0}{2}\sqrt{1+\frac{U}{2U_0}\left(1+7{\delta }^2-(1-{\delta }^2)\frac{{\epsilon }_{\mathbf{k}}}{{\epsilon }_0}\right)}.$$

Full text

Hubbard model, slave boson, collective mode

I Introduction

Our understanding of excitations in correlated electron systems has been strongly influenced by the seminal works of Hubbard Hub63 , Landau Lan56 , and Pines PinBoh . In his study Hub63 of the model Hamiltonian that is now associated with his name, Hubbard put forward one very important feature of strongly correlated electrons: the splitting of the bands into the so-called upper Hubbard band (UHB) and the lower Hubbard band (LHB). Its origin can be traced back to the atomic limit where a gap of the order of the interaction strength separates two sets of states, one at $\omega\approx 0$ and the other one at $\omega\approx U$. One then expects that by ramping up the hopping between the sites, the hybridization of the atomic orbitals progressively results in the delocalized states forming the dispersive LHB and UHB. However Hubbard’s treatment fails to produce the predicted Fermi liquid for weak coupling. Indeed, in this regime the self-energy he postulated neither reduces to the perturbation theory result, nor does it yield the correct Fermi surface for the metallic phase Edw68 . In the latter limit Landau’s theory of the Fermi liquid Lan56 ; PinNoz ; Vol84 has proved to be a successful paradigm for understanding a large variety of fermion systems at low temperature such as normal liquid $\,{}^{3}$He, metals or semimetals, and nuclear matter. This phenomenological approach is based on the physical intuition that the low-energy properties of interacting particles can be modeled from a gas of elementary excitations, referred to as quasiparticles, which are formed with a life time that is infinite on the Fermi surface but rapidly decays away from it. Alternatively to Landau’s original formulation this result can be obtained within perturbation theory. Using the latter to compute correlation functions within the random phase approximation (RPA), Pines and Bohm PinBoh showed that the response functions are composed of a continuum generated by the incoherent response of quasiparticles supplemented by peaks signaling collective excitations, that arise as dynamical fluctuations of the groundstate. The dispersion and attenuation of the collective modes are then indicative of the nature of the phase. Prominent examples are the Goldstone modes which appear when a continuous symmetry is spontaneously broken at a phase transition, such as phonons for rotational and translational symmetries, or magnons for the spin-rotational symmetry.

Reconciling the Fermi liquid with Hubbard’s local physics remains an important and largely unsolved problem for correlated electrons. Indeed numerical approaches generically face finite-size effects Bor91 because required computing resources blow up exponentially with increasing system size. Yet results can be obtained in limiting cases such as the infinite-coordination lattice where the dynamical mean-field theory Geo96 ; Bul01 catches the Hubbard gap in the vicinity of half filling. However in the doped system a clear picture is still missing, possibly because of the formation of incommensurate phases with large unit cells which cannot be captured by the method and its cluster extensions Shr90 ; Fre91 ; SeiSi ; Sei02 ; Rac06 ; Cha10 ; stripes ; Cor14 ; Lep15 . Furthermore, the above mentioned approaches are mainly focused on a self-consistent calculation of one-particle correlation functions, which can be directly related to experimental observations such as photoemission. However other experimental techniques require the knowledge of two-particle quantities such as the charge and spin response functions that are probed in neutron scattering experiments. Computing two-particle correlations is more challenging because of the need to include vertex corrections Geo96 ; Haf14 .

The purpose of this work is to compute the charge response function of the Hubbard model using an extension of the Kotliar and Ruckenstein slave-boson representation. One of our main results is that it reduces to the RPA susceptibility for weak coupling. The obtained charge excitation spectrum generically consists of i) a continuum the width of which decreases with increasing interaction strength and density, ii) a collective mode with anisotropic zero-sound (ZS) velocity, and iii) a high frequency mode at $\omega\approx U$ which is the signature of the UHB. Hence our scheme reconciles the Fermi liquid physics — including collective modes — with Hubbard’s local physics embedded in the split bands. The calculation is carried out in a paramagnetic phase, free of symmetry breaking, in the thermodynamical limit. It allows us to resolve the full momentum dependence of the spectra. At first glance, neglecting magnetic instabilities puts severe constraints on the parameter range where the calculation may be meaningfully performed. However, as shown below, the incommensurate magnetic instabilities are strongly suppressed with increasing temperature, so that they essentially disappear for $T\approx t/6$.

Since Mott insulating groundstates arise at large $U$ and at half filling, we perform our investigations in a framework which is able to capture interaction effects beyond the physics of Slater determinants. We use an extension of the Kotliar and Ruckenstein slave-boson representation that reproduces the Gutzwiller approximation on the saddle-point level Kot86 . It entails the interaction driven Brinkman-Rice metal-to-insulator transition Bri70 . A whole range of valuable results have been obtained with Kotliar and Ruckenstein Kot86 and related slave-boson representations Li89 ; FW , which motivates the present study. In particular they have been used to describe anti-ferromagnetic Lil90 , spiral Fre91 ; Igo13 ; Fre92 ; Doll2 , and striped SeiSi ; Fle01 ; Sei02 ; Rac06 phases. Furthermore, the competition between the latter two has been addressed as well RaEPL . Besides, it has been obtained that the spiral order continuously evolves to the ferromagnetic order in the large $U$ regime ($U\gtrsim 60t$) Doll2 so that it is unlikely to be realized experimentally. Consistently, in the two-band model, ferromagnetism was found as a possible groundstate only in the doped Mott insulating regime Fre02 . Yet, adding a ferromagnetic exchange coupling was shown to bring the ferromagnetic instability line into the intermediate coupling regime lhoutellier15 . A similar effect has been obtained with a sufficiently large next-nearest-neighbor hopping amplitude FW98 or going to the fcc lattice Igo15 .

The influence of the lattice geometry on the metal-to-insulator transition was discussed, too Kot00 . For instance, a very good agreement with Quantum Monte Carlo simulations on the location of the metal-to-insulator transition for the honeycomb lattice has been demonstrated Doll3 . Also, strongly inhomogeneous polaronic states that have been found in correlated heterostructures have also been addressed using this formalism applied to the Hubbard model extended with inter-site Coulomb interactions Pav06 . Most recently the approach has been used to address possible capacitance enhancement in a capacitor consisting of strongly correlated plates separated by a dielectric Ste17 . Furthermore, comparison of groundstate energies to existing numerical solutions have been carried out for the square lattice, too. For instance, for $U=4t$ it could be shown that the slave-boson groundstate energy is larger than its counterpart by less than 3% Fre91 . For larger values of $U$, it has been obtained that the slave-boson groundstate energy exceeds the exact diagonalization data by less than 4% (7%) for $U=8t$ ($20t$) and doping larger than 15%. The discrepancy increases when the doping is lowered Fre92 . It should also be emphasized that quantitative agreement to quantum Monte Carlo charge structure factors was established Zim97 .

The paper is organized as follows. In Sec. II we give a brief presentation of the spin-rotation-invariant (SRI) Kotliar and Ruckenstein slave-boson representation of the Hubbard model and the method used to calculate dynamical response functions (more details can be found in, e.g., review fresard12 ). Sec. III presents the paramagnetic saddle-point solution and discusses its temperature dependence. In addition, phase diagrams summarizing the temperature dependence of magnetic and charge instabilities are established. We evaluate the spin and charge susceptibilities from fluctuations captured within the one-loop approximation in Sec. IV and investigate the dispersion of their collective modes in Sec. V. Our results are discussed in comparison with the perturbation Hartree-Fock (HF)+ RPA prediction, as well as with available numerical investigations (exact diagonalization and Quantum Monte Carlo method) in Sec. VI. And we summarize the paper in the conclusion.

II Model and method

The Hubbard Hamiltonian in the SRI Kotliar and Ruckenstein slave-boson representation Kot86 ; fresard12 is expressed with auxiliary boson operators $e_{i}$, $p_{i\mu}$, $d_{i}$ (for atomic states with respectively zero, single and double occupancy) and pseudo-fermion operators $f_{i\sigma}$ as

[TABLE]

Here hopping occurs between nearest neighbor sites with amplitude $t_{ij}=-t$. A key feature of the representation lies in the reduction of the on-site Coulomb interaction into a term bilinear in bosonic operators, at the expense of a more complicated hopping term. In order to preserve spin rotation symmetry Li89 ; FW the canonical operator $p_{i\mu}$ build a $2\times 2$ matrix in spin space that is expanded into the identity matrix $\underline{\tau}^{0}$ and the Pauli matrices as ${\underline{p}}_{i}=\frac{1}{2}\sum_{\mu=0}^{3}p_{i\mu}{\underline{\tau}}^{\mu}$. In this space the occupancy-change operator $\underline{z}_{i}$ in the hopping term is also a matrix defined as

[TABLE]

with

[TABLE]

where $\tilde{\underline{p}}_{i}=\frac{1}{2}(p_{i0}{\underline{\tau}}^{0}-{\bf p}_{i}\cdot\boldsymbol{\underline{\tau}})$.

In the augmented Fock space generated by the auxiliary boson operators, the subspace of physical states is the intersection of the kernels of operators

[TABLE]

i.e. in this subspace ${\cali A}_{i}=0$ that is the constraint of one atomic state per site, and ${\cali B}_{i\mu}=0$ which equates the number of fermions to the number of $p$ and $d$ bosons.

The partition function is calculated as a functional integral li91 ; Zim97 with the effective Lagrangian ${\cal L}={\cal L}^{\rm B}+{\cal L}^{\rm F}$ where the purely bosonic part is

[TABLE]

with ${\cali B}_{i\mu}^{\rm B}$ being the bosonic part of the operator ${\cali B}_{i\mu}$, and the mixed-fermion-boson part can be written as

[TABLE]

after the fermion fields have been integrated (here $\mu$ is the chemical potential). The constraints that define the physical states are enforced with Lagrange multipliers $\alpha_{i}$ and $\beta_{i\mu}$. The internal gauge symmetry group of the representation allows to simplify the problem. The phases of $e$ and $p_{\mu}$ can be gauged away by promoting the Lagrange multipliers to time-dependent fields FW , leaving us with radial slave-boson fields Fre01 . Their values obtained at the saddle-point level may be viewed as an approximation to their exact expectation values that are generically non-vanishing Kop07 . The slave-boson field corresponding to double occupancy $d_{i}=d^{\prime}_{i}+{\rm i}d^{\prime\prime}_{i}$ however has to remain complex as emphasized by several authors Jol91 ; Kot92 ; FW . Since $e_{i}$ and $p_{i\mu}$ are now real, their kinetic terms drop out of ${\cal L}^{\rm B}$ due to the periodic boundary conditions on boson fields.

Within the approximation of Gaussian fluctuations, the action is expanded to second order in field fluctuations

[TABLE]

around the paramagnetic saddle-point solution

[TABLE]

as

[TABLE]

(the matrix $S$ is given in Appendix A). We have introduced the notation $k=({\bf k},\nu_{n})$, where $\nu_{n}=2\pi nT$, and $\sum_{k}=T\sum_{\nu_{n}}L^{-1}\sum_{{\bf k}}$ with $L$ the number of lattice sites. The correlation functions of boson fields are then Gaussian integrals which can be obtained from the inverse of the fluctuation matrix $S$ as $\langle\psi_{\mu}(-k)\psi_{\nu}(k)\rangle=\frac{1}{2}S^{-1}_{\mu\nu}(k)$. For instance the slave-boson representation of the spin fluctuation $\delta{\cali S}_{z}=\delta(p_{0}^{\dagger}p_{3}+p_{3}^{\dagger}p_{0})$ yields the spin susceptibility

[TABLE]

Similarly, using the density fluctuation $\delta{\cali N}=\delta(d^{\dagger}d-e^{\dagger}e)$, the charge susceptibility is

[TABLE]

Dynamical response functions are eventually evaluated within analytical continuation ${\rm i}\nu_{n}\rightarrow\omega+{\rm i}0^{+}$.

The saddle-point approximation is exact in the large degeneracy limit, and the Gaussian fluctuations provide the $1/N$ corrections FW . Moreover it obeys a variational principle in the limit of large spatial dimensions where the Gutzwiller approximation (GA) becomes exact for the Gutzwiller wave function Met89 .

III Paramagnetic saddle-point solution

III.1 Characterization of the paramagnetic phase

At the paramagnetic saddle-point, the field $\underline{z}_{i}$ reduces to $z_{0}\underline{\tau}^{0}$ with

[TABLE]

where $\delta=1-\langle{\cali N}\rangle$ is the hole doping from half-filling. The factor $z_{0}^{2}$ plays the role of a quasiparticle residue, and it also renormalizes the quasiparticle dispersion as

[TABLE]

with the bare dispersion $t_{{\bf k}}=-2t(\cos k_{x}+\cos k_{y})$ for the square lattice.

The boson saddle-point values can be expressed with the doping and the variable $x=e+d$ as

[TABLE]

(the expressions result from the constraints on physical states $e^{2}+p_{0}^{2}+d^{2}=1$ and $p_{0}^{2}+2d^{2}=1-\delta$, and saddle-point conditions). Here the coupling scale

[TABLE]

has been introduced in terms of the semi-renormalized kinetic energy

[TABLE]

and the Fermi function $n_{F}(\epsilon)=1/(\exp(\epsilon/T)+1)$.

As discussed in Ref. Vol87 ; FW , the paramagnetic solution for fixed values of doping $\delta$ and coupling $U$ is found by determining the chemical potential via the filling condition

[TABLE]

and the solution of the saddle-point equation

[TABLE]

The procedure is carried out self-consistently with the evaluation of $z_{0}$ since the latter renormalizes the dispersion. It is however simplified at $T=0$ because then, for a fixed filling, $\varepsilon_{0}$ and $U_{0}$ have the same values for all finite $z_{0}$. This implies they do not vary with the coupling, except at $\delta=0$ where they vanish discontinuously above a critical coupling $U_{c}$. Fig. 1 displays their variations with the doping. Increasing the temperature from zero reduces their amplitudes and it smooths out the discontinuity at half filling while enlarging the collapse around it, as shown by the curves plotted at temperature $T=t/10$.

For most values of coupling and doping, the saddle-point equation possesses one finite solution $x>0$ corresponding to a metallic state. As shown in Fig. 2 saddle-point values converge in the infinite coupling limit where $x=\sqrt{|\delta|}$ and $z_{0}^{2}=2|\delta|/(1+|\delta|)$. A remarkable phenomenon occurs at half filling where $x$ vanishes above the critical coupling that is $U_{c}=-8\varepsilon_{0}=2(8/\pi)^{2}t\approx 12.97t$ at $T=0$. This solution corresponds to an insulating state since $z_{0}^{2}=0$ results in a diverging quasiparticle mass and a vanishing quasiparticle residue: This is the Brinkman-Rice mechanism Bri70 for the Mott metal-to-insulator transition. Note that at finite temperature, for small doping and $U<U_{c}$, the equation can have up to three positive solutions Doll3 , among which the groundstate is determined by minimizing the free energy

[TABLE]

The degeneracy of solutions gives rise to a first-order transition when increasing the coupling from a metallic state into either an insulating state at half filling FK97 , or a bad metallic one characterized by a small quasiparticle residue $z_{0}^{2}$ for finite doping (see Fig. 2). Fig. 3 shows the transition line in the $(\delta,U)$-phase diagram at different temperatures. It is terminated by a critical endpoint at finite doping, and as well as the parameter region with multiple solutions, it shrinks with lowering temperature and vanishes at $T=0$ Doll3 ; camjayi07 . Contrary to an ordinary band insulator where thermal excitations of quasiparticle enhance the conductivity, increasing the temperature in the strongly correlated Hubbard model can induce a transition from a low-temperature metal to a high-temperature insulator as thermal fluctuations destroy the poor coherence of the small-$z_{0}$ metallic state in a fashion similar to the transition observed in V2O3 McWhan73 .

III.2 Instabilities

Let us now look for the parameter range in which the above solution is stable. In Ref. doll93 it was found that the zero-temperature slave-boson paramagnetic phase is stable at low density, even at large coupling, while incommensurate magnetic instabilities develop at large densities. One then may ask what is the picture at finite temperature, especially since early estimates at half-filling yield a temperature at which magnetic instabilities are destroyed to be of order $t^{2}/U$ Doll3 . Furthermore, while there is a regain of interest in charge instabilities at $T=0$ lhoutellier15 ; Ste17 , little attention has been paid to them at finite temperature. Hence the robustness of the saddle-point solution against spin and charge fluctuations is investigated by looking for a divergence of the respective static response functions (see Eq. (37) and Eq. (20)). The instabilities of the paramagnetic phase at different temperatures are mapped in Fig. 4. The static spin susceptibility $\chi_{s}({\bf k},\omega=0)$, given by Eq. (37), has no pole at high temperature but a magnetic instability appears below $T\approx t$, around half filling and for a finite but not large coupling. Its domain in the $(\delta,U)$-phase diagram then grows with lowering temperature, with a significant variation between $T=t/6$ and $t/8$. Earlier studies Fre92 ; doll93 have found that the instability boundary in the phase diagram at $T=0$ signals a magnetic ordering into a spiral groundstate. The doping range of the magnetic phase increases with the coupling up to the maximum doping $\delta\approx 0.63$ reached at $U\gtrsim 60t$. Contrary to the magnetic behavior, the domain of the charge instability shrinks with lowering temperature. It is limited to small doping and occurs at all coupling above a moderate threshold value which increases with lowering temperature. The charge instability is related to a tendency towards a phase separation Fre92 or towards the more complicated stripe phases SeiSi ; Sei02 ; Rac06 .

IV Expressions of the dynamical response functions

The evaluation of correlation functions is simplified in the paramagnetic state because the Gaussian fluctuations decouple into spin and charge channels. This results into a matrix $S$ that is block diagonal with a charge $6\times 6$ submatrix and three identical $2\times 2$ blocks for the components of the spin. As discussed in, e.g., Refs. Lav90 ; Zim97 ; li91 ; li94 , the blocks can be independently inverted to yield the spin (see Appendix B) and the charge dynamical response function

[TABLE]

with

[TABLE]

The susceptibilities are particle-hole symmetric as expected for the Hubbard model on the square lattice.

The expression of $\chi_{c}(k)$ given in Ref. Zim97 is valid only at zero frequency because the matrix elements omitted in the previous work vanish in the static limit. We have checked that the numerical discrepancies between the charge structure factors evaluated in Ref. Zim97 and using Eq. (20) are minor. They do not alter the previous conclusion that slave-boson results are in very good agreement with Quantum Monte Carlo calculations Zim97 ; DziUnpub . However the missing matrix elements are crucial in the investigation of charge collective modes. Without them, the poles of the dynamical response function (or their residues) would not vanish in the free-particle limit. Furthermore their dispersions would depend on the sign of the doping, which is in conflict with the particle-hole symmetry expected for the Hubbard model on a bipartite lattice.

In the weak-coupling limit slave-boson expressions yield the textbook results derived from perturbation methods. This is obtained by writing the charge susceptibility (20) as

[TABLE]

then by expanding the function $f^{s}(k)$ in powers of the coupling. Here the Lindhard function $\chi_{0}(k)$, given by Eq. (24), solely differs from the charge response function of a Fermi gas through the quasiparticle mass renormalization $z_{0}^{2}$. The function $f^{s}(k)$ is related to the Landau parameter of Fermi-liquid theory lhoutellier15 by $F_{0}^{s}=\chi_{0}(0)f^{s}(0)=N_{\rm F}f^{s}(0)$ where $N_{\rm F}$ is the density of states at the Fermi level. Its expansion to first order $f^{s}(k)=U/2$ is in perfect correspondence with the expected RPA result. This generalizes Li et al.’s results li91 ; li94 to arbitrary momentum and frequency. Including the next order in $U$ yields

[TABLE]

with the ratio $\gamma(k)=\chi_{1}(k)/\varepsilon_{0}\chi_{0}(k)$ and

[TABLE]

The ratio $\gamma(k)$ in the second order expansion has a complex value. Hence the function $f^{s}(k)$ actually possesses an imaginary part. Its real part becomes negative just below a critical energy at which it diverges (see Fig. 5). The domain with ${\rm Re}f^{s}<0$ is largest around M. Its size increases with the doping and the coupling. As can be inferred from the structure of $f^{s}(k)$, we show in the next section that the charge susceptibility (20) has a rich spectrum that cannot be captured within the conventional HF + RPA framework.

A theory going beyond the Landau Fermi liquid model and the RPA approximation has been developed by Pines and coworkers Pin81 for the excitations and transport properties of quantum liquids. The so-called polarization potential (PP) theory is a semi-phenomenological approach that describes the collective action of the particles by an averaged self-consistent field which can be polarized by particle-hole excitations via an effective screened potential. Using parameters obtained from static measurements and sum rule considerations, it attempts to describe both liquid 4He and 3He within a unified formalism. In particular, the theory can reproduce the experimental dispersion of the ZS collective mode, beyond the Landau Fermi liquid regime. They obtained a density response of the form

[TABLE]

within the linear response theory. Two contributions enter the PP. The first function $f^{s}_{\rm pp}({\bf k})$ is the Fourier transform of the potential of an effective static particle interaction. The second term corresponds to the effect of the so-called backflow, that is the additional screening caused by longitudinal current fluctuations accompanying the density fluctuations. In the long wavelength limit these quantities are related to the Landau parameters by $f^{s}_{\rm pp}(0)=F_{0}^{s}/N_{\rm F}$ and $g^{s}_{\rm pp}(0)=mF_{1}^{s}/3\langle{\cali N}\rangle$ where $m$ is the particle mass. A reasonable description of the neutron-scattering data on 4He and 3He can be obtained by assuming the PP to be essentially the same for both liquids. The influence of the statistics is mainly present in the screened density response function $\chi^{\rm sc}(k)$. Using a sum rule argument, the latter is defined as the weighted sum of the expression for a free Bose or Fermi gas of particles with an effective mass $m^{*}$, and a structureless multiparticle contribution that is fitted to the experimental data.

Comparing expression (22) of the density response function with equation (25), one can note two distinguishing features. First the PP used in $\chi_{\rm pp}(k)$ appears to be an expansion of the function $f^{s}(k)$ to second order in the frequency. With such a frequency dependence, the PP is not singular and $\chi_{\rm pp}(k)$ possesses one single pole corresponding to the ZS mode. It cannot then produce a second collective excitation, contrary to our result. As shown by Eq. (23), $f^{s}(k)$ can diverge and it can then give rise to another collective mode. As shown below, for strong coupling, the latter disperses around $\omega\approx U$ and we therefore call it the UHB mode. However the PP theory includes a phenomenological multiparticle contribution in the screened density response function $\chi_{\rm sc}(k)$, which is absent from the approximation level used in the present work. Multiparticle processes may have a significant influence on the collective modes as, for instance, within the PP theory they soften the ZS mode at large wavevectors. Including them in our approach could, in principle, be achieved with an expansion of the action going beyond the Gaussian fluctuation approximation.

V Charge collective modes

The charge susceptibility possesses two collective modes that appear at finite coupling. These excitations form narrow peaks at well defined energies in the spectrum of the inelastic response ${\rm Im}\chi_{c}(k)$. As shown in Fig. 6 the spectrum is composed of a broad continuum that results from incoherent single-particle excitations. Beyond its upper boundary $\omega_{\rm cont}({\bf k})$ lie the peaks of the two modes. The typical evolution of the charge response function with the coupling is plotted in Fig. 7 and the effect of doping is shown in Fig. 8. The continuum contribution to $\chi_{c}(k)$ is roughly reduced by a factor $\sim(1+UN_{\rm F}/2)$ while its energy width shrinks as $z_{0}^{2}$. The mode at lower energy $\omega_{\rm ZS}({\bf k})$ is the zero-sound mode. It has a linear dispersion at long wavelength, that is around the k-point $\Gamma$. It appears as a resonance at the upper edge $\omega_{\rm cont}({\bf k})$ and it changes into a well defined peak that departs from the continuum when increasing the coupling. The second mode is the upper-Hubbard-band mode which occurs at higher energy $\omega_{\rm UHB}({\bf k})$. It appears at small coupling with no dispersion at $\omega=U_{0}/2$ and it then develops with a gap at ${\bf k}=\Gamma$ that grows as $U$ in the strong-coupling limit.

The dispersions of the collective modes are presented below in more details. Since our results are best understood at $T=0$ we postpone the discussion of temperature effects to the end of the section.

V.1 Zero-sound mode

The conditions under which the collective modes develop can be discussed with the weak-coupling expressions (22) and (23) for the susceptibility. To first order in the coupling, $f^{s}(k)\approx U/2$ so the denominator of $\chi_{c}(k)$ can vanish only if $\chi_{0}(k)$ is real and negative. As shown in Fig. 6, these conditions are met beyond the upper edge $\omega_{\rm cont}({\bf k})$ of the response continuum, which corresponds to the largest energy of the particle-hole excitations with momentum ${\bf k}$. Generally ${\rm Re}\chi_{0}(k)$ has a deep minimum at $\omega_{\rm cont}({\bf k})$ that can even diverge if ${\rm Im}\chi_{0}(k)$ varies discontinuously. Hence, the charge susceptibility can develop a pole in the vicinity of the upper edge, which results in the onset of the ZS mode even for a small coupling. In this regard the $(1,0)$ and $(0,1)$-direction are special. The particle-hole susceptibility $\chi_{0}$ takes the form of a 1D-response for $\bf k$ along $\Gamma$-X. On a large range of doping around half filling, this results in a square-root singularity at $\omega_{\rm cont}({\bf k})$ which ensures the existence of the ZS mode along the symmetry axis and around the ${\bf k}$-point X. Note however that the mode is suppressed just below the UHB mode energy because ${\rm Re}f^{s}(k)$ becomes negative (see Fig. 5).

Close to half filling the ZS mode exists for nearly all momenta. As shown in Figs. 7 and 8, at strong coupling, the intensity of the charge response is largely transfered from the single-particle processes to the collective modes. Increasing the doping results in the softening of the ZS mode, while the response continuum grows as the quasiparticle mass is less renormalized. Eventually, at large doping, the ZS pole is suppressed for nearly all wavevectors as the singularity of ${\rm Re}\chi_{0}(k)$ at the continuum boundary is smoothed out.

At long wavelength, that is in the vicinity of $\Gamma$, the dispersion of the pole is proportional to $|{\bf k}|$ and one can define the ZS velocity as

[TABLE]

For the Hubbard model on the square lattice the sound velocity is anisotropic. The maximum is along the M-direction while the minimum is along the X-direction (see Fig. 9). However the anisotropy vanishes in two limiting cases: at large doping $|\delta|\approx 1$ as the quasiparticle dispersion around the Fermi energy tends to a parabolic dispersion and, more surprisingly, close to half filling for strong coupling. In the latter case, the isotropy is approached because the ZS pole is located far above the strongly renormalized edge $\omega_{\rm cont}({\bf k})$, at an energy where the functions $\chi_{m}(k)$ at long wavelength are dominated by their $s$-wave component.

The sound velocity along the two high-symmetry directions is plotted in Fig. 10 for different values of coupling and doping. At small coupling the collective mode appears close to the continuum upper boundary which is $\omega_{\rm cont}({\bf k})=z_{0}^{2}{\rm max}_{{\bf v}_{\rm F}^{0}}({\bf v}_{\rm F}^{0}\cdot{\bf k})$ for small momentum. Here ${\bf v}_{\rm F}^{0}=\left.\frac{\partial t_{\bf q}}{\partial{\bf q}}\right|_{{\bf q}_{\rm F}}$ is the bare Fermi velocity. Hence the weak-coupling approximation yields the velocity in the M-direction c_{s}(\frac{\pi}{4})\approx v_{\rm F}^{0}(\frac{\pi}{4})\big{[}1+\big{(}\frac{1}{(18t)^{2}}-\frac{1-\delta^{2}}{U_{0}^{2}}\big{)}U^{2}\big{]} where v_{\rm F}^{0}(\frac{\pi}{4})=2t\sqrt{2\big{(}1-({\mu}/{4t})^{2}\big{)}}. In the X-direction, for doping $|\delta|\gtrsim 0.63$, the expression remains the same though with v_{\rm F}^{0}(0)=2t\sqrt{\big{(}1-(1-{|\mu|}/{2t})^{2}\big{)}}. Otherwise for smaller doping we find c_{s}(0)\approx 2t\big{[}1+\big{(}\frac{1}{(9t)^{2}}-\frac{1-\delta^{2}}{U_{0}^{2}}\big{)}U^{2}\big{]}. The evolution of the velocity with the coupling is complicated since it is governed by two opposite trends. On the one hand the increase of the quasiparticle mass reduces it. On the other hand the increase of $f^{s}$ with $U$ moves the ZS pole to higher energy. As a result, at large doping $|\delta|\approx 1$ where the mass renormalization can be neglected, the velocity increases with increasing coupling. Then at a smaller doping the renormalization is more important and the velocity variation depends on the propagation angle: $c_{s}(\frac{\pi}{4})$ decreases while $c_{s}(0)$ increases before eventually decreasing at strong coupling. Lastly, in the vicinity of half filling the variation of $c_{s}$ is non-monotonic (see Fig. 11). The velocity reaches a maximum at a coupling below $U_{c}$ before collapsing to c_{s}\sim 2|\varepsilon_{0}|\sqrt{|\delta|\big{(}1+\frac{U_{0}}{U}\big{)}} in the bad-metal state. The behavior at $\delta=0$ is even discontinuous: $c_{s}$ abruptly falls at $U_{c}$ from its maximum value $\approx 3.2t$ to zero. As previously noted the velocity around half filling becomes isotropic at large coupling.

V.2 Upper-Hubbard-band mode

A charge excitation with an energy of the order of $U$ has been predicted as the result of strong correlation effects since the early days of the Hubbard model. Indeed in the vanishing hopping limit $t=0$, all particles rest localized at the energy of the atomic levels $\omega=0$ or $\omega=U$. A perturbative inclusion of the hopping, as done by Hubbard and extended by Pairault et al. Pai00 , results in the broadening of the atomic levels and the formation of dispersive bands around each one, the lower and upper Hubbard band. Hence excitations resting on the UHB are expected from this physical picture.

The slave-boson approach yields such a collective excitation, below denoted the UHB mode, which occurs at an energy $\omega_{\rm UHB}({\bf k})$ that grows as $U$ for strong coupling. Like the ZS mode, the UHB mode has an energy dispersion with a minimum at $\Gamma$ and a maximum at M, but pushed to a higher energy (see Figs. 7 and 8). Actually the excitation energy does not vanish at $\Gamma$, even at small coupling. Numerical evaluations find that the peak weight is zero at $\Gamma$ and maximum at M. These features are illustrated in Figs. 12 and 13 where the dispersion with momenta along $\Gamma-{\rm M}$ is plotted for different dopings and couplings. The mode appears at weak coupling around $\omega=U_{0}/2$ which is the frequency where $f^{s}(k)$ diverges (see the second-order expansion (23)). No dispersion is observed at the onset of the mode. Although its pole exists at any finite coupling, the mode disappears in the uncorrelated limit as its residue vanishes at $U=0$. A shift to higher energy can be observed with increasing doping or coupling. Their influences on the dispersion width are opposite. A widening is obtained by increasing the coupling while the effect of doping is to narrow the dispersion to the point that it vanishes at $|\delta|=1$. On the whole the mode has its maximum weight at M and it is most clearly observed for a moderately large coupling $U\sim 5t$ at small but finite doping $|\delta|\sim 0.1$. Indeed its weight decreases in the close vicinity of half filling, and vanishes at $\delta=0$. We found that it is also vanishingly small at $|\delta|=1$.

The UHB mode can be distinguished from the ZS mode and the response continuum because its energy is generally larger than $\omega_{\rm ZS}({\bf k})$. However this is not necessarily the case at weak coupling. Spectra of ${\rm Im}\chi_{c}(k)$ in Fig. 12 show that it enters the quasiparticle continuum for momenta around M. This results in the damping of the excitation by quasiparticle scattering and the mode peak is replaced by a depletion around $\omega\approx U_{0}/2$ in the charge response continuum. At moderate coupling ($U\sim 3t$) the ZS mode that appears just beyond the continuum edge hybridizes with the UHB mode around M, and there is only one single peak around M that continuously becomes the UHB peak as ${\bf k}$ goes to $\Gamma$. The depletion associated to the UHB mode moves to higher energy with increasing coupling. After it exits the continuum the ZS mode can extend until M where it forms a second well-defined peak below the UHB one.

Analytical expressions for the dispersion of the UHB mode can be obtained at weak and strong coupling. The mode mostly occurs far beyond the continuum where $\chi_{m}(k)\sim 1/\omega^{2}$ and in particular $\chi_{0}(k)\approx 2z_{0}^{2}(\varepsilon_{0}-\varepsilon_{\bf k})/\omega^{2}$ with

[TABLE]

To first order in the high-energy expansion the denominator (20) behaves as $\omega^{2}-\omega_{\rm HB}^{2}$. The charge response function then possesses two poles, one at negative energy $\omega_{\rm LHB}=-\omega_{\rm HB}$ and one at positive energy $\omega_{\rm UHB}=\omega_{\rm HB}$.

At small coupling $U\ll U_{0}$ , the saddle-point solution can be approximated with $x^{2}\approx 1-(1-\delta^{2})U/U_{0}$ which yields

[TABLE]

The weak-coupling expression highlights several features of the UHB mode dispersion. Firstly, the collective mode appears around the energy $U_{0}/2$ with a dispersion that is vanishingly small. The expression also shows that doping results into a narrower dispersion that is shifted to a higher energy, as seen in Figs. 12 and 13. The dispersion width is approximately equal to $(1-\delta^{2})U/4$ and it vanishes for $|\delta|=1$.

The approximation in the strong-coupling limit is obtained with x^{2}\approx|\delta|/{\sqrt{1-\frac{U_{0}}{U}\Big{(}1-{|\delta|}/{\sqrt{1-U_{0}/U}}\Big{)}}} which gives

[TABLE]

The dispersion thus has its minimum $\approx U-U_{0}(\frac{1}{2}-|\delta|)$ at $\Gamma$, and its width is approximately $(1-|\delta|)U_{0}/2$. Hence at large coupling the energy of the mode grows linearly as the on-site Coulomb interaction $U$. This genuine strong correlation effect is one of the most important results of this work. Being of order $U$ the mode follows from the UHB, that is not captured by the conventional HF + RPA approach. It should also be emphasized that Eqs. (28) and (29) hold for arbitrary lattices with one atom in the unit cell, irrespective of the dimensionality.

V.3 Effect of temperature

Within our theory the impact of temperature on the collective modes manifests itself in two different ways. Firstly the collective mode dispersion shrinks with increasing temperature. This results from the decrease of the saddle-point values, most notably for doping $|\delta|\lesssim 0.1$ and strong coupling. As shown in Fig. 1 and Fig. 14 the averaged kinetic energy $\varepsilon_{0}$, the coupling scale $U_{0}$, and the inverse-mass renormalization factor $z_{0}^{2}$ vary significantly with temperature for this regime of parameters. However in this region of the phase diagram (Fig. 4) the paramagnetic solution is unstable toward phase separation or incommensurate magnetic ordering. Outside this regime, where our investigation is of better relevance, the effect of temperature is a mild reduction of the amplitudes of the saddle-point values. Thus increasing the temperature up to $T=t/3$ slightly scales down the spectrum along the energy axis. We will not discuss the regime of high temperature where the approximation of Gaussian fluctuations certainly becomes insufficient. For instance we expect that incoherent multi-particle processes, which are not taken into account here, get more prominent and modify significantly the charge response of the system, as exemplified by the physics of liquid helium.

The second notable effect of temperature is the broadening of the collective mode peak, which results from scattering of thermally excited quasiparticles. Let us first remark that in the absence of incoherent multi-particle processes, the peak is not damped above the energy $\Delta E_{\rm max}({\bf k})=4tz_{0}^{2}(|\sin\frac{k_{x}}{2}|+|\sin\frac{k_{y}}{2}|)$ of the most energetic one-particle transition with momentum transfer $\bf k$. The UHB peak generally lies above it so its shape is hardly affected by increasing the temperature. This is not the case of the ZS mode for small wavevectors at large doping, and in the vicinity of $\Gamma$ at any finite doping. The ZS peak is broadened because the charge response continuum does not extend up to $\Delta E_{\rm max}({\bf k})$ for small wavevectors. The reason comes from the Fermi statistics which, at zero temperature, excludes some one-particle transitions, among which can be found the most energetic one that occurs between the states of momenta $(\frac{\pi-k_{x}}{2},\frac{\pi-k_{y}}{2})$ and $(\frac{\pi+k_{x}}{2},\frac{\pi+k_{y}}{2})$. As a result the ZS peak can be located between the continuum upper edge $\omega^{T=0}_{\rm cont}({\bf k})$ at $T=0$ and $\Delta E_{\rm max}({\bf k})$. Increasing the temperature then smears the Fermi distribution, which populates the response continuum in this energy range, and eventually broadens the ZS peak.

VI Comparison with other approaches

VI.1 Comparison with HF + RPA result

In the weak-coupling limit the charge response obtained within the slave-boson method is mostly similar to the standard HF + RPA result. But, as stated earlier, the former possesses a supplementary collective mode at high energy, the UHB mode. And, although the perturbation method also produces a ZS mode, it fails to account for the correlation effects which strongly renormalize the quasiparticle mass around half filling and for the dynamical screening of the electron interaction. This is shown in Fig. 15 where the slave-boson charge response is compared with the HF + RPA response

[TABLE]

Here $\chi_{0}^{(0)}(k)$ is the charge response function of a Fermi gas, i.e. with no mass renormalization. At moderate coupling $U=t$ the only observable difference between the two responses is the dispersionless UHB mode. The contribution of the latter is small and the weight of its peak actually vanishes in the limit $U=0$. However at large coupling $U=8t$, the two spectra are quite different. The slave-boson response has two well separated collective modes while the perturbation method only yields the ZS mode. Furthermore, the continuum width and the ZS dispersion shrink due to the quasiparticle mass enhancement, whereas such a correlation effect is not captured by $\chi_{\rm RPA}(k)$. The mass renormalization is not the only effect of correlations. At large doping, the ZS peak in $\chi_{c}(k)$ disappears around ${\bf k}={\rm M}$, in contrast to the HF + RPA prediction. This is because the bare electron interaction $U/2$ of the perturbation result is replaced by the complex function $f^{s}(k)$ within the slave-boson approach. The latter depends on frequency and momentum, and it can have a negative real part near $\omega_{\rm cont}({\bf k})$ (see Fig. 5) which thus suppresses the ZS pole.

VI.2 Comparison with time-dependent GA

The Kotliar and Ruckenstein slave-boson approach has historically been designed to reproduce the Gutzwiller approximation at the saddle-point level Kot86 , thereby strongly tiding both schemes. Later on, a method to calculate excitations at zero temperature has been built on the GA and the RPA Sei01 . It takes the form of the above RPA series with an effective interaction, therefore missing the physics of the Hubbard split bands in the charge response function. Yet a refined treatment has been proposed in Ref. Bun13 , which we now compare to the slave-boson result.

We restrict the analysis to the double-occupancy excitations for which the comparison is simplified. We note that all three terms in Eq. (11) contribute to the particle-hole continuum, implying a damping in the double-occupancy excitation spectra that is absent from the time-dependent Gutzwiller approximation (TDGA) Bun13 . From a quantitative point of view, one can observe that the pole of the double-occupancy propagator found by the TDGA Eq. (100) in Bun13 is located at an energy smaller that the slave-boson one. The discrepancy is largest at ${\bf k}={\rm M}$, for small doping, and strong coupling. For instance the TDGA (slave-boson) pole disperses from $\omega/t=7.1$ to 8.8 (7.2 to 10) for $U=8t$ and $\delta=0.2$, and from $\omega/t=14.4$ to 17.8 (14.6 to 20.1) for $U=20t$ and $\delta=0.1$. Hence the excitations computed within TDGA show both qualitative and quantitative differences to our results which are controlled by the $1/N$ expansion FW .

VI.3 Comparison with numerical methods

We have compared the charge response function evaluated with the slave-boson method to exact diagonalization (ED) data Bun13 ; Jia12 ; Wan14 and quantum Monte Carlo (QMC) simulations Zim97 ; DziUnpub ; Che94 ; Buh99 ; Bec00 ; Pre97 ; Gro00 ; Koh04 ; Kun15 available in the literature. The low-temperature phase found by the numerical methods at half filling is an antiferromagnetic insulator. But, as confirmed by our investigation of instabilities, the paramagnetic solution becomes predominant with increasing doping and temperature. Keeping this in mind when comparing our evaluation of the charge response, we note that the spectra computed at finite doping by both numerical approaches show salient features that can be naturally explained by the two collective modes found in the present work. In particular the variations of their dispersions with the coupling and the doping qualitatively agree with the behavior we have described.

ED are performed on finite clusters and the small size of the system enhances the energy level separation. As a result the obtained spectrum is a set of peaks rather than a continuous function of the frequency. The energy quantization is visible in the spectra of the charge susceptibility calculated in Ref. Bun13 at small density $\langle\cali{N}\rangle\approx 0.03$. They remarkably show two distinctive peaks at the energies where we have found the ZS peak at $\bf{k}=$X and the UHB peak. Confirming our results, the first peak is exactly in the middle of the main contribution to the momentum-integrated response, that corresponds to the continuum of single-particle excitations. As for the second peak at higher energy, we note that the dispersion of the ED computation is narrow and the peak weight is vanishingly small, which can be explained by the UHB mode found by our theory close to doping $|\delta|=1$. The charge response function has also been computed around half filling, but for the Hubbard model including hopping between next-nearest-neighbor sites Jia12 ; Wan14 . The latter is known to break particle-hole symmetry. So the comparison with our results for the simple Hubbard model should be taken with caution. One can nevertheless remark a encouraging agreement for hole doping. For the large value of coupling $U=10t$ the ZS-like structure at the boundary of the continuum is found to decrease in energy with increasing hole doping, which is also predicted by the slave-boson method. Besides the high-energy feature moves up in energy as the UHB mode.

Early QMC simulations of the Hubbard model have been mainly focused on the static spin and charge structure factors Che94 ; Buh99 ; Bec00 . As previously discussed in Ref. Zim97 , the SRI slave-boson approach is in a very good quantitative agreement with the numerical evaluations of these quantities. Concerning the dynamical response functions, and in particular the search for collective modes, the analysis of the QMC results meets two hurdles. Indeed statistical averages computed by QMC simulations yield the values of the correlation functions on the imaginary-frequency axis. Their values on the real-frequency axis are then approximated by different numerical schemes, such as the maximum entropy method, which limits the obtained frequency definition. Furthermore the simulations of a doped system are restricted to the high-temperature regime $T\gtrsim t/3$ by the sign problem. As a consequence QMC spectra may lack the necessary energy resolution to distinguish fine structures, such as several collective-mode peaks close to one another, or a peak with a small weight which is the case of the UHB mode for a large set of parameters.

The charge response function of a doped system is computed in the QMC simulations Pre97 ; Gro00 ; Koh04 ; Kun15 at temperature $T=t/3$ for coupling values $U=4t$ and $8t$. At this temperature a sensible comparison with our theory may be made for doping $|\delta|\gtrsim 0.1$ at which the paramagnetic phase should prevail. The slave-boson results are consistent with the obtained QMC spectra. The latter show that the continuum response is strongly reduced at low doping and the intensity is mainly located beyond it, around the energies of the ZS and the UHB modes. For instance for ${\bf k}=M$ and $U=8t$, the intensity mainly spreads from $\omega\approx 4t$ to $\omega\approx 12t$. This can be interpreted as the response of the collective modes that interact with the background of incoherent multi-particle processes. By increasing the doping, the UHB mode energy increases and because of its small weight, its signature can no longer be distinguished from the structureless background in the QMC spectra. Meanwhile the continuum response is less renormalized away from half filling and the ZS mode energy decreases. The most satisfying comparison is found with the QMC simulations of Ref. Koh04 performed for $U=4t$. The spectra show two clear structures, one similar to the ZS peak at the edge of the continuum response, and the other one around $\omega\approx 8t$ which possesses a slight dispersion as the UHB mode.

VII Conclusion

We have derived the expression of the charge susceptibility of the Hubbard model in its Kotliar and Ruckenstein slave-boson representation. We have shown that it reduces to the conventional HF + RPA result when expanded to lowest order in $U$. They markedly depart from one another already to next order in $U$. We then investigated spin and charge instabilities as well as charge collective modes of the 2D Hubbard model in the thermodynamical limit. To that aim we used the spin rotation invariant formulation of the above representation. Extending previous work, our calculations showed that magnetic instabilities of the paramagnetic phase essentially disappear for temperature $T\geq t/6$, which lays ground for the computation of the charge susceptibility in this regime. In the strong coupling regime, the charge excitation spectrum splits into a low frequency branch, and a high frequency collective mode. En passant, an approximated analytical form of the latter has been derived. It applies to arbitrary lattices containing one site in the unit cell. This mode, that may not be accounted for within the conventional HF + RPA framework or self-consistent perturbative schemes such as FLEX, disperses around $\omega\simeq U$ and therefore follows from the upper Hubbard band.

At low energy the charge excitations form a continuum, which width scales with the quasiparticle residue $z_{0}^{2}$, again in contrast to the conventional HF + RPA framework result. A collective mode lies above its upper boundary. The velocity of this zero-sound mode is anisotropic both off half-filling and away from the low density limit. We did not find a universal behavior in its dependence on the coupling strength, because it results from two opposite trends: on the one hand the increase of the effective mass reduces it while on the other hand the zero-sound excitation is shifted to higher energy. Nevertheless some trends could be identified; for instance it shows very small dependence on $U$ in the small density regime. Furthermore, for small to intermediate doping, the zero-sound velocity decreases, once $U$ exceeds the band width. To some extend, our results could be interpreted within Pines’ polarization potential theory. Indeed striking similarities are found at low frequency when the ZS and UHB modes are well split. Yet the polarization potential theory does not entail a UHB mode and it therefore fails to describe the regime where the ZS and the UHB modes strongly hybridize. We also studied the temperature dependence of the charge excitation spectrum. We found the small wavevector zero-sound excitation to broaden with increasing temperature, while the other features show little temperature dependence.

Acknowledgments

We gratefully thank D. Braak, T. Kopp, M. Raczkowski, and A.-M. S. Tremblay for several stimulating discussions. The authors acknowledge the financial support of the French Agence Nationale de la Recherche (ANR), through the program Investissements d’Avenir (ANR-10-LABX-09-01), LabEx EMC3, the Région Basse-Normandie, the Région Normandie, and the Ministère de la Recherche.

Appendix A Elements of the fluctuation matrix $S_{ij}$

The fluctuation matrix is symmetric except for off-diagonal elements $S_{\mu 3}(k)=-S_{3\mu}(k)$. It is composed of four blocks, one for the charge fluctuations and three for the spin fluctuations.

As emphasized by Kot92 ; FW , it is essential to notice the absence of a full radial gauge in order to describe the UHB mode. Indeed, in the early calculations Ras88 ; Lav90 ; Lil90 ; li94 the erroneous conclusion that the phase of all slave-boson fields could be gauged away resulted in a $5\times 5$ matrix for the fluctuation matrix in the charge channel. However following the observation that one slave-boson field has to be complex yields a $6\times 6$ matrix that possesses the supplementary dynamics introduced by the time derivative of this boson field Jol91 . As a result the charge susceptibility acquires an $\omega^{2}$-dependence in addition to the frequency dependence contained in the fermionic bubbles $\chi_{m}(k)$, and a second pole describing the UHB mode. The charge fluctuation matrix has thus been obtained in the limit ${\bf q}=0$ and $\delta=0$ where the softening of the UHB mode has been found at the Mott-Hubbard transition Kot92 . Later a general expression of $\chi_{c}(k)$ for arbitrary momentum and density has been derived within the SRI representation Zim97 . But, as stated earlier, it does not include several matrix elements that are present in the correct expression (20). It turns out that the missing terms do not contribute to the correlation functions in the limits $\omega=0$ or ${\bf q}=0$, which may explain why they have been overlooked until now. However they are crucial to reproduce the RPA result at weak coupling.

The non-zero terms of the charge part are

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The spin blocks are given by

[TABLE]

[TABLE]

[TABLE]

We have used

[TABLE]

[TABLE]

with the fermionic bubbles

[TABLE]

and $\varepsilon_{\bf k}$ is given by Eq. (27).

The expressions of the derivatives of $z$ may be gathered from Ref. li91 ; Zim97 . Note however a misprint in Zim97 which should be corrected as

[TABLE]

Appendix B Spin susceptibility $\chi_{s}(k)$

As shown by Lav90 ; li91 ; Zim97 , inverting the fluctuation matrix $S$ yields the spin dynamical response function

[TABLE]

where

[TABLE]

Similarly to the charge dynamical response function, we have found that in the weak-coupling limit the expression can be simplified as

[TABLE]

where $f^{a}(k)$ can be expanded to the second order in $U$ as

[TABLE]

with $\gamma(k)=\chi_{1}(k)/\varepsilon_{0}\chi_{0}(k)$. Reducing this result to first order in $U$ yields an exact agreement with the perturbation theory, and the Landau’s Fermi-liquid spin parameter is obtained with $F_{0}^{a}=N_{F}f^{a}(0)$.