Solving the Bethe-Salpeter Equation in Real Frequencies at Finite Temperature

The Bethe-Salpeter equation (BSE) was initially introduced to study two-particle Green’s function and the corresponding bound states in nuclear physics [1]. Later on it was formulated for studies of optical spectra in solids [2, 3] and molecules [4, 5]. In material science/quantum chemistry methods one often starts by establishing the Hartree-Fock (HF) basis and proceeds with the electronic structure calculations and solutions of the BSE for two-particle Green’s function and related properties [6, 7, 8]. In strongly correlated models and quantum field theories BSE is used to calculate the effect of vertex corrections [9, 10]. This is especially important in the case of polarization when looking at charge and spin responses because neglecting vertex corrections leads to wrong results [11].

The standard finite-temperature approach is to first solve the problem in the Matsubara domain [12] and then apply a numeric analytic continuation procedure to obtain real frequency results for comparison with experimental observations. However, even if original data is known with high accuracy, the final step distorts, or even fails to resolve, important spectral features unless one imposes constrains on how the result is supposed to behave [13]. Until recently, this problem standing in the way of accurate theoretical descriptions of experimentally relevant observables was considered unavoidable.

The breakthrough development done in the context of diagrammatic Monte Carlo (diagMC) for Fermi-Hubbard model in Refs. [14, 15, 16] (see also Refs. [17, 18]) was that one can automatically express real-frequency answers for an arbitrary diagram using the so-called algorithmic Matsubara integration (AMI) protocol and, thus, completely eliminate the need for numeric analytic continuation. It works best for expansions in terms of “bare” (non-interacting) Green’s functions and frequency-independent interactions. However, for strongly correlated regimes, when diagrammatic expansions in terms of original potentials diverge—a typical example would be the fundamental homogeneous electron gas (HEG), or jellium, model—expansions should be performed in terms of “dressed” propagators and screened potentials, which include some interactions effects non-perturbatively. The corresponding generalization of the AMI protocol was described in Ref. [11] but its implementation for high-order diagrammatic simulations poses significant technical challenges.

Current applications of the AMI protocol for jellium are based on expansion in terms of the Yukawa potential, $Y=4\pi e^{2}/(q^{2}+\kappa^{2})$, with optimally chosen screening momentum $\kappa$ and self-consistent HF single particle propagators, $G$, see Refs. [19, 20, 21] and the text below Eq. (1). In brief, the AMI technique lists all diagrams of order $N$ for an observable of interest, eliminates all sums over internal Matsubara frequencies for every listed diagram by an exact analytic transformation, and stochastically samples the remaining momentum integrals. However, it remains biased because all poles are regularized by adding small but finite imaginary part to the external frequency, $\Omega\to\Omega+i\eta$. The exact answer is recovered by taking the $\eta\to+0$ limit which requires extremely long computation times for high-order diagrams [19].

In this work, we apply the AMI and series resummation techniques for solving the ladder-type finite-temperature BSE equation for the $3$-point vertex function $\Gamma^{(3)}$ convoluted with two Green’s functions to obtain system’s polarization directly on the real-frequency axis. We take the crucial new step of implementing the limit of $\eta\to+0$ explicitly; the same procedure may be used for any model with instantaneous interaction. The resulting scheme is not only numerically exact, it is also very efficient and allows one to reach expansion orders high-enough for performing reliable resummations of divergent series. The same expansion orders cannot be reached if simulations are performed with appropriately small finite-$\eta$ regularization.

For illustration, we compute the polarization of the HEG in the HF basis, and compare it to the random-phase-approximation (RPA) predictions (also in the HF basis) to see the effects of multiple particle-hole scattering processes. In particular, we find that the Landau damping coefficient is significantly increased relative to the RPA prediction.

Hamiltonian of the HEG is defined by

$$H=\sum_{i}\frac{k_{i}^{2}}{2m}+\sum_{i<j}\frac{e^{2}}{|\mathbf{r}_{i}-\mathbf{r}_{j}|}-\mu N,$$

(1)

with $m$ the electron mass, $\mu$ the chemical potential, and $e$ the electron charge. We use the inverse Fermi momentum, $1/k_{F}$, and Fermi energy, $\varepsilon_{F}=k_{F}^{2}/2m$, as units of length and energy, respectively; the definition of the Coulomb parameter $r_{s}$ in terms of the particle number density, $\rho=k_{F}^{3}/3\pi^{2}$, and Bohr radius, $a_{B}=1/me^{2}$, is standard: $4\pi r_{s}^{3}/3=1/\rho a_{B}^{3}$.

For a meaningful diagrammatic expansion one has to work with screened interactions and counter-terms. In this work we follow the setup formulated in Ref. [20] and expand in the Yukawa potential on top of the self-consistent HF solution for the Green’s function: $G^{-1}=G^{-1}_{0}-\Sigma_{HF}[G]$, where $G_{0}$ is the bare Green’s function and $\Sigma_{HF}$ is the Fock proper self-energy shown in Fig. 1(a) (Hartree diagram is absent by charge neutrality). In the converged HF solution, $G=(i\omega_{s}-k^{2}/2m-\Sigma_{HF}(\mathbf{k})+\mu)^{-1}\equiv(i\omega_{s}-\epsilon_{k}+\mu)^{-1}$, the chemical potential is also self-consistently adjusted to keep the electron density fixed, see upper part of Fig. 1. By incorporating Fock diagrams into $G$ one simplifies the expansion by omitting all diagrams with Fock-type self-energy insertions.

Formally, not only the screening momentum $\kappa$ in the Yukawa potential, but the entire potential itself is arbitrary as long as it does not feature divergent behavior at $q=0$. Any difference between the Coulomb, $V(q)=4\pi e^{2}/q^{2}$, and expansion potentials is taken care of by a counter-term $C(q)=1/V(q)-1/Y(q)$ which keeps the field theoretical representation of the model exact [22, 23, 20]. For example, the physical screened interaction is given by $1/W=1/Y-[\Pi-C(q)]\equiv 1/V-\Pi$ where $\Pi$ is the system polarization. Ultimately, final results are supposed to show little dependence on the choice of $Y(q)$. Flexibility in reformulating the diagrammatic expansion can, and should be, used for optimization of series convergence properties. This is precisely how recent work reported in Refs. [20, 21] achieved unprecedented accuracy for single-particle and static linear response properties of HEG.

In this work we apply the above field-theoretical setup to explore how the polarization is modified by multiple scattering of the particle-hole pair at finite temperature. Since the calculation is performed directly on the real frequency axis, our results can be immediately used to obtain the Landau damping parameter, see Ref. [24], or describe the so-call charge loss function, $\mathrm{Im}W$. The set of diagrams involved is presented in Figs. 1 (b), (c), and (d) where the ladder-diagram series for polarization are first expressed in terms of the BSE for the dressed three-point vertex function $\Gamma^{(3)}$, see 1 (c), which is then convoluted with two Green’s functions, see 1 (b). As noticed in Ref. [11], if the calculation is restricted to have no more than one ladder rung, the result for $\Pi=\Pi_{0}+\Pi_{1}$, where $\Pi_{i}$ is the contribution of the diagram with $i$ rungs Fig. 1(d), may violate causality. Also, an observation that $\Pi_{1}$ is not small compared to $\Pi_{0}$ raises the question of importance of contributions from high-order ladder diagrams. Summation of the entire ladder series is supposed to solve both problems.

Before proceeding, let us elaborate on the importance of considering very small values of $\eta$ when using it for regularization of poles. Even the lowest order diagram demonstrates significant finite-$\eta$ correction at low temperature, see Fig. 2. At $T/\varepsilon_{F}=0.02$ the finite-$\eta$ effect in $\Pi_{0}$ is as large as $20\%$ for $\eta/\varepsilon_{F}=0.005$. One may attempt to extrapolate results towards $\eta=0$, but this procedure is not guaranteed to work for all higher-order contributions because the $\eta$-dependence is not necessarily monotonic.

The $i$-th order contribution to the ladder diagram series for $\Pi$, see Figs. 1(b) and 1(c), is given by (in what follows we treat $\Omega$ and ${\bf Q}$ as external parameters and do not mention them explicitly as function variables)

$$\Pi_{i}=-2\left(\prod_{j=1}^{i+1}\int\frac{d{\bf p}_{j}}{(2\pi)^{3}}\right)\prod_{j=1}^{i+1}{\cal F}_{{\bf p}_{j}}\prod_{j=1}^{i}Y({\bf p}_{j+1}-{\bf p}_{j})\,,$$

(2)

where

$${\cal F}_{{\bf p}}=\frac{n_{{\bf p}+{\bf Q}}-n_{\mathbf{p}}}{\Omega-\epsilon_{{\bf p}+{\bf Q}}+\epsilon_{\mathbf{p}}+i\eta},$$

(3)

with $n_{\mathbf{p}}=(e^{\epsilon_{p}-\mu}+1)^{-1}$. While this setup is immediately suitable for Monte Carlo simulations with small finite $\eta$, it quickly becomes too expensive as $i$ increases. The reason is the severe sign problem coming from Monte Carlo sampling of the principle part integrals around poles of $F$ functions.

This problem can be eliminated by taking the limit of $\eta\to+0$ explicitly with the help of the $\lim_{\eta\to+0}\frac{1}{x+i\eta}=-i\pi\delta(x)+{\cal P}\frac{1}{x}$ formula were ${\cal P}$ stands for the principal part value. In practice it means that integrals featuring simple poles can be transformed identically as

$$\int_{a}^{b}dx\frac{h(x)}{\Omega-g(x)+i\eta}=\int_{a}^{b}dx\left[\alpha(x,x_{0})+i\gamma(x,x_{0})\right]\,,$$

(4)

where $x_{0}$ is the solution of the $g(x_{0})=\Omega$ equation and

	$\displaystyle\gamma(x,x_{0})$	$\displaystyle=$	$\displaystyle-\pi\frac{h(x_{0})}{|g^{\prime}(x_{0})|}N(x,x_{0})\;\;\mbox{if}\;x_{0}\in(a,b),$
	$\displaystyle\gamma(x,x_{0})$	$\displaystyle=$	$\displaystyle 0\qquad\qquad\qquad\qquad\;\;\;\mbox{otherwise},$		(5)

	$\displaystyle\alpha(x,x_{0})$	$\displaystyle=$	$\displaystyle\frac{1}{2}\left[\frac{h(x)}{\Omega-g(x)}+\frac{h(x_{r})}{\Omega-g(x_{r})}\right]\;\;\mbox{if}\;(x_{0},x_{r})\in(a,b)$
	$\displaystyle\alpha(x,x_{0})$	$\displaystyle=$	$\displaystyle\frac{h(x)}{\Omega-g(x)}\qquad\qquad\qquad\qquad\;\;\mbox{otherwise},$		(6)

with $x_{r}=2x_{0}-x$. Here $N(x,x_{0})$ is an arbitrary function normalized to unity on the $(a,b)$ interval–it may explicitly depend on $x_{0}$. By construction, the $\alpha+i\gamma$ function is not singular at the pole location. In our case, we choose variables $x=\cos\widehat{{\bf p}{\bf Q}}$ to regularize poles. For a given spherically symmetric dispersion relation $\epsilon_{\mathbf{p}}$ the location of the pole $x_{0}(p,Q,\Omega)$ as a function of $p$ is tabulated on a fine grid at the start of the simulation and subsequently is treated as a known function.

However, taking the limit $\eta\to+0$ also involves “splitting” the contribution of the diagram with variable $x$ into up to three contributions from variables $x$, $x_{0}$, and $x_{r}$ if $(x_{0},x_{r})\in(a,b)$. As the ladder diagram order increases, one now faces the problem of computing up to $3^{i+1}$ contributions for some multi-dimensional momentum points $\{{\bf p}_{j}\}$ because momentum variables are linked by potentials $Y$, see (2). This remaining problem is solved by “decoupling” momentum integrals using Fourier transforms of Yukawa potentials

$$Y({\bf p}_{j+1}-{\bf p}_{j})=\int d{\bf r}_{j}Y(r_{j})\,e^{i({\bf p}_{j+1}-{\bf p}_{j})\cdot{\bf r}_{j}}\,,$$

(7)

where $Y(r)=(e^{2}/r)\,e^{-\kappa r}$. Now

$$\Pi_{i}=-2\left(\prod_{j=1}^{i}\int d{\bf r}_{j}Y(r_{j})\right)\prod_{j=1}^{i+1}\int\frac{d{\bf p}_{j}}{(2\pi)^{3}}\tilde{\cal F}_{{\bf p}_{j},{\bf r}_{j-1}-{\bf r}_{j}}$$

(8)

with

$$\tilde{\cal F}_{{\bf p},{\bf R}}=e^{i{\bf R}\cdot{\bf p}}\frac{n_{{\bf p}+{\bf Q}}-n_{\mathbf{p}}}{\Omega-\epsilon_{{\bf p}+{\bf Q}}+\epsilon_{\mathbf{p}}+i\eta}.$$

(9)

Here we assume that ${\bf r}_{0}={\bf r}_{i+1}=0$ for convenience of compact notations. According to (8), for a given set of integration variables the integrand is a product of functions which depend only on one momentum variable. When the $\eta\to+0$ limit is taken, each $\tilde{\cal F}$ function may be split into up to three contributions, but the computational cost of evaluating the total contribution is now reduced to $3(i+1)$ instead of $3^{i+1}$. Additional real space integrals do not pose any technical or efficiency problems for Monte Carlo simulations, which all can be performed on a single laptop.

The convergence of the ladder series strongly depends on $R=\Omega/Qv_{F}$ ratio with $v_{F}$ the Fermi velocity, see panels (a) and (b) in Fig. 3 where we plot partial contributions to polarization $\Pi$ from ladder diagrams in $\Gamma^{(3)}$. For large and small $R$ we observe excellent convergence for all values of $\kappa$ considered in this work. However, in the $R\sim 1$ region the series start diverging for $\kappa\leq 1$, and the problem is getting progressively more severe as $\kappa$ is reduced. At $Q/k_{F}\approx 0.1$, the series for $\kappa=0.8$ are already visibly beyond the convergence radius. The dimensionless parameter controlling series convergence can be estimated from

$$g\sim\rho_{F}\frac{4\pi e^{2}}{\kappa^{2}}\sim\frac{\kappa_{TF}^{2}}{2\kappa^{2}}$$

(10)

where $\rho_{F}$ is the Fermi-energy density of states per spin. When the screening momentum is increased to $\kappa=1.2$, the series start to converge for any value of $R$ (similarly to what is observed for larger momentum and higher temperature in panels (a) and (b) in Fig. 4).

If we formally consider ladder diagrams as Taylor series expansion in powers of parameter $\xi$ with the physical result corresponding to $\xi_{0}=1$, see the lower-right block in Fig. 1, then series resummation can be performed with the help of the conformal map transformation, $z=\xi f(\xi)\rightleftarrows\xi=zs(z)\equiv\sum_{m=1}^{\infty}c_{m}z^{m}$. By substituting $\xi(z)$ into the expression for $\Pi$ in Fig. 1(e), one generates Taylor series in powers of $z$ with the physical result corresponding to $z_{0}=\xi_{0}f(\xi_{0})$. The map is designed to have $z_{0}$ within the convergence radius so that the new Taylor series in powers of $z_{0}$ converge. We find that all our data can be successfully resummed by using a simple-pole map, $z=\xi/(\xi+\xi_{p}(R))$, with $\xi_{p}(R=1)\sim 1$ and increasing away from the $R=1$ point. The final results for three values of the Yukawa screening are shown in panels (c) and (d) in Figs. 3 and 4. Remarkably, despite substantial differences from the zeroth-order RPA-type results, there is extremely close agreement between the final curves for different values of $\kappa$ (note that $\kappa^{2}$ was changed by more than a factor of two). Resummation of Taylor series, as opposed to self-consistent and skeleton sequences, does not lead to multi-valued solutions of the type observed in Ref. [25], not to mention that resummed and converges ”as is” curves end up on top of each other.

At large momentum $Q/k_{F}=1$ and higher temperature $T/\varepsilon_{F}$, the series converge well for all values of $\kappa/k_{F}$ considered and can be summed up “as is”. The final converged results are shown in panels (c) and (d) in Fig. 4. The dependence on the Yukawa screening is barely detectable in this parameter regime. While Figs. 3 and 4 present results for $r_{s}=2$, we did not see any notable differences in efficiency of calculations neither for $r_{s}=1$ nor for $r_{s}=4$, in agreement with Eq. (10) which predicts that for $\kappa\sim\kappa_{TF}$ one has $g\sim 1$.

Within the RPA, the slope of $\textrm{Im}\Pi\propto\Omega/Q$ dependence is determined by the electron dispersion near the Fermi surface. Multiple scattering changes this result by nearly $100\%$, which has direct implications for accuracy of Landau damping (and other energy looses in metals) estimates based on the RPA-type calculations (for more details on Landau damping see Ref. [24]).

We presented an efficient technique for computing ladder-diagram series for polarization on the real frequency axis at finite temperature, or, equivalently, for solving the Bethe-Salpeter equation for the $3$-point vertex function in the ladder approximation and convoluting it with two Green’s functions in the Hartree-Fock basis. Our method is numerically exact once all approximations at the level of diagrams involved are specified, i.e. it does not involve the ill-conditioned numerical analytic continuation, finite regularization of the poles, or momentum grids for evaluating multi-dimensional integrals.

Remarkably, final results are nearly independent of the Yukawa screening used to setup the expansion, indicating small bias coming from this auxiliary parameter. We find that multiple scattering of the particle-hole pair results in significant changes of the polarization dependence on the $\Omega/Q$ ratio and, correspondingly, of the Landau damping. We expect that our scheme can be adapted for analogous simulations of realistic materials, and further advanced to include retardation effects for effective interaction potentials.

This work was supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, under Award DE-SC0023141.