Fermi blockade of the electron-phonon interaction: why strong coupling effects may not be seen in optimally doped high temperature superconductors

We study how manifestations of strong electron-phonon interaction (EPI) depend on the carrier concentration by solving the two-dimensional Holstein model for the spin-polarized fermions using an approximation free bold-line diagrammatic Monte Carlo (BDMC) method. We show that the strong EPI, obviously present at very small Fermion concentration, is masked by the Fermi blockade effects and Migdal's theorem to the extent that it manifests itself as moderate one at large carriers densities. Suppression of strong EPI fingerprints is in agreement with experimental observations in doped high temperature superconductors

Discussions on the role of the EPI in the physics of cuprate compounds with high superconducting transition temperatures (high T c ) have been going for decades [1][2][3][4][5][6] without resulting in a consensus opinion. While the role of EPI in superconductivity is still under debate, its strong manifestations were clearly observed in numerous other phenomena in high T c materials [5][6][7][8][9][10][11][12][13][14]. The apparent puzzle is that strong EPI effects seen in spectroscopic data of undoped and weakly doped compounds become much less pronounced with hole doping [15][16][17][18]. Hence, having a clear picture of how the EPI effects change with the carrier concentration is of seminal importance for understanding the nature of unconventional superconductors where rigorous studies are hindered by the complexity of many-body fermion problem. Accurate results on the EPI in many-fermion systems may provide the way to reconcile the observed fingerprints of the strong EPI in the underdoped regime with successful descriptions of the strongly doped high T c materials by models based on direct electron-electron interactions alone.
More generally, it is a long standing fundamental problem to reveal how the Migdal's theorem [19,20] emerges at the large fermion concentration and eliminates the need for vertex corrections even for strong EPI, provided the Fermi-liquid state remains stable. The crossover between the two regimes is expected to take place at ω ph ∼ ε F , where ω ph is the phonon frequency and ε F is the Fermi energy, and it can be addressed by the approximation diagrammatic Monte Carlo methods [8,[21][22][23]. To this end, we consider a spin polarized (SP) twodimensional (2D) lattice system in order to avoid system instabilities that would be triggered by the strong EPI in continuous and spin-balanced systems, such as structural transitions or a singlet on-site bipolaron formation at λ ≈ 0.5 (in 2D) [24] with the concomitant superconducting state. An essential feature of the SP model resembling that of the t − J model near half-filling [25,26] (which is prototypical for description of high T c superconductors) is that in both cases one can only create one hole per site.
In this work we employ the BDMC technique developed for many-body systems with EPI in Ref. [23]. For the same system parameters the determinant Monte Carlo [27,28] method would suffer from a severe sign problem. The dynamical mean-filed method (DMFT) [29,30], would be inadequate because the EPI self-energy is strongly momentum dependent at low carrier concentration [23], in violation of the key DMFT assumption. The BDMC technique is based on the expansion of irreducible free-energy Feynman diagrams in terms of exact electron, G, and bare, D (0) , phonon propagators [31] and is free from the above limitations. In more detail, see Ref. [23], the electron self-energy Σ (m) is expanded into the series of irreducible skeleton graphs up to the largest order m defined by the number of D (0) propagators, with self-consistency implemented by a feedback loop involving the solution of the algebraic Dyson equation, [G(k, ω )] −1 = [G (0) (k, ω )] −1 −Σ (m) (k, ω ), in momentum, k, and Matsubara frequency, ω = 2πT ( +1/2), representation ( is an integer). The 2D Holstein model on a square lattice reads: where c † i /b † i are standard notations for electron/phonon creation operators, t is the nearest neighbor hopping amplitude, ω ph = 0.5t is the energy of the local optical mode, and g is the EPI coupling. The electron gas is spin-polarized and, hence, any site can be occupied by no more than one electron. It is standard to characterize the strength of the EPI by a dimensionless coupling constant λ = g 2 /(4ω ph t). The lattice constant a, amplitude t, and Planck's constanth are used to set units of length, energy, and time, respectively. In this study we chose λ = 1.07 beyond the crossover from weakto strong-coupling regimes for single polarons and the threshold for the singlet bipolaron formation. For convenient systematic error-free handling of the data in momentum space we perform simulations for finite systems arXiv:2007.09888v1 [cond-mat.supr-con] 20 Jul 2020

M i g d a l t h e o r e m s t e p b y s t e p v i o l a t i o n
FIG. 1. Quasi-particle residue at the Fermi (FS) as a function of ratio between the Fermi energy and phonon frequency without (m = 1) and with vertex corrections (m > 1). Symbols and dashed lines represent data obtained by skeleton expansions truncated at some finite order m. The solid red line with stars is obtained by extrapolation to the infinite diagramorder limit m → ∞. The errorbars, if not visible, are smaller than the symbol sizes. with 16 × 16 sites, large enough to reproduce the infinite system results with high accuracy (see Supplemental Material [32]). The temperature is set to T = t/20, which is an order of magnitude smaller than all energy scales of the model parameters. In the zero-density limit an alternative exact (numerically) diagrammatic Monte Carlo (DMC) approach for single polarons [21,22]  Our main results are presented in Figs. 6, 7, and 8. Figure 6 shows the dependence of the QP residue on the adiabaticity ratio γ = ε F /ω ph . One can see in Fig. 6 that at large γ ≥ 3 the Migdal's theorem ensures that vertex corrections are small and the lowest-order m = 1 skeleton diagram for self-energy (also known, depending on the context, as the non-crossing, the self-consistent Born, and the Eliashberg approximations) well describes the EPI renormalization even at strong coupling. In contrast, for smaller values of γ one has to account for highorder vertex corrections; up to order m = 7 at γ = 1 and all the way to m > 20 for γ → 0 with extrapolation to the infinite diagram-order limit. An immediate conclusion is that EPI strongly suppresses the QP residue to values smaller that 0.5 (indicative of strong coupling) only at a rather small filling factor when γ < 1. In Fig. 7(a) we further quantify the role of vertex corrections at low and high carrier density (or occupation number per site), δ, in both adiabatic and anti-adiabatic regimes. Vertex corrections become important at γ < 3, and at low values of γ and δ it is not sufficient to take into account just m = 2, or even m = 3 contributions; in this parameter regime the convergence is reached only for m ≥ 16 in the skeleton expansion, see Fig. 7(a)). Signatures of strong EPI are observed at δ < 0.1 that roughly corresponds to γ ≈ 1. The key conclusion that clear manifestations of strong EPI are limited to small doping is consistent with experimental findings for high T c superconductors [15][16][17][18].
One evidence for Fermi blockade of the EPI with doping comes from angle resolved photoemission experiments [15]. It was shown that the kink angle, related to the ratio, v high /v low , between the phase velocities of the dispersion relation above and below the Debye frequency, decreases with doping. Our simulations reveal a similar trend, see Fig. 8. The QP dispersion relation ω(k) was obtained from the energy of the lowest peak in the Lehmann spectral function, see Fig. 9, extracted from the imaginary time Matsubara Green function G(τ ) by the stochastic optimization with consistent constraints method of analytic continuation [22,33].
All data for the QP residues at the FS, also denoted as Z F S , were deduced from the Fermi-liquid relation, , for: δ = 3.8 × 10 −4 (squares connected by the solid line), δ ≈ 0.006 (circles connected by the dashed line), and δ ≈ 0.015 (diamonds connected by the dotted line). In the inset we present the effective coupling constant λ e deduced from the scaling relation (2) using experimental data for LSCO [18] (squares connected by a dashed line) and locations of theoretical spectral density maxima in Fig. 9 (circles connected by a solid line). We also re-plot the same theoretical data by using 4δ for the horizontal axis (diamonds connected by a dotted line). Spectral densities were computed for self-energies evaluated up to order m = 16 (δ = 3.8 × 10 −4 ), m = 7 (δ = 0.006), and m = 5 (δ = 0.015). These expansion orders are enough to have converged results for the corresponding carrier density (see Supplemental Material [32], Table I).
In the low-temperature limit, the self-energy derivative at zero frequency is accurately obtained from the ratio −Im[Σ(k F , )]/ω at the lowest Matsubara frequencies.
As expected, this procedure works perfectly at large carrier concentration. However, in the zero density limit the Fermi surface shrinks to a point at zero momentum, and the entire protocol becomes questionable. Spectral density offers an alternative way of computing the QP residue from the integrated weight of the lowest frequency peak (we denote it as Z GF ), see Fig. 9. Somewhat surprisingly, we find that even in the zero-density limit both procedures produce consistent results at any expansion order m, see inset in Fig. 9. At small, but finite concentration δ = 0.01526 (or γ = 0.334), with Fermi-momentum k F ≈ π/8 the agreement is even more precise: at order m = 5 we find that Z F S = 0.605 and Z GF (k F = π/8) = 0.611.
Calculations of the frequency dependent optical conductivity [16] and angle resolved photoemission spectra [18] in the low-concentration limit (one hole) of the t − J-Holstein model revealed that the experimental dependence of both quantities on δ can be reproduced theoretically if one introduces effective EPI coupling constant λ e (δ) that decreases with doping. It can be deduced from the photoemission spectra using scaling relation [18] derived from nonperturbative calculations for the t − J-Holstein model, where v low (v high ) is the velocity above (below) the kink energy ω ph . Note, the doubling of the spectral peak around the kink energy ω ph is a general feature of theoretical calculations [7,18,34,35]. These two peaks merge into a customary experimental picture of a single-peak kink at ω = ω ph when the theoretical spectra are broadened by experimental resolution or additional damping processes [18,34]. We compare λ e (δ) deduced from experimental data of Ref. [18] with our theoretical analysis in the inset of Fig. 8, dashed versus solid line. To have a meaningful quantitative comparison we also need to account for the difference between the non-degenerate spectrum of the spin-polarized Holstein model and fourfold degenerate ground state minimum of the experimental system. To this end we re-plot theoretical data using 4δ for the carrier concentration (dotted line). We observe semi-quantitative agreement between the theory and experiment despite a number of signifi- cant differences between the two cases at the microscopic level.
As already mentioned in connection with Figs. 6 and Fig. 7(a), at small doping the skeleton expansion needs to go beyond m = 16 in order to obtain correct results for the QP residue. However, both Z and the polaron energy E at the FS accurately follow an empirical scaling relation, a + b/ √ m, at any carrier concentration δ, see Fig. 5. This allows us to perform an extrapolation to the infinite-order limit to eliminate the remaining systematic error as shown in Figs. 6-7. The extrapolation procedure is validated by an excellent agreement between the BDMC result for the ground state energy of singlepolarons, E(m → ∞) = −4.89 and the DMC benchmark E 1 = −4.891. The single polaron zero temperature residue Z 1 = 0.238 is renormalized to Z 1 (β = 20) ≈ 0.31 due to finite temperature projection of the low energy self-trapping states [36,37] (see Supplemental Material [32]) which is also consistent with extrapolated value Z(m → ∞) ≈ 0.33.
The violation of Migdal's theorem for T = 0 is apparent in Fig. 6 for all filling factors except the two largest ones. At the lowest carrier concentrations the condition ε F T does not hold any more, but this fact is barely relevant for the discussion because the theorem is severely violated well before that, at ε F ∼ ω ph T . Thus, our finite temperature results are still valid for interpretation of the EPI suppression in high T c materials, which is observed from low to room temperatures [15][16][17][18].
Conclusions. We computed approximation-free results for the concentration dependence of the quasiparticle residue Z and kink angle caused by the strong electronphonon interaction in the spin-polarized two-dimensional Holstein model on the square lattice. We demonstrated that clear signatures of strong electron-phonon coupling at small carrier concentration are quickly suppressed for Fermi energies exceeding the phonon frequency. Our results provide detailed account for importance of highorder vertex corrections across the adiabatic crossover and demonstrate that Fermi blockade of the electronphonon interaction and irrelevance of vertex corrections both proceed in agreement with the Migdal's theorem. This picture explains experimental results reporting radical weakening of the electron-phonon coupling effects in lightly doped high temperature superconductors.
Supplementary material for "Fermi blocade of the electron-phonon interaction: why strong coupling effects may not be seen in optimally doped high temperature superconductors".
Here we provide additional details on calculations performed for the 2D Holstein model described in the main text. The lattice constant a, hopping amplitude t, and Planck's constanth are used to set units of length, energy, and time, respectively. The phonon frequency ω ph = 0.5t is nearly an order of magnitude smaller than the particle bandwidth, and the dimensionless coupling constant λ = 1.07 corresponds to the strong coupling regime (see also below).

Size dependence
To check whether the system size N × N = 16 × 16 is sufficient to reproduce properties of the Holstein model for single polarons when the largest finite-size effects are expected, we calculated various characteristics of the polaron by the diagrammatic Monte Carlo [21,22] and compared them with known infinite lattice results. In the simulations of finite lattice all momenta in the reciprocal space also form a lattice k x,y = (2π/N )j, −N/2 ≤ j < N/2 .
In Fig. 6 we show how the polaron energy, E, and quasiparticle residue, Z, depend on the lattice size for N = 4, 8, 16, 32, 64, 128, ∞, and conclude that N = 16 results reproduce the infinite system limit with accuracy of three to four significant digits.  Convergence properties of the skeleton expansion strongly depend on the fermion density δ (or chemical potential, µ, in the grand canonical ensemble). In Fig. 7. we present our BLDMC data for density dependence on the expansion order at low temperature T = t/20 and different values of µ. At low density one needs to account for vertex corrections up to order 16 to obtain converged results. Note that the chemical potential µ is not directly related to the Fermi energy counted counted from the bottom of the dispersion relation which is strongly renormalized by interactions. Table I provides final relations between all quantities, including the required expansion order. In Fig. 8 we present the spectral function of a single polaron in the infinite in finite N 2 = 16 2 systems. Nearly perfect agreement (well within the analytic continuation procedure uncertainties) proves that finite-size effects in this case are negligible not only for ground state energies but also for excited states.
To determine the quasiparticle residue and interaction induced energy shift, ∆E = E − (−4t), we rely on the standard reliable method: at large imaginary time the asymptotic decay of the Green's function is given by see [21,22], allowing one to extract Z and ∆E from a simple exponential fit. The leading correction decays with exponent controlled by the lowest excited state (the second polaron state according to the spectral density analysis). In Fig. 9 we show how Z and E estimates change when we move the fitting interval [0.95τ max , 1.05τ max ] to larger values of τ max . It is clear from Fig. 8 that the energy dependence on τ max within the range τ max ∈ [14, 80] is very weak (about 2%). This is in sharp contrast, with the quasiparticle residue estimates: Z increases by nearly 40% when τ max decreases from τ max = 80 to τ max = 14. This sensitivity explains the discrepancy between the calculations performed at finite temperature T = t/20 and at T = 0. We attribute it to the presence of the second polaron state with comparable Z factor and relatively small excitation energy E 2 − E G ≈ 0.17t.
Relation of the single polaron parameters and results of extrapolation procedure for BDMC data The extrapolation procedure is validated by an excellent agreement between the BDMC result for the ground state energy of single-polarons, E(m → ∞) = −4.89 and the DMC benchmark E 1 = −4.891 . In the same limit, the extrapolated result for the QP residue Z(m → ∞) ≈ 0.33 turns out to be larger than that for single polarons, Z 1 = 0.238. The reason for the discrepancy is a combination of the finite temperature effect and self-trapping phenomenon [36,37], manifesting itself as a second, low energy, E 2 − E 1 ≈ 0.17 < ω ph , excited polaron state with rather large spectral weight, Z 2 ≈ 0.3, clearly observed in the spectrum of single polarons at T = 0, see Fig. 8. Because of this soft excitation, the standard procedure of extracting Z from the large-τ asymptotic behavior of the imaginary time Green function G(τ ) [21,22] turns out to be sensitive to the choice of the large imaginary time used to fit the data (for τ < 40), whereas the estimate for energy remains accurate even for τ < 20, see Fig. 9 . Therefore, at T = t/20 we detect the QP weight that overestimates Z 1 of single polarons in the ground state. Semi-quantitatively, the finite temperature BDMC result can be understood from the relation Z β/2 = Z 1 +Z 2 exp[−(β/2)(E 2 −E 1 )] ≈ 0.31, which accounts for the activated second polaron state contribution at τ max = β/2.