Phase diagram of the two-dimensional Hubbard-Holstein model: enhancement of $s$-wave pairing between charge and magnetic orders

We investigate the role of electron-electron and electron-phonon interactions in strongly correlated systems by performing unbiased quantum Monte Carlo simulations in the square lattice Hubbard-Holstein model at half-filling. We study the competition and interplay between antiferromagnetism (AFM) and charge-density wave (CDW), establishing its very rich phase diagram. In the region between AFM and CDW phases, we have found an enhancement of superconducting pairing correlations, favoring (nonlocal) $s$-wave pairs. Our study sheds light over past inconsistencies in the literature, in particular the emergence of CDW in the pure Holstein model case.

Introduction: The electron-phonon (e-ph) interaction is a central issue in condensed matter, in particular when discussing properties of conventional superconductivity (SC) and charge ordering [1]. While Bardeen, Cooper and Schrieffer used this interaction in their seminal work to explain pairing [2], Peierls took it into account to provide a mechanism, based on Fermi surface nesting (FSN), that leads to charge-density wave (CDW) [3].
Recently, the debate about the role of the e-ph coupling has been intensified due to the occurrence of unconventional (non Peierls-like) CDW phases, and their competition with SC, in some classes of materials, such as transition-metal dichalcogenides [4][5][6][7][8][9]. Even for cuprates, materials known by their strong electronelectron (e-e) interactions, recent findings provided evidence for the occurrence of CDW in the doped region, with competing effects with SC [9][10][11][12][13][14], e.g., on doped LBCO and YBCO [15][16][17][18][19][20][21]. These results have suggested that the phase diagram of high-T c superconductors [22] is far more complex than previously supposed, and have raised issues about the relevance of the e-ph coupling for correlated materials, rather than just e-e interactions.
From a theoretical point of view, a simplified Hamiltonian that captures the interplay between antiferromagnetism (AFM), CDW, and SC is the singleband Hubbard-Holstein model (HHM) [23]. It exhibits Coulomb repulsion between electrons, leading to spin fluctuations; and also electron-ion interactions, which enhance charge/pairing correlations. The emergence of long-range order depends on the competition between these tendencies. This model was vastly studied in onedimensional systems, with well-known phase diagrams presenting spin-density wave, bond-order-wave, CDW, and also metallic or phase separation behavior [24][25][26][27][28][29][30][31][32][33]. A remarkable feature in 1D systems is the occurrence of a quantum phase transition from a metallic Luther-Emery liquid phase to a CDW insulator at a finite critical e-ph coupling, in the limit case of the pure Holstein model (U = 0), despite of the FSN [34,35]. By contrast, the properties of the HHM in two-dimensional systems are not entirely clear, even for simple geometries, such as the square lattice. For instance, the existence of such a metal-CDW quantum critical point (QCP) is matter of controversies for the pure Holstein model in 2D lattices [36][37][38][39][40][41]. The scenario is much less clear in presence of a repulsive Hubbard term (U = 0), in spite of the large effort to characterize the model [23,[39][40][41][42][43][44][45][46][47][48][49][50][51], since no unbiased results are available for quantum AFM/CDW transitions, to the best of our knowledge. In view of these open issues, and as a step towards a better understanding of the role of e-ph interactions in strongly interacting systems, we investigate in this Letter the competition between AFM and CDW in the square lattice HHM at half-filling, as well as its pairing response, using unbiased quantum Monte Carlo (QMC) simulations. We determine precise critical points for the HHM at intermediate interaction strengths, presenting benchmarks for lattices with linear size up to L = 64 (i.e., 4096 sites) in some cases. Our main results are summarized in the ground state phase diagram displayed in Fig. 1. Here we highlight [i] the absence of a finite critical e-ph coupling for the pure Holstein model, i.e., the CDW phase sets in for any λ > 0 (and U = 0); [ii] the existence of a finite AFM critical point on the line U = λ, which is strongly dependent on the phonon frequency; and [iii] an enhancement of nonlocal s-wave pairing in the region between the AFM and CDW phases. These results are also compared with other methodological approaches, such as variational QMC. Methodology: The Hubbard-Holstein Hamiltonian reads where d † iσ (d iσ ) is a creation (annihilation) operator of electrons with spin σ (=↑, ↓) at a given site i on a two-dimensional square lattice under periodic boundary conditions, with i, j denoting nearest-neighbors, and n iσ ≡ d † iσ d iσ being number operators. The first two terms on the right hand side of Eq. (1) correspond to the kinetic energy of electrons, and their chemical potential (µ) term, respectively, while the on-site Coulomb repulsion between electrons is included by the third term. The ions' phonon modes are described in the fourth term, in whichP i andX i are momentum and position operators, respectively, of local quantum harmonic oscillators with frequency ω 0 .
The last term corresponds to local electron-ion interactions, with strength g. Hereafter, we define the mass of the ions (M ) and the lattice constant as unity.
It is also worth to introduce additional parameters, due to the effects of the phonon fields to the electronic interactions. From a second order perturbation theory on the e-ph term [23], one obtains an effective dynamic e-e interaction, U eff (ω) = U − g 2 /ω 2 0 1−(ω/ω0) 2 , with g 2 /ω 2 0 ≡ λ being the energy scale for polaron formation. The appearance of such a retarded attractive interaction, depending on the phonon frequency, leads us to define ω 0 /t as the adiabaticity ratio, and λ/t as the strength of the e-ph interaction. To facilitate the following discussion, we also define U eff ≡ U − λ, which gives us information about the local effective e-e interaction, and is also an important parameter to our methodological approaches. Furthermore, we set the electron density at half-filling, i.e., n iσ = 1/2.
We investigate the properties of Eq. (1) by performing two different unbiased auxiliary-field QMC approaches: the projective ground state auxiliary-field (AFQMC) [52,53], and the finite temperature determinant quantum Monte Carlo (DQMC) methods [53][54][55][56]. Briefly, the DQMC (AFQMC) approach is based on the decoupling of the non-commuting terms of the Hamiltonian in the partition function (projection operator) by Trotter-Suzuki decomposition, which discretizes the imaginarytime coordinate τ in small intervals ∆τ , with the inverse of temperature T (projection time) being β = M ∆τ . Throughout this Letter, we choose ∆τ t = 0.1, with β in unit of t. Detailed introduction for these methodologies can be found, e.g., in Refs. 57-59.
Following the procedures described in Ref. 60, we implemented a sign-free AFQMC approach to the halffilling of the HHM, allowing us to analyze large lattice sizes, but, conversely, being restricted to the U eff ≥ 0 region. The properties of the U eff < 0 region, forbidden to our AFQMC method, are investigated by DQMC simulations. We recall that the DQMC method may exhibit sign problem for the HHM, depending on the strength of parameters. However, the average sign is less affected when U < λ, in particular to intermediate interaction strengths, allowing us to obtain the physical quantities of interest, in some cases up to L = 14 and β = 28. In fact, the DQMC average sign is strongly suppressed just around U ≈ λ, where our signfree AFQMC approach works. Thus, our AFQMC and DQMC simulations are used complementarily.
The charge and magnetic responses are quantified by their respective structure factors, i.e., and N = L × L being the number of sites. This allows us to probe their critical behavior by means of the correlation ratio with |δq| = 2π/L, q = (π, π), and ν ≡ cdw or afm. According to well established finite size scaling analysis, the critical region is determined by the crossing points of and next-nearest-neighbors (NNN) spin-singlet pairing operators for the s-wave symmetry, which are denoted by α ≡ s, s * , and s * * , respectively [66]; and also consider the d x 2 −y 2 -wave symmetry, α ≡ d (see, e.g., Ref. 67). Results: We first discuss the limit case of U = 0 in Eq. (1), i.e., the pure Holstein model. While the Peierls' argument [68] suggests an insulating charge ordered ground state for any λ > 0, one-dimensional systems exhibit a metal-CDW transition for a finite critical λ c , despite the perfect FSN. In two dimensions, the occurrence of a finite critical λ c in the square lattice is controversial: while variational QMC approaches provide evidence for λ c /t ≈ 0.8 [40,41], unbiased QMC results suggest the nonexistence of a finite critical point [38,39].
Here, we address this controversy by analyzing the critical behavior given by the CDW correlation ratio, Eq. (2), by means of DQMC simulations [69]. The quantum critical region may be estimated by the crossing of R cdw (L) as function of λ/t, as displayed in Fig. 2 for several lattice sizes. Here we investigate the ground state behavior by assuming β ∼ L z , with z = 1 being the dynamic critical exponent. Notice that, the correlation ratio increases monotonically as a function of the lattice size starting from λ/t = 0.8, clearly supporting long-range CDW ordering at this interaction strength and above. Furthermore, as showed in the inset of Fig. 2, the λ c (L, L − ∆L), i.e. the values of λ/t for the crossing points between R cdw (L) and R cdw (L − ∆L), are reduced when L increases, suggesting that, whether a metal-insulator transition exists, it should occur at smaller coupling strengths. A thorough determination of the existence of a critical point is given by a finite size scaling analysis of λ c (L, L − ∆L). Following the procedure adopted in Ref. 39, which hereafter is used to determine the CDW transitions, we perform a power law scaling [f (L) = a + bL c ] of the crossing points, as displayed in the inset of Fig. 2. Within this scaling, λ c (U = 0) is consistent with a vanishing or very small value, even when we adopt the less accurate polynomial fit, indicating that a finite critical e-ph coupling is not plausible for the square lattice Holstein model. The difference between the square lattice and one-dimensional systems may stem on the larger electronic susceptibility of the former, which diverges with the square logarithm of temperature. We now discuss the behavior of the HHM for U = 0, investigating initially the particular case of U = λ, by means of AFQMC simulations. We recall that, the HHM in the antiadiabatic limit (ω 0 → ∞) should exhibit a metallic behavior along the line U = λ. However, the occurrence of a retarded interaction in the adiabatic limit leads to a more complex ground state. Our QMC results for ω 0 /t ≤ 1 (not shown) exhibit an enhancement of the spin-spin correlations as a function of U/t, on the line U = λ, while the charge-charge response remains weak for any interaction strength. This result suggests an AFM ground state, with long-range order being probed by its correlation ratio. Similarly to CDW analysis, Fig. 3 displays the U c (L, αL), i.e. the values of U/t for the crossing points of R afm (L) and R afm (αL), for different phonon frequencies, and fixing β = 2L. Due to the large lattice sizes achieved in our AFQMC simulations, here we adopt a linear finite size scaling for the AFM transitions. As expected, the pure Hubbard model (black square symbols in Fig. 3) is AFM for any U > 0, i.e., U c = 0. However, in presence of e-ph coupling (along the line U = λ), a quantum phase transition occurs, changing from a correlated metallic-like ground state to an ordered AFM one, for a given coupling strength. For instance, for ω 0 /t = 1, one finds U c /t = λ c /t = 0.73 (6), that is, the ground state is AFM for any U = λ > 0.73t. The position  of this QCP strongly depends on the choice of ω 0 /t, as showed in the inset of Fig. 3. Such increasing behavior is consistent with the expectation of an emergent metallic behavior in the antiadiabatic limit. The properties of this correlated metallic-like state are discussed later. We proceed by investigating the quantum critical behavior for the general case, U = λ. To this end, one may analyze the correlation ratio as function of U/t, for fixed λ/t and ω 0 /t. For instance, Fig. 4 (a) displays AFQMC results for R afm (L) as function of U/t, for λ/t = 2 and ω 0 /t = 1. The crossing points U c (L, L−∆L) between R afm (L) and R afm (L − ∆L), as well as their finite size scaling, are displayed in Fig. 4 (c), leading to an AFM quantum phase transition at U AFM c /t = 1.88 (2). Similarly, the CDW QCP may be obtained by DQMC simulations of R cdw (L), as presented in Fig. 4 (b), with crossing points and finite size scaling shown in Fig. 4 (c), leading to U CDW c /t = 1.63(1). When the above analysis is repeated for other values of λ/t, we obtain the phase diagram presented in Fig. 1. It is worth mentioning that, for the range of parameters analyzed, we obtain continuous transitions for both AFM and CDW phases, without coexistence, and with a metallic-like (or SC) region between them. First order transitions may occur for stronger coupling, as suggested in Refs. 40 and 41. Finally, it is instructive to discuss the properties of such a region between AFM and CDW phases, from which it is expected the emergence of SC. Although establishing long-range SC order is challenging, its properties may be investigated from the tendency of  the pairing susceptibility as a function of temperature. For instance, by using the DQMC method, and fixing λ/t = 2, and U/t = 1.7, we observe an enhancement of χ α(sc) at low temperature, with the dominant symmetry being the d-wave [51], as presented in Fig. 5 (a). However, a more appropriate quantity for SC is given by extracting the particle-particle contribution of χ α(sc) , which defines the effective pair (vertex) susceptibility, i.e., χ eff α(sc) = χ α(sc) −χ α(sc) , withχ α(sc) being the noninteracting susceptibility [67]. A positive (negative) response of χ eff α(sc) corresponds to an enhancement (weakening) of pair correlations for the αth symmetry. Figure 5 (b) exhibits χ eff α(sc) for the data of panel (a), showing negative tendency towards all of the examined channels. In particular, the on-site s-wave has the largest negative response, which shows the harmfulness of the Hubbardlike term for local pairs formation. Since T sc ∼ ω 0 from the BCS theory [2], further insights about the nature of this region may be given by increasing ω 0 , while keeping U eff fixed, as displayed in Figs. 5(c) and 5(d), for ω 0 /t = 4, λ/t = 2, and U/t = 1.7. For these parameters, despite the increasing dominant behavior of χ α(sc) for the d-wave, only the NNN s-wave exhibits a positive effective susceptibility.
It is important to mention that, despite the dominant character of the d-wave in χ α(sc) , for both adiabatic and antiadiabatic limits, its negative tendencies for χ eff α(sc) suggest that such a channel may not lead to long-range SC order in the ground state of the half-filled HHM. Conversely, whether SC emerges, the results of Fig. 5 provide evidence for (nonlocal) s-wave. Indeed, since this metallic-like region is mostly in the negative U eff side of the phase diagram in Fig. 1, the s-wave symmetry seems more plausible than d-wave. Interestingly, such an enhancement of nonlocal s-wave response suggests that short-range charge/spin correlations may suppress the formation of local (s-wave) and NN (s * -wave) Cooper pairs, making the NNN ones (s * * -wave) the main channel for pairing. Conclusions: We have presented results for the HHM in the square lattice, using unbiased AFQMC and DQMC methods complementarily, which provide a broader picture about the physical responses of this model. In particular, we have shown that, different from onedimensional systems, the emergence of the CDW phase in the square lattice occurs for any λ > 0, for U = 0. However, these CDW correlations are strongly affected by a Coulomb interaction (U = 0), with the occurrence of AFM even at U eff < 0. We also observed the existence of a correlated metallic-like region between CDW and AFM phases, with an enhancement of nonlocal s-wave pairing, rather than d-wave. Despite the difficulty to establish long-range order for SC, one may expect that s-wave SC sets in at zero temperature. These findings constitute a significant step towards a better understanding of this model, by providing precise QCPs, and shedding lights over past theoretical inconsistencies. Furthermore, these results emphasize the role of the e-ph interaction in strongly correlated systems, which is crucial to charge order and pairing, and may be relevant for the physics of cuprates and transition-metal dichalcogenides.