Dominant charge-density-wave correlations in the Holstein model on the half-filled square lattice

TL;DR

This study uses quantum Monte Carlo simulations to investigate charge-density-wave correlations in the Holstein model on a half-filled square lattice, finding strong CDW tendencies across various phonon frequencies.

Contribution

It provides the first comprehensive analysis across all phonon frequencies, demonstrating persistent CDW correlations in the Holstein model at zero temperature.

Findings

CDW correlations dominate over pairing correlations at low temperatures

Results are consistent with a CDW ground state

Analysis spans adiabatic to antiadiabatic limits

Abstract

We use an unbiased, continuous-time quantum Monte Carlo method to address the possibility of a zero-temperature phase without charge-density-wave (CDW) order in the Holstein and, by extension, the Holstein-Hubbard model on the half-filled square lattice. In particular, we present results spanning the whole range of phonon frequencies, allowing us to use the well understood adiabatic and antiadiabatic limits as reference points. For all parameters considered, our data suggest that CDW correlations are stronger than pairing correlations even at very low temperatures. These findings are compatible with a CDW ground state that is also suggested by theoretical arguments.

Equations

$$\displaystyle\hat{H}=-t\sum_{\langle i,j\rangle\sigma}\hat{c}^{\dagger}_{i\sigma}\hat{c}^{\phantom{{\dagger}}}_{j\sigma}+\sum_{i}\left[\mbox{$\frac{1}{2M}$}\hat{P}^{2}_{i}+\mbox{$\frac{K}{2}$}\hat{Q}_{i}^{2}\right]-g\sum_{i}\hat{Q}_{i}\hat{\rho}_{i}\,.$$

$${\chi }_{\text{c}}(\mathbf{q})$$

$${\chi }_{\text{p}}(\mathbf{q})$$

$$R_{\text{CDW}}^{\chi }$$

$${\chi }_{\text{CDW}}=L^{2-\eta }f[L^z/(T-T_c^{\text{CDW}})],$$

$${\chi }_{\text{SC}}=L^{2-\eta }$$

$${\chi }_{\text{c}}$$

$${\chi }_{\text{s}}$$

Full text

I Introduction

Charge-density-wave (CDW) and superconducting (SC) phases are ubiquitous in quasi-two-dimensional (quasi-2D) materials and often arise from electron-phonon coupling. Holstein's molecular-crystal model [1] of electrons coupled to quantum phonons has played a central role for the investigation of such phenomena. However, even after decades of research, fundamental questions are still unanswered. Quantum Monte Carlo (QMC) approaches have played a key role in the study of this problem. Despite important recent methodological advances [2; 3; 4; 5; 6; 7], simulations of electron-phonon problems remain significantly more challenging than, for example, those of purely fermionic Hubbard models.

The Holstein-Hubbard model captures the interplay of electron-phonon coupling ($\sim\lambda$) and electron-electron repulsion ($\sim U$). For $U=0$, it reduces to the Holstein model simulated in the following. For the much studied case of a half-filled square lattice, earlier work [8; 9; 10; 11; 12; 13; 14] agreed on either long-range CDW or AFM order at $T=0$ depending on $\lambda$ and $U$, consistent with theoretically expected instabilities of the Fermi liquid. In contrast, relying on variational QMC simulations, two recent papers [14; 15] reported the existence of an intermediate metallic phase with neither CDW nor AFM order, see Fig. 1. Instead, it was characterized as either SC or paramagnetic [14; 15]. The prospect of a metallic ground state has to be distinguished from metallic behavior emerging at finite temperatures [13] simply via the thermal destruction of AFM order [2]. The predicted existence of the metallic state even at $U=0$ (see Fig. 1), i.e., the Holstein model, appears to rule out competing interactions as its origin.

Apart from the challenges due to small gaps and order parameters at weak coupling, the methods used in Refs. [14; 15] are variational and involve an ansatz for the ground-state wave function that may bias the results. A similar controversy regarding metallic behavior in the 1D Holstein-Hubbard model was recently resolved. For the latter, approximate strong-coupling results in combination with unfounded conclusions from numerical simulations [16] as well as insufficiently accurate renormalization-group (RG) approaches were contradicted by unbiased numerical simulations and functional RG calculations. For the 1D case, a disordered phase has been firmly established [17], although claims of dominant pairing correlations [18] in this regime were refuted [19]. Whereas the Fermi liquid is expected to be unstable at $T=0$ in the particular 2D setting considered, Refs. [14; 15] also suggest that the non-CDW region could have SC order. SC correlations have been found to be enhanced in Holstein models with next-nearest-neighbor hopping [20], dispersive phonons [21], frustration [22], or finite doping [23]. An extended semimetallic phase is supported by theory and numerics in the Holstein model on the honeycomb lattice [24; 25].

Here, to provide further insight into this problem in the limit $U=0$, we exploit the properties of the continuous-time interaction-expansion (CT-INT) QMC method [26]. Compared to other approaches, it can in principle access rather low temperatures in the weak-coupling regime. Although simulations are partially restricted by a sign problem, we obtain evidence for long-range CDW order at very low temperatures for parameters where a non-CDW phase was predicted [14; 15].

The paper is organized as follows. In Sec. II, we define the model and summarize previous work and theoretical arguments. Section III provides the necessary details about the CT-INT simulations. Our results are discussed in Sec. IV, followed by our conclusions in Sec. V.

II Model

References [14; 15] presented phase diagrams for the Holstein-Hubbard model on the half-filled square lattice. Selected results are reproduced in Fig. 1. Because the purported non-CDW region is most extended for a vanishing Hubbard repulsion ($U=0$), we focus on the simpler Holstein Hamiltonian [1]

[TABLE]

Here, $\hat{c}^{\dagger}_{i\sigma}$ creates an electron with spin $\sigma$ at lattice site $i$ and the first term describes nearest-neighbor hopping with amplitude $t$. Lattice vibrations are modeled in terms of independent harmonic oscillators with frequency $\omega_{0}=\sqrt{K/M}$, displacements $\hat{Q}_{i}$, and momenta $\hat{P}_{i}$. The electron-phonon interaction takes the form of a density-displacement coupling to local fluctuations of the electron number; we have $\hat{\rho}_{i}=\hat{n}_{i}-1$ with $\hat{n}_{i}=\sum_{\sigma}\hat{n}_{i\sigma}$ and $\hat{n}_{i\sigma}=\hat{c}^{\dagger}_{i\sigma}\hat{c}^{\phantom{{\dagger}}}_{i\sigma}$. All simulations were done on $L\times L$ square lattices with periodic boundary conditions and for a half-filled band ($\langle\hat{n}_{i}\rangle=1$, chemical potential $\mu=0$). We define a dimensionless coupling constant $\lambda=g^{2}/(WK)$ ($W=8t$ is the free bandwidth), set $\hbar$, $k_{\text{B}}$, and the lattice constant to one and use $t$ as the energy unit.

Hamiltonian (1) has been the subject of numerous QMC investigations [23; 27; 28; 29; 30; 20; 31; 32; 33]. With the exception of Refs. [14; 15], a CDW ground state was assumed to exist for any $\lambda>0$, as suggested by several theoretical arguments. First, for classical phonons, corresponding to $\omega_{0}=0$, mean-field theory is exact at $T=0$ and reveals a gap and CDW order for any $\lambda>0$ [29]. The origin of this weak-coupling instability is the combination of perfect nesting on a half-filled square lattice with nearest-neighbor hopping and a zero-energy Van Hove singularity in the density of states [34; 27]. These features give rise to a noninteracting charge susceptibility [defined in Eq. (2) below] that diverges as $\chi^{(0)}_{\text{CDW}}\sim\ln^{2}\beta t$ [34; 27], where $\beta=1/T$. Both, at the mean-field level and in numerical simulations, such a divergence produces CDW order for any $U<0$ in the attractive Hubbard model [35]. The latter is an exact limit of the Holstein model for $\omega_{0}\to\infty$ ($M\to 0$), with $U=\lambda W$. Hence, long-range CDW order at $T=0$ is established for the Holstein model both for $\omega_{0}=0$ and $\omega_{0}=\infty$. In contrast, the s-wave pairing susceptibility [Eq. (3)] has a weaker divergence, $\chi^{(0)}_{\text{SC}}\sim\ln\beta t$, because nesting plays no role. This is consistent with the observation that for $\omega_{0}<\infty$ SC correlations are weaker than CDW correlations at half-filling [29; 23] but not with an SC phase [14; 15]. However, earlier work did not consider the weak-coupling regime and unbiased, high-precision finite-size scaling analyses only appeared recently [2; 5; 4; 7].

In Refs. [14; 15], the challenging problem of determining the ground-state phase diagram was approached using zero-temperature variational QMC methods. In contrast, most other work (for exceptions see Refs. [8; 22]) infer ground-state properties from simulations at low but finite temperatures. Whereas the AFM phase of the Holstein-Hubbard model shown in Fig. 1 exists only at $T=0$, long-range CDW order is associated with an Ising order parameter and persists up to a critical temperature $T^{\text{CDW}}_{c}$ [29; 2]. Similarly, the U(1) SC order parameter also permits a nonzero transition temperature $T^{\text{SC}}_{c}$. In both cases, given an ordered ground state, we therefore expect a finite-temperature phase transition. An important exception is the limit $\omega_{0}=\infty$, corresponding to the attractive Hubbard model. The latter has an enhanced symmetry that combines the CDW and SC order parameters into an SU(2) vector [35]. According to the Mermin-Wagner theorem [36], long-range order is therefore confined to $T=0$.

Let us address the purported intermediate phase at weak coupling and $\omega_{0}>0$ reported in Refs. [14; 15] in the light of these arguments. The overall size of this phase increases with increasing $\omega_{0}/t$ in Refs. [14; 15], similar to the case of the 1D Holstein-Hubbard model [17]. In 1D, quantum lattice fluctuations promote the proliferation of domain walls in the Ising CDW order parameter. There, the ground state is metallic up to $\lambda=\lambda_{c}$ with $\lambda_{c}\to\infty$ for $\omega_{0}\to\infty$ (attractive Hubbard model) [17]. In contrast, the 2D Holstein model is CDW-ordered in the antiadiabatic limit $\omega_{0}\to\infty$. An explicit comparison of data for $\omega_{0}=\infty$ and $\omega_{0}<\infty$ will be made in Sec. IV. No theoretical arguments were given in Refs. [14; 15] against the weak-coupling instability expected from the divergence of $\chi^{(0)}_{\text{CDW}}$. Interestingly, although $\chi_{\text{CDW}}^{(0)}$ and $\chi_{\text{AFM}}^{(0)}$ diverge in the same way, the methods of Refs. [14; 15] successfully detect the weak-coupling AFM instability at $\lambda=0$ but not the CDW instability at $U=0$ (see Fig. 1). Another apparent inconsistency is that the non-CDW region at $U=0$ is significantly larger for $\omega_{0}/t=8$ than for $\omega_{0}/t=1$ in Ref. [14], whereas it remains virtually unchanged between $\omega_{0}/t=1$ and $\omega_{0}/t=15$ in Ref. [15]. For the value $\omega_{0}/t=1$ analyzed in both works and shown in Fig. 1, Ref. [14] predicts a non-CDW ground state up to $\lambda\approx 0.11$ at $U=0$, whereas Ref. [15] reports a critical value of $\lambda\approx 0.125$ (using our definition of $\lambda$).

III Method

The application of the CT-INT method [26] to electron-phonon models goes back to the work by Assaad and Lang [37]. For investigations of 2D Holstein and Holstein-Hubbard models, see Refs. [2; 24]. Its general, action-based formulation makes it suitable for retarded fermion-fermion interactions that arise naturally from electron-phonon problems after integrating out the phonons in the path-integral representation of the partition function [37]. The weak-coupling expansion can be shown to converge for fermionic systems in a finite space-time volume [26], so that the method is exact apart from statistical errors. General reviews have been given in Refs. [38; 39].

The numerical effort scales cubically with the average expansion order $n$, where $n\approx{O}(\beta\lambda L^{2})$ for the Holstein model. While other methods formally scale linearly in $\beta$ [40], CT-INT is typically less limited by autocorrelation times. For the present work, its use is motivated by a significant speedup at small $\lambda$ that permits us to study reasonably large system sizes up to $L\leq 12$ at inverse temperatures $\beta t\leq 96$. For intermediate phonon frequencies, the method is ultimately limited by a sign problem that arises from the absence of an exact symmetry between the spin-$\uparrow$ and spin-$\downarrow$ sectors [2]. The results for $\omega_{0}=\infty$ were obtained by directly simulating the attractive Hubbard model with the CT-INT method. We used 1000 single-vertex updates and 8 Ising spin flips per sweep for all simulations. Although our method is entirely unbiased, as opposed to the algorithms of Refs. [14; 15], limitations arise regarding model parameters, temperatures, and system sizes.

IV Results

To detect CDW and/or s-wave SC order, we carried out a finite-size scaling analysis based on the charge and pairing susceptibilities (with $\hat{\Delta}_{i}=\hat{c}_{i\uparrow}\hat{c}_{i\downarrow}$)

[TABLE]

We define $\chi_{\text{CDW}}\equiv\chi_{\text{c}}(\bm{Q}_{\text{CDW}})$ with $\bm{Q}_{\text{CDW}}=(\pi,\pi)$ and $\chi_{\text{SC}}\equiv\chi_{\text{p}}(\bm{Q}_{\text{SC}})$ with $\bm{Q}_{\text{SC}}=(0,0)$. The factor $2$ in the definition of $\chi_{\text{p}}$ ensures $\chi_{\text{CDW}}\equiv\chi_{\text{SC}}$ for the attractive Hubbard model ($\omega_{0}\to\infty$) and at $\lambda=0$.

Long-range CDW order can be detected from by the renormalization-group invariant correlation ratio

[TABLE]

At a fixed $\lambda$, $R^{\chi}_{\text{CDW}}$ depends only on $L^{z}/(T-T^{\text{CDW}}_{c})$, so that data for different $L$ are expected to intersect (up to corrections to scaling) at the transition temperature $T^{\text{CDW}}_{c}$. By definition, $R^{\chi}_{\text{CDW}}\to 0$ as $L\to\infty$ in the absence of long-range CDW order, whereas $R^{\chi}_{\text{CDW}}\to 1$ as $L\to\infty$ if $\chi_{\text{CDW}}$ diverges with $L$. The ratio $R^{\chi}_{\text{CDW}}$ can be expected to have smaller scaling corrections than $\chi_{\text{CDW}}$ itself [41]. The use of susceptibilities rather than static structure factors suppresses background contributions to critical fluctuations. We will also consider the finite-size scaling of the susceptibility itself, which is described near the Ising critical point by the scaling form

[TABLE]

with $\eta=0.25$ known from the exact solution of the 2D classical Ising model.

A potential transition to an SC phase should be in the Berezinskii-Kosterlitz-Thouless universality class with power-law correlations below the critical temperature $T^{\text{SC}}_{c}$. In the absence of long-range order and hence a divergence of $\chi_{\text{SC}}$, we exploit the finite-size scaling form exactly at the critical temperature

[TABLE]

with $\eta=0.25$ [42]. Equation (6) again implies a crossing point of results for different $L$ at $T=T^{\text{SC}}_{c}$, as recently observed for the half-filled Holstein model on the frustrated triangular lattice [22].

To detect gaps for long-wavelength charge and spin fluctuations, we also consider the static (uniform) charge and spin susceptibilities

[TABLE]

Here, $\hat{S}^{x}_{i}=\hat{c}^{\dagger}_{i\uparrow}\hat{c}^{\phantom{{\dagger}}}_{i\downarrow}+\hat{c}^{\dagger}_{i\downarrow}\hat{c}^{\phantom{{\dagger}}}_{i\uparrow}$, giving a maximum magnetization per site $\langle M\rangle/L^{2}=1$.

Based on the arguments in Sec. II, we expect CDW order rather than SC order at half-filling. Mean-field theory predicts a CDW transition temperature $T^{\text{CDW}}_{c}\sim e^{-1/\sqrt{\lambda}}$ that appears consistent with simulations for $\omega_{0}=0$ [2] and renders the weak-coupling regime challenging. Moreover, $T^{\text{CDW}}_{c}$ decreases with increasing $\omega_{0}$, and vanishes for $\omega_{0}=\infty$ [35]. Before discussing the case of $\omega_{0}/t=1$ depicted in Fig. 1, we consider $\omega_{0}/t=0.1$ (close to the mean-field limit) and $\omega_{0}/t=\infty$ (the attractive Hubbard model) as useful reference points. To address the findings of Refs. [14; 15] for $\omega_{0}/t=1$, we consider the couplings $\lambda=0.075$ and $\lambda=0.025$, both inside the purported non-CDW phase in Fig. 1.

Results for $\omega_{0}/t=0.1$ and $\lambda=0.25$ are shown in Fig. 2. The CDW transition for these parameters was previously investigated using equal-time correlation functions [2], so that we can benchmark our diagnostics. The data reveal a strong increase of the CDW susceptibility at low temperatures [Fig. 2(a)], whereas the SC susceptibility is virtually independent of $L$ [Fig. 2(b)]. The rescaled CDW susceptibility in Fig. 2(c) exhibits a clean crossing point at $T^{\text{CDW}}_{c}/t\approx 0.2$, in agreement with previous findings [2]. This crossing is consistent with that observed in the results for the correlation ratio in Fig. 2(e), a close-up of which is shown in Fig. 2(f). In contrast, the rescaled SC susceptibility [Fig. 2(d)] is strongly suppressed for $T\lesssim T_{c}^{\text{CDW}}$. Finally, the uniform charge and spin susceptibilities in Figs. 2(g) and (h), respectively, reveal a gap in both sectors at sufficiently low temperatures, as expected for a CDW insulator.

In the opposite, antiadiabatic regime $\omega_{0}=\infty$, we can rely on previous results for the ground state of the attractive Hubbard model [35; 43] to interpret our finite-temperature data. To gain insight into the behavior expected in the weak-coupling regime, we consider $\lambda=0.075$ and focus on low temperatures $T/t\leq 0.1$. The CT-INT data in Fig. 3 exhibit significant differences compared to Fig. 2. The CDW and SC susceptibilities in Fig. 3(a) and Fig. 3(b) are identical due to the SO(4) symmetry of the Hubbard model. Whereas a crossing point is not visible for the rescaled susceptibilities in Figs. 3(c) and (d) at the temperatures considered, the CDW correlation ratio [Figs. 3(e),(f)] again approaches 1 for $T\to 0$. Such behavior is consistent with long-range CDW order at $T=0$. Finally, the charge susceptibility in Fig. 3(g) is consistent with metallic behavior (due to the coexistence of CDW and SC order), whereas Fig. 3(h) reveals the expected spin gap.

Having established the physics, but also the limitations of our simulations, in the undisputed adiabatic and antiadiabatic limits, we turn to parameters that are directly relevant for Fig. 1, specifically $\omega_{0}/t=1$ and $\lambda=0.075$, where Refs. [14; 15] predict a paramagnetic or SC state. The corresponding results are shown in Fig. 4. Comparing the CDW and SC susceptibilities in Figs. 4(a) and (b) reveals that CDW correlations are significantly stronger than SC correlations at a given $T$. Whereas a crossing point in the rescaled CDW susceptibility [Fig. 4(c)] at temperatures below the accessible range is plausible, we are unable to reach temperatures comparable to $T_{c}^{\text{CDW}}$. On the other hand, a critical point signaling SC order is not expected based on the results of Fig. 4(d), especially upon comparison with Fig. 3(d) for $\omega_{0}/t=\infty$. The latter case has stronger pairing correlations than observed for $\omega_{0}/t=1$ even though $T^{\text{SC}}_{c}=0$. The CDW correlation ratio [Figs. 4(e),(f)] is also consistent with CDW order at $T=0$. The uniform susceptibilities in Figs. 4(g) and (h) are not entirely conclusive but consistent with a gap for charge and spin excitations at $T=0$.

We also simulated a weaker coupling $\lambda=0.025$, deep inside the predicted intermediate phase in Fig. 1. Of course, any type of order will be extremely delicate to detect at such weak interactions on finite systems. Moreover, CDW and SC correlations are necessarily degenerate at $\lambda=0$ (free fermions). Nevertheless, Fig. 5 does indicate somewhat stronger CDW than SC correlations, which again seems to contradict the claims of Refs. [14; 15]. At the same time, the expected spin gap is only visible in Fig. 5(h) at the lowest temperatures, whereas the expected charge gap is beyond the accessible temperature range in Fig. 5(g).

The dependence of the CDW and SC susceptibilities on the coupling strength $\lambda$ at $T=1/24$ can more clearly be seen in Fig. 6. Starting from identical values at $\lambda=0$, $\chi_{\text{CDW}}$ increases significantly with $\lambda$, whereas $\chi_{\text{SC}}$ flattens after a weak initial increase.

Finally, Fig. 7 compares the temperature-dependent CDW and SC susceptibilities at different phonon frequencies. The CDW susceptibility in Figs. 7(a),(b) evolves continuously, with values for intermediate $\omega_{0}$ falling between those for $\omega_{0}/t=0.1$ and $\omega_{0}/t=\infty$. For the SC susceptibility, Figs. 7(c),(d), the data suggest the possibility of non-monotonic behavior: $\chi_{\text{SC}}$ for $\omega_{0}/t=1$ in Fig. 7(d) is equal to that for $\omega_{0}/t=\infty$ at intermediate temperatures yet still smaller than $\chi_{\text{CDW}}$.

V Conclusions

Although limitations regarding lattice size and temperature preclude definitive conclusions regarding the ground state, we believe that our unbiased results point rather strongly toward long-range CDW order in the half-filled Holstein model on the square lattice. By extension, it seems reasonable to expect only CDW and AFM ground states in the Holstein-Hubbard model.

Our main arguments are as follows.

(i) For the parameters considered, including those where Refs. [14; 15] predict no CDW order, we find that CDW correlations are stronger than SC correlations, consistent with long-range CDW order at $T=0$.

(ii) CDW (SC) correlations are stronger (weaker) than for the attractive Hubbard model with the same effective interaction $U=\lambda W$. The latter corresponds to the Holstein model in the antiadiabatic limit $\omega_{0}\to\infty$. Because the Hubbard model is known to have long-range CDW order at $T=0$, this suggests long-range CDW order also for the Holstein model with $\omega_{0}<\infty$. Weaker SC correlations do not rule out SC order at $T=0$. However, the coexistence of CDW and SC order in the attractive Hubbard model is linked to an enhanced symmetry that is absent in the Holstein case for $\omega_{0}<\infty$ [16]. Even if SC order exists at $T=0$, the stronger CDW order conflicts with the claims of Refs. [14; 15].

(iii) Since we infer the nature of the ground state from simulations at $T>0$, there is in principle a possibility of a non-monotonic temperature dependence, with a phase transition to an SC phase at even lower temperatures. However, we do not observe any signatures or precursor effects of this scenario, such as a decrease of the CDW susceptibility at low temperatures.

(iv) Our results are consistent with the theoretical arguments for a weak-coupling CDW instability due to nesting and a Van Hove singularity, which should apply to the weak-coupling regime where a non-CDW region was reported in Refs. [14; 15].

It is beyond the scope of this work to determine the origin of the different findings in Refs. [14; 15]. However, the necessity of choosing a variational wave function seems the most likely source for different physics in trying to distinguish a paired Fermi liquid from either a pair crystal (CDW state) or a pair condensate (SC state). In the weak-coupling regime that is of interest here, different states are expected to be close in energy. Moreover, the deviations between the critical values estimated with the help of the same QMC method in Refs. [14; 15] and visible in Fig. 1, also with respect to the strong-coupling phase boundary $U=\lambda W$, suggest uncertainties that significantly exceed the reported error bars.

We expect the $T=0$ phase diagram in the ($\lambda$,$U$) plane to contain a single line of critical points that emanates from the point $\lambda=U=0$ and separates CDW and AFM phases. For further progress on this problem, functional RG calculations with a suitable treatment of the energy and momentum dependence of the interaction appear promising [44] to detect CDW order at weak coupling. The combination of projective QMC simulations with improved updates based on recent ideas [4; 5; 7] could yield $T=0$ results without the variational approximation of Refs. [14; 15]. Finally, the use of pinning fields together with an extrapolation to the thermodynamic limit may also prove useful [45].

Acknowledgements.

We thank F. F. Assaad for helpful discussions and the DFG for support via SFB 1170. We gratefully acknowledge the computing time granted by the John von Neumann Institute for Computing (NIC) and provided on the supercomputer JURECA [46] at the Jülich Supercomputing Centre.