Competition between antiferromagnetic and charge density wave fluctuations in the extended Hubbard model

TL;DR

This study uses an advanced parquet approach to analyze the competition between charge and spin fluctuations in the extended Hubbard model, revealing that antiferromagnetic fluctuations dominate except near the charge density wave transition.

Contribution

It extends the parquet method to include non-local Coulomb interactions, providing a comprehensive analysis of fluctuation competition in the extended Hubbard model.

Findings

Charge and spin fluctuation competition is narrow near the phase transition.

Antiferromagnetic fluctuations dominate due to lattice geometry.

Multiple observables consistently support the conclusions.

Abstract

By extending our {\it victory} implementation of the parquet approach to include non-local Coulomb interactions, we study the extended Hubbard model on the two-dimensional square lattice with a particular focus on the competition of the non-local charge and spin fluctuations. Surprisingly, we find that their competition, as the mechanism driving the phase transition towards the charge density wave, dominates only in a very narrow parameter regime in the immediate vicinity of the phase transition. Due to the special geometry and the Fermi surface topology of the square lattice, antiferromagnetic fluctuations dominate even for sizable next-nearest neighbor interactions. Our conclusions are based on the consistent observations in both the single- and two- particle quantities, including the self-energy, the single-particle spectral function, the two-particle susceptibility, the density-density vertex function and the optical conductivity. Our work unbiasedly establishes the connection of these quantities to the charge fluctuations, and the way of interpretation can be readily applied to any many-body method with access to the two-particle vertex.

Equations

$$H=$$

$$+U\sum _i{n}_{i\uparrow }{n}_{i\downarrow }+\frac{V}{2}\sum _{⟨ij⟩,\sigma {\sigma }^{\prime }}{n}_{i\sigma }{n}_{j{\sigma }^{\prime }},$$

$$U_d^{\mathbf{k}{\mathbf{k}}^{\prime }\mathbf{q}}$$

$$U_m^{\mathbf{k}{\mathbf{k}}^{\prime }\mathbf{q}}$$

$$U_s^{\mathbf{k}{\mathbf{k}}^{\prime }\mathbf{q}}$$

$$U_t^{\mathbf{k}{\mathbf{k}}^{\prime }\mathbf{q}}$$

$${\Phi }_{d/m}^{k{k}^{\prime }q}$$

$${\Phi }_{s/t}^{k{k}^{\prime }q}$$

$${\Sigma }_k=$$

$$\displaystyle\bigg{.}+V_{\mathbf{q}}{F}^{kk^{\prime}q}_{d}\bigg{]}+\frac{1}{N\beta}\sum_{k^{\prime}}G_{k}\left[U+2V_{\mathbf{q}=0}-V_{\mathbf{k}^{\prime}-\mathbf{k}}\right]\;,$$

$$\chi (q)$$

$${\lambda }_{d/m}^q{\varphi }_{d/m}^{kq}$$

$${\lambda }_{pp}^q{\varphi }_{pp}^{kq}$$

$$\chi =$$

$$\underset{\equiv \chi }{\underbrace{-\frac{2}{N\beta }\sum _k{\left[{\gamma }_P^{\mathbf{k}\frac{\mathbf{q}}{2}}\right]}^2{G}_k{G}_{q+k}}}$$

$${\gamma }_P^{\mathbf{kq}}$$

$${\mathbf{j}}_{\mathbf{q}=0}(\omega )$$

Full text

I Introduction

The study of strongly correlated electron systems is notoriously difficult and represents one of the greatest challenges in contemporary condensed matter physics. The competition of the various entangled degrees of freedom, such as the charge, spin and orbital contributions, is often the source of emerging exotic physics. On the other hand, it imposes a great challenge to theory. The interest in correlated electron models, especially the Hubbard model with an interaction strengths comparable to the bandwidth, is highly motivated by the unconventional physics beyond the simple itinerant or localized picture discovered in high-temperature superconductors 1; 2; 3; 4; 5; 6. The half-filled, two-dimension Cu-O plane shows magnetic long-range order, which can be gradually removed by electron- or hole-doping, necessary for the emergence of superconductivity.

Besides magnetic order and fluctuations also the relevance of charge density wave fluctuations and order is discussed—already in the early days of high-temperature superconductivity 7 but more intensively after experimental evidence thereof 8; 9; 10. These charge-density wave fluctuations are fostered by non-local Coulomb interactions, leading to the so-called extended Hubbard model (EHM) 11; 12; 13; 14; 15; 16 as a minimal model. The EHM Hamiltonian reads

[TABLE]

where $c_{i\sigma}^{(\dagger)}$ is the annihilation (creation) of electrons with spin $\sigma\in\{\uparrow,\downarrow\}$ at lattice site $i$. The particle density is hence given by $n_{i\sigma}=c^{\dagger}_{i\sigma}c^{\phantom{\dagger}}_{i\sigma}$. The kinetic term of the Hamiltonian (1) describes the hopping of electrons with amplitude $t$ to neighboring sites only, via $\langle ij\rangle$; The screened long-range Coulomb potential is described by an on-site density-density interaction with strength $U$ and by a $V$-term including interaction between electrons on neighboring sites. In this paper, we restrict ourselves to the case of only positive $U$ and $V$ since the Coulomb interaction is repulsive, albeit low-energy effective models may also be attractive and show interesting physics 17.

As already mentioned, the presence of non-local interaction term introduces, in addition to the antiferromagnetic correlation caused by the on-site interaction, the intersite charge-density correlations. The competition of the two types of correlations leads to a phase transition between the spin-density wave (SDW) and the charge-density wave (CDW) in the low-temperature regime. From the low-energy effective theory point of view, a weak non-local interaction can be viewed as an additional screening of the local interaction. An effective local interaction, after integrating out the non-local interaction, can then be obtained. However, the intersite fluctuations (or the non-local charge fluctuations) are absent in such effective theory. That is, in the effective theory the bare interaction is renormalized by the frequency/energy dependent dielectric function determined from the charge-charge correlation function and the effective local interaction becomes dynamic, i.e. $U(\omega)$ 18; 19; 20; 21 With the increase of the non-local interactions, the ”screening” picture is, however, no longer appropriate due to the underestimation of the non-local interaction effect and a more complete description of the two interactions in the EHM is required.

It is, however, difficult to treat the two types of interaction in Eq. (1) simultaneously. Let us start with the exact solution in the atomic limit, i.e. for $t_{ij}=0$. In one dimension at half-filing, the ground states of the atomic EHM is in a CDW ordered phase when $U<0$ irrespective the value of $V$. This phase is characterized by the alternative appearance of the empty and the doubly occupied states, forming the CDW with two-site periodicity. It extends to the positive value of $U$ as long as $U/2<V$, otherwise the Mott phase is realized where only singly occupied states are allowed 22. Similar analysis can be carried over to higher dimensions, and the transition to the CDW phase was found to occur at $U=2n_{dim}V$, where $n_{dim}$ is the system dimensionality 23; 24.

Wherever the electronic hopping process is recovered, the EHM loses its exactly solvable limit. Either approximations or numerical methods have to be consulted 25; 26; 27; 28; 29. Figure 1 displays the phase diagram obtained in the Hartree-Fock approximation and the dynamical mean-field theory (DMFT) for the two-dimensional EHM 30. A strong coupling expansion around the atomic limit predicts the transition to occur at values of $V>U/4$, 23 i.e. for $V$’s larger than in the atomic limit. In DMFT, the non-local interaction $V$ reduces to its Hartree contribution 31 but the difference to the Hartree-Fock solution of the EHM is that the local interaction $U$ remains dynamic. This DMFT solution 30 predicts the transition at a larger $V$ value than in the Hartree-Fock or atomic limit solution and it increases with the increase of $U$ as shown in Fig. 1. The light/dark blue triangles in Fig. 1 denote the parameters at which our calculations were performed, which will be discussed later in Section III. The shift of the phase transition line towards larger value of $V$ was confirmed by exact numerical simulations of DQMC 32.

In contrast to the DMFT, the extended DMFT (EDMFT) 33 contains the non-local $V$-term in its construction and treats it self-consistently in the effective polarization function. In EDMFT the transition was also found at $4V/U>1$ in the strong coupling regime, which is similar to the DMFT solution. The similarity is partially due to the fact that the EDMFT approximation also neglects the non-local fluctuations which become crucial in low-dimensional systems. The diagrammatic extensions of EDMFT, i.e. the dual boson approach 34, captures the non-local fluctuations and hints to a transition closer to $4V/U=1$ (still however with $4V/U>1$) in the intermediate coupling regime as well. The transition boundary in the dual boson approach lies in between the solid and dashed lines in Fig. 1.

To go beyond the DMFT/EDMFT local approximations, either a cluster 36 or a diagrammatic extension 37 of DMFT is required. In a recent dynamical cluster approach (DCA) study 38; 39, a similar tendency was confirmed. In this study, a momentum cluster of 20-sites was considered where the short-range fluctuations within the cluster were fully taken into account. Note, however, that for the larger values of $U$ used in this DCA study, the system without or with weak non-local interaction $V$ is in the ordered antiferromagnetic (AFM) phase, which is a finite-size effect. This AFM ordering is suppressed within the parquet approach used here already for the finite lattices considered.

The aim of this paper is to provide an unbiased study of the EHM with both local and short-range fluctuations included in a way that both the single- and two-particle Green’s functions are self-consistently determined. We study, for a wide range of non-local interactions, the competition between magnetic and charge fluctuations before the phase transition. Both the one- and two-particle quantities indicate that the charge fluctuations dominate only in the vicinity of the transition, leading to a suppression of the AFM fluctuations. But otherwise, even for sizable non-local interactions, AFM fluctuations prevail. The precise transition boundary is not of the prime interest of this work, instead we want to show how the transition is approached and how it is probed in the parquet approach, which is also of interest to many other cluster methods that have the access to the two-particle vertex functions. This work also represents a methodological progress which further extends the victory solver of the parquet equations 40; 41; 42 to include the non-local Coulomb interactions in a lattice model.

Throughout the paper, $t\equiv 1$, $\hbar\equiv 1$, and $k_{B}\equiv 1$ are set as the units of energy, frequencies and temperature. The square lattice with average number of electrons per site $n\equiv 1$ will be studied at temperature $T=t/6$ and onsite interactions are taken as $U=2t$ or $U=4t$. The nearest-neighbor interaction $V$, which reads $V_{\mathbf{q}}=2V[\cos(q_{x})+\cos(q_{y})]$ in momentum space, will be varied as fractions of $U/4$ for a given $U$ to gradually approach the transition.

Section II recapitulates the parquet equations and explains the methodological development of the present paper. Section III presents the results obtained, except for the optical conductivity which is discussed in Section IV.

II Method: parquet approximation to the EHM

First, we will explain the methodological development of the victory package 42 with respect to the non-local interaction. The extension is straightforward in the sense that the kernel approximation, which is the key concept of the victory implementation of the parquet approach, can still be applied. However, due to the non-local interaction the various two-particle vertex functions become more strongly momentum-dependent. Thus, one has to introduce the two-level kernel approximations 41 in each momentum patch, as will be explained in this section.

In the parquet approach, a set of exact equations that couple the one- and two-particle vertex functions is solved iteratively until a self-consistency in both levels is achieved. This is conceptually different from any single-particle theory, where only the self-consistency imposed by the Dyson equation, such as the DMFT. In the parquet approach, all the single-particle quantities including the Green’s function and the self-energy are determined by the two-particle vertex functions in all fluctuating channels. As a result, it allows for an unbiased treatment of all fluctuations simultaneously. This is also in sharp contrast to any ladder approximation where only a set of a priori chosen ladder diagrams are considered in favor of only certain type of fluctuations. More details about the parquet approach can be found in references 43; 44; 45; 46; 47; 48; 49; 40; 50; 51; 41; 42. Here we would like to concentrate on the extension of it to the EHM. Thus, only the most relevant formulas will be given.

In the parquet approach, the only input needed is the fully irreducible two-particle vertex $\Lambda^{kk^{\prime}q}$. Here, the four-vector notation with $k=(\mathbf{k},\nu_{n})$ and $q=(\mathbf{q},\omega_{n})$ is used with the momenta $\mathbf{k}$ and $\mathbf{q}$, the discrete Matsubara frequencies are $\nu_{n}=\frac{\pi}{\beta}(2n+1)$ (fermionic) and $\omega_{n}=\frac{\pi}{\beta}2n$ (bosonic) with $n\in\mathbb{Z}$ and the inverse temperature is $\beta=1/T$. In this work we employ the parquet approximation (PA),40 in which the fully irreducible two-particle vertex $\Lambda^{kk^{\prime}q}$ is approximated by its frequency independent lowest-order contribution, namely the bare interaction. Compared to the Hubbard model, the inclusion of the nearest-neighboring Coulomb interaction $V$ leads to additional terms in the lowest-order vertex functions:

[TABLE]

in the respective spin channels: density ($d$), magnetic ($m$), singlet ($s$) and triplet ($t$) (For details of the spin-diagonalized notation used throughout this manuscript see e.g. Ref. 37). More complicated than the case of the Hubbard model, all three bare vertex functions become momentum-dependent. As a result, the reducible vertex functions $\Phi_{d/m/s/t}$ depend not only on $q$ but also $k$ and $k^{\prime}$ even in the first iteration:

[TABLE]

The crucial observation of the victory implementation is the simpler structure of the reducible vertex functions which are used for extrapolation to reach the frequency-asymptotic of various two-particle vertex functions. In the Hubbard model, different momentum patches have essentially the same frequency-asymptotics 41 (known as kernel-I approximation in our notation) as can be seen from the only $q$ dependence of $\Phi_{d/m/s/t}$ calculated in the first iteration of Eqs. (6) and (7) if $U^{\mathbf{k}\mathbf{k}_{1}\mathbf{q}}$ is momentum independent for the Hubbard model, cf.  Eqs. (2)-(5). With the presence of the $V$-interaction, this changes and $\Phi_{d/m/s/t}$ depend on all three momenta. As a result, the reducible vertex functions display a $\mathbf{k},\mathbf{k^{\prime}}$-dependence even in the kernel-I function. We, therefore, modified the kernel-I function to faithfully incorporate such difference in different momentum patches. Instead of just taking one background value of the reducible vertex function, in each momentum patch (characterized by the different combinations of $\mathbf{k}$ and $\mathbf{k^{\prime}}$) we take a different background value to represent the kernel-I function of this momentum patch.

The change of the bare vertex function also affects the evaluation of the self-energy $\Sigma_{k}$ in the Dyson-Schwinger equation

[TABLE]

with $F$ being the full vertex in $d/m$ channels. The other parquet equations formally remain unchanged.

Following Ref. 52, we also improve the evaluation of Bethe-Salpeter equations by including, in the frequency part of the $k_{1}$ sum, the exactly known high frequency asymptotics of the full vertex $F$ and the irreducible vertex $\Gamma$ 37; 53. A similar regularization (extension to high frequencies) is employed for the Dyson-Schwinger equation (8). Furthermore, victory also benefits from the point group symmetry of the square lattice, where only the transfer momentum $\mathbf{q}$ in the irreducible Brillouin Zone (BZ) need to be stored in memory. With these changes *victory *can be used efficiently to study single-band correlated many-body systems with non-local interactions.

III Results

The evaluation of the parquet equations is numerically very demanding 41; 42; 54 and, at the moment, only possible for small clusters. All the results presented in this work were obtained on a $6\times 6$ momentum cluster (corresponding to a $6\times 6$ cluster in real-space with periodic boundary conditions), which, nevertheless, represents the state of the art in the solution of the parquet equation. Whenever possible, we used a finer grid for the Green’s function, as explained in Ref. 54. Since the bottleneck of the calculations is the memory consumption, only relatively small number of Matsubara frequencies can be taken, which restricts the computations to rather high temperatures. Here the temperature was set to $T=t/6$ and $N_{f}=96$ Matsubara frequencies were used.

III.1 Single-particle self-energy and spectral function

In Fig. 3 we first show the single-particle self-energy $\Sigma(k)$ as a function of the Matsubara frequency. Two different values of $U$ were considered in these calculations, i.e. $U=4t$ in (a-c) and $U=2t$ in (d-f). Each row of Fig. 3 corresponds to a different $\mathbf{k}$-point in the BZ, i.e. $(0,0)$, $(0,\pi)$ and $\left(\frac{\pi}{3},\frac{2\pi}{3}\right)$. The latter two $\mathbf{k}$-points, i.e. $(0,\pi)$ and $\left(\frac{\pi}{3},\frac{2\pi}{3}\right)$, reside on the Fermi surface at half-filling. We vary the values of $V$ to probe the competition between the local and non-local Coulomb interactions in the EHM occurring around the transition boundary. We denote the results corresponding to different $V$ by different symbols and colors in Fig. 3.

As can be clearly seen in Fig. 3(a) and (d), away from the Fermi surface, the self-energy becomes smaller (in absolute terms) with the increase of $V$, in nice agreement with the expectation of the non-local interaction effect. As we know, if interaction between any two electrons becomes a constant regardless of their distance, the system is only determined by the single-particle hopping. The self-energy becomes frequency independent, i.e. a constant for the occupied and unoccupied $\mathbf{k}$-points. The EHM with only nearest-neighbor $V$ interpolates the Hubbard model and the above-mentioned infinite-range interaction model. Consequently, it will display a finite but less correlated self-energy compared to that of the Hubbard model, cf. Ref. 52.

A similar reduction is observed for the points on the Fermi surface, i.e. for $\mathbf{k}=(0,\pi)$ and $\mathbf{k}=\left(\frac{\pi}{3},\frac{2\pi}{3}\right)$, as long as $4V/U<1$. Calculations with parameters $4V/U=0.2$ and $4V/U=0.4$ at $U=4t$ and $U=2t$ (not shown here) also confirm this trend.

However, close to the transition boundary at $4V/U\sim 1$ and $U=4t$ (and at a slightly smaller $V$ for $U=2t$) the self-energy changes dramatically, i.e. the first few lowest-frequency points of $\text{Im}\Sigma$ clearly deviate from those for smaller $V$ values. They even start to increase in absolute values for $U=2t$ and $4V/U=1.2$, see Fig. 3(e) and (f). At $4V/U=1.2$ and $U=2t$ the curvature of the self-energy changes completely, leading to a large negative $\text{Im}\Sigma$ at $\nu_{n}\to 0$. This indicates a pseudogap or at least bad metallic paramagnetic phase close to the CDW transition. As it will become clear later, this is due to the enhanced charge fluctuations caused by the non-local Coulomb interaction $V$ (cf. the eigenvalues and susceptibilities in Fig. 4). We hence interpret the large negative $\text{Im}\Sigma$ in Fig. 3(e) and (f) as a precursor of the insulating CDW phase, which would be further enhanced at lower temperatures. That is, we do not yet have CDW order yet but already long-range CDW correlations, which leads to a (pseudo)gap similar as in case of antiferromagnetic fluctuations in the Hubbard model 55.

Fig. 3 shows the local spectral function $A_{loc}$ obtained from the $\mathbf{k}$-integrated spectral function $A_{k}=-\frac{1}{\pi}\text{Im}G_{\mathbf{k}}(\omega)$ for real frequencies $\omega$. The analytical continuation was done with the maximum entropy method (MEM)56; 57. In accordance with the self-energy calculations, $A_{loc}(\omega)$ is shown for the same parameters as in Fig. 3. For $4V/U=1.2t$ at $U=2t$ the Padé interpolation58 is also shown as it results in a slightly different spectrum.

Independent of the values of $V$, at $U=4t$ $A_{loc}(\omega)$ consists of three contributions stemming from the dominant quasiparticle peak at zero frequency and the upper/lower Hubbard bands at $\omega\sim\pm U/2$ indicating the overall metallic nature of the solutions. In case of $U=2t$, the peak positions of the Hubbard bands are slightly displaced at $\omega\sim\pm 2t$. Increasing $V$ suppresses and broadens the quasiparticle peak. This suppression can be assigned to the increase of charge fluctuations (cf. Fig. 4 below). The effect is however small for the case of $4V/U\lesssim 1$.

For a better comparison, the Fermi liquid parameters, namely the renormalization factor $Z_{\mathbf{k}}=\left[1-\left.\text{Im}\Sigma_{k}/\nu_{n}\right|_{\nu_{n}\to 0}\right]^{-1}$ and the scattering rate $\Gamma_{\mathbf{k}}=\left|\text{Im}\Sigma_{k}\right|_{\nu_{n}\to 0}$, are extracted from the self-energy and are shown for $\mathbf{k}=\left(\frac{\pi}{3},\frac{2\pi}{3}\right)$ in the upper and lower insets of Fig. 3, respectively. To do this, the imaginary part of self-energy was fitted with a second order polynomial for the first three Matsubara frequencies and the error estimate was obtained from the comparison to a fit with a third order polynomial. The quasiparticle renormalization $Z_{\mathbf{k}=\left(\frac{\pi}{3},\frac{2\pi}{3}\right)}$ slightly grows with increasing $V$ up to $4V/U\sim 1$ above which the growth becomes much sharper. Similarly, the scattering rate $\Gamma_{\mathbf{k}=\left(\frac{\pi}{3},\frac{2\pi}{3}\right)}$ slightly decreases with increasing $V$ until $4V/U\sim 1$. Afterwards, it substantially grows. The particular behavior of the Fermi liquid parameters is another evidence of the charge fluctuations that only appear in the close vicinity of the CDW phase transition. For all other values of $V$ with $4V/U\lesssim 1$, the CDW fluctuations have very little impact on the one-particle properties.

For $4V/U=1.2$ and $U=2t$, the spectral function is suppressed at small frequencies. For these parameters the analytic continuation by means of Padé interpolation resulted in a dip at the Fermi energy, cf. Fig. 3, indicating a tendency to opening of a gap. Despite the different solutions from the MEM and the Padé for this parameter, both qualitatively agree with the expectation of the charge fluctuations effect, and they are consistent with the large negative $\text{Im}\Sigma$ at the Fermi surface (see Fig. 3). We believe that with the further reduction of temperature a full gap can be obtained which smoothly evolves into the CDW gap in the ordered phase.

III.2 Two-particle susceptibilities and eigenvalues

A more direct measure of the fluctuations is the susceptibility in various channels. As the parquet approach determines the single- and two-particle quantities simultaneously, the signals of the CDW observed in the self-energy (Fig. 3) and the spectral function (Fig. 3) will also leave their fingerprints in the two-particle vertex functions. Actually, as it will be seen soon the characterization of the phase transition can be more conveniently done in the two-particle level, as the two-particle susceptibilities are the fluctuations of the corresponding single-particle quantities. These fluctuations become more significant in the vicinity of the phase transition, and are not directly accessible in any single-particle self-consistent theory. In the parquet approach, the susceptibilities in the magnetic, density (charge) and particle-particle (including both singlet $s$ and triplet $t$) channels are given as

[TABLE]

with $F_{p\overline{\uparrow\downarrow}}^{kk^{\prime}q}=\frac{1}{2}[F_{s}-F_{t}]^{kk^{\prime}q}$. The divergence of the susceptibility in a given channel signals the phase transition towards the breaking of the corresponding symmetry of this channel. Similarly, Eqs. (9) and (10) also indicate that the phase transition can be equivalently characterized by the leading eigenvalue $\lambda$ of a corresponding Bethe-Salpeter equation which reaches the value of $1$. Therefore the eigenvalues themselves constitute another measure of the fluctuation strength in a given channel close to a phase transition. The corresponding eigenvalue equations read

[TABLE]

where $\Gamma$ denotes the irreducible vertex in a given channel and $\phi$ is the eigenvector. Since the Bethe-Salpeter equations are diagonal in transfer momentum $\mathbf{q}$ and frequency $\omega_{n}$, the above eigenvalues $\lambda^{q}_{d/m/pp}$ also inherit this dependence. In the following we discuss eigenvalues only for the zeroth bosonic frequency ($\omega_{n}=0$), $\mathbf{q}=(\pi,\pi)$ for $\lambda_{d/m}$, and $\mathbf{q}=(0,0)$ for $\lambda_{pp}$, where the eigenvalues in the respective channels are dominant.

In Fig. 4 the eigenvalues (top) and (inverse) susceptibilities (bottom) are shown as functions of $V$. Again results for two values of the local interaction are shown: $U=4t$ (solid lines) and $U=2t$ (dashed lines). In the overview phase diagram Fig. 1, the ($U$, $V$) points of Fig. 4 are marked as upper and lower triangles.

Changes are visible in the charge susceptibility \raisebox{2.0pt}{\chi}_{d} which increases when $V$ is enhanced, but not surpassing the magnetic susceptibility \raisebox{2.0pt}{\chi}_{m}. The inverse of \raisebox{2.0pt}{\chi}_{m} is close to zero and likewise, its corresponding eigenvalue is close to one, indicating very strong antiferromagnetic fluctuations ($\mathbf{q}=(\pi,\pi)$) in this parameter regime. For the smaller value of $U$ both AFM and CDW fluctuations are significantly weakened with the AFM fluctuations being still dominant in the system.

It is only at a rather large value of $V$, $4V/U>0.8$ that the density eigenvalue starts to rapidly increase (see the top panel of Fig. 4). At the same time the magnetic eigenvalue slightly decreases with increasing $V$. When $4V/U$ becomes slightly greater than $1$, the magnetic eigenvalue is eventually surpassed by the density eigenvalue. Only then the CDW fluctuations become more important, which can also be seen in the charge susceptibility that becomes bigger than the magnetic one when $4V/U\sim 1$ (see the bottom panel of Fig. 4). For the smaller local interaction, $U=2t$, the increase of the charge susceptibility and eigenvalue is less steep as compared to the case of $U=4t$ (cf. Fig. 4). As a result, the dominance of CDW fluctuations occurs at a slightly larger $V$ value. Nevertheless, for both cases ($U/t=2$ and $U/t=4$) the fluctuations in the charge channel become stronger than in the magnetic channel but only very close to the transition boundary. The leading eigenvalue and the near-divergence of the charge susceptibility occur at $\mathbf{q}=(\pi,\pi)$ which corresponds to a new periodicity in real-space of $\sqrt{2}\times\sqrt{2}$. This is in agreement with the strong coupling expansion for this model23; 24 and was also seen in forth-order perturbation theory.23

The $pp$ eigenvalue and the corresponding susceptibility (not shown here) remain small for all values of $V$, which leads to a conclusion that, for the parameters considered, pairing fluctuations do not play any significant role. Note, however, that for $U=4t$ the $pp$ eigenvalue is almost twice as big as that for $U=2t$, i.e. a relatively strong local interaction is important for pairing fluctuations, as it is also for the AFM and CDW ones (cf. top panel of Fig. 4).

Although it is not the prime aim of this paper to precisely determine the transition boundary of the EHM, the above analysis on the leading eigenvalue and the susceptibility provide an efficient way to estimate the critical value of $V$. The value of $\lambda_{d}$ at $\mathbf{q}=(\pi,\pi)$ reaches almost $1$ for $4V/U\to 1.06$ at $U=4t$ and for $4V/U\to 1.2$ at $U=2t$, respectively. A more reliable estimate of the phase transition point within the parquet method is very difficult. It would require calculations for different and large cluster sizes to properly resolve the deviations from mean-field critical exponents. Further, we are restricted to second order phase transitions because the parquet equations are only solved in the paramagnetic phase. If there is a first order transition, as might well be the case in some parameter range 59; 39, the thus calculated $V$ values are too large.

We note that very similar critical values of $4V/U$ for the CDW phase transition were obtained in many different works 34; 38; 30; 33; 32. DQMC predicted a phase transition in the regime of $4V/U\in\{1,1.25\}$ for $U=4t$ and $4V/U\in\{1,1.4\}$ for $U=2t$ 32. In the DCA calculations 38, they are $4V/U=1.04$ for $U=4t$ and $4V/U=1.216$ for $U=2t$ . At a lower temperature ($T=t/12.5$), the dual boson approach 34 identifies the phase transition at $4V/U=1.08$ for $U=4t$ and $4V/U=1.04$ for $U=2t$. Our estimation based on the parquet approximation is in agreement with these studies, confirming the feasibility of the parquet approach in the study of phase transitions.

III.3 Full two-particle vertex function

In addition to the susceptibility and the leading eigenvalue, the competing fluctuations can also be seen in the full two-particle vertex function $F^{kk^{\prime}q}$. In Fig. 5 we show the fermionic frequency dependence (i.e. as a function of $\nu_{n}$ and $\nu^{\prime}_{n}$) of the full vertex in the density, magnetic and particle-particle channels (different columns) for different values of $V$ (different rows) at $U=4t$. The momenta are set to $\mathbf{k}=\mathbf{k^{\prime}}=(0,0)$ (please note that the dependence on the fermionic momenta ${\bf k}$ and ${\bf k^{\prime}}$ is very weak), $\mathbf{q}=(\pi,\pi)$ and the bosonic frequency $\omega_{n}=0$.

Increasing the non-local interaction $V$ affects mostly the density channel. Comparing however the three rows of Fig. 5 we can conclude that the onset of strong charge fluctuations is visible only shortly before the phase transition, i.e. at $4V/U\sim 1$. For $V=0$, the full vertex function is mainly described by the main diagonal with $\nu_{n}=\nu^{\prime}_{n}$. As $V$ is increased to $4V/U=0.8$, the overall amplitude of the full vertex function increases and becomes negative except for the main diagonal components. This background, however, is nicely approximated by the first order contribution Eq. (2), $2V_{\mathbf{q}=(\pi,\pi)}-V_{\mathbf{k}^{\prime}-\mathbf{k}=(0,0)}=-12V$, which for $4V/U=0.8$ is equal to $-9.6$ in our units. When the charge fluctuations become stronger ($4V/U=1$, bottom row of Fig. 5), these contributions significantly increase far beyond the first order approximation. In the magnetic and the particle-particle channel (cf. the second and the third column of Fig. 5) the structure of $F$ stays nearly unchanged compared to the density channel. At $4V/U=1$ we see however a strong broadening of the plus-structure, e.g. along $\nu_{0}=\frac{\pi}{\beta}$ and $\nu^{\prime}_{0}=\frac{\pi}{\beta}$ in the magnetic vertex $F_{m}$. This indicates a suppression of magnetic contributions due to the influence of the now very strong CDW fluctuations. In this regime, likewise, the magnetic eigenvalue and magnetic susceptibility decrease (see Fig. 4). Analyzing the momentum dependence of the full vertex function (not shown here) leads to the same conclusion. Predominant effects of the charge fluctuations are seen in $F_{d}$ only close to the phase transition for $4V/U\gtrsim 1$. If the phase transition was of first-order and hence at a smaller $V$, AFM fluctuations might even prevail throughout the entire paramagnetic phase up to the CDW phase transition line. Turning back to the magnetic contributions: we observe a global suppression for $4V/U=1$, whereas for smaller non-local interactions this effect is very small. The contributions of the particle-particle channel stay almost unchanged with the increase of the non-local interaction.

IV Response to electric field

With the two-particle vertex function available in the parquet approach, we are also able to study the response to external perturbations easily. Like in the two-particle susceptibility and the density-density vertex function, the non-local charge fluctuations are also encoded in response functions which are experimentally accessible. As an example, we apply an electric field along the $x$-direction to the EHM, i.e. $\mathbf{E}=E\,\mathbf{a}_{x}$, and study the response to it at different value of $V$. The response function (or paramagnetic current-current correlation function) within the Peierls approximation reads

[TABLE]

for a coupling along the field direction, $x$, with

[TABLE]

Because of the long optical wavelength, the variation of electric field in space can be neglected. Thus, it is sufficient to restrict the calculations to $\mathbf{q}=0$. The optical conductivity which is obtained generally from linear response theory as

[TABLE]

Following Eq. (15), one can first transform the current-current correlation function from Matsubara to real frequency and then divide it with $\omega+i\delta$ to get the optical conductivity. Alternatively, one can directly transform the \text{Im}\raisebox{2.0pt}{\chi}_{jj,\mathbf{q}=0}(i\Omega_{m})/i\Omega_{m} to real frequency with MEM, which was the strategy taken in our calculations 57.

The response fuctions to the applied electric field (in the direction of the lattice vector $\mathbf{a}_{x}$) are displayed on the real-frequency axis in Fig. 6. The response functions \raisebox{2.0pt}{\chi}_{P,\mathbf{q}=0} and the optical conductivity $\sigma_{P}$, are displayed for various values of $V$ at $U/t=4$: from $4V/U\lesssim 1$ (weak $V$-limit) to $4V/U>1$ (strong $V$-limit) in which either antiferromagnetic or charge fluctuations prevail (cf. Sec. (III.2)). Additionally, the bubble and vertex contributions, defined in Eq. (13), are displayed in Fig. 6 by light green dashed/dotted lines and red dotted lines, respectively.

The correlation function exhibits a complex two peak structure with the maximum peak value increasing with $V$ up to $4V/U\sim 1$ and decreasing afterwards. However, there is only one peak observed in the optical conductivity. The small high-frequency feature observed in $\mbox{Im}\chi_{P,q=0}$ is suppressed by the large frequency denominator leaving only the Drude-like peak visible in Fig. 6. We expect a more visible high-frequency peak when the temperature becomes lower or $U$ is increased so that pronounced and well separated Hubbard bands start to form.

We notice that the zero-frequency peak in the optical conductivity increases when increasing $V$ from $V=0$ up to $4V/U=1$. Both the bubble and vertex contributions become larger. In particular, the vertex contribution at $\omega=0$ changes its sign and continuously increases with the increase of $V$, indicating the correlation effect is always enhanced with increasing $V$. On the other hand, at $4V/U=1.06$ the zero-frequency peak in the optical conductivity clearly drops down, which mainly results from the decrease of the bubble contributions. The vertex part, however, still strongly contributes. The bubble susceptibility in turn contributes less and less due to the suppression of the single-particle density of states at the Fermi level, see Fig. 3. As already discussed in the context of the self-energy and one-particle spectrum, this is because of the emergence of long-range CDW fluctuations at $4V/U>1$ as is also reflected in the dramatic increase of the CDW susceptibility in Fig. 4. For the optical conductivity calculations, we observe the same effect as discovered in the self-energy and two-particle vertex functions: only in the close vicinity of the phase transition (around $4V/U\sim 1.06$ at $U=4t$), the effect of the charge fluctuations becomes dominant. Away from the phase transition, for this half-filled two-dimensional square lattice, spin fluctuations are always dominant.

V Conclusion

By extending the parquet approach to include the non-local Coulomb interactions, we studied the extended Hubbard model on a two-dimensional square lattice using the parquet approximation. Through the calculations of the single-particle self-energy (Fig. 3), spectral function (Fig. 3) and the two-particle susceptibility (Fig. 4), the vertex function (Fig. 5), as well as the optical conductivity (Fig. 6), we are able to consistently identify the strong competition of the non-local charge fluctuations and the spin fluctuations. Quite surprisingly, the spin fluctuations are dominating over a wide range of parameters and only in the immediate vicinity of the phase transition to the CDW order charge fluctuations prevail. In case of a first order transition 59; 39, AFM fluctuations might even dominate in the entire paramagnetic phase up to the CDW phase.

In the narrow regime, where the charge susceptibility surpasses the antiferromagnetic susceptibility, the self-energy in the paramagnetic phase turns from a metallic to an insulating-like behavior — the vanguard of the CDW insulator. Likewise the optical conductivity is suppressed by the strong CDW fluctuations. The parquet approach alows us to look at this competition and suppression also on a deeper level, unseen directly in the measured quantities. The otherwise very sharp frequency structure of the two-particle magnetic vertex becomes broadened close to the transition and the magnetic fluctuations are suppressed by the very large vertex in the charge channel. For the pairing (particle-particle) fluctuations, on the other hand, this has very little influence.

Due to the fact that the self-consistency at both, one- and two-particle, levels is simultaneously satisfied in the parquet approach, it is not surprising to obtain a consistent picture for all the above-mentioned quantities. Although it is not the prime interest of this work to locate the precise transition boundary, our rough estimation agrees with the published results. Our work provides a valuable explanation of how the single- and two-particle quantities are internally related. The transparent structure of the parquet equation provides a predicting power in the sense that the change in one quantity can be naturally evolved to interpret the change of other quantities.

VI acknowledgment

We would like to thank Josef Kaufmann for the help with the analytical continuation. P.P. has been supported by the Austrian Science Fund (FWF) through the Doctoral School “Building Solids for Function”, A.K. and K.H. through FWF project P 30997 and P 32044. G.L. acknowledges the starting grant of ShanghaiTech University, the Program for Professor of Special Appointment (Shanghai Eastern Scholar) and the support from the National Natural Science Foundation of China (Grant No. 11874263). Calculations have been done in part on the Vienna Scientific Cluster (VSC).