# Pseudogaps in the 2D half-filled Hubbard model

## Introduction

For over a decade it has been recognized that the normal state properties of high-$`T_c`$ superconductors are unusual and appear to have non-Fermi liquid characteristics. One of the most remarkable features of the normal state is a suppression of the density of states at the Fermi energy in a temperature regime above $`T_c`$ in underdoped samples. Angular resolved photoemission experiments show that this pseudogap in the spectral function has a d-wave anisotropy, the same symmetry as the superconducting order parameter in these materials. This, along with theories that short-ranged spin fluctuations mediate pairing in the high-$`T_c`$ cuprates, emphasizes the importance of understanding the normal state, insulating phase.

It is thought by many that the two-dimensional Hubbard model, or closely related models, should capture the essential physics of the high-$`T_c`$ cuprates.  Yet, despite years of effort, neither the precursor pseudogap nor d-wave superconducting order have been conclusively seen in the Hubbard model.

Intuitively, one may expect that the Hubbard model should show pseudogap behavior. At half-filling, the ground state of the 2D Hubbard model is an antiferromagnetic insulator and the spectrum is therefore gapped. However, the Mermin-Wagner theorem precludes any transition at finite $`T`$, so as the temperature is lowered one may anticipate that a pseudogap will develop. This question has been previously addressed in the 2D Hubbard insulator by finite-size lattice Quantum Monte Carlo (QMC) and approximate many-body techniques . The results have been contradictory and inconclusive as to the existence of a pseudogap at low temperatures, due to limitations of these techniques.

Using the recently developed Dynamical Cluster Approximation (DCA) we find that at sufficiently low temperatures a pseudogap opens in the single particle spectral weight $`A(𝐤,\omega )`$ of the 2D Hubbard model with a simultaneous destruction of the Fermi liquid state due to critical fluctuations above the $`T=0`$ transition temperature. This occurs in the weak-to-intermediate coupling regime $`U<W`$, where $`U`$ is the on-site Coulomb energy and $`W`$ the non-interacting band width.

Using finite-sized techniques, it is difficult to determine if a gap persists in the thermodynamic limit. At half filling, finite-size QMC calculations display a gap in their spectra as soon as the correlation length exceeds the lattice size, so they tend to overestimate the pseudogap as it would appear in the thermodynamic limit. Finite-size scaling is complicated by the lack of an exact scaling ansatz for the gap and the cost of performing simulations of large systems. Calculations employing Dynamical Mean Field Approximation (DMFA) in the paramagnetic phase do not display this behavior since they take place in the thermodynamic limit rather than on a finite-size lattice. However, the DMFA lacks the non-local spin fluctuations often believed to be responsible for the pseudogap. The Dynamical Cluster Approximation (DCA) is a fully causal approach which systematically incorporates non-local corrections to the DMFA by mapping the problem onto an embedded impurity cluster of size $`N_c`$. $`N_c`$ determines the order of the approximation and provides a systematic expansion parameter $`1/N_c`$. While the DCA becomes exact in the limit of large $`N_c`$ it reduces to the DMFA for $`N_c=1`$. Thus, the DCA differs from the usual finite size lattice calculations in that it is a reasonable approximation to the lattice problem even for a “cluster” of a single site. Like the DMFA, the DCA solution remains in the thermodynamic limit, but the dynamical correlation length is restricted to the size of the embedded cluster. Thus the DCA tends to underestimate the pseudogap.

## Method

The DCA is based on the assumption that the lattice self energy is weakly momentum dependent. This is equivalent to assuming that the dynamical intersite correlations have a short spatial range $`bL/2`$ where $`L`$ is the linear dimension of the cluster. Then, according to Nyquist’s sampling theorem, to reproduce these correlations in the self energy, we only need to sample the reciprocal space at intervals of $`\mathrm{\Delta }k2\pi /L`$. Therefore, we could approximate $`G(𝐊+\stackrel{~}{𝐤})`$ by $`G(𝐊)`$ within the cell of size $`\left(\pi /L\right)^D`$ (see, Fig. 1) centered on the cluster momentum $`𝐊`$ (wherever feasible, we suppress the frequency labels) and use this Green function to calculate the self energy. Knowledge of these Green functions on a finer scale in momentum is unnecessary, and may be discarded to reduce the complexity of the problem. Thus the cluster self energy can be constructed from the coarse-grained average of the single-particle Green function within the cell centered on the cluster momenta:

$$\overline{G}(𝐊)\frac{N_c}{N}\underset{\stackrel{~}{𝐤}}{}G(𝐊+\stackrel{~}{𝐤}),$$
(1)   
where $`N`$ is the number of points of the lattice, $`N_c`$ is the number of cells in the cluster, and the $`\stackrel{~}{𝐤}`$ summation runs over the momenta of the cell about the cluster momentum $`𝐊`$ (see, Fig. 1). For short distances $`rL/2`$ the Fourier transform of the Green function $`\overline{G}(r)G(r)+𝒪((r\mathrm{\Delta }k)^2)`$, so that short ranged correlations are reflected in the irreducible quantities constructed from $`\overline{G}`$; whereas, longer ranged correlations $`r>L/2`$ are cut off by the finite size of the cluster.

This coarse graining procedure and the relationship of the DCA to the DMFA is illustrated by a microscopic diagrammatic derivation of the DCA. For Hubbard-like models, the properties of the bare vertex are completely characterized by the Laue function $`\mathrm{\Delta }`$ which expresses the momentum conservation at each vertex. In a conventional diagrammatic approach $`\mathrm{\Delta }(𝐤_1,𝐤_2,𝐤_3,𝐤_4)=_𝐫\mathrm{exp}\left[i𝐫(𝐤_1𝐤_2+𝐤_3𝐤_4)\right]=N\delta _{𝐤_1+𝐤_2,𝐤_3+𝐤_4}`$ where $`𝐤_1`$ and $`𝐤_2`$ ($`𝐤_3`$ and $`𝐤_4`$) are the momenta entering (leaving) each vertex through its legs of $`G`$. However as $`D\mathrm{}`$ Müller-Hartmann showed that the Laue function reduces to

$`\mathrm{\Delta }_D\mathrm{}(𝐤_1,𝐤_2,𝐤_3,𝐤_4)=1+𝒪(1/D)\text{.}`$ (2)   
The DMFA assumes the same Laue function, $`\mathrm{\Delta }_{DMFA}(𝐤_1,𝐤_2,𝐤_3,𝐤_4)=1`$, even in the context of finite dimensions. Thus, the conservation of momentum at internal vertices is neglected. Therefore we may freely sum over the internal momenta at each vertex in the generating functional $`\mathrm{\Phi }_{DMFA}`$. This leads to a collapse of the momentum dependent contributions to the functional $`\mathrm{\Phi }_{DMFA}`$ and only local terms remain.

The DCA systematically restores the momentum conservation at internal vertices. As discussed above the Brillouin-zone is divided into $`N_c=L^D`$ cells of size $`(2\pi /L)^D`$. Each cell is represented by a cluster momentum $`𝐊`$ in the center of the cell. We require that momentum conservation is (partially) observed for momentum transfers between cells, i.e. for momentum transfers larger than $`\mathrm{\Delta }k=2\pi /L`$, but neglected for momentum transfers within a cell, i.e less than $`\mathrm{\Delta }k`$. This requirement can be established by using the Laue function 

$$\mathrm{\Delta }_{DCA}(𝐤_1,𝐤_2,𝐤_3,𝐤_4)=N_c\delta _{𝐌(𝐤_1)+𝐌(𝐤_3),𝐌(𝐤_2)+𝐌(𝐤_4)}$$
(3)   
where $`𝐌(𝐤)`$ is a function which maps $`𝐤`$ onto the momentum label $`𝐊`$ of the cell containing $`𝐤`$ (see, Fig. 1).

With this choice of the Laue function the momenta of each internal leg may be freely summed over the cell. This is illustrated for the second-order term in the generating functional in Fig. 2. Thus, each internal leg $`G(𝐤_1)`$ in a diagram is replaced by the coarse–grained Green function $`\overline{G}(𝐌(𝐤_1))`$, defined by Eq. 1

The diagrammatic sequences for the generating functional and its derivatives are unchanged; however, the complexity of the problem is greatly reduced since $`N_cN`$. We showed previously that the DCA estimate of the lattice free-energy is minimized by the approximation $`\mathrm{\Sigma }(𝐤)\overline{\mathrm{\Sigma }}(𝐌(𝐤))`$, where $`\delta \mathrm{\Phi }_{DCA}/\delta \overline{G}=\overline{\mathrm{\Sigma }}`$.

The cluster problem is then solved by usual techniques such as QMC, the non-crossing approximation or the Fluctuation-Exchange approximation. Here we employ a generalization of the Hirsh-Fye QMC algorithm to solve the cluster problem. The initial Green function for this procedure is the bare cluster Green function $`𝒢(𝐊)^1=\overline{G}(𝐊)^1+\overline{\mathrm{\Sigma }}(𝐊)`$ which must be introduced to avoid over-counting diagrams. The QMC estimate of the cluster self energy is then used to calculate a new estimate of $`\overline{G}(𝐊)`$ using Eq. 1. The corresponding $`𝒢(𝐊)`$ is used to reinitialize the procedure which continues until the self energy converges to the desired accuracy.

## Results

We study the 2D Hubbard Hamiltonian:

$`H`$ $`=`$ $`t{\displaystyle \underset{i,j,\sigma }{}}(c_{i\sigma }^{}c_{j\sigma }+c_{j\sigma }^{}c_{i\sigma })`$ (5)   
$`+U{\displaystyle \underset{i}{}}(n_i{\displaystyle \frac{1}{2}})(n_i{\displaystyle \frac{1}{2}})\mu {\displaystyle \underset{i,\sigma }{}}n_{i\sigma }.`$   
where $`c_{i\sigma }^{}(c_{i\sigma })`$ creates (destroys) an electron at site $`i`$ with spin $`\sigma `$, $`U`$ is the on-site Coulomb potential, and $`n_{i\sigma }=c_{i\sigma }^{}c_{i\sigma }`$ is the number operator. We set the overlap integral $`t=1`$ and measure all energies in terms of $`t`$. We work at $`\mu =0`$ where the system is half-filled ($`n=1`$). We choose $`U=5.2`$, which is well below the value $`UW`$ believed to be necessary to open a Mott-Hubbard gap.

We also calculate the angle integrated dynamical spin susceptibility shown in the inset. It does not have a pseudogap, as expected for the half-filled model since the spin-wave spectrum is gapless. Since a (spin) charge gap is generally defined as one which appears in the (spin) charge dynamics or thermodynamics, we conclude that the pseudogap is only in the charge response and is due to short-ranged antiferromagnetic spin correlations.

Fig. 3 shows the spectral density $`A(𝐤,\omega )`$, and the real $`\mathrm{Re}\mathrm{\Sigma }(𝐤,\omega )`$ and imaginary $`\mathrm{Im}\mathrm{\Sigma }(𝐤,\omega )`$ parts of the self-energy for the 2D Hubbard model via the DCA with a paramagnetic host at the Fermi surface $`X`$ point $`𝐤=(\pi ,0)`$ for a 64-site cluster ($`N_c=64`$) at various temperatures. We obtain the spectral function $`A(𝐤,\omega )`$ via the Maximum Entropy Method (MEM). As the temperature is lowered, the system first builds a Fermi-liquid-like peak in $`A(𝐤,\omega )`$. By $`\beta =2.6`$, a pseudogap begins to develop in $`A(𝐤,\omega )`$. The pseudogap builds as the temperature is further lowered.

Fig. 4 shows the spectral function $`A(𝐤,\omega )`$ at the half-filled Fermi surface point $`𝐤=(\pi /2,\pi /2)`$. The qualitative features are similar to Fig. 3, but the pseudogap opens at a lower temperature and the distance between the peaks is less than that seen at the $`X`$ point. This behavior is reminiscent of the anisotropy of the pseudogap observed experimentally in the insulating and in the superconducting cuprates, but is not large enough to be comparable with that seen experimentally. We speculate that the anisotropy seen here may be due to a difference in the number of states near the Fermi energy at these two points in the zone.

The DCA self-energy spectra Figs. 3 & 4 support the spectral evidence. At the $`X`$ point, the slope of the real part $`\mathrm{Re}\mathrm{\Sigma }(𝐤,\omega )`$ becomes positive below $`\beta =2.6`$, the temperature at which we observed the opening of a pseudogap. This signals the appearance of two new solutions in the quasiparticle equation $`\mathrm{Re}(\omega ϵ_𝐤\mathrm{\Sigma }(𝐤,\omega ))=0`$ in addition to the strongly damped solution at $`\omega =0`$ which is also present in the noninteracting system. These two new quasiparticle solutions for the same $`𝐤`$-vector indicate precursor effects of the onset of antiferromagnetic ordering which entails a doubling of the unit cell. They are referred as shadow states and are caused by antiferromagnetic spin fluctuations in the paramagnetic state. At these temperatures, the imaginary part $`\mathrm{Im}\mathrm{\Sigma }(𝐤,\omega )`$ displays a local minimum at $`\omega =0`$ indicating the breakdown of the Fermi liquid behavior. We note that a different conclusion was previously reached in a FLEX study, which found that $`\mathrm{Im}\mathrm{\Sigma }(𝐤,\omega )`$ has a local minimum at $`\omega =0`$ which was not accompanied by an opening of a pseudogap. Since the pseudogap is due to short-range spin correlations, we conclude that FLEX underestimates these correlations.

It is instructive to compare the DCA results with those obtained by finite-size QMC calculations. Fig. 5 shows the spectral density $`A(\pi ,0,\omega )`$ obtained by analytically continuing both finite-size and DCA QMC data. In spite of the difference in the two methods, the information they provide is complimentary. In the finite-size results, (b), we see a similar opening of a pseudogap. However, as the length of the antiferromagnetic (AF) correlations reach the longest length on the finite-size lattice, the system develops a full gap. Thus, the finite-size QMC overestimates the size of the gap. In the DCA results, (a), the pseudogap emerges as soon as $`N_c>1`$. The temperature $`T^{}`$ at which the pseudogap first becomes apparent in the spectra, as well as the full width $`\mathrm{\Delta }`$ measured from peak to peak is plotted in the inset. Both $`T^{}`$ and $`\mathrm{\Delta }`$ increase with $`N_c`$. Since the DCA calculation remains in the thermodynamic limit, a full gap due to antiferromagnetic correlations alone cannot open until their correlation length diverges. However, since these correlations are restricted to the size of the cluster, the DCA systematically underestimates the size of the gap. Thus, if a pseudogap exists in the DCA for finite $`N_c`$, it should persist in the limit as $`N_c\mathrm{}`$.

In summary, we have employed the recently developed DCA to study the long-open question of whether the half-filled Hubbard model has a pseudogap due to AF spin fluctuations. We find conclusive evidence of a pseudogap in the charge dynamics and have shown unambiguously that the $`T=0`$ phase transition of the half-filled model is preceded by an opening of a pseudogap in $`A(𝐤_F,\omega )`$ accompanied by pronounced non-Fermi liquid behavior in $`\mathrm{\Sigma }(𝐤_F,\omega )`$.

## Acknowledgments

We would like to acknowledge useful conversations with P. van Dongen, B. Gyorffy, M. Hettler, H.R. Krishnamurthy R.R.P. Singh and J. Zaanen. This work was supported by the National Science Foundation grants DMR-9704021, DMR-9357199, and the Ohio Supercomputing Center.