Diagrammatic Monte Carlo Method for Many-Polaron Problems

We introduce the first bold diagrammatic Monte Carlo approach to deal with polaron problems at a finite electron density nonperturbatively, i.e., by including vertex corrections to high orders. Using the Holstein model on a square lattice as a prototypical example, we demonstrate that our method is capable of providing accurate results in the thermodynamic limit in all regimes from a renormalized Fermi liquid to a single polaron, across the nonadiabatic region where Fermi and Debye energies are of the same order of magnitude. By accounting for vertex corrections, the accuracy of the theoretical description is increased by orders of magnitude relative to the lowest-order self-consistent Born approximation employed in most studies. We also find that for the electron-phonon coupling typical for real materials, the quasiparticle effective mass increases and the quasiparticle residue decreases with increasing the electron density at constant electron-phonon coupling strength.

Figure 1

Green’s function of a single polaron with zero momentum at $\lambda =0.45$, $\mu /t=-4.5$ and zero temperature obtained by the BDMC technique (blue dots). It is compared to the asymptotic behavior $G_{\tau \to \infty }(k=0)=Z_s\exp [-E_s\tau ]$ (dashed line) with $Z_s=0.826$ and $E_s=-0.302$ calculated by the conventional single-polaron diagrammatic MC approach [9]. All error bars are smaller than symbol sizes. Inset: Ratio of the polaron self-energies $\mathrm{\Sigma }(\mathbf{k},\tau )$ at $\mathbf{k}=(\pi ,0)$ and $\mathbf{k}=0$ as a function of imaginary time $\tau$.

Figure 2

Electron density computed from diagrams up to order $N$ without (open blue squares) and with (filled black circles) the self-consistent renormalization of phonons at $\lambda =0.45$, $\gamma =1.008$, and $T=0.01$. We subtracted the static local polarization of the ideal electron gas at the same density from the phonon self-energy to mimic the standard material science protocol [6] where this contribution is already accounted for in the value of ${\omega }_0$. Inset: Momentum dependence of the static self-energy $\mathrm{\Sigma }(\mathbf{k},0)$ along the ${\widehat{k}}_x$ axis.

Figure 3

Momentum distributions for various electron densities per site (from right to left: $n=0.6498$, 0.4055, 0.3026, 0.2086, 0.1222, 0.08169, 0.04323, 0.02846, 0.01413, 0.001416) at $\lambda =0.45$. For $n>0.1$ we consider $T/t=0.025$ (thin lines); for $n<0.1$ the temperature is reduced down to $T/t=0.01$ (thick lines). The dashed line corresponds to the ideal gas case at $T/t=0.025$ and $\mu /t=-1$. With dash-dotted lines we show (approximately) the locations of the $T=0$ jumps in the distribution function.

Figure 4

Effective mass along $\widehat{x}$ (black open squares) and ($\widehat{x}$, $\widehat{y}$) (black filled squares) directions and the quasiparticle residue (blue circles) as functions of the adiabatic parameter at $\lambda =0.45$ deduced from the same set of simulations as in Fig. 3. The upper horizontal axis provides an approximate density scale. The error bars for $Z$ at large ${\gamma }^{-1}$ are not statistical; they indicate the anisotropic spread of $Z$ values on the FS. The anisotropy of the effective mass is significant only for the largest density, while for all other points it is unmeasurably small.