Short-range charge fluctuations in the two-dimensional Hubbard model
Xinyang Dong Emanuel Gull
cond-mat.str-el

apsrev4-1




Department of Physics, University of Michigan, Ann Arbor, MI 48109, USA

Department of Physics, University of Michigan, Ann Arbor, MI 48109, USA


We investigate charge fluctuations in the two-dimensional Hubbard model as a function of doping, interaction strength, next-nearest-neighbor hopping, and temperature within the eight-site dynamical cluster approximation. In the regime of intermediate interaction strengths, we find that d density wave fluctuations, which had previously been postulated using analytical arguments, are present and strong but cannot be interpreted as the cause of the pseudogap in the model, due to their evolution with doping and interaction strength. For all parameters away from half filling, the charge fluctuations investigated, including d density wave fluctuations, are weaker than d-wave superconducting fluctuations.



Short-range charge fluctuations in the two-dimensional Hubbard model
    Emanuel Gull
    December 30, 2023
====================================================================




§ INTRODUCTION


The presence of large competing fluctuations of several kinds is one of the defining aspects of correlated electron systems. These fluctuations may then condense into phases that exhibit remarkable properties, including unusually high superconducting transition temperatures and interesting magnetism.

The cuprate superconductors are a paradigmatic example for such a competition.  Antiferromagnetic fluctuations are strongest in the `undoped' parent compounds but present over a large part of phase space. Superconducting fluctuations lead to a superconducting dome for dopings smaller than ∼ 20%. Charge phenomena  <cit.>, such as the famous stripes at 1/8th doping <cit.> or the features observed in scanning-tunneling experiments <cit.>, are present in several parts of phase space. Charge order with d-wave symmetry has also been found in RXS <cit.>.

A minimal model that describes many of the salient features of these materials is the single-band Hubbard model <cit.>. While unable to describe excitations involving high-lying orbitals, the model reproduces much of the observed low-energy phenomenology, including a pseudogap (PG) <cit.>, superconductivity <cit.>, and the response functions of Raman spectroscopy <cit.>, optical conductivity <cit.>, nuclear magnetic resonance <cit.>, and neutron spectroscopy <cit.>. It is therefore interesting to examine properties of the model in the context of cuprate physics. Due to the non-perturbative parameter regime relevant to the materials, reliable predictions have to resort to numerics, and a wide range of efficient numerical methods are able to describe the relevant parameter regime with consistent results <cit.>.

Hubbard model calculations find spin, charge, and superconducting fluctuations. Spin fluctuations are well understood and mainly dominant near half filling <cit.>. Calculations also find that strong short-wavelength spin fluctuations are  primarily responsible for the formation of the pseudogap <cit.>, i.e., the suppression of the density of states near the antinode but not near the node. Superconductivity is found unambiguously in the weak coupling regime <cit.>, and strong indications from dynamical cluster calculations show that superconductivity does persist to larger couplings <cit.>. In contrast, results from some newer methods find that in the absence of a next-nearest-neighbor hopping, it is charge (rather than superconducting) order that dominates the ground state <cit.>. However, all orders are in very close competition.
The precision to which the energetics of these phases is known is much better than the uncertainty in the model parameters, indicating that phenomena beyond simple Hubbard model physics may well force the system to choose one order over the other. 

At finite temperature, charge fluctuations, in contrast to antiferromagnetic and superconducting fluctuations, are less well investigated for the model without additional nearest-neighbor interactions, whereas the “extended” model has been studied extensively in recent years <cit.>. This is despite the fact that theoretical approaches have proposed unusual charge phenomena, such as the d density wave (DDW) order <cit.>, as candidates responsible for pseudogap physics <cit.>. It is therefore interesting to investigate the extent to which charge fluctuations are present in the model, and the extent to which they correspond to the proposed d density wave fluctuations, using numerical methods that generate these fluctuations dynamically from an underlying Hamiltonian. As we shall show below, the results are unexpected. First, we do find substantial d density wave fluctuations. However, while we find no area in parameter space where those fluctuations are dominant, they are comparable in magnitude to superconducting fluctuations (DSC). Also, while fluctuations are large in the general area of the pseudogap, the behavior with doping, interaction strength, and next-nearest neighbor hopping is not consistent with d density wave fluctuations as a mechanism for the suppression of the density of states.

The remainder of this paper is organized as follows. In Sec. <ref>, we describe the method used in this paper. In Sec. <ref>, we show the comparison between leading fluctuations in the 2D Hubbard model. In Sec. <ref>, we show the temperature evolution of DDW and DSC. In Sec. <ref>, we show the next-nearest-neighbor hopping evolution of DDW. Sections <ref> and <ref> discuss our results and present conclusions. 





§ METHOD
 
We investigate the two-dimensional Hubbard model on a square lattice with on-site interaction U and chemical potential μ,


    H=∑_kσ(ε_k-μ)c^†_kσc_kσ+U∑_in_i↑n_i↓.


Here i labels the lattice site, k the momentum, c^(†) annihilation (creation) operators, and n the density. ε_k=-2t(cos k_x + cos k_y) -4t' cos k_x cos k_y is the dispersion with hopping t and next-nearest-neighbor hopping t'. 

We define the single particle Green's function as 
G_σ(k_1τ_1,k_2τ_2) = ⟨ T_τ (c_k_1σ^† (τ_1) c_k_2σ(τ_2))⟩
and the two-particle Green's function as 
G_2, σ_1σ_2σ_3σ_4(k_1τ_1,...,k_4τ_4) = ⟨ T_τ (c_k_1σ_1^† (τ_1) c_k_2σ_2(τ_2) c_k_3σ_3^† (τ_3) c_k_4σ_4(τ_4))⟩, , with conservation of momentum k_1 + k_3 = k_2 + k_4, and  τ_i denote imaginary time points.


The generalized susceptibility is defined as <cit.> χ_σσ'(k_1τ_1,k_2τ_2,k_3τ_3,k_4 0) = G_2,σσσ'σ'(k_1τ_1,k_2τ_2,k_3τ_3,k_4 0) - 
G_σ(k_1τ_1,k_2τ_2)G_σ'(k_3τ_3,k_4 0).
Its Fourier transform in the particle-hole notation is

    χ_ph, σσ'^ωω'Ω(k,k',q) = ∫_0^β e^-iωτ_1e^i(ω+Ω) τ_2e^-i(ω'+Ω) τ_3
       ×χ_σσ'(kτ_1,(k+q)τ_2,(k'+q)τ_3,k' 0) dτ_1 dτ_2 dτ_3.

The “density channel" susceptibility is then defined as

    χ_d^ωω'Ω(k, k', q)=χ_ph, ↑↑^ωω'Ω(k, k', q) +χ_ph, ↑↓^ωω'Ω(k, k', q),

and can be decomposed using the Bethe-Salpeter equation into

    χ_d^ωω'Ω   (k, k', q) = χ_0^ωω'Ω(k, k', q) - 1/β^2 N^2χ_0^ωω_1Ω(k, k_1, q)
       ×Γ_d^ω_1ω_2Ω(k_1, k_2, q)χ_d^ω_2ω'Ω(k_2, k', q),

with χ_0^ωω'Ω(k, k', q) = -β N G_σ (iω, k)G_σ (iω+iΩ, k+q)δ_ωω'δ_k k' the bare susceptibility and Γ_d^ω_1ω_2Ω(k_1, k_2, q) the irreducible vertex in the density channel. N denotes the number of k points for the summation over k_1, k_2.











Linear response theory relates χ_d to a generating field Λ(k) as <cit.>

    -1/N_0∫_0^β dτ' ∑_k k'δ G_σ(k+q,0,k,0; Λ)/δΛ(k', τ') g(k)g(k')|_Λ=0   
    = 1/β^2 N_0∑_ωω',k k'χ_d^ωω'0(k, k', q)g(k)g(k')

with N_0 = ∑_k k'|g(k)g(k')| a normalization and g(k) a symmetry factor. The related order parameter is D = i ∑_k, σ g(k) c_k+q,σ^† c_k,σ <cit.>.

Superconducting fluctuations are related to the susceptibility in the particle-particle channel which follows from an analogous derivation <cit.>, 

    1/N_0∫_0^β dτ' ∑_k k'δ F(k,0; η)/δη(k', τ') g(k)g(k')|_η=0   
    = 1/β^2 N_0∑_ωω',k k'χ_pp ↑̅↓̅^ωω'0(k, k', q = 0)g(k)g(k')

with F(k, τ)=-⟨ T_τ c_k, ↑(τ) c_-k, ↓(0) ⟩ the anomalous Green's function, and η(k) the generating field. The related order parameter is P = ∑_k g(k) c_k,↑ c_-k,↓ <cit.>.

We define the right hand side quantity in both Eq.<ref> and Eq. <ref> as χ_g,

    χ_g = 1/β^2 N_0∑_ωω',k k'χ^ωω'0(k, k', q)g(k)g(k').

This is the central quantity investigated in this paper.

Phase transitions are indicated by a divergence of the susceptibility χ^ω,ω',Ω. Since we expect large values of χ^ω,ω',Ω to be caused by vertex contributions, we define

    V_g :=    (χ - χ_0)_g = 1/β^2 N_0∑_ωω',k k' g(k)g(k')  ×
       (χ^ωω'0(k, k', q) - χ_0^ωω'0(k, k', q))

in analogy to P_g in the superconducting case <cit.>.
This quantity highlights vertex contributions by subtracting band-structure and single-particle effects contained within the bare susceptibility.

We will show that charge fluctuations with momentum transfer q = (π, π) and d-wave symmetry are large.
These fluctuations are known as d density waves <cit.>.
Ref. Nayak2000 defines density order in analogy with superconductivity, with a singlet order parameter in the form ⟨ψ_α ^†(k+Q) ψ_β(k) ⟩ = Φ_Q f(k) δ_α,β for symmetry factor f(k).
Two of the possible orders are f(k) = sin k_x for p_x symmetry and f(k) = cos k_x - cos k_y for d_x^2-y^2 symmetry.

The d density wave, also called the “staggered flux state", occurs at Q = (π, π). The p_x density wave (PDW), also called the “bond order wave", happens at Q = (π, π) and (π, 0). For a real space representation of DDW, see Fig. 2 of Ref. Nayak2000 and Fig. 2 of Ref. Chakravarty2001.

The exact solution of the Hamiltonian Eq. <ref> is unknown. Here we use the dynamical cluster approximation (DCA) on a cluster with N_c = 8 sites, which approximates the self-energy by N_c “coarse-grained” patches in which the self-energy is momentum independent, but retains the full frequency dependence <cit.>.
The method is approximate but controlled in the sense that the thermodynamic limit for local quantities are approached ∼ N_c^-2 as cluster size N_c →∞ <cit.>.
The eight-site cluster used here is a compromise chosen large enough to accommodate a clear nodal-antinodal differentiation <cit.>, pseudogap regime <cit.>, and superconductivity <cit.>, while remaining cheap enough for simulation of a wide range of parameters.

Figure <ref> shows the phase diagram of the model obtained within the 8-site DCA with parameters U/t=7 and t'=-0.15t. 

These parameters are chosen to represent the overall phase diagram common to several cuprates <cit.>.
The pseudogap regime is obtained by observing a suppression in the single-particle spectral function. The superconducting phase is computed in a Nambu formulation and defined as the area where the anomalous Green's function F(k,τ) = -⟨ T_τc_k,↑(τ)c_-k,↓(0)⟩ at k = (0, π), (π, 0) becomes nonzero.








DCA only yields the cluster Green's functions and corresponding cluster susceptibilities χ_c,d^ωω'Ω(K, K', Q) at the cluster momenta K, K', and Q. 
The corresponding approximation to Eq. <ref> is obtained by interpolating χ_d(k, k', Q) = (χ_0^-1(k, k', Q) + χ_c,d^-1(K, K', Q) - χ_c,0^-1(K, K', Q))^-1. This is analogous to identifying the DCA cluster vertex with the lattice vertex Γ_d^ω_1ω_2Ω(k_1, k_2, q) <cit.>.


Within the eight-site cluster, the symmetry factors we used below corresponding to s-, p_x, d_xy, and d_x^2-y^2 symmetry are defined as
g_s(K) = 1, 
g_p_x(K) = sin(K_x), 
g_d_xy(K) =  sin(K_x)sin(K_y), and
g_d_x^2-y^2(K) = cos(K_x)-cos(K_y).
The DCA approximation generates strong antiferromagnetic (AFM) fluctuations with a correlation length comparable to the cluster size. If the establishment of long-range AFM order is allowed, the system chooses an ordered state at a temperature above the onset of the PG or superconductivity. This is a finite size effect <cit.>. For this reason, we suppress magnetic long-range order and only show results obtained in the paramagnetic state, which have the correlation length of AFM fluctuations restricted to the cluster size <cit.>.

Our results are obtained with a continuous-time auxiliary field quantum Monte Carlo impurity solver <cit.> based on the ALPS <cit.> libraries. The summation over fermionic frequencies in Eq. <ref> and Eq. <ref> goes over all frequencies from -∞ to ∞. In our calculations, only a finite number of frequencies are available, but the asymptotic behavior of χ_0 is known analytically <cit.>, . In the results presented here, we use 36 fermionic frequencies for β t= 5, 10; 50 fermionic frequencies for β t= 15, 20; and 80 fermionic frequencies for β t= 30 on both positive and negative sides to compute the vertices. The relative change for omitting the last eight frequencies on each side is on the order of 10^-3. The asymptotic behavior of vertex Γ and F is also analyzed in Ref. wentzell2016highfrequency which provides an alternative way of treating the high frequency behavior.

χ is computed with the number of frequencies listed above for the vertex correction part, plus χ_0 computed with 1024 fermionic frequencies (both positive and negative) and supplemented with an analytically known asymptotic correction.













§ CHARGE AND SUPERCONDUCTING FLUCTUATIONS
 
Figure. <ref> shows six panels for the leading fluctuations at β t = 20 and t'/t = -0.15. The three rows represent p_x and d_x^2 - y^2 density fluctuations, and d_x^2 - y^2 superconducting fluctuations. The left columns show |χ_g| (see Eq. <ref>), the right columns  |(χ - χ_0)_g| (see Eq. <ref>). Each panel displays data as a function of U and doping δ, with δ = 0 corresponding to half filling. For superconductivity, only fluctuations with d_x^2 - y^2 symmetry (DSC) are large <cit.>. 


Fig. <ref> (b), (d) and (f) illustrate that χ itself is not a good measure for the correlation contribution that may eventually drive the system to an ordered state, as most of χ stems from χ_0 = G G.
Fig. <ref> (c)-(f) show that the amplitudes of DDW and DSC fluctuations are comparable (at the same order), implying competing fluctuations. However, we find numerically that DSC fluctuations are always larger than DDW for the parameters examined. 
In a small regime of parameter space, where U ∼ 6.5t with slight electron doping, p-density wave fluctuations are the dominant charge fluctuation (panel (a) and (b)). 

The maximum of DSC is on the electron-doped side, while the maximum of DDW is on the hole-doped side, both at intermediate interaction strength U. In addition, DSC fluctuations are suppressed in the psedogap regime starting from U ∼ 6t (see Refs. Chen2015, Gull09_8site), while DDW fluctuation starts to show suppression around half filling for U ∼ 6.5t, which corresponds to the onset of the Mott insulator <cit.> in this approximation. PDW does not show any suppression by either the pseudogap or the Mott insulating state; its maximum is near U ∼ 6.5t, which is the same interaction strength where DDW shows a suppression near half filling.



§ TEMPERATURE EVOLUTION
 
In order to investigate the competition between d density wave and d wave superconducting fluctuations in more detail, we explore their temperature evolution with different dopings in Fig. <ref>. 
We show the results at U = 7 t and t' = -0.15t.
Panel (b) shows that away from half filling, the vertex part of DSC increases as temperature decreases for all doping levels investigated (see also Ref. Chen2015), whereas panel (a) shows that the vertex correction part of DDW increases as temperature decreases in the underdoped regime away from half filling but rapidly decays to zero for large doping. The maximum of the DSC fluctuations are near the maximum T_c <cit.>, while the corresponding maximum DDW fluctuations occur at slightly lower doping. The transition to superconductivity on the hole-doped side will take place near β t = 35 at optimal doping in this model, i.e., at a temperature about twice below where these results have been obtained.

The amplitude of the vertex correction part of DDW and DSC fluctuations as a function of temperature at the doping level corresponding to largest DDW fluctuation and largest DSC fluctuation is shown in Fig. <ref>. These results are obtained in the paramagnetic state but reach temperatures just above the superconducting transition. At the doping level where DDW fluctuations are strongest (corresponding to μ/t=-1.4, δ∼ 0.09), DSC fluctuations are substantially larger than DDW, and increase faster as temperature decreases. We have been unable to find a region of parameter space where d density wave order prevails over superconductivity around optimal doping. At the doping level corresponding to the largest DSC fluctuation(μ/t=1, δ∼ 0.09) we could find, DDW fluctuations first increase as temperature decreases, then start to decrease at β t∼ 20. This result is consistent with the findings of Refs. Macridin2004, Honerkamp2002.




























§ PARTICLE-HOLE ASYMMETRY
 
Figure. <ref> shows the dependence of PDW and DDW on interaction and doping for three values of t'. t' is to shift the Van Hove singularity of the density of states toward hole doping and destroys particle hole symmetry.

We find that as t' / t is changed from 0 to -0.2 and correspondingly the Van Hove singularity more toward the hole-doped side, PDW fluctuations spread out over a larger area of the parameter space at the electron-doped side, but their maximum does not change significantly.
On the other hand, the maximum of DDW fluctuations shifts toward hole doping and their intensity decreases substantially.

The particle-hole asymmetry of DSC shows a different trend from DDW <cit.>. As -t' / t increases, the maximum of DSC fluctuations moves to the electron-doped side. This is consistent with a scenario where the  establishment of a pseudogap on the hole-doped side suppresses d-wave superconducting fluctuations <cit.>. 

Panel (a) of Fig. <ref> shows the amplitude of the vertex corrections of DDW and DSC with different t' at U = 7t, β t= 15, and doping levels corresponding to the largest DDW fluctuations. As -t' increases, the maximum of DDW shifts toward the hole-doped side. Since the overall intensity of the fluctuations decreases, the maximum value of DDW fluctuations saturates and then drops. At the same time, for all the t' values we explore, DDW fluctuations are still weaker than DSC fluctuations around the optimal doping for DDW fluctuations.
Panel (b) shows the doping levels where the largest DDW fluctuations are found.

Eight-site DCA shows a narrowing of the pseudogap on the electron-doped side, and for t' < -0.15t has a first-order transition between the Mott insulator and a momentum-dependent Fermi liquid state <cit.>. This rapid change in the single-particle quantities does not show an analog in the DDW or PDW fluctuations.



§ DISCUSSION
 
The causes and consequences of DDW have been much debated. Early studies <cit.> considered such fluctuations as a candidate for the mechanism of the pseudogap. Numerical calculations of the Hubbard model <cit.> argued that this type of fluctuation is not strong enough to form an ordered state in the strongly correlated region, and is always dominated by DSC fluctuations, while renormalization group studies find that for SU(N), when N > 6, DDW becomes the leading instability <cit.>.
Other works have studied the coexistence of DDW and DSC orders <cit.>.
The d-wave fluctuations at Q = (0, 0) or d-nematic fluctuations have  been studied in Ref .Honerkamp2002, where it is shown that they grow together with DDW, DSC, and antiferromagnetic fluctuations.

More recent variational Monte Carlo studies <cit.> found staggered flux states (DDW states) in the strongly correlated underdoped regime of the Hubbard model, for interactions smaller than the Mott transition. They proposed this state as a candidate for an anomalous “normal state” competing with DSC since its properties are similar to those of the pseudogap. In these more recent works, the state does not coexist with DSC.

Our results clarify some of these arguments. 
Our interaction, doping, and t' evolution of DDW correlations show that the evolution of the pseudogap and that of DDW correlations do not track each other. 
This can be seen from the fact that in the area around half filling, where there is a pseudogap, DDW correlations show a suppression. However, the interaction strength where DDW starts to show a suppression is at U ∼ 6.5t, not at U ∼ 6t, where the pseudogap opens.
As -t'/t increases, the pseudogap regime moves to the hole-doped side (shown in Fig. 3 and Fig. 4 of Ref. Gull09_8site), while the suppression of DDW moves to the electron doped side. These trends show little overall correlation between areas with the largest DDW fluctuations and the appearance of a pseudogap. 

Without entering an ordered phase, it is difficult to make statements about a potential coexistence of DDW and DSC. However, comparison between the temperature evolution of DDW and DSC at U = 7t (Fig. <ref> and Fig. <ref>) clearly shows that around the optimal doping for DDW, DSC correlations are dominant over DDW and will order first. Thus, if there is a coexistence regime, it is likely fully contained inside the DSC dome.

Previous works <cit.> investigated the competition between DSC and DDW in the Hubbard model with DCA on 4-site clusters at the two doping levels δ = -0.05 and δ = -0.25. The first doping level is in the pseudogap regime, the second far in the overdoped.
The main finding, namely that the susceptibility corresponding to the d density wave does not diverge, indicating the absence of a possible transition to the DDW state, is reproduced by our calculations. However, the claim that both DDW and DSC correlation functions are enhanced in the pseudogap regime is inconsistent with our more detailed calculations.
The reason is that the intermediate maximum, which we find around δ = -0.09, is missed by the coarse doping resolution employed in Ref. Macridin2004.

DDW, the staggered flux (SF) state, and its competition with DSC, have been studied extensively in the t-J model <cit.>, which maps to the Hubbard model in certain limits. In the t-J model, a well-defined SF state exists in a small doping region in the phase diagram <cit.>.
Large-N calculations on the model <cit.> related the DDW self-energy to the pseudogap via the formation of Fermi arcs. This is different from the results presented here, but whether the differences should be attributed to the different models, parameter regimes, approximations, or post-processing procedures is a question for future study. 


Fluctuation diagnostics <cit.> attributes contributions to the self-energy to fluctuations of different types. Reference. Dong2019 shows that in the parameter regime we explore in this paper, the dominant fluctuation in the normal state is always magnetic. It is shown in Ref. Gunnarsson15 that these magnetic fluctuations are responsible for the major contribution to the single-particle self-energy, and thereby for the suppression of the density of states. The pseudogap in the Hubbard model can therefore clearly be attributed to magnetic fluctuations. 
DDW and DSC fluctuations do not directly contribute to a single-particle self-energy, since the summation of a fluctuation with a d-wave symmetry tends to cancel out in the calculation of the self-energy <cit.>.


The DCA simulations performed here are insensitive to the stripe order widely found in experiment and in numerical ground-state calculations. In order to find a transition to an ordered state in DCA, the ordering vector typically needs to be commensurate with the cluster geometry. However, the stripe orders, e.g., found in Ref. Zheng17_stripe are too large to fit into the DCA cluster studied here. Thus, while the method is sensitive to charge fluctuations on a length scale smaller than the cluster size, DCA is not expected to find the period 4 and period 5 stripes of Ref. Zheng17_stripe. Appropriately chosen larger clusters may find these stripes, but finite size effects would likely overestimate their contribution. The unbiased detection of such orders with DCA or lattice methods is an important open problem.

Our conclusion on the relation between PG and DDW, however, is not influenced by the omission of the stripe order. With an eight-site cluster simulation, we can see a clear signal of a pseudogap from the single particle spectrum, while DDW fluctuations have been shown to be unrelated to it. 






§ CONCLUSIONS
 
In conclusion, we have investigated the physics of short-range charge fluctuations in the strongly correlated regime of the 2D Hubbard model within the 8-site DCA approximation.
We have found that the dominant charge fluctuations have d-wave symmetry, apart from a small regime near the Mott transition, where we find p-wave charge fluctuation.
For all doping, interaction, and next-nearest-neighbor hopping parameters investigated, we showed that superconducting d-wave fluctuations are always stronger than charge fluctuations away from half filling. At U = 7t, at the doping level that is most favorable for DDW, DSC will order first as temperature decreases, showing that if there were a coexistence of DDW and DSC orders, the coexistence area would likely be fully contained inside the DSC dome.

Our parameter scans show that DDW fluctuations cannot be viewed as the cause of the pseudogap in the single-particle density of states, as the change of DDW fluctuations does not match the evolution of the pseudogap. This is consistent with the results of several recent works, including Refs. Gunnarsson15,Dong2019, that show convincingly that the pseudogap can be attributed to strong short-wavelength AFM correlations.





This work was supported by NSF DMR-1606348 and the Simons Collaboration on the Many-Electron Problem. We thank Kai Sun for insightful discussions.





apsrev4-1

