Magnetic Correlations in the Two-dimensional Repulsive Fermi Hubbard Model

TL;DR

This paper investigates magnetic fluctuations in the two-dimensional repulsive Fermi Hubbard model near half-filling, identifying the temperature scales and interaction strengths where magnetic correlations are most prominent, relevant for cold atom experiments.

Contribution

The study employs a self-consistent diagrammatic approach to quantify magnetic fluctuation temperature scales and their dependence on doping and interaction strength in the Hubbard model.

Findings

Maximum magnetic fluctuations occur at U/t ≈ 4-6.

The characteristic temperature T_M(n) varies with doping levels.

Magnetic fluctuations are strongest near U/t ≈ 4-6, facilitating experimental observation.

Abstract

The repulsive Fermi Hubbard model on the square lattice has a rich phase diagram near half-filling (corresponding to the particle density per lattice site $n=1$): for $n=1$ the ground state is an antiferromagnetic insulator, at $0.6 < n \lesssim 0.8$, it is a $d_{x^2-y^2}$-wave superfluid (at least for moderately strong interactions $U \lesssim 4t$ in terms of the hopping $t$), and the region $1-n \ll 1$ is most likely subject to phase separation. Much of this physics is preempted at finite temperatures and to an extent driven by strong magnetic fluctuations, their quantitative characteristics and how they change with the doping level being much less understood. Experiments on ultra-cold atoms have recently gained access to this interesting fluctuation regime, which is now under extensive investigation. In this work we employ a self-consistent skeleton diagrammatic approach to quantify the characteristic temperature scale $T_{M}(n)$ for the onset of magnetic fluctuations with a large correlation length and identify their nature. Our results suggest that the strongest fluctuations---and hence highest $T_{M}$ and easiest experimental access to this regime---are observed at $U/t \approx 4-6$.

Figures

Equations

$$H=-t\sum _{⟨i,j⟩\sigma }{\hat{c}}_{i,\sigma }^{†}{\hat{c}}_{j,\sigma }+U\sum _i{\hat{n}}_{i,\uparrow }{\hat{n}}_{i\downarrow }-\mu \sum _{i,\sigma }{\hat{n}}_{i,\sigma }$$

$$G_{\sigma }=\frac{G_{\sigma }^{(0)}}{1-G_{\sigma }^{(0)}{\Sigma }_{\sigma }},W=\frac{U}{1-U\Pi }-U,$$

$$G_{pp}=\frac{U}{1-U{\Sigma }_{pp}}-U,G_{ph}=\frac{U}{1-U{\Sigma }_{ph}}-U,$$

$$\chi ({\omega }_n,\vec{q})=\mathrm{Tr}{S}_z\frac{\Pi ({\omega }_n,\vec{q})}{1-\Pi ({\omega }_n,\vec{q})U}{S}_z,$$

Full text

The fermionic Hubbard model Hubbard (1963); Anderson (1963); Anderson et al. (1997); Lee et al. (2006), defined by the square lattice Hamiltonian

[TABLE]

has for years played a crucial role in studies of correlated electrons in solids; it is regarded as one of the “standard models” of condensed matter physics to introduce and discuss Mott insulating phases, antiferromagnetic (AFM) correlations, novel mechanisms of superconducting pairing, non-Fermi-liquid behavior, etc. In Eq. (1) and in what follows, the nearest-neighbor hopping amplitude $t$ is set to be the energy and temperature unit (distances are measured in units of the lattice constant); $U$ is the on-site repulsive coupling constant; $\mu$ is the chemical potential; $\hat{c}^{\dagger}_{i,\sigma}$ and $\hat{c}_{i,\sigma}$ create and annihilate (respectively) a fermion of the spin component $\sigma\in\{\uparrow,\downarrow\}$ on the site $i$, and $\hat{n}_{i,\sigma}$ counts the number of fermions of a particular spin on a given lattice site.

On the one hand, Eq. (1) involves a number of crucial simplifications that make it qualitatively different from real materials, such as high-$T_{c}$ superconductors. The most important ones include (i) strictly two dimensional (2D) as opposed to the strongly anisotropic 3D geometry, (ii) neglect of long-range Coulomb interactions, (iii) suppression of hopping matrix elements beyond the nearest-neighbor ones ($t^{\prime}=0$), (iv) single-band approximation, and (v) absence of electron-phonon coupling. Correspondingly, the model (1) cannot feature an ordered AFM phase at a non-zero temperature, but is allowed to have a first-order transition between phases with different electron densities, not to mention that $t^{\prime}=0$ leads to the Fermi surface nesting and particle-hole symmetry at $n=\left<n_{i\uparrow}+n_{i\downarrow}\right>=1$. As a result, the schematic phase diagram of (1) in the doping-temperature plane shown in Fig. 1 (see discussion below) is distinct from the “canonical” picture of high-$T_{c}$-type materials Lee et al. (2006). On the other hand, advances in ultra-cold atomic experiments have made it possible to accurately emulate the model (1) on optical lattices Jaksch et al. (1998); Lewenstein et al. (2007); Köhl et al. (2005); Jördens et al. (2008); Schneider et al. (2008); Hulet et al. (2015); Greif et al. (2015); Parsons et al. (2016); Mazurenko et al. , bringing ultra-cold atom experiments into the region of the phase diagram (Fig. 1), where they can cooperate with the state-of-the-art numerical methods to reveal the underlying physics. Numerical results can now be directly compared to experiments and vice versa dramatically increasing the importance of producing reliable data sets.

Recent years have seen a remarkable progress in unveiling the $T=0$ phase diagram of the Hubbard model (1). Well-understood regions include the limit of vanishing densities $n\rightarrow 0$ Kagan and Chubukov (1989); Baranov and Kagan (1992); Chubukov and Lu (1992); Chubukov (1993); Zanchi and Schulz (1996); Halboth and Metzner (2000); Fukazawa and Yamada (2002) and vanishing interaction strength Hlubina (1999); Raghu et al. (2010); Rømer et al. (2015); Šimkovic IV et al. (2016). For densities $n\leq 0.7$ and coupling strength $U\leq 4$, the ground state is a BCS superfluid (with the $d_{x^{2}-y^{2}}$-wave symmetry at density $n>0.6$) Deng et al. (2015). At half-filling $n=1$, the ground state is an AFM insulator for any $U$ Hirsch (1985); Hirsch and Tang (1989); White et al. (1989); Furukawa and Imada (1992); Varney et al. (2009). Being a qualitative property, AFM order can only disappear (with doping) by a quantum phase transition, with the simplest scenario being that of phase separation (PS). The first PS state was proposed to be a mixture of AFM and ferromagnetic (FM) orders in the region of small doping $\delta=1-n\ll 1$ and large $U$ Penn (1966); Nagaoka (1966); Tasaki (1998); this conjecture has been later supported numerically for $U>25$ Zitzler et al. (2002). The instability of the model towards incommensurate AFM and domain wall formation was also reported in Refs. Schulz (1990); Lin (1991). Recently, PS for small values of $U$ was observed in Auxiliary-field Quantum Monte Carlo Chang and Zhang (2008); Sorella (2015) and variational Tocchio et al. (2016) studies.

Much less is known conclusively about the finite-temperature behavior. Given that the correlation length $\xi$ for AFM correlations at $n=1$ diverges exponentially fast when $T\to 0$, there exists a relatively high temperature $T_{M}$ below which magnetic correlations extend over many lattice sites and electronic degrees of freedom are getting locked in collective modes. The characteristic temperature scale $T_{M}(n)$ is supposed to decrease with increasing the doping level. Since magnetic correlations and fluctuations are considered to be the prime reason for PS and BCS phenomena near $n=1$, quantifying their behavior is paramount to understanding the finite-$T$ phase diagram of the Fermi Hubbard model in 2D.

In this Letter, we employ a skeleton diagrammatic approach to quantify the $T_{M}(n)$ scale, see Fig. 2, and shed light on the structure of dominant fluctuations. Experiments with ultra-cold ${}^{6}\text{Li}$ fermions have now reached temperatures low enough to directly observe magnetic correlations in the model (1) Hulet et al. (2015); Greif et al. (2015); Parsons et al. (2016), with the lowest temperatures attained on the scale of $T\sim 0.25$ Mazurenko et al. .

Method. The imaginary-time spin correlation function $\chi(\tau,i)=\langle S_{z}(\tau,i)S_{z}(0,0)\rangle$ and its Fourier transform $\chi(\omega_{m},\vec{q})$ at bosonic Matsubara frequencies $\omega_{m}=2\pi Tm$ ($m$ is an integer) was computed within the skeleton diagrammatic framework based on self-consistently renormalised (“dressed”) elements in four separate channels: particle-particle and particle-hole pair propagators, screened interaction, and single-particle propagator.

To define the framework, we first note that in the standard weak-coupling expansion in powers of $U$ (for introduction see, e.g., Abrikosov et al. (1975)) every instance of the bare interaction vertex can be replaced by either of the three types of infinite sums shown in Fig. 3. These sums originate from three possible ways of connecting bare vertexes by non-interacting Green’s functions $G^{(0)}_{\sigma}$ to form a geometric series and are commonly referred to as particle-particle ($G_{pp}^{(0)}$), particle-hole ($G_{ph}^{(0)}$), and bubble (screened interaction $W_{\sigma\sigma^{\prime}}^{(0)}$) ladders. Note that the functions $G_{pp}^{(0)}$, $G_{ph}^{(0)}$, and $W_{\sigma\sigma^{\prime}}^{(0)}$, have the same structure as the single-particle propagators $G_{\sigma}^{(0)}$; i.e., thanks to the local nature of the Hubbard interaction $U$, they depend only on one lattice coordinate and time. Therefore, they are represented diagrammatically as lines, whereas the bare vertex is a point 111Bare interaction needs to be treated separately to define precise rules for avoiding double counting.; the complete set of elements is shown schematically in Fig. 4 (left panel). An arbitrary diagram can now be constructed by taking any number of these elements and connecting their incoming and outgoing ends with propagator lines, as exemplified in the the right panel of Fig. 4. The resulting series contains significantly fewer terms because a large fraction of diagram topologies in the weak-coupling expansion are accounted for by the ladder sums.

The standard step that allows further reduction of the diagrammatic space is replacing all the bare lines with the ones dressed by an infinite geometric series of all sorts of diagrammatic insertions, referred to as self-energies, and dropping a substantial fraction of diagram topologies that are double-counted as a result. This leads to the self-consistent so-called skeleton technique Abrikosov et al. (1975), where one computes the diagrammatic series for the self-energies $\Sigma_{\sigma}$, $\Sigma_{pp}$, $\Sigma_{ph}$ and $\Pi_{\sigma\sigma^{\prime}}$ constructed from some approximation to the dressed lines $G_{\sigma}$, $G_{pp}$, $G_{ph}$, $W_{\sigma\sigma^{\prime}}$ and obtains the next approximation by solving the corresponding Dyson equations:

[TABLE]

[TABLE]

which are simple algebraic relations in the momentum/frequency representation. The diagrammatic calculation of self-energies built on the solutions of Eqs. (2), (3) for all the diagram lines is repeated iteratively until convergence. The spin-spin correlation function is then directly related to the polarization $\Pi$ by

[TABLE]

with the trace taken over the spin index. This scheme can be abbreviated as GGGW to emphasise the four renormalization channels.

To avoid double counting, diagrams for the self-energies must be constrained to the skeleton set in all the dressed channels: they must remain connected after cutting any two lines of the same kind; an example is shown in Fig. 4. In addition, one has to enforce two rules concerning dots: (1) no two dots can be connected directly by two Green’s functions, and (2) a dot cannot be connected by two Green’s functions to the same end of $G_{pp}$ or $G_{ph}$ or $W$. This is necessary because, by construction, such diagrams are already accounted for in the corresponding ladder sums. Finally, there is one exception to the rule: to avoid triple counting of the same diagram contributing to the lowest-order $\Sigma_{\sigma}$, one has to perform subtraction of the diagram based on two points, see Fig. 5.

All our results are based on the lowest-order GGGW scheme, in which the self-energies are given by the skeleton diagrams shown in Fig. 5 and the lines are computed self-consistently by Eqs. (2), (3). In principle, higher-order skeleton graphs can be summed by the standard diagrammatic Monte Carlo (DiagMC) scheme Van Houcke et al. (2010); Kozik et al. (2010); Kulagin et al. (2013); Deng et al. (2015) with obvious modifications required to handle a larger set of diagrammatic elements and self energies. We have implemented the DiagMC scheme and used it to assess systematic errors of the lowest-order approximation.

Formally, the exact answer follows from summing the whole infinite series of all skeleton diagrams for the self-energies. However, it is known that in strongly correlated regimes close to half-filling ($U\gtrsim 4$, $T\lesssim 0.5$, $n\sim 1$) skeleton sums for the Hubbard model cannot converge to the correct answer being attracted to an unphysical branch of the Luttinger-Ward functional Kozik et al. (2015). On the other hand, at weaker interactions skeleton series quickly converge to the exact solution. Therefore, by continuity it is natural to expect that a low-order skeleton theory produces qualitatively and even quantitatively accurate results somewhat into the strongly correlated regime, provided the self-consistent dressing adequately captures fluctuations in relevant channels. We observe that it is indeed the case for the first-order GGGW approach employed here by benchmarking the results at $n=1$, which is notoriously the most difficult regime for skeleton schemes Kozik et al. (2015), against the numerically exact determinant diagrammatic Monte Carlo method Burovski et al. (2006); Kozik et al. (2013). In particular, we observe that in the range or interaction strengths considered here, our approach produces accurate thermodynamic observables at the level of a few percent. The functional form of $\Sigma_{\sigma}(\omega_{n},\vec{q})$ and $\chi(\omega_{n},\vec{q})$ displays all the features and the deviation of the overall amplitude from the exact answer is at most $30\%$. For the purposes of the qualitative analysis carried out below this level of accuracy is sufficient. This makes the lowest-order GGGW scheme a practical computationally-inexpensive tool that captures complex correlation effects in the regime of parameters accessible in current cold-atom experiments.

Results. The very notion of $T_{M}$ as the crossover temperature between the high-temperature regime, in which the magnetic correlation length $\xi$ is less than or comparable to the lattice constant $a$, and the regime of strong long-range correlations with $\xi\gg a$, implies that its definition is not unique. To characterise the onset of magnetic correlations we examine the momentum dependence of the static magnetic susceptibility $\chi(\omega=0,\vec{q})$ and monitor development of a narrow peak structure. We define $T_{M}$ at a given value of interaction $U$ and density $n$ as the highest temperature at which the amplitude of the peak in $\chi(\omega=0,\vec{q})$ is an order of magnitude larger than its minimum value over the Brillouin zone, $\chi_{\text{max}}/\chi_{\text{min}}=10$. Our results are summarized in Fig. 2.

As expected, the largest values of $T_{M}$ (at fixed $U$) are observed at half-filling where the crossover temperature can be as high as $T\sim 0.25$ (or about $1000$ K in units representative of the $CuO_{2}$ superconductors with hopping amplitute $t\sim 0.3\,eV$ Gull et al. (2013)). As a function of interaction, $T_{M}$ goes to zero at both large and small values of $U$, and features a smooth maximum around $U\sim 4$. This appears to be the optimal spot for experimental studies of magnetism in the Hubbard model (1) where reaching low temperatures remains challenging. The magnetic crossover scale eventually goes to zero with doping but remains relatively high for intermediate values of $U$ even at doping levels $\delta\in(0.15,0.25)$. We did not see evidence for PS at $T_{M}$, meaning that the PS dome takes place within the magnetic region, see Fig. 1.

The character of spin correlations undergoes a dramatic transformation with doping. A mismatch between the largest momentum transfer at the Fermi surface and the reciprocal lattice vector results in the incommensurate spin-wave fluctuations that take the form of AFM domains with diagonal domain walls (or “diagonal carpet”, for brevity). In Fig. 6, we show a typical example of the emerging spin structure (for $T=0.05$, $U=4$ and $n=0.8085$). In the left panel, we see that the otherwise dominant peak around the commensurate vector $\left(\pi,\pi\right)$ is split and features a minimum at $\left(\pi,\pi\right)$ surrounded by two maxima at the incommensurate vectors. The real-space spin texture behind this split-peak signal is shown in the right panel with different colors corresponding to the sign of $\chi(0,\vec{r})$. It is plausible that in the PS region, see Fig. 1, the AFM order is intermixed with the diagonal carpet, and the mechanism for the $d_{x^{2}-y^{2}}$-wave pairing is based on coupling to these spin fluctuations.

Conclusions. We discussed the finite-temperature phase diagram of the repulsive Fermi Hubbard model on a square lattice and identified the overarching dome defining the onset of strong magnetic correlations that change their structure from commensurate antiferromagnetism to incommensurate diagonal texture as the doping level is increased. Given relatively high values of $T_{M}$ that extend well into the doping region where optimal values for transition temperatures to the $d_{x^{2}-y^{2}}$ superfluidity are expected to take place, magnetic correlations appear to be the key ingredient behind both the PS and superfluidity near half-filling.

Further development of the Diagrammatic Monte Carlo approach is required to obtain controllable results at temperatures below $T_{M}$ where convergence of the diagrammatic expansion becomes problematic. Large magnetic correlation length and the possibility of phase separation should be treated with extreme care by any numerical method based on finite-cluster calculations because this physics imposes restrictions on the minimal acceptable cluster size and questions homogeneity of the solution. In particular, the superfluid states proposed in Gull et al. (2013); Chen et al. (2015); Zheng and Chan (2016) in the region of parameter space where one can expect phase separation Chang and Zhang (2008); Sorella (2015); Tocchio et al. (2016) could result from the intrinsic bias of computational technique. To find high-$T_{c}$ regions one has to avoid PS near half-filling by doping or modify the model to include non-zero values of the next-nearest-neighbor hopping $t^{\prime}>0$ Chen et al. (2015).

Ultra-cold atom experiments are expected to have a major impact on revealing the finite-temperature phase diagram. Current experiments have already reached temperatures $T\sim T_{M}$ Mazurenko et al. , and are well positioned to explore the structure of strong magnetic correlations. Detection and characterization of the PS state requires reaching lower temperature scales. However, there is no a priori reason for the PS dome to take place at $T\ll T_{M}$ given that known correlations saturate quickly below $T_{M}$.

Quantifying magnetic correlations is also of significant interest in relation with copper oxides Kastner et al. (1998) as neutron scattering experiments have revealed the coexistence of commensurate and incommensurate magnetic structures at finite doping. For $\text{La}_{2-p}\text{Sr}_{p}\text{Cu}_{4}$ an incommensurate state with a magnetic structure wave vector was found at small dopings; for $\text{YBa}_{2}\text{Cu}_{3}\text{O}_{6+y}$ a wide doping window exists where commensurate AFM fluctuations are observed at low temperatures Matsuda et al. (2002).

The authors are grateful to Shiwei Zhang, Sandro Sorella, Luca F. Tocchio, and Markus Greiner for useful discussions of their results. Fruitful exchanges with Ulrich Schollwöck about the possibility of phase separation at larger $U$’s are acknowledged. Fedor Šimkovic would like to thank USTC Hefei for generous hospitality during a period when parts of this paper were written. This work was supported by the Simons Collaboration on the Many Electron Problem, National Science Foundation under the grant PHY-1314735, the MURI Program “New Quantum Phases of Matter” from AFOSR, and the Swiss National Science Foundation, NSFC Grant No. 11625522, CAS, and NKBRSFC Grant No. 2016YFA0301600.