Extended crossover from Fermi liquid to quasi-antiferromagnet in the half-filled 2D Hubbard model

TL;DR

This paper investigates how increasing interactions in the half-filled 2D Hubbard model lead to a crossover from Fermi liquid to an insulating state driven by antiferromagnetic correlations, revealing a non-Fermi liquid intermediate phase.

Contribution

It demonstrates an extended crossover driven by AFM correlations using advanced Monte Carlo and cluster approximation methods, without a sharp metal-insulator transition.

Findings

Long-range AFM correlations induce a crossover from Fermi liquid to insulator.

An intermediate non-Fermi liquid state with a partially gapped Fermi surface is identified.

No sharp metal-insulator transition occurs in the studied parameter regime.

Abstract

The ground state of the Hubbard model with nearest-neighbor hopping on the square lattice at half filling is known to be that of an antiferromagnetic (AFM) band insulator for any on-site repulsion. At finite temperature, the absence of long-range order makes the question of how the interaction-driven insulator is realized nontrivial. We address this problem with controlled accuracy in the thermodynamic limit using self-energy diagrammatic determinant Monte Carlo and dynamical cluster approximation methods and show that development of long-range AFM correlations drives an extended crossover from Fermi liquid to insulating behavior in the parameter regime that precludes a metal-to-insulator transition. The intermediate crossover state is best described as a non-Fermi liquid with a partially gapped Fermi surface.

Equations

$$H=\sum _{\mathbf{k},\sigma }\left({ϵ}_{\mathbf{k}}-\mu \right)c_{\mathbf{k}\sigma }^{†}{c}_{\mathbf{k}\sigma }+U\sum _i{n}_{i\uparrow }{n}_{i\downarrow },$$

$${\Sigma }_{\mathbf{k}\sigma }(i\omega ,T,\mu =U/2,U)=\sum _{n=1}^{\infty }{a}_{2n,\mathbf{k}\sigma }(i\omega ,T)U^{2n}$$

Full text

The interaction driven metal-to-insulator (MIT) transition has been for many years a problem of central focus for the field of strongly correlated electron systems (see, e.g., Refs. Georges et al. (1996); Schäfer et al. (2015) and references therein). Particularly challenging has been the quantitative, and even qualitative, understanding of the MIT in two-dimensional systems. Here, the basic model—the single-band Hubbard model with nearest-neighbor hopping on the square lattice—can nowadays be accurately emulated and probed with ultracold atoms in optical lattices Bloch (2005); Lewenstein et al. (2007); Greif et al. (2015); Parsons et al. (2016); Mazurenko et al. (2017); Nichols et al. (2019) at ever decreasing temperatures, putting controlled experimental studies of the this problem within reach. The model is given by the Hamiltonian

[TABLE]

where $\mu$ is the chemical potential, $\mathbf{k}$ (quasi)momentum with the lattice constant set to unity, $n_{i\sigma}$ the number operator of fermions with spin $\sigma$ on the square lattice site $i$, $U$ the on-site repulsion, and $\epsilon_{\mathbf{k}}=-2t\left[\cos(k_{x})+\cos(k_{y})\right]$. For the description of the MIT, nonperturbative numerical methods, such as the dynamical mean-field theory (DMFT) and related cluster and diagrammatic extensions Georges et al. (1996); Maier et al. (2005); Park et al. (2008); Werner et al. (2009); Gull et al. (2013); Schäfer et al. (2015); Rohringer and Toschi (2016); Fratino et al. (2017); van Loon et al. (2018), have played a central role. In the single-site paramagnetic DMFT Georges et al. (1996) solution, which becomes exact in the limit of infinite dimensions, the metallic phase at half filling (the average density per site $\langle n\rangle=1$) persists down to zero temperature at weak interactions. It is separated from the Mott insulator by a first-order MIT at a sufficiently large value of $U=U_{c}$, ending at a finite temperature with an Ising critical point. Extensions of DMFT to small (up to 16 sites) real-space clusters Park et al. (2008); Werner et al. (2009); Gull et al. (2013) have shown that the inclusion of short-range spin fluctuations changes this picture substantially—a non-Fermi-liquid (nFL) state with a Fermi surface (FS) gap in certain momentum sectors continuously develops at a finite $U$ before the transition, the value of $U_{c}$ is reduced, and the slope of the first-order line is inverted.

It is, however, well known that the ground state of the model (1) at $\langle n\rangle=1$ is an antiferromagnet (AFM) at any $U>0$. As revealed by Slater Slater (1951), the FS nesting, i.e. the existence of a single wave vector $\mathbf{Q}=(\pi,\pi)$ that connects any point on the FS to another FS point, makes the interacting Fermi gas unstable against formation of the spin density wave with the wave vector $\mathbf{Q}$ already at infinitesimally small $U$. The corresponding unit-cell doubling makes the ground state a band insulator. While the Mermin Wagner theorem Mermin and Wagner (1966) forbids the long-range order at $T>0$, the AFM correlation length $\xi$ is exponentially large at low temperature, $\log\xi\propto t/T$; i.e., for practical purposes the system is best described as a quasi-AFM. Indeed, experiments with ultracold atoms Mazurenko et al. (2017) observed a perfect AFM state in model (1) on a $\sim 10\times 10$ lattice at temperatures as high as $T\sim 0.25t$. Such AFM correlations are explicitly truncated (and typically suppressed) in cluster DMFT calculations unless linear cluster sizes are comparable to $\xi$ LeBlanc and Gull (2013), which becomes computationally prohibitive at low $T$ even at half filling LeBlanc et al. (2015). Recent work Schäfer et al. (2015) based on the dynamical vertex approximation (D$\Gamma$A), a diagrammatic extension of DMFT capable of capturing long-range correlations approximately, and determinant quantum Monte Carlo (DQMC) simulations indicates that the low-temperature crossover from the Fermi liquid (FL) to the quasi-AFM insulator preempts and precludes the MIT.

In this Letter, we aim at establishing the picture of developing an insulating AFM state with controlled accuracy using the recently introduced $\Sigma$DDMC approach Šimkovic and Kozik (2019) (a similar approach was developed in Ref. Moutenet et al. (2018)). The method deterministically sums all topologies of Feynman diagrams for self-energy (for introduction, see, e.g., Ref. Abrikosov et al. (1975)) by means of determinants Rubtsov (2003); Rubtsov et al. (2005) with a recursive scheme in the spirit of Rossi’s algorithm Rossi (2017) to extract only one-particle irreducible diagrams. Integration over internal variables is performed by Monte Carlo sampling, and the thermodynamic limit (TDL) is taken explicitly. Compared to the standard DiagMC approach Van Houcke et al. (2010); Kozik et al. (2010), where diagram topologies are sampled stochastically, $\Sigma$DDMC enables access to substantially higher ($\sim 10-12$) expansion orders and more accurate determination of the self-energy.

Our main result is summarized in Fig. 1. At sufficiently low temperature, $T\lesssim 0.25t$, we observe the crossover from the low-$U$ metallic FL state with a well-defined FS (the red region) to the quasi-AFM insulator with temperature-activated quasiparticles (the blue region). In between, there is a transitional nFL regime (the gray region) bound by two lines $T_{\text{an}}(U)$ (red points) and $T_{\text{n}}(U)$ (blue points): as $T$ is lowered, the quasiparticle gap continuously develops along the FS, first appearing at the antinodal point $\mathbf{k}_{\text{an}}=(\pi,0)$ at $T_{\text{an}}(U)$ and proliferating to the nodal point $\mathbf{k}_{\text{n}}=(\pi/2,\pi/2)$ at $T_{\text{n}}(U)$. Thus, the crossover involves a regime with a partially gapped FS and damped gapless quasiparticles elsewhere on the FS, the so-called pseudogap driven by extending AFM fluctuations Tremblay et al. (2006); Jin et al. (2011); Sénéchal and Tremblay (2004); Macridin et al. (2006); Jin et al. (2011); Gunnarsson et al. (2015); Wu et al. (2017, 2018), similar to that found at the MIT for the model (1) in eight-site cluster-DMFT Werner et al. (2009). It extends over an appreciable range of parameters at larger $U$ (or $T$). When quasiparticle properties could be meaningfully defined ($T\lesssim 0.25t$), we find that already at $U\gtrsim 4t$ the FL is lost and the self-energy reveals a charge gap. This value is significantly smaller than the critical $U_{c}\sim 5-6t$ found for the MIT in small-cluster DMFT results Park et al. (2008); Werner et al. (2009); Gull et al. (2013). This leaves no room for the MIT in the Hubbard model (1) without additional frustration of antiferromagnetic correlations—the FL-quasi-AFM crossover destroys the FL before it can undergo a first-order transition everywhere where the FL can be defined. $U=4t$ is the upper bound on the interaction strength beyond which the low-$T$ behavior is not qualitatively different from that of the Heisenberg model. As the crossover is expected to become increasingly mean-fieldlike at smaller $U<2.5t$, driven by magnetic correlations with large $\xi$, it is rather instructive that in this regime the crossover temperatures $T_{\text{an}}(U)$ and $T_{\text{n}}(U)$ approximately coincide and both are captured by the mean-field Néel temperature ansatz $a\exp(-b/\sqrt{U})$ with empirical parameters that agree with estimates found in Ref. Šimkovic IV et al. (2017).

We verify our results in Fig. 1 against large-scale dynamical cluster approximation (DCA) calculations at higher temperatures. DCA is a nonperturbative momentum space variant of cluster DMFT with which we utilize an auxiliary-field cluster impurity solver Hettler et al. (2000); Jarrell et al. (2001); Maier et al. (2005). Results for cluster sizes up to $144$ sites reveal very slow convergence of the self-energy with cluster size. We further illustrate the significance of finite-size errors by comparing $\Sigma$DDMC results in the TDL to $\Sigma$DDMC calculations on finite lattices.

In the FL theory, the quasiparticle residue at the chemical potential is a positive number less than unity. It is defined through $\lim\limits_{\omega\to 0}\left[\partial\rm{Re}\Sigma_{\mathbf{k}}(\omega)/\partial\omega\right]$ (with $\Sigma_{\mathbf{k}}(\omega)$ the self-energy at the momentum $\mathbf{k}$ and real frequency $\omega$). On the Matsubara (imaginary-frequency) axis, in the low-temperature limit, this is equivalent through a Wick rotation to $\lim\limits_{i\omega_{n}\to 0}\left[\partial\rm{Im}\Sigma_{\mathbf{k}}(i\omega_{n})/\partial\omega_{n}\right]$. In contrast, $\Sigma_{\mathbf{k}}(\omega)$ in an insulator exhibits a pole at $\omega=0$. At sufficiently low temperature, when discrete values of $\omega_{n}=2\pi T(n+1/2)$ remain closely spaced to perform the limit $\omega_{n}\to 0$, this qualitative difference provides a metric to define the state: if $\rm{Im}\Sigma_{\mathbf{k}}(i\omega_{0})>\rm{Im}\Sigma_{\mathbf{k}}(i\omega_{1})$ for all $\mathbf{k}$ on the FS the state is metallic (FL-like, referred to as FL below); when the reverse is true, $\rm{Im}\Sigma_{\mathbf{k}}(i\omega_{0})<\rm{Im}\Sigma_{\mathbf{k}}(i\omega_{1})$ on the whole FS, the system exhibits an insulating behavior by opening a quasigap at finite temperature. Throughout we use a shorthand notation, $\Delta\Sigma_{\mathbf{k}}=\rm{Im}\Sigma_{\mathbf{k}}(i\omega_{0})-\rm{Im}\Sigma_{\mathbf{k}}(i\omega_{1})$, positive (negative) values of which imply FL (insulator) states. This characterization loses its meaning in a thermal state when the first frequency $\omega_{1}=3\pi T$ is of the order of the Fermi energy.

$\Sigma$DDMC performs a numerically exact evaluation of the coefficients $a_{n}$ in the Taylor-series expansion of $\Sigma$ in powers of $U$—which at half filling takes on the form

[TABLE]

—up to a maximal order $n_{*}$, the truncation of the series being the only approximation. In general, reconstructing $\Sigma_{\mathbf{k}\sigma}(i\omega)$ from its series is a problem by itself Šimkovic and Kozik (2019). In the regime of interest, however, the crossover lines are always within the series convergence radius, and, in principle, the result can be obtained by taking the sum up to a sufficiently high $n_{*}$ to ensure that the truncation error is negligible compared to the statistical error. We are able to compute $\{a_{n}\}$ with statistical errors $\lesssim 10\%$ up to order $n_{*}=12$ for temperatures $T\geq 0.1$ and up to order $n_{*}=10$ for $0.025\leq T\leq 0.1$. We further extrapolate the series by using the standard Dlog-Padé-type approximants Baker Jr (1961); Hunter and Baker Jr (1979), and verify that the systematic error of the extrapolation procedure is small compared to the statistical error (see the Supplemental Material Sup ).

Figure 2 illustrates how the crossover diagram, Fig. 1, was obtained. It shows the variation with $U$ of $\rm{Im}\Sigma(i\omega_{0})$ and $\rm{Im}\Sigma(i\omega_{1})$ at the antinodal and nodal points ($\rm{Re}\Sigma=0$ at the FS). At $U=2t$ the values at $\omega_{0}$ are higher than at $\omega_{1}$, which is typical for a FL. As $U$ is increased, $\Delta\Sigma_{\mathrm{k}}$ shows a trend toward nFL behavior by first changing its sign at $\mathbf{k}_{\text{an}}$ (we take it as the onset of the nFL behavior). Following our measure, we mark the region of $U$ and $T$ where $\Delta\Sigma_{\mathbf{k}}>0$ for all momenta on the FS as FL (red shading). Similarly, the insulating region (blue shading) corresponds to $\Delta\Sigma_{\mathbf{k}}<0$ for all $\mathbf{k}$ on the FS. The nFL pseudogap regime (gray) falls in between the two: it has $\Delta\Sigma_{\mathbf{k}}>0$ at some momenta on the FS and $\Delta\Sigma_{\mathbf{k}}<0$ at others. Correspondingly, the nFL state is bounded by the temperature scales $T_{\text{an}}(U)$ and $T_{\text{n}}(U)$ where $\Delta\Sigma_{\mathbf{k}}=0$ at momenta $\mathbf{k}_{\text{an}}$ and $\mathbf{k}_{\text{n}}$, respectively. As the gap proliferates along the FS between $T_{\text{an}}(U)$ and $T_{\text{n}}(U)$, we expect that the heat capacity $C_{V}$ to $T$ ratio decreases before growing due to the AFM quasiorder near the insulator boundary and eventually reaching the asymptotic $C_{V}/T\propto T$ law due to spin waves deep in the insulating regime.

At small $U$, both crossover temperatures scale exponentially according to the BCS solution for the mean-field AFM transition, $T_{\text{an}}=a\exp(-b/\sqrt{U})$ and $T_{\text{n}}=a^{\prime}\exp(-b^{\prime}/\sqrt{U})$, with the fit parameters $a=6.99$, $b=6.51$, $a^{\prime}=4.7$, $b^{\prime}=6.08$ (see Fig. 1), suggesting the crossover is being driven by extended AFM correlations.

At high temperature $T\gtrsim 0.25t$, the data points corresponding to $\Delta\Sigma_{\mathbf{k}_{\text{an},\text{n}}}=0$ lose their meaning as boundaries between FL, nFL, and insulator regimes (white region in Fig. 1). Given that $\rm{Im}\Sigma(i\omega_{n})$ is a negative-valued function approaching zero at large frequencies, it is clear that by increasing $T$ in the FL regime $\Delta\Sigma$ will change the sign, marking a crossover to the thermal gas, not the nFL, state. Nonetheless, we plot the results for $\Delta\Sigma_{\mathbf{k}_{\text{an},\text{n}}}=0$ at high temperature as TDL benchmarks for other numerical techniques.

Providing the controlled extrapolation of finite-system numerical results to the TDL has been long recognized as important. Thus far it has been accomplished, in particular, by large-system-size studies of the 3D Hubbard model near the AFM Néel transition Kent et al. (2005); Fuchs et al. (2011); Kozik et al. (2013). In 2D, the extremely slow finite-size scaling due to the exponential growth of the AFM correlation length near the crossover makes Fig. 1 challenging to reproduce by finite-size methods, even if they do not suffer from the fermionic sign problem at half filling (cf. DQMC results for this problem in Ref. Schäfer et al. (2015)). In the Supplemental Material Sup , we further illustrate the difficulty of obtaining $T_{\text{an}}(U)$, $T_{\text{n}}(U)$ through extrapolation to the TDL by the example of $\Sigma$DDMC calculations on finite-size lattices of dimensions $L\times L$, where any value of $L$ is accessible at the same computational cost.

To verify our results by an independent method we resort to the DCA, which produces unbiased results after extrapolation to the TDL. At low temperatures the extrapolation is extremely challenging. Figure 3 shows DCA results for $\rm{Im}\Sigma_{\mathbf{k}}(i\omega_{n})$ at $U=3t$, $T=0.12t$, and various cluster sizes $16\leq N_{c}\leq 144$. At $\mathbf{k}=\mathbf{k}_{\text{n}}$ (left panel), the results show the FL behavior for all accessible system sizes, while the $\Sigma$DDMC results are insulatorlike. At $\mathbf{k}=\mathbf{k}_{\text{an}}$ (right panel), the character changes from FL to insulating as a function of $N_{c}$ when cluster sizes exceed $100$, reaching qualitative agreement with the $\Sigma$DDMC data. Bigger clusters are required for a quantitatively accurate extrapolation. Note that if the data for only $N_{c}\leq 72$ were available, one would be led to conclude that the state at $U=3t$, $T=0.12t$ is a FL. [Similarly, in the DQMC study of Ref. Schäfer et al. (2015), small-system-size data in what is actually a FL regime show insulatorlike behavior, resulting in a large error bar of the extrapolation to the TDL.] It is not surprising then that past work limited to substantially smaller cluster sizes Park et al. (2008); Werner et al. (2009); Gull et al. (2013) observed the MIT at $U_{c}>5t$, while the TDL system is, in fact, already insulating at smaller $U$. In general, finite-size effects are less severe at higher $T$, and a controlled extrapolation of DCA data to the TDL is feasible: extrapolated DCA results for the crossover temperatures are in quantitative agreement with $\Sigma$DDMC for $T\gtrsim 0.2t$ (see the Supplemental Material Sup ).

Arbitrary momentum resolution of $\Sigma$DDMC allows us to observe the structure correlations in the momentum space. Figure 4 shows the Brillouin zone map of $\rm{Im}\Sigma_{\mathbf{k}}(\omega_{0})$ and $\Delta\Sigma_{\mathbf{k}}$ and relates them to the measurable with ultracold atoms (quasi)momentum distribution $n(\mathbf{k})$ at the experimentally relevant $T=0.2t$ and three values of $U$ that are seen in Fig. 1 to be: in the metallic ($U=2t$), nFL ($U=3.6t$), and insulator ($U=4t$) regimes, respectively. As $U$ is increased, $|\rm{Im}\Sigma_{\mathbf{k}}(\omega_{0})|$ exhibits relatively small change away from the FS and grows substantially along the FS, with the largest value at $\mathbf{k}_{\text{an}}$ and the smallest within the FS at $\mathbf{k}_{\text{n}}$. At the same time, $\Delta\Sigma_{\mathbf{k}}$ remains positive for all $\mathbf{k}$ in the metallic regime, but islands of negative values appear around $\mathbf{k}_{\text{an}}$ at $U=3.6t$ [below $T_{\text{an}}(U)$]. At $U=4t$ [below $T_{\text{n}}(U)$], $\Delta\Sigma_{\mathbf{k}}$ is negative along the whole FS, with the lowest (most insulatorlike) values around $\mathbf{k}_{\text{an}}$. Interestingly, away from the FS, the change in $\Delta\Sigma_{\mathbf{k}}$ is opposite: it grows, becoming more FL-like, upon increasing $U$. The corresponding $n(\mathbf{k})$ shows gradual smearing of the step at the FS as $U$ increases and the system crosses from the FL to insulating regime. This smearing is more pronounced around $\mathbf{k}_{\text{an}}$: the diamond-shaped region of occupied states in the metallic regime shrinks and evolves toward a circle across the crossover, providing an observable signature of developing correlations.

In conclusion, we have revealed the scenario of the metal-to-quasiantiferromagnetic-insulator crossover in the 2D Hubbard model (1) qualitatively different from the MIT previously suggested for this system Park et al. (2008); Werner et al. (2009); Gull et al. (2013) and in qualitative agreement with the recent D$\Gamma$A Schäfer et al. (2015); Rohringer and Toschi (2016) and DQMC Schäfer et al. (2015) results. The crossover could not be captured by the small-cluster DMFT restricted to the paramagnetic solution, predicting the MIT instead. The insulating regime sets in at all values of $U$ due to extended AFM correlations that transform the system into the quasi-AFM after an intermediate nFL pseudogap regime. The quantitative shape of the crossover is different from that reported in Refs. Schäfer et al. (2015); Rohringer and Toschi (2016): it is described by the mean-field AFM transition at small $U$ and features a nFL regime that transforms to insulating behavior below $U\approx 4t$. All our results are obtained with controlled accuracy and offer guidance for precision experiments with ultracold atoms in optical lattices, as well as unbiased numerical techniques, in the ongoing effort to describe the phase diagram of the Hubbard model. In particular, the most nontrivial correlated regime is realized for $2.5t<U<4t$ and temperatures $T\lesssim 0.25t$. At weaker coupling, $U<2.5t$, the low-temperature behavior is governed by the mean-field BCS-type physics, while at $U>4t$ the low-temperature state is expected to qualitatively resemble that of the Heisenberg model Kim et al. (2019). By continuity of the key mechanism, the long-range AFM correlations (quantified, e.g., in Ref. Šimkovic IV et al. (2017)), this qualitative crossover picture is valid in a range of nonzero next-to-nearest-neighbor hopping $t^{\prime}$ and doping $\delta$. The question of whether the conventional MIT scenario is realized at certain (large-enough) $t^{\prime}$ requires further systematic investigation.