Relative importance of nonlinear electron-phonon coupling and vertex corrections in the Holstein model

Determining the range of validity of Migdal’s approximation for electron-phonon (e-ph) coupled systems is a long-standing problem. Many attempts to answer this question employ the Holstein Hamiltonian, where the electron density couples linearly to local lattice displacements. When these displacements are large, however, nonlinear corrections to the interaction must also be included, which can significantly alter the physical picture obtained from this model. Using determinant quantum Monte Carlo and the self-consistent Migdal approximation, we compared superconducting and charge-density-wave correlations in the Holstein model with and without second-order nonlinear interactions. We find a disagreement between the two cases, even for relatively small values of the e-ph coupling strength, and, importantly, that this can occur in the same parameter regions where Migdal’s approximation holds. Our results demonstrate that questions regarding the validity of Migdal’s approximation go hand in hand with questions of the validity of a linear e-ph interaction. Phonon mediated superconductivity is understood using Migdal’s approximation and models with linear electron-phonon coupling, but there are limitations to the information these approximations can provide. Here, the authors investigate the influence of nonlinear electron-phonon interactions on the charge density wave and superconducting phases to determine their relevance on the validity of Migdal’s approximation.

O ur modern understanding of phonon-mediated superconductors is largely based on results from ab initio approaches 1 coupled with Migdal's approximation [2][3][4] . Migdal's approximation 2 neglects corrections to the electron-phonon (e-ph) interaction vertex, which scale as Oðλ _Ω E F Þ, where λ is a dimensionless measure of the e-ph coupling strength, ℏΩ is the typical phonon energy, and E F is the Fermi energy. Physically, this approximation neglects processes leading to polaron formation, and determining precisely when these processes become important and their impact on transport properties is a long-standing problem [5][6][7][8][9][10][11][12] .
Many attempts to address this question have utilized nonpertubative simulations of simplified effective models like the Holstein 13 or Fröhlich 14 Hamiltonians, where the electron density couples linearly with phonon fields. For example, owing to it's relative simplicity, the Holstein model and its extensions have been studied extensively using quantum Monte Carlo (QMC) 5,[15][16][17][18][19][20][21][22][23][24][25][26][27][28][29] , and serves as a prototype for studying different polaronic regimes. Recently, it was shown that even if _Ω E F < 1, one can find instances where vertex corrections (i.e., polaron formation) become important for λ ≈ 0.4-0.5 8,9 . It is generally understood that small (large) polarons form when the polaron binding energy is larger (smaller) than the hopping energy of the carriers 30 . However, small polarons are often accompanied by sizable lattice distortions and a tendency toward localization and charge order. This observation has motivated some work to include higher-order nonlinear e-ph coupling terms to study changes in polaron formation 31 and on charge-density-wave (CDW) and superconducting (SC) pairing correlations 32,33 . These studies found that small positive (negative) nonlinear terms decrease (increase) the effective mass of the carriers and contracts (enlarges) the local lattice distortions surrounding the carriers 31 . Furthermore, mean-field treatments aiming to recover a linear model via effective model parameters fail to capture the quantitative nature of the true nonlinear model [31][32][33] , indicating that nonlinearities cannot be renormalized out of the problem. For example, one can tune the parameters of an effective linear model to capture either the electronic or phononic properties of the nonlinear model but not both simultaneously 33 . This failure is important to note in the context of polaron formation, where the electrons and phonons become highly intertwined. To capture this physics accurately, an effective model must describe both degrees of freedom on an equal footing, and an effective linear description of a nonlinear e-ph model will not do this.
These results raise an important question about the priority of investigations into the validity of the aforementioned approximations. Are there scenarios where the breakdown of the linear approximation supersedes the breakdown of Migdal's approximation? In this work, we show that this is indeed the case. Specifically, by comparing QMC simulations of the (non)linear Holstein model with results obtained with the Migdal approximation's, we show that nonlinear corrections can be more important than vertex corrections, and that this can occur even when Migdal's approximation appears to be valid. Our results have consequences for any conclusions drawn about the validity of Migdal's approximation from model Hamiltonians and highlight a critical need to move beyond such models for a complete understanding of strong e-ph interactions.

Model.
We study an extension of the Holstein model that includes nonlinear e-ph interaction terms and defined on a twodimensional (2D) square lattice. The Hamiltonian isĤ ¼Ĥ el þĤ lat þĤ int , wherê describe the noninteracting electronic and phononic parts, respectively, and describes the e-ph interaction to kth order in the atomic displacement. Here,ĉ y i;σ (ĉ i;σ ) create (annihilate) spin σ(=↑, ↓) electrons on site i,n i ¼ P σĉ y i;σĉi;σ is the number operator, μ is the chemical potential, t is the nearest-neighbor hopping integral, and 〈i, j〉 restricts the summation to nearest-neighbors only. Each ion has a mass M with position and momentum operators denoted bŷ X i andP i , respectively. Quantizing the lattice vibrations leads to Einstein phonons created (annihilated) by the operatorb y i (b i ) with associated phonon frequency Ω. Lastly, α k and g k are the e-ph interaction strengths in the two representations, which are Following previous works on the nonlinear Holstein model 32,33 , we truncate the series inĤ int to second order and introduce the ratio ξ = g 2 /g 1 to quantify the relative size of the two e-ph couplings. (The standard Holstein model is recovered by setting g 2 = 0.) This simplification is sufficient to assess the relative importance of nonlinear interactions relative to Migdal's approximation. Additional orders up to k = 4 have been studied in the single carrier limit 31 , where they produce the same qualitative picture.
To facilitate comparison with previous work, we set k B = ℏ = t = M = 1. The scale of atomic displacements is set by the oscillation amplitude of the free harmonic oscillator A ¼ ffiffiffiffiffiffiffiffiffiffiffi 1=2Ω p ( ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi _=2MΩ p with the physical units restored). When reporting expectation values of X and its fluctuations, we explicitly divide by A in model units, thereby making the results dimensionless. To get an idea for what these values mean in reality, one can simply multiply the results by A in physical units. Later, we will consider FeSe to estimate the strength of the nonlinear interactions. In that case, the prefactor (in physical units) is A~0.036 Å, which is obtained after adopting a selenium mass M = 1.31 × 10 −25 kg and the experimental phonon energy ℏΩ = 20.8 meV of the A 1g mode (see Supplementary Note 1 for more details). Alternatively, we obtain a comparable scale of A~0.051 Å for the transition metal oxides, where ℏΩ = 50 meV and M = 2.66 × 10 −26 kg are typical for the optical oxygen phonons. Finally, we adopt the standard definition for the dimensionless linear e-ph interaction strength λ ¼ 2g 2 1 =WΩ, where W = 8t is the bandwidth.
In what follows, we compare results obtained using determinant quantum Monte Carlo (DQMC) 34 and the self-consistent Migdal approximation (SCMA) 16 (see "Methods" section and Supplementary Note 2). To determine the relative importance of nonlinear e-ph interactions against vertex corrections to Migdal's approximation, we juxtapose results obtained using these methods for the linear and nonlinear models. For example, comparing results obtained from DQMC and the SCMA for the linear model reveals the importance of vertex corrections. Analogously, comparing DQMC results for the linear and nonlinear models provides a measure for the importance of nonlinear interactions while treating the two models exactly. This methodology will allow us to isolate the source of any observed discrepancies. Deviations between SCMA and DQMC for the linear model must be due to vertex corrections, while disagreement between DQMC results for the linear and nonlinear models must arise from the additional quadratic interaction.
Comparison of susceptibilities at half-filling. We begin by comparing the susceptibilities for charge-density-wave (CDW) χ CDW (π, π) and pairing χ SC correlations for a few illustrative cases (see "Methods" section). The first comparison takes place at halffilling n hn i i ¼ 1, where both CDW correlations and lattice displacements are significant. For example, at λ = 0.2, Ω = 0.5t, and T = 0.1t, we obtain |〈X i,l 〉|/A~1.97 and 2.53 for ξ = 0.05 and ξ = 0, respectively. Taking A ≈ 0.036 Å for FeSe, these values correspond to~2.4 and 3.1% of the 2.95 Å Fe-Fe bond length. Similarly, taking A ≈ 0.051 Å translates to 5.1 and 6.6% of the typical 1.96 Å Cu-O bond-distance in a high-T c superconducting cuprate. These displacements are not negligible (as we will show) when the nonlineararities are included, particularly given the weak values of the coupling we consider here. Later, we also discuss the size of the corresponding vibrational fluctuations. Figure 1 presents results for the temperature dependence of χ CDW (π, π) and χ SC using Ω = 0.5t, N = 8 × 8, and λ = 0.1 and 0.2. This parameter set corresponds to weak coupling and satisfies the adiabatic criterion Ω E F < 1, where we expect the SCMA to hold. Indeed, when λ = 0.1 (Fig. 1a, b), there is fair agreement between the DQMC results obtained from both the linear (ξ = 0) and nonlinear (ξ = 0.05) e-ph models (symbols with solid curve), as well as the SCMA results for the linear model (dash-dot curve). When λ = 0.2 (Fig. 1c, d), however, we find significant disagreement between the results for ξ = 0 and ξ = 0.05 in both susceptibilities, especially at lower temperatures. In Fig. 1d, the DQMC and SCMA results mostly agree for the linear Holstein model (ξ = 0), but a small nonlinear correction of ξ = 0.05 yields a marked suppression the CDW correlations. The rapid onset of CDW order in the ξ = 0 case (Fig. 1d) coincides with a sharp downturn in χ SC (Fig. 1c), a feature which is not captured by the SCMA result.
The suppression of CDW correlations (Fig. 1d) in the presence of nonlinear e-ph coupling demonstrates the importance of higher-order interactions over vertex corrections in this case. As we show later, the need for nonlinear e-ph coupling is greatest near half-filling, where the CDW correlations are strongest. Of course, the downturn of the pairing susceptibility obtained from DQMC (ξ = 0) at lower temperatures (Fig. 1c) appears to indicate that vertex corrections are also important for capturing the low temperature behavior of χ SC at λ~0.2. This value of λ is smaller than the breakdown values reported in Esterlis et al. 9 , however, our models differ slightly. For one, they suppress the effects of Fermi-surface nesting by situating the electron density away from half-filling and also include hopping between next nearestneighbors. Second, they use an alternate definition for λ = α 2 N (E F )/MΩ 2 , where N(E F ) is the density of states evaluated at the Fermi energy. Nevertheless, the results in Fig. 1c, d reveal that the nonlinear corrections to the linear model are non-negligible at high temperature, even before the breakdown of Migdal's approximation becomes apparent. Pairing susceptibilities for large phonon frequency. Now we consider a counter comparison in the antiadiabatic regime with intermediate coupling by setting λ = 0.4, Ω = 4t, N = 10 × 10, and n = 0.55 (Fig. 2). Away from half-filling, the pairing correlations grow more rapidly in part due to the larger Ω, but also because of less competition with (incommensurate) CDW correlations (Supplementary Note 3). Each of the curves in Fig. 2 shows that the system has strong pairing correlations, but they would yield very different estimates for T c . The SCMA significantly overestimates χ SC , which is not surprising because Migdal's approximation is ill justified in this case (i.e., λ Ω E F $ 1). Interestingly, the nonlinear corrections become important at low temperature despite the presence of smaller lattice displacements (e.g., |〈X i,l 〉|/A ≈ 0.22).
Comparison over doping. Fig. 3 shows results for three combinations of λ and Ω/t over a wide range of electronic filling and at a fixed temperature T = 0.25t. The DQMC results for ξ = 0 (ξ = 0.05) are represented by open (closed) symbols in all three panels, while the SCMA results are shown as dashed or dotted lines in Fig. 3a, b. For reference, Fig. 3c shows the corresponding the average lattice displacement, obtained by averaging over all spacetime points hX i;l i ¼ 1 i;l X i;l . We caution that 〈X i,l 〉 provides a rough measure of the typical lattice displacements, and not a complete picture of the ionic subsystem. We will return to this subtle issue later, when we discuss the displacement fluctuations. Additional results away from half-filling are also provided in Supplementary Fig. 2.
Case (1), λ = 0.15, Ω = 0.5t: These parameters are nearly identical to those used in Fig. 1. The deviations in χ CDW (π, π) for ξ = 0 and ξ = 0.05 (Fig. 3a) become apparent near n ≥ 0.6, whereas the SCMA result starts to deviate from DQMC for n ≥ 0.8. At this temperature, the results for χ SC essentially agree (Fig. 3b), but the nonlinear model yields a smaller average displacement 〈X i,l 〉/A (Fig. 3c). These results further reinforce our prior observation that Migdal's approximation and the linear model can break down in different parameter regimes (in this case doping).
Case (2), λ = 0.3, Ω = 0.5t: Now we double λ while keeping the Ω fixed. The increase in λ produces larger average displacements (Fig. 3c) and more pronounced nonlinear corrections. It also induces a stronger CDW (Fig. 3a) for the linear (ξ = 0) model. The SCMA qualitatively captures the CDW correlations of the linear model in panel (a), but underestimates their strength, which can be attributed to solely to the vertex corrections, consistent with the conclusions of Esterlis et al. 9 . Due to the large CDW correlations, there is a suppression 35,36 in χ SC for ξ = 0, which is mostly captured by the SCMA (Fig. 3b). The introduction of the nonlinear interaction significantly reduces the CDW correlations and their competition with SC, which enhances χ SC at larger values of n.
Case (3), λ = 0.15, Ω = 4.0t: Now we look at the large phonon frequency results for DQMC (green and crimson triangles) and SCMA (green dotted line). The larger Ω boosts pairing  Error bars on the DQMC data were estimated using jackknife resampling; however, all one-sigma error bars are smaller than the symbol size and have been suppressed for clarity. Lines connecting DQMC data are used only to guide the eye. correlations (Fig. 3b) across the entire doping range and all of the χ SC 's are in fair agreement. However, we know from Fig. 2 that larger differences between each curve will emerge at lower temperatures. In fact, the SCMA already overestimates χ CDW (π, π) near half-filling at this temperature, a feature we may attribute to antiadiabaticity. The increased value of Ω means that the lattice vibrations are characterized by stiffer spring constants. We obtain smaller average lattice displacements at all n as a result (Fig. 3c), which reduces the importance of the nonlinear interaction and produces better agreement between ξ = 0 and ξ = 0.05 DQMC results.
We should be careful in interpreting the results at lower filling in each of the cases above. On the one hand, our examples suggest that corrections to the e-ph interaction are most important for describing the CDW phase transition near half-filling, which appears at higher temperatures. On the other hand, corrections could become important in the dilute carrier region at much lower temperatures. Nonetheless, our results suggest that the linear Holstein model is sensitive to nonlinear corrections over a large parameter space and that regions of this space overlap with regions where Migdal's approximation is not valid. But perhaps more importantly, there are regions where the linear approximation breaks down before Migdal's approximation does.
Average lattice displacement and its fluctuations. In Fig. 3c, we showed that the magnitude of the mean displacement 〈X i,l 〉/A had an approximately linear dependence on the filling n. These displacements become larger when the dimensionless e-ph coupling λ is increased or when the phonon energy Ω is decreased (or, equivalently, when the spring constants are softer). The behavior of 〈X i,l 〉 as a function of doping can be loosely understood by considering the atomic limit. In this case, the effect of the linear e-ph interaction is to shift the equilibrium position of the lattice to X 0 % Àn α 1 MΩ 2 37 . Indeed, the results shown in Fig. 3c for the linear model are well described by this function. Based on this observation, one might then be tempted to try to eliminate the nonlinear interactions by defining new lattice operatorsX 0 1 X À X 0 ¼X in hopes that the displacements ofX 0 remain small.
Unfortunately, this procedure is not viable for several reasons. The first reason is that global shift of the equilibrium position will only be effective in the case of a uniform charge distribution. This certainly will not be the case when the CDW correlations are significant. For example, in the (π, π) CDW phase, half of the sites are doubly occupied with an average displacement of ≈2X 0 while the remaining sites are unoccupied with an average displacement of zero. In this instance, 〈X i,l 〉 = X 0 , consistent with our results in Fig. 3c, but shifting the origin to X = X 0 will not eliminate the large lattice displacements at each site.
The second reason why redefining the origin will not work is that such transformations do not affect the displacement fluctuations, which are also significant for the linear Holstein model. To show this, we examine the root-mean-square (rms) displacement of X i,l in our system with a linear e-ph coupling strength λ, which is defined as and is formally identical to a standard deviation. The value of σ X in the limit n → 0 approaches the (thermal) rms displacement for the free harmonic oscillator, which we denote as σ X (0) and is given by where n B ðΩÞ ¼ ½e βΩ À 1 À1 is the Bose occupation function. Figure 4 shows results for [σ X (λ) − σ X (0)]/A, as a function of filling for the same parameters used in Fig. 3. (For reference, for Ω/t = 0.5 and 4.0, we obtain σ X (0)/A = 1.146 and 0.354, respectively.) Here, we see that the fluctuations of the linear Holstein model are quite sensitive to the size of the dimensionless linear e-ph interaction λ. Moreover, the magnitude of the fluctuations generally grows monotonically with filling until reaching a maximum at half-filling. There, the largest fluctuations shown correspond to σ X (λ = 0.3)/A ≈ 2.43, which is more than double the size of captured by σ X (0)/A. Again, taking FeSe or a typical cuprate as references, these fluctuations correspond tõ 3.0 and~6.3% of the respective lattice constants. It is important to note that σ X (0)/A rises sharply for even smaller (and more realistic) model values of Ω/t ≈ 0.02-0.1. Such values, however, are typically inaccessible to DQMC due to prohibitively long autocorrelation times.
The overall effect of the nonlinear coupling is to suppress the rms displacements relative to the linear case, especially near halffilling. Only when (ξ, λ, Ω/t) = (0, 0.15, 4.0) and (0.05, 0.15, 4.0) do we find close agreement between the linear and nonlinear models, and typical lattice displacements that are a small fraction of the lattice spacing.
How big are nonlinear interactions in materials? Throughout this work, we used ξ = 0.05, but how representative is this value for a realistic system? To address this question, we considered the case of bulk FeSe, a quasi-2D material where the position of the Se atoms influences the on-site energies of the Fe 3d orbitals 38 , somewhat akin to the 2D Holstein model. To determine the strength of the linear and nonlinear e-ph coupling, we constructed a Wannier function basis from density functional theory (DFT) calculations to determine the on-site energy of the Fe 3d orbitals ϵ diag (z Se ) as a function of the Se atom's static displacement z Se along the c-axis (Supplementary Fig. 1). We then applied a polynomial fit of the form f ðz Se Þ ¼ a 0 þ a 1 ðz Se À z 0 Þ þ a 2 ðz Se À z 0 Þ 2 to the site energy ϵ diag (z Se ) for each orbital and computed ξ ¼ A a 2 a 1 . Here, the oscillation amplitude % 0:036 Å, adopting a selenium mass The results are summarized in Table 1, where ξ ranges from −0.1640 to 0.0148, with the strongest nonlinearity appearing for the d xy orbital. We do not investigate ξ < 0 in our model calculations because others have shown that it leads to increased softening of the phonon dispersion and larger CDW correlations 33 . Nevertheless, our results show |ξ| ≈ 0.05 is certainly not out of the question for a real material.

Discussion
We have demonstrated that the linear approximation to the e-ph interaction in the Holstein model breaks down in commonly studied parameter regimes. Importantly, this breakdown regime overlaps with ones where Migdal's approximation captures the DQMC result, even if only qualitatively. This observation indicates that nonlinear corrections to the underlying linear lattice model may be important even when vertex corrections are not. We also studied the example of bulk FeSe from first principles and found that nonlinear e-ph interactions in a real materials can be quite significant and on par with, or even larger than our model choice of |ξ| = 0.05.
It is natural to wonder which parameter regimes might be best for ensuring lattice displacements remain small enough justify the use of a linear interaction. We have found that tuning λ to smaller values suppresses the lattice displacements and their fluctuations, but also pushes the growth of correlations to lower temperatures, making computations more expensive. (Some groups 29 have recently managed to access such temperatures in QMC, however.) Alternatively, one could also shrink the displacements by choosing antiadiabatic parameters (i.e., Ω > E F ). But even for a strongly antiadiabatic choice of (λ, Ω/t, n) = (0.4, 4.0, 0.55), nonlinear corrections to the e-ph interaction produced considerable differences in the resulting temperature dependence of the superconducting susceptibility. Unfortunately, focusing on smaller phonon energies, which are relevant for real materials, will also produce larger lattice displacements and fluctuations that are inconsistent with a linear interaction. While our results are not comprehensive across the entire parameter space of the Holstein model, we are forced to conclude that they do call large portions of this space into question. For instance, our results imply that combinations of λ ≳ 10 −1 and Ω ≲ 4t yield sizable displacements and displacement fluctuations, which would necessitate additional nonlinear interactions and/or anharmonic lattice potentials 24 . Our results indicate a clear and present need for more work extending beyond the simplest effective models, especially when one is trying to describe the physics of a real system. The Holstein model and Migdal's approximation have long served as cornerstones in the study of electron-phonon interactions. Their relative simplicity has helped shape our intuition about superconductivity, its competition with charge order, and polaron formation, and studying the Holstein model can address the essential physics of these processes. While it is clear that these models are built on the assumption of small lattice displacements, it is not always clear how large these displacements will be in practice or whether additional nonlinear interactions will modify the physics of the model. One must, therefore, be careful when extrapolating results from effective models to real materials when they are driven outside their range of validity. For example, we have shown that the Holstein model can produce displacements that begin to approach the Lindemann criteria for melting (particularly as Ω is reduced), but the model cannot describe such a transition. Instead, it over predicts various tendencies towards ordered phases in these cases. Similarly, it is unclear how one should map critical λ values derived for the breakdown of Migdal's approximation onto real materials.

Methods
Determinant quantum Monte Carlo. DQMC is an auxiliary field, imaginary-time technique that computes expectation values within the grand canonical ensemble. It is inherently nonperturbative and includes all Feynman diagrams. While QMC simulations of the (non)linear Holstein model face long autocorrelation times 39 , they are free of a sign problem. We refer the reader to Supplementary Note 2 for more details on our DQMC implementation and more generally to White et al. 34 , Scalettar et al. 15 , and Johnston et al. 37 .
The self-consistent Migdal approximation. The self-consistent Migdal approximation (SCMA) is a diagrammatic approach that neglects higher-order corrections to the e-ph interaction vertex. Here, we use a recently developed SCMA code that treats both the electron and phonon self-energies on an equal footing and captures the competition between the CDW and SC instabilities 40 .
Susceptibilities. The effects of nonlinear e-ph coupling or the omission of vertex corrections will manifest uniquely in different observables. Here, we focus on twoparticle correlation functions.
The CDW correlation at momentum q is measured by the charge susceptibility whereρ q ðτÞ P i;σ e ÀiqÁR in i;σ ðτÞ and hÂBi c ¼ hÂBi À hÂihBi denotes a connected correlation function. Near half-filling, the Fermi surface is well nested and the charge susceptibility has a single commensurate peak at q max ¼ ðπ; πÞ. Moving away from half-filling removes the nesting condition and eventually redistributes the weight of the single peak in χ CDW (q) into four incommensurate peaks ( Supplementary Fig. 3).
The e-ph coupling is also responsible for spin-singlet s-wave pairing, resulting in superconducting correlations measured by the pair-field susceptibility whereΔðτÞ ¼ P iĉi;" ðτÞĉ i;# ðτÞ.

Data availability
Data are available upon request. Results from fits of the on-site energy ϵdiag(zSe) for each d-orbital of Fe as a function of the height of the Se atom zSe measured with respect to the Fe-plane. Fitting ϵdiag(zSe) with a simple polynomial of the form fðz Se Þ ¼ a 0 þ a 1 ðz Se À z 0 Þ þ a 2 ðz Se À z 0 Þ 2 , we estimate the nonlinear coupling ratio ξ from the fitting parameters a2/a1.