Numerically exact generalized Green's function cluster expansions for electron-phonon problems

We generalize the family of approximate momentum average methods to formulate a numerically exact, convergent hierarchy of equations whose solution provides an efficient algorithm to compute the Green's function of a particle dressed by bosons suitable in the entire parameter regime. We use this approach to extract ground-state properties and spectral functions. Our approximation-free framework, dubbed the generalized Green's function cluster expansion (GGCE), allows access to exact numerical results in the extreme adiabatic limit, where many standard methods struggle or completely fail. We showcase the performance of the method, specializing three important models of charge-boson coupling in solids and molecular complexes: the molecular Holstein model, which describes coupling between charge density and local distortions, the Peierls model, which describes modulation of charge hopping due to intersite distortions, and a more complex Holstein$+$Peierls system with couplings to two different phonon modes, paradigmatic of charge-lattice interactions in organic crystals. The GGCE serves as an efficient approach that can be systematically extended to different physical scenarios, thus providing a tool to model the frequency dependence of dressed particles in realistic settings.

Figure 1

Image of a $L=4,$ $\mathbf{n}=[2,0,0,3]$ boson cloud (blue), such that $n_{\mathrm{T}}=5,$ contained within a variational space specified by a maximum cloud extent $M=5$ and maximum number of bosons $N=6.$ In this example, the constraint $M=5$ spans sites $i-1$ to $i+4$ so bosons can be created only on these sites (e.g., green circles) and are not allowed outside of the $M$-site cloud (e.g., red circle). $N=6$ implies that states with two more bosons than those depicted in the figure (blue circles) are omitted from the variational space. Note that the carrier (not shown) is allowed to be anywhere on the chain.

Figure 2

Ground-state energy $E_{\mathrm{GS}}/t$ as a function of coupling strength $\lambda ={\alpha }^2/2\mathrm{\Omega }t$ for the Holstein model in adiabaticity limits extending from anti $(\mathrm{\Lambda }\gg 1)$ to extreme adiabatic $(\mathrm{\Lambda }\ll 1).$ Values for $\mathrm{\Omega }/t$ are shown in grey and GGCE results for $\mathrm{\Omega }/t=0.1$ and 0.5 are compared to DMC data (symbols) obtained from Ref. [61]. Ground-state peak locations are converged with respect to $N$ and generally require $\approx 10$ bosons at small couplings but up to $\approx 30$ at large couplings.

Figure 3

Spectral function $A(k,\omega )$ at $k=0,\pi /2$ and $\pi$ for the Holstein model in (a) the mildly adiatabic ($\mathrm{\Omega }/t=0.5$) and (b) adiabatic limits ($\mathrm{\Omega }/t=0.1$), for dimensionless couplings $\lambda =0.2,0.5$, and 0.8. In this figure, we use $\eta =0.005$. Various values of $(M,N)$ are shown in the legend. Specifically, the gradient from red to black shows convergence with respect to $M$ for fixed $N.$ We demonstrate convergence with respect to $N$ via the blue line, which shows results that use the largest used $M,$ but with one less boson than the largest-used $N.$ The onset of the continuum is shown in shaded gray, and is defined as $E_{\mathrm{GS}}+\mathrm{\Omega },$ where $E_{\mathrm{GS}}$ is the polaron ground-state energy.

Figure 4

Spectral function $A(k,\omega )$ (scaled to a maximum of 1) of the Peierls model for $\mathrm{\Omega }/t=1,\eta =0.05$, and various values of the dimensionless coupling strength, $\lambda .M=5$ and $N=10$ used here are sufficient for convergence of the bands on the scale of the plot. It is worth noting that fine structures in states above the lowest energy band can be resolved on a finer grid and smaller value of $\eta ,$ although the intensity of these states is at most roughly an order of magnitude smaller than that of the lowest energy band.

Figure 5

Spectral function $A(k,\omega )$ (scaled to a maximum of 1) of the mixed-boson mode Holstein$+$Peierls model for various values of ${\mathrm{\Omega }}_{\mathrm{H}}$ and ${\mathrm{\Omega }}_{\mathrm{P}},{\lambda }_{\mathrm{H}}={\lambda }_{\mathrm{P}}=1,t=1$, and $\eta =0.05.$ For these calculations, we use $M_{\mathrm{H}}=M_{\mathrm{P}}=3$, and a maximum total cloud length, or absolute extent, $A=3$ (see Appendix pp3), and $N_{\mathrm{H}}=N_{\mathrm{P}}=5,$ for which semiquantitive convergence is achieved.

Figure 6

Convergence of the energy of the ground-state polaron band, $E_{\mathrm{P}}\left(k\right)$, for parameters shown in Fig. 5 against combinations of $M_{\mathrm{H}},M_{\mathrm{P}},A,N_{\mathrm{H}}$, and $N_{\mathrm{P}}.$

Figure 7

Example of a configuration of HP bosons corresponding to $n_{\mathrm{H}}=[1,0,2,1,0]$ and $n_{\mathrm{P}}=[0,0,1,3,1].$ Similar to the single boson models, we require that ${\sum }_j{n}_{\mathrm{H}}^{\left(j\right)}\le N_{\mathrm{H}},{\sum }_j{n}_{\mathrm{P}}^{\left(j\right)}\le N_{\mathrm{P}},L_{\mathrm{H}}\le M_{\mathrm{H}},L_{\mathrm{P}}\le M_{\mathrm{P}}$, and $L\le A.$