Collective charge excitations and the metal-insulator transition in the square lattice Hubbard-Coulomb model

In this article, we discuss the nontrivial collective charge excitations (plasmons) of the extended square lattice Hubbard model. Using a fully nonperturbative approach, we employ the hybrid Monte Carlo algorithm to simulate the system at half-filling. A modified Backus-Gilbert method is introduced to obtain the spectral functions via numerical analytic continuation. We directly compute the single-particle density of states which demonstrates the formation of Hubbard bands in the strongly correlated phase. The momentum-resolved charge susceptibility also is computed on the basis of the Euclidean charge-density-density correlator. In agreement with previous extended dynamical mean-field theory studies, we find that, at high strength of the electron-electron interaction, the plasmon dispersion develops two branches.

Figure 1

The effect of regularization on the spectral functions. We demonstrate the solution of (19) for the DOS using standard Tikhonov regularization. Here we plot the estimator for the spectral function obtained according to (30). The frequency $\omega$ corresponds to the center of the resolution function, and the filled area shows the statistical error. The example data are taken for the interaction strength $U/\kappa =1.66,V/\kappa =0.62$, and temperature $T/\kappa =0.046$. The lattice with spatial size $N_s=20$ and $N_{\tau }=160$ Euclidean time slices is used. Resolution functions for $\lambda ={10}^{-7}$ can be found in Fig. 2.

Figure 2

The resolution functions $\delta ({\omega }_0,\omega )$, obtained with standard Tikhonov regularization [see (41)] during the solution of the Green-Kubo relation for DOS (19). The parameter ${\omega }_0$ labels the center of the resolution function. The value of the regularization parameter is $\lambda ={10}^{-7}$. This setup will be used in all cases where we compute the DOS from the Monte Carlo data.

Figure 3

The half-peak width of the resolution function centered around ${\omega }_0=0$ versus the statistical error. The calculation is performed with different regularization algorithms for the correlator $G\left(\tau \right)$ calculated on the lattice with spatial sizes $N_s=20$ and $N_{\tau }=160$ for the interaction strength $U/\kappa =1.66,V/\kappa =0.62$, and temperature $T/\kappa =0.046$. The statistical error is controlled by the regularization parameter $\lambda$. The resolution functions are computed in the process of inverting the Green-Kubo relation for the DOS (19). The width is given in units of temperature, and one can see that, as $\lambda$ vanishes, each curve converges to $\sim 2$. This is as expected as the resolution in frequency space ultimately is limited by the temperature.

Figure 4

Two variants of correlator averaging. The charge-density-density correlator $C(\mathbf{q},\tau )$ is calculated for $U/\kappa =3.33,V/\kappa =1.26$, and $T/\kappa =0.046$ using a lattice with spatial size $N_s=20$ and $N_{\tau }=160$ Euclidean time slices. The momentum $q$ belongs to the $X\text{−}M$ line in the Brillouin zone, which is the most physically interesting case (see below in Sec. 5).

Figure 5

Spectral functions $\mathrm{Im}{\chi }_{\rho }(\omega ,q)$ computed using the relations (28) and (29) for two different variants of the averaged correlator shown in Fig. 4. In both cases the ordinary Tikhonov regularization is used with $\lambda =5.0\times {10}^{-6}$. We plot the estimator for the spectral function obtained according to (30). The frequency $\omega$ corresponds to the center of the resolution function, and the filled area shows the statistical error.

Figure 6

Plot of the DOS for the lattice ensemble with electron-electron interaction (a) $U/\kappa =0.83,V/\kappa =0.5$, (b) $U/\kappa =1.66,V/\kappa =0.62$, and (c) $U/\kappa =3.33,V/\kappa =1.26$. In all cases we have used a spatial lattice size of $N_s=20$ with $N_{\tau }=160$ steps in Euclidean time and a temperature of $T/\kappa =0.046$. The standard Tikhonov regularization with $\lambda ={10}^{-7}$ is employed during the analytical continuation. Thus, in all cases, we have the same resolution in frequency. No additional averaging in Euclidean time is applied. The filled areas show the statistical error computed with the data binning procedure, and the frequency $\omega$ corresponds to the center of the resolution function. The average value represents the estimator for the DOS computed according to (30). One can see the formation of the Hubbard bands, characterized by the peak at $E/\kappa \approx 0.5$, indicating that one is in the strongly coupled phase for the strongest interaction strength. For weaker interaction strengths we obtain a strong peak at zero energy indicating that the system is in the metallic phase.

Figure 7

Plot of the DOS for the lattice ensemble with electron-electron interaction $U/\kappa =4.4,V/\kappa =1.26$, which corresponds to the weakest interaction strength studied in Ref. [13] ($U^{*}=1.1$ in their notation). We have used a spatial lattice size of $N_s=20$ with $N_{\tau }=160$ steps in Euclidean time and two times higher temperature $T/\kappa =0.092$, which is almost equal to the one used in Ref. [13]. The setup for the analytic continuation is exactly the same as for Fig. 6. One can see that, in contrast to the results displayed in Ref. [13], the Hubbard bands are formed already, which suggests that the metal-insulator phase transition is shifted to a weaker interaction strength even at the same temperature.

Figure 8

Charge susceptibility for the square lattice Hubbard-Coulomb model for three different interaction strengths. The color scale corresponds to $V\left(q\right)\mathrm{Im}{\chi }_{\rho }(\omega ,q)$. For the spectral function $\mathrm{Im}{\chi }_{\rho }(\omega ,q)$ we use the estimated value after analytic continuation [see (30)]. The same lattice ensembles were used as in Fig. 6. The standard Tikhonov regularization with $\lambda =5\times {10}^{-6}$ and additional time averaging according to Fig. 4 is used in all cases.

Figure 9

Comparison of the spectral function for the charge-charge correlator at the $X$ point for different interaction strengths. The lattice setup and the parameters of the analytic continuation procedure are identical to the ones used for Fig. 8. The filled areas show the statistical error computed with the data binning procedure, and the frequency $\omega$ corresponds to the center of the resolution function. The spectral functions for weaker interaction strengths are rescaled by factors of 2 and 4 in order to fit to the same scale.