Predominance of non-adiabatic effects in zero-point renormalization of the electronic band gap

Electronic and optical properties of materials are affected by atomic motion through the electron-phonon interaction: not only band gaps change with temperature, but even at absolute zero temperature, zero-point motion causes band-gap renormalization. We present a large-scale first-principles evaluation of the zero-point renormalization of band edges beyond the adiabatic approximation. For materials with light elements, the band gap renormalization is often larger than 0.3 eV, and up to 0.7 eV. This effect cannot be ignored if accurate band gaps are sought. For infrared-active materials, global agreement with available experimental data is obtained only when non-adiabatic effects are taken into account. They even dominate zero-point renormalization for many materials, as shown by a generalized Fr\"ohlich model that includes multiple phonon branches, anisotropic and degenerate electronic extrema, whose range of validity is established by comparison with first-principles results.

The electronic band gap is arguably the most important characteristic of semiconductors and insulators. It determines optical and luminescent thresholds, but is also a prerequisite for characterizing band offsets at interfaces and deep electronic levels created by defects 1 . However, accurate band gap computation is a challenging task. Indeed, the vast majority of first-principles calculations relies on Kohn-Sham Density-Functional Theory (KS-DFT), valid for ground state properties 2 , that delivers a theoretically unjustified value of the band gap in the standard approach, even with exact KS potential [3][4][5] .
The breakthrough came from many-body perturbation theory, with the so-called GW approximation, first non-selfconsistent (G 0 W 0 ) by Hybertsen and Louie in 1986 6 , then twenty years later self-consistent (GW) 7 and further improved by accurate vertex corrections from electron-hole excitations (GWeh) 8 . The latter methodology, at the forefront for bandgap computations, delivers a 2%-10% accuracy, usually overestimating the experimental band gap. GWeh calculations are computationally very demanding, typically about two orders of magnitude more than G 0 W 0 , itself two orders of magnitude more time-consuming than KS-DFT calculations, roughly speaking. a) Anna Miglio and Véronique Brousseau-Couture contributed equally to this work. b) Electronic mail: xavier.gonze@uclouvain.be Despite being state-of-the-art, such studies ignored completely the electron-phonon interaction. The electron-phonon interaction drives most of the temperature dependence of the electronic structure of semiconductors and insulators, but also yields a zero-point motion gap modification at T = 0 K, often termed zero-point renormalization of the gap (ZPR g ) for historical reasons.
The ZPR g had been examined 40 years ago, by Allen, Heine and Cardona (AHC) 9,10 who clarified the early theories by Fan 11 and Antoncik 12 . Their approach is, like the GW approximation, rooted in many-body perturbation theory, where, at the lowest order, two diagrams contribute, see Fig. 1, the so-called "Fan" diagram, with two 1 st -order electron-phonon vertices and the "Debye-Waller" diagram, with one 2 nd -order electron-phonon vertex. In the context of semi-empirical calculations, the AHC method was applied to Si and Ge, introducing the adiabatic approximation, in which the phonon frequencies are neglected with respect to electronic eigenenergy differences and replaced by a small but non-vanishing imaginary broadening 10 . It was later extended without caution to GaAs and a few other III-V semiconductors 13,14 .
In this work, we present first-principles AHC ZPR g calculations beyond the adiabatic approximation, for 30 materials. Comparing with experimental band gaps, we show that adding ZPR g improves the GWeh first-principles band gap, and moreover, that the ZPR g has the same order of magnitude as the G 0 W 0 to GWeh correction for half of the materials (typically materials with light atoms, e.g. O, N ...) on which GWeh has been tested. Hence, the GWeh level of sophistication misses its target for many materials with light atoms, if the ZPR g is not taken into account. With it, the theoretical agreement with direct measurements of experimental ZPR g is improved. This also demonstrate the crucial importance of phonon dynamics to reach this level of accuracy.
Indeed, first-principles calculations of the ZPR g using the AHC theory are very challenging, and only started one decade ago, on a case-by-case basis (see Refs [15][16][17][18][19][20][21][22][23][24][25][26] , and the SI II.D), usually relying on the adiabatic approximation, and without comparison with experimental data. An approach to the ZPR g , alternative to the AHC one, relies on computations of the band gap at fixed, distorted geometries, for large supercells (see Refs.) 18,[27][28][29][30][31][32][33][34] and the SI 26 , Sec.II.D). As the adiabatic approximation is inherent in this approach, we denote it as ASC, for "adiabatic supercell". A recent publication by Karsai and coworkers 35 presents ASC ZPR g based on DFT values as well as based on G 0 W 0 values for 18 semiconductors, with experimental comparison for 9 of them. Both AHC and ASC methodologies have been recently reviewed 36 .
Although the adiabatic approximation had already been criticized by Poncé et al. in 2015 24 (SI Sec.III.A) as causing an unphysical divergence of the adiabatic AHC expression for infrared-active materials with vanishing imaginary broadening, thus invalidating all adiabatic AHC calculations except for non-infrared-active materials, the full consequences of the adiabatic approximation have not yet been recognized for ASC. We will show that the non-adiabatic AHC approach outperforms the ASC approach, so that the predictions arising from the mechanism by which the ASC approach bypass the adiabatic AHC divergence problem for infrared-active materials are questionable. This is made clear by a generalized Fröhlich model with a few physical parameters, that can be determined either from first principles or from experimental data.
Although the full electronic and phonon band structures do not enter in this model, and the Debye-Waller diagram is ignored, for many materials it accounts for more than half the ZPR g of the full first-principles non-adiabatic AHC ZPR g . As this model depends crucially on non-adiabatic effects, it demonstrates the failure of the adiabatic hypothesis be it for the AHC or the ASC approach.
By the same token, we also show the domain of validity and accuracy of model Fröhlich large-polaron calculations based on the continuum hypothesis, that have been the subject of decades of research [37][38][39][40][41][42][43] . Such model Fröhlich Hamiltonian capture well the ZPR for about half of the materials in our list, characterized by their strong infrared activity, while it becomes less and less adequate for decreasing ionicity. In the present context, the Fröhlich large-polaron model provides an intuitive picture of the physics of the ZPR g .

RESULTS
Zero-point renormalization: experiment versus first principles. Fig. 2 compares first-principles ZPR g with experimental values. As described in the METHODS section, and in SI Sec.II.C, the correction due to zero-point motion effect on the lattice parameter, ZPR lat g , has been added to fixed volume results from both non-adiabatic AHC (present calculations) and ASC methodologies 35 . While for a few materials experimental ZPR g values are well established, within 5-10%, globally, experimental uncertainty is larger, and can hardly be claimed to be better than 25% for the majority of materials, see Sec.II.A in the SI. This will be our tolerance.
Let us focus first on the ASC-based results. For the 16 materials present in both Ref. 35 and the experimental set described in the SI Table S1, the ASC vs experimental discrepancy is more than 25% for more than half of the materials 44 . There is a global trend to underestimation by ASC, although CdTe is overestimated.
By contrast, the non-adiabatic AHC ZPR g (blue full circles) and experimental ZPR g agree with each other within 25% for 16 out of the 18 materials. The outliers are CdTe with a 43% overestimation by theory and GaP with a 33% underestimation. For none of these the discrepancy is a factor of two or larger. On the contrary, from the ASC approaches, several materials show underestimation of the ZPR g by more than a factor of 2. The materials showing such large underestimation (CdS, ZnO, SiC) are all quite ionic, while more covalent materials (C, Si, Ge, AlSb, AlAs) are better described.
Therefore, Fig. 2 clearly shows that the non-adiabatic AHC approach performs significantly better than the ASC approach. AHC ZPR g and ASC ZPR g also differ by more than a factor of two for TiO 2 and MgO (see SI Sec.II.D), although no experimental ZPR g is available for these materials to our knowledge.
We now examine band gaps. Fig. 3 presents the ratio between first-principles band gaps and corresponding experimental values, for 12 materials. The best first-principles values at fixed equilibrium atomic position, from GWeh 8 , are represented, as well as their non-adiabatic AHC ZPR corrected values.
For GWeh without ZPR g , a 4% agreement is obtained only for two materials (CdS and GaN). There is indeed a clear, albeit small, tendency of GWeh to overestimate the band gap value, except for the 3 materials containing shallow core delectrons (ZnO, ZnS and CdS) that are underestimated. By FIG. 2. Absolute values of first-principles band gap renormalization ZPR g compared with experimental ones. Blue full circles: present calculations, using non-adiabatic AHC, based on KS-DFT ingredients; red empty triangles: adiabatic supercell KS-DFT Dashed lines: limits at which the smallest of both ZPR is 25% smaller than the largest one (note that the scales are logarithmic). For consistency, the computed lattice ZPR lat g were added to both ASC and AHC. The majority of non-adiabatic AHC-based results fall See numerical values in Table S2 of the SI. contrast, if the non-adiabatic AHC ZPR g is added to the GWeh data (blue dots), a 4% agreement is obtained for 9 out of the 12 materials (8 if ZPR lat g is not included). For ZnO and ZnS, with a final 10-12% underestimation, and CdS with a 5% underestimation, we question the GWeh ability to produce accurate fixed-geometry band gaps at the level obtained for the other materials, due to the presence of rather localized 3d electrons in Zn and 4d electrons in CdS.
As a final lesson from Fig. 3, we note that for 4 out of the 12 materials (SiC, AlP, C, and BN), the ZPR g is of similar size than the G 0 W 0 to GWeh correction, and it is a significant fraction of it also for Si, GaN and MgO. As mentioned earlier, GWeh calculations are much more time-consuming than G 0 W 0 calculations, possibly even more time-consuming than ZPR calculations (although we have not attempted to make a fair comparison). Thus, for materials containing light elements, first row and likely second row (e.g. AlP) in the periodic table, GWeh calculations miss their target if not accompanied by ZPR g calculations. A review of the variance and accuracy of G 0 W 0 calculations for these materials is discussed in Ref. 45. Table S5 of the SI gathers our full set of ZPR results, beyond those present in the ASC or experimental sets. It also includes 10 oxides, while the experimental set only includes three materials containing oxygen (ZnO, MgO and SrTiO 3 ), and there are none in Karsai's ASC set. The ZPR g of the band gap for materials containing light elements (O or lighter) is FIG. 3. Ratio between first-principles band gaps and experimental ones 46 . First-principles results without electron-phonon interation 8 are based on G 0 W 0 (diamonds) or on GWeh (squares). Blue full circles are the GWeh results to which AHC ZPR g and lattice ZPR lat g have been added, while empty red circles are the GWeh results to with only AHC ZPR g was added. GWeh usually overshoots, while adding the electron-phonon interaction gives better global agreement. See numerical values in Table S3 of the SI. between -157 meV (ZnO) and -699 meV (BeO), while, relatively to the experimental band gap, it ranges from -4.6% (ZnO) to -10.8% (TiO 2 -t). This can hardly be ignored in accurate calculations of the gap.

Generalized Fröhlich model.
We now come back to the physics from which the ZPR g originates, and explain our earlier observation that the ASC describes reasonably well the more covalent materials, but can fail badly for ionic materials. We argue that, for many polar materials, the ZPR g is dominated by the diverging electron-phonon interaction between zone-center LO phonons and electrons close to the band edges, and the slow (nonadiabatic) response of the latter: the effects due to comparatively fast phonons are crucial. This was already the message from Fröhlich 37 and Feynman 38 , decades ago, initiating large-polaron studies. However, large-scale assessment of the adequacy of Fröhlich model for real materials is lacking. Indeed, the available flavors of Fröhlich model, based on continuum (macroscopic) electrostatic interactions, do not cover altogether degenerate and anisotropic electronic band extrema or are restricted to only one phonon branch, unlike most real materials [47][48][49][50] .
In this respect, we introduce now a generalized Fröhlich model, gFr, whose form is deduced from first-principles equations in the long-wavelength limit (continuum approximation). Such model covers all situations and still uses as input only long-wavelength parameters, that can be determined either from first principles (SI Sec.IV.D) or from experiment (SI Sec.IV.E). We then find the corresponding ZPR g from perturbation theory at the lowest order, evaluate it for our set of materials, and compares its results to the ZPR g obtained from full first-principles computations. The corresponding Hamiltonian writes (see the SI Sec.IV.B for detailed explanations): with (i) an electronic part that includes direction-dependent effective masses m * n (k), governed by so-called Luttinger parameters in case of degeneracy, electronic creation and annihilation operators,ĉ + kn and c kn , with k the electron wavevector and n the band index, (ii) the multi-branch phonon part, with direction-dependent phonon frequencies ω j0 (q), phonon creation and annihilation operators,â + q j andâ q j , with q the phonon wavevector and j the branch index, and finally (iii) the electron-phonon interaction part The k and q sums run over the Brillouin zone. The sum over n and n runs only over the bands that connect to the degenerate extremum, renumbered from 1 to n deg . The generalized Fröhlich electron-phonon interaction 51,52 is This expression depends on the directionsk,q, andk (k = k + q), but not explicitly on their norm (hence only longwavelength parameters are used), except for the 1 q factor. The electron-phonon part also depends only on few quantities: the Born effective charges (entering the mode-polarity vectors p j ), the macroscopic dielectric tensor ε ∞ , and the phonon frequencies ω j0 , the primitive cell volume Ω 0 , the Born-von Karman normalisation volume V BvK corresponding to the k and q samplings. The s tensors are symmetry-dependent unitary matrices, similar to spherical harmonics. Eqs. (1-29) define our generalized Fröhlich Hamiltonian. Although we will focus on its T = 0 properties within perturbation theory, such Hamiltonian could be studied for many different purposes (non-zero T, mobility, optical responses ...), like the original Fröhlich model, for representative materials using first-principles or experimental parameters.
The corresponding ZPR gFr c can be obtained with a perturbation treatment (SI Sec.IV.C), giving ZPR gFr A similar expression exists for the valence ZPR The markers identify materials of similar ionicity: oxides (purple circles), other II-VI materials containing S, Se or Te (yellow squares), nitrides (green diamonds), other III-V materials, containing P, As, or Sb (dark blue triangles), and group IV material SiC (light blue hexagon). Dashed lines: limits at which the smallest of both ZPR is 25% smaller than the largest one. Dotted lines: limits at which the Fröhlich model ZPR misses respectively 50% and 66% of the AHC values (note that the scales are logarithmic). See numerical values in Table S5 of the SI. from experimental measurements, but are most easily computed from first principles, using density-functional perturbation theory with calculations only at q = Γ (e.g. no phonon band structure calculation). Eq. (39) can be evaluated for all band extrema in our set of materials, irrespective of whether the extrema are located at Γ or other points in the Brillouin Zone (e.g. X for the valence band of many oxides, with anisotropic effective mass), whether they are degenerate (e.g. the three-fold degeneracy of the top of the valence band of many III-V or II-VI compounds), and irrespective of the number of phonon branches (e.g. 3 different LO frequencies for TiO 2 , moreover varying with the direction along which q → 0). Fig. 4 compares the band gap ZPR from the first-principles non-adiabatic AHC methodology and from the generalized Fröhlich model. The 30 materials can be grouped into five sets, based on their ionicity: 11 materials containing oxygen, rather ionic, for which the Born effective charges and the ZPR are quite large, 6 materials containing chalcogenides, also rather ionic, 4 materials containing nitrogen and 5 III-V materials, less ionic, and 4 materials from group-IV elements, non-ionic, except SiC.
For oxygen-based materials, the ZPR ranges from 150 meV to 700 meV, and the gFr model captures this very well, with less than 25% error, with only one exception, BeO. The chalcogenide materials are also reasonably well described by the gFr model, capturing at least two-third of the ZPR. Globally their ZPR is smaller (note the logarithmic scale).
For the nitride materials and for SiC, the gFr captures about 50% of the quite large ZPR (between 176 meV and 406 meV). The adequacy of the gFr model decreases still with the III-V materials and the three non-ionic IV materials. In the latter case, the vanishing Born effective charges result in a null ZPR within the gFr model. (These three materials are omitted from Fig. 4).
For the oxydes and chalcogenides, the ZPR is thus dominated by the zone-center parameters (including the phonon frequencies), and the physics corresponds to the one of the large-polaron picture 38 , namely, the slow electron motion is correlated to a phonon cloud that dynamically adjusts to it. This physics is completely absent from the ASC approach. Even for nitrides, the gFr describes a significant fraction of the ZPR.
A perfect agreement between the non-adiabatic AHC firstprinciples ZPR and the generalized Fröhlich model ZPR is not expected. Indeed, differences can arise from different effects: lack of dominance of the Fröhlich electron-phonon interaction in some regions of the Brillouin Zone, departure from parabolicity of the electronic structure (obviously, the electronic structure must be periodic so that the parabolic behavior does not extend to infinity), interband contributions, phonon band dispersion, incomplete cancellation between the Debye-Waller and the acoustic phonon mode contribution.
It is actually surprising to see that for so many materials, the generalized Fröhlich model matches largely the firstprinciples AHC results. Anyhow, as a conclusion for this section, for a large number of materials, we have validated, a posteriori and from first principles, the relevance of large-polaron research based on Fröhlich model despite the numerous approximations on which it relies.

DISCUSSION
We focus on the mechanism by which the AHC divergence of the ZPR in the adiabatic case for infrared-active materials 24 is avoided, either using the ASC methodology or using the non-adiabatic AHC methodology. As Fröhlich and Feynman have cautioned us 37,38 , and already mentioned briefly in previous sections, the dynamics of the "slow" electron is crucial in this electron-phonon problem.
In the ASC approach, the bypass of this divergence can be understood as follows, see conduction band has a Bloch type wavefunction, with an envelope phase factor characterized by the wavevector k c multiplying a lattice-periodic function. Its density is lattice periodic. With such long-wavelength potential, as a function of the amplitude of the atomic displacements, the corresponding electronic eigenenergy changes first quadratically (as the average of the lowering and increase of potential for this Bloch wavefunction forbids a linear behaviour except in case of degeneracy), but for larger amplitudes, it behaves linearly, as the electron localizes in the lowered potential region and the minimum of the potential is linear in the amplitude of the atomic displacements. This is referred to as "nonquadratic coupling" in Ref. 33. A wavepacket is formed, by combining Bloch wavefunctions with similar lattice periodic functions but slightly different wavevectors (k c , k c + q, k c − q, etc), coming from a small interval of energy ∆ε 0 ∝ q 2 ∝ 1/(∆L) 2 , see By contrast, in the time-dependent case, as illustrated in Fig. 5(b), the wavepacket will require a characteristic time ∆τ, to form or to displace. This will be given by the Heisenberg uncertainty relation, ∆ε 0 ∆τ h, hence ∆τ h/∆ε 0 ∝ ∆L 2 . For long wavelengths, the characteristic time diverges. As soon as ∆τ is larger than the phonon characteristic time τ ph ∼ 1/ω LO , the "slow" electron will lag behind the phonon, and the static or adiabatic picture described above is no longer valid.
In all adiabatic approaches, either AHC or ASC, the elec-tron is always supposed to have the time to adjust to the change of potential, in contradiction with the time-energy uncertainty principle. Furthermore, the adiabatic AHC approach only considers the quadratic region for the above-mentioned dependence of eigenvalues with respect to amplitudes of displacements. This results in a diverging term 24 . At variance with the AHC case, the ASC approach samples a whole set of amplitudes, including the onset of the asymptotic linear regime, in which case the divergence does not build up. However, this ASC picture does not capture the real physical mechanism that prevent the divergence to occur, the impossibility for the electron to follow the phonon dynamics, that we have highlighted above. By contrast, such physical mechanism is present both in the non-adiabatic AHC approach and in the (generalized) Fröhlich model: the "slow" electron does not follow adiabatically (instantaneously) the atomic motion. The divergence of the adiabatic AHC is indeed avoided in the nonadiabatic picture by taking into account the non-zero phonon frequencies.
Thus the ASC avoids the adiabatic AHC divergence for the wrong reason, which explains its poor predictive capability for the more ionic materials emphasized by Fig. 2. To be clear, we do not pretend the nonquadratic effects are all absent, but the non-adiabatic effects have precedence, at least for materials with significant infrared activity, and the nonquadratic localization effects will be observed only if the electrons have the time to physically react. The shortcomings of the ASC approach are further developed in the SI Sec.V, and as a consequence of such understanding, all the results obtained for strongly infrared-active materials using the adiabatic frozenphonon supercell methodology should be questioned.
For non-infrared-active materials, the physical picture that we have outlined, namely the inability of slow electrons to follow the dynamics of fast phonons, is still present, but does not play such a crucial role: the electron-phonon interaction by itself does not diverge in the long-wavelength limit as compared to the infrared-active electron-phonon interaction, see Eq. (29), only the denominator of the Fan self-energy diverges, which nevertheless results in an integrable ZPR 24 . In such case, neglecting non-adiabatic effects, as in the ASC approach, is only one among many approximations done to obtain the ZPR.
Beyond the discovery of the predominance of non-adiabatic effects in the zero-point renormalization of the band gap for many materials, in the present large-scale first-principles study of this effect, we have established that electron-phonon interaction diminishes the band gap by 5% to 10% for materials containing light atoms like N or O (up to 0.7 eV for BeO), a decrease that cannot be ignored in accurate calculations of the gap. Our methodology, the non-adiabatic Allen-Heine-Cardona approach, has been validated by showing that, for nearly all materials for which experimental data exists, it achieves quantitative agreement (within 25%) for this property.
We have also shown that most of the discrepancies with respect to experimental data of the (arguably) best available methodology for the first-principles band-gap computation, denoted GWeh, originate from the first-principles zero-point renormalization: after including it, the average overestimation from GWeh nearly vanishes. There are some exceptions, materials in which transition metals are present, for which the addition of zero-point renormalization worsens the agreement of the band gap. For the latter materials, we believe that the GWeh approach is not accurate enough.

First-principles electronic and phonon band structures.
Calculations have been performed using ABINIT 53 with norm-conserving pseudopotentials and a plane-wave basis set. Table S1 of the SI provides calculation parameters: planewave kinetic cut-off energy, electronic wavevector sampling in the BZ, and largest phonon wavevector sampling in the BZ. For most of the materials, the GGA-PBE exchangecorrelation functional 54 has been used and the pseudopotentials have been taken from the PseudoDojo project 55 . For diamond, BN-zb, and AlN-wz, results reported here come from Ref. 24, where the LDA has been used, with other types of pseudopotentials.
The calculations have been performed at the theoretical optimized lattice parameter, except for Ge for which the gap closes at such parameter, for GaP, as at such parameter the conduction band presents unphysical quasi-degenerate valleys, and for TiO 2 , as the GGA-PBE predicted structure is unstable 56 . For these, we have used the experimental lattice parameter. The case of SrTiO 3 is specific and will be explained later.
Density-functional perturbation theory [57][58][59] has been used for the phonon frequencies, dielectric tensors, Born effective charges, effective masses, and electron-phonon matrix elements. First-principles calculations of zero-point renormalization. We first detail the method used for the AHC calculations. In the many-body perturbation theory approach, an electronic self-energy Σ appears due to the electron-phonon interaction, with Fan and Debye-Waller contributions at the lowest order of perturbation 36 : The Hartree atomic unit system is used throughout (h = m e = e = 1). An electronic state is characterized by k, its wavevector, and n, its band index, ω being the frequency. These two contributions correspond to the two diagrams presented in Fig. 1.
Approximating the electronic Green's function by its noninteracting KS-DFT counterpart without electron-phonon interaction, gives the standard result for the T = 0 K retarded Fan self-energy 36 : In this expression, contributions from phonon modes with harmonic phonon energy ω q j are summed for all branches j, and wavevectors q, in the entire Brillouin Zone (BZ). The limit for infinite number N q of wavevectors (homogeneous sampling) is implied. Contributions from transitions to electronic states k + qn | with KS-DFT electron energy ε k+qn and occupation number f k+qn (1 for valence, 0 for conduction, at T = 0 K) are summed for all bands n (valence and conduction). The H (1) q j is the self-consistent change of potential due to the q j-phonon 36 . Limit of this expression for vanishing positive η is implied. For the Debye-Waller selfenergy, Σ DW kn , we refer to Refs. 21 and 36. In the non-adiabatic AHC approach, the ZPR is obtained directly from the real part of the self-energy, Eq. (7), evaluated at ω = ε kn 36 : If the adiabatic approximation is made, the phonon frequencies ω q j are considered small with respect to eigenenergy differences in the denominator of Eq. (8) and are simply dropped, while a finite η, usually 0.1 eV, is kept. With a vanishing η, the adiabatic AHC ZPR at band edges diverges for infraredactive materials, see SI Sec.IV.A. Summing the Fan and Debye-Waller self-energies, and working also with the rigid-ion approximation for the Debye-Waller contribution delivers the non-adiabatic AHC ZPR, given explicitly in Eqs. (16) and (17)  The imaginary smearing of the denominator in the ZPR computation is 0.01eV, except for SiC, where 0.001 eV is used. Other technical details are similar to previous studies by some of ours, e.g. Refs. 24 and 25.
The dependence of the electronic structure on zero-point lattice parameter corrections is computed from where the lattice parameters {R fix i } minimize the Born-Oppenheimer energy without phonon contribution, while {R T =0 i } minimizes the free energy that includes zero-point phonon contributions, see SI Sec.II.C.
At variance with the AHC approach, in the ASC the temperature-dependent average band edges (here written for the bottom of the conduction band) are obtained from where β = k B T (with k B the Boltzmann constant), T the temperature, Z I the canonical partition function among the quantum nuclear states m with energies E m (Z I = ∑ m exp(−β E m )), and ε c m the band edge average taken over the corresponding many-body nuclear wavefunction. At zero Kelvin, this gives an instantaneous average of the band edge value over zeropoint atomic displacements, computed while the electron is NOT present in the conduction band (or hole in the valence band) thus suppressing all correlations between the phonons and the added (or removed) electron. Convergence of the calculation. As previously noted 24,25 , the sampling of phonon wavevectors in the Brillouin zone is a delicate issue, and has been thoroughly analyzed in Sec.IV.B.2 of Ref. 24. In particular, for infrared-active materials treated with the non-adiabatic effects, at the band structure extrema, a N −1/3 q convergence of the value is obtained. We have taken advantage of the knowledge acquired in Ref. 24 to accelerate the convergence by three different methodologies. In the first one, we fit the N −1/3 q behavior for grids of different sizes, and extrapolate to infinite N q . In the second one, we estimate the missing contribution to the integral around q = 0, at the lowest order, see SI Sec.I.B, using ingredients similar to those needed for the generalized Fröhlich model, except the effective masses. For ZnO and SrTiO 3 , a third correcting scheme, further refining the region around the band edge with an extremely fine grid, is used. The special case of SrTiO 3 . While phonons in most materials in the present study are well addressed within the harmonic approximation, this is not the case for SrTiO 3 . This material is found in the cubic perovskite structure at room temperature, and undergoes a transition to a tetragonal phase below 150 K, characterized by tilting of the TiO 6 octahedra 60 . Within our first-principles scheme in the adiabatic approximation, the cubic phase remains unstable with respect to tilting of the octahedra. Quantum fluctuations of the atomic positions actually play a critical role in stabilizing the cubic phase at high temperature 61 , as well as suppressing the ferroelectric phase at low temperature 62,63 .
Since SrTiO 3 is however a material for which large polaron effects are clearly identified, see e.g. Ref. 49 and references therein, we decided to study it as well. We addressed the anharmonic stabilization problem using the state-of-theart TDep methodology 64 . We used VASP molecular dynamics to generate 40 configurations in a 2x2x2 cubic cell of STO at 300 K, producing 20,000 steps with 2fs per step and sampling the 40 configurations out of the last 5000 steps. Then we computed the forces with ABINIT and performed TDep calculations with the ALAMODE code 65 . Our calculation stabilizes the acoustic phonon branches and yields a phonon band structure in good agreement with experimental data. Sources of discrepancies between experiment and theory. The anharmonic corrections to phonon frequencies are not the only reasons for potential differences between the experimental ZPR g and our non-adiabatic AHC ZPR g values. The following phenomena may also play a role: (1) the rigidion approximation 18,21 ; (2) the nonquadratic behaviour of the eigenenergies with collective displacements of the nuclei, in reference to the ASC, especially emphasized in Ref. 33; (3) the reliance on GGA-PBE eigenenergies and eigenfunctions, instead of more accurate (e.g. GW) ones, as in Ref. 30, also discussed in Ref. 35; (4) self-trapping effects, overcoming the quantum fluctuations, yielding small polarons, see e.g. Ref. 66. There is still little knowledge about each of these effects when correctly combined to predict the ZPR g beyond the AHC picture.
As an example, in Ref. 35, the difference between the ASC-PBE and the ASC-GW was argued to be only a few meV, but a more careful look at their values show that it is often bigger than 10% of the ASC-PBE. Unfortunately, the convergence of the ASC-GW results with respect to supercell size could not be convincingly achieved in Ref. 35. It remains to be seen whether a non-adiabatic AHC treatment based on GW matrix elements would differ by such relative ratio, see Sec.V.C of the SI. Altogether, it would be hard to claim more than 25% accuracy with respect to experimental data, from our non-adiabatic AHC ZPR g calculations. Together with the experimental uncertainties, this explains our choice for the 25% accuracy comparative limit used in Fig. 2.

ACKNOWLEDGMENTS
We acknowledge fruitful discussions with Y. Gillet and S. Poncé

COMPETING INTERESTS
The authors declare no competing financial or non-financial interests.

DATA AVAILABILITY
The numerical data used to create all the figures in the main text have been collected in the the Supplementary Information.

CODE AVAILABILITY
ABINIT is available under GNU General Public Licence from the ABINIT web site (http://www.abinit.org).

SUPPLEMENTARY INFORMATION
In the Supplementary Information, the tables with the numerical data used to create the figures in the main text have been collected, as well as information and comments about such data (including the corresponding bibliographical references for results that have not been computed in the present work). Experimental data for the ZPR are critically examined and discussed. The gap ZPR values are also discussed with respect to previously published values, highlighting the importance of the control of phonon wavevector grid convergence and imaginary broadening. The non-quadratic effects and GW approximation versus DFT approach are further discussed in the light of previously published data. We have evaluated Eqs. (8) and (9) of the main text for thirty materials. The list of these materials, with Materials Project IDs 67 , is presented in Table S1. It includes three non-polar materials (C, Si, Ge), and one of their combinations (SiC), nine III-V compounds (GaAs, AlAs, AlP, GaN-w, GaN-zb, AlN, BN, AlSb, GaP) among which four nitrides, ten oxides (TiO 2 , ZnO, SnO 2 , BaO, SrO, CaO, Li 2 O, MgO, SiO 2 -stishovite, BeO) and SrTiO 3 , as well as six non-oxide II-VI compounds (CdTe, CdSe, CdS, ZnS, ZnSe, ZnTe).
With this set of materials, we have (i) a large overlap with well established experimental ZPR g values, as discussed in Sec. II A, and complete coverage of the set of materials for which ASC ZPR g calculations have been done by Karsai et al. 35 ; (ii) nearly entire coverage of the set of materials for which GWeh computations have been performed, from Shishkin et al. 68 ; and we include also (iii) a dozen of binary oxyde infrared-active materials for which the trends in ZPR can be analyzed.
Table S1 also lists the parameters used in the first-principles computations. Within Density-Functional Perturbation Theory (see METHODS section), the computation of phonons at arbitrary q-points is done without any supercell calculation, thus allowing us to rely on the very fine q-point sampling grids mentioned in Table S1 at affordable CPU time cost. Still the convergence with respect to such sampling is extremely slow, and is adressed below.
B. Accuracy of the Brillouin Zone sampling in the AHC case As mentioned in the METHODS section of the main text, the convergence with respect to the q-point sampling can be improved by three different techniques: either by a linear extrapolation based on the expected scaling (N q ) −1/3 , or by the computation of the coefficient of such scaling behaviour, as explained below, or by explicitly refining the region around Γ. For ZnO and SrTiO 3 , we have used the latter, with a 192x192x192 fine grid. The second methodology is the following. According to Ref. 24, the (N q ) −1/3 behavior is associated with a q = 0 discarded (missing) part of the Brillouin zone in the Fan contribution, Eq. (8) of the main text. We denote Ω q=0 this region of the BZ, and estimate its contribution using the generalized Fröhlich electron-phonon interaction, however neglecting the electronic dispersion in this small region (unlike in the generalized Fröhlich model of the main text). The estimation of the dominant correction to the ZPR AHC c thus reads Note that the sum over degenerate states has been cancelled by the n deg denominator. The volume of the missing region is actually the BZ volume (Ω BZ = (2π) 3 /Ω 0 ) divided by N q , the number of points sampling the BZ. The equivalent sphere has a cut-off radius We then replace the integral over Ω q=0 by an integral in the sphere with radius q c , evaluate the radial part of the integral, and obtain where K av is the angular average of theq-dependent quantities in Eq. (12), namely, The correction is of opposite sign for the top of the valence band, due to the change of occupation number. This technique and the linear extrapolation technique have been tested for most materials in our list. Later (Tables S2 and  S5), we will report the values obtained with the correction from Eq. (14), except for the following materials: the noninfrared-active materials (C, Si and Ge), for which there is no such correction, and also AlN-w and BN-zb, for which the convergence study had been done in Ref. 24, based on the linear extrapolation 69 .  14) does not exactly removes the 1/N MP behaviour, although it decreases it by more than a factor of ten. The coefficient to the (N q ) −1/3 factor computed from this equation is only approximate.
Still another technique to speed up the convergence, based on Fröhlich model, but not restricted to lowest order, has been sketched in Ref. 71. However, in this reference, it has only fully been elaborated for materials with isotropic characteristics. It should be possible to develop it further using the generalized Fröhlich model introduced in the present work.
C. Comparing the AHC Brillouin Zone sampling and the ASC supercell size One might wonder why even with the above-mentioned corrections, the AHC approach needs a grid with typically N MP = 20 to obtain reasonably converged values, while for the ASC approach, results with 5x5x5 or 6x6x6 supercells, corresponding to grid N MP = 5 or N MP = 6 are apparently converged at the same level 35 . We argue that the needs of both methods are not identical, due to different physics, and the rate of convergence is quite different.
In the ASC methodology, for which the divergence of the harmonic approximation is removed by the inclusion of the nonquadratic coupling, the convergence is much faster. This is partly because, as argued in the DISCUSSION section of the main text, "for larger amplitudes, the eigenenergy behaves linearly, as the electron localizes in the lowered potential region and the minimum of the potential is linear in the amplitude of the atomic displacements.". This is at variance with the quadratic behaviour found from the AHC perturbative theory.
Let us examine a numerical proof of this difference. In Ref. 35, Karsai and coauthors quote their detailed convergence study for diamond (their Table 2) from 3x3x3 to 6x6x6 supercells, with apparent convergence within 6 meV for a ZPR around 330 meV (values: −0.278, −0.365, −0.315, and −0.321 eV). With the AHC, we checked that fluctuations of ZPRs between corresponding q-point grids are on the order of 100 meV or larger, which is in line with the results from Ref. 22 for random wavevector grids. Fluctuations becomes less than 100 meV only above 8x8x8. Thus, converged values are much easier to reach with the ASC methodology than with the AHC methodology. Still, as emphasized in the main text, the ASC avoids the adiabatic AHC divergence for the wrong physical reason, and the better convergence rate of ASC is misleading, as the converged value is not the right one.  Table S2 compares ZPR g data from experiment and from different computations. These data are used in Fig. 2 of the main text. Further computed ZPR g data are provided in Table  S3, although the latter focuses on band gaps, to be analyzed in the next section (Sec. III). Gathering experimental ZPR g data needed some care, as explained below. Also, theoretical data ought to be discussed in view of previous works. For the in-depth discussion of the theoretical data among themselves, we refer to Sec. V. We also complete the list of references in which ZPR g have been computed, that had been mentioned in the main text.
A. Experimental data: discussion As explained by Cardona and Thewalt, in relation to the Table 3 of their review Ref. 46, the experimental ZPR are obtained by two different techniques. The first one starts from measurements of the band gap energy for different samples of a material with varying isotopic content. The derivative(s) of the band gap with respect to the mass(es) are extracted, followed by extrapolation to infinite mass(es). The second one relies on measurements of the band gap energy in a large temperature range, including the low-temperature regime, but also the high-temperature regime, where a linear asymptote must be reached. In this case, the ZPR is the difference between the measured zero-temperature limit and the linear extrapolation of the high-temperature asymptote to zero temperature, as shown for Ge in Fig. 21 of Ref. 46. When both techniques agree, the experimental result can be reasonably trusted.
For the latter technique (extrapolation of temperaturedependent data), we rely on the analysis of Pässler 74,75 , that superceeds his earlier work 86 on which Ref. 46 relied. The two publications by Pässler 74,75 yield essentially equivalent results, except for AlN. For this material, we have taken the mean of the two reported values as reference. However, also in the latter case, our results are within the 25% limit. We have surveyed the literature after 2003, but were unable to find more reliable data.
For the isotopic extrapolation, to our surprise, a large fraction of the data provided in the related column of Table 3 of Ref. 46 do NOT report actual isotopic extrapolation analysis. We have thus extensively examined the primary literature, and have obtained that such isotopic extrapolation is available reliably for the following materials in our list: Ge, Si, C, CdS, ZnO and ZnS. They are clearly identified as such in table S2. Moreover, we have also realized that a non-negligible number of data available through other secondary references are misleading. We will not elaborate further on this problem, but generally speaking, we warn the reader about the unreliability of several publications in the field. The above-mentioned works of Pässler are apparently free of such problem.
Anyhow, the two techniques agree very well for Ge or C, or reasonably (within 25%) for Si and ZnS, but clearly dif-  [75] fer for CdS (about a factor of 2). Pässler 74 does not provide data for ZnO. In Fig.2 of the main text, we have taken the isotopic extrapolation result as reference, since the extrapolation of temperature-dependent data is often plagued with insufficient sampling of the high-temperature region, which is crucial for the correct extrapolation, as discussed at length in Ref. 86. For the above-mentioned CdS discrepancy, the global expected trend due to the change of mass, comparing with CdSe and CdTe, or with the Zn series, appears more reasonable from the isotopic technique than from extrapolation of temperature-dependent data. Also, the comparison of On this basis, one can infer that, globally speaking, the reliability of the experimental data is on the order of 25%, which is the reason of our use of this value in the main text. Note that Karsai et al. 35 quotes experimental data (most from secondary sources) for about half the materials they compute, but miss many available data, especially those that do not compare well with their results.

B. Theoretical AHC data: discussion
The values quoted in Tables S2 and S3 for the band gap ZPR of C, BN, and MgO, respectively -330 meV, -406 meV, and -524 meV, differ from those previously published in Ref. 23, authored by some of us, giving respectively -366 meV, -370 meV, and -341 meV. We emphasize that the same first-principles method as in our current manuscript (non-adiabatic AHC) was used to obtain these data. The differences with the present data are entirely due to the sampling of the Brillouin Zone and associated imaginary broadening factor.
In Ref. 23, the sampling of the Brillouin Zone was already emphasized as an important technical issue, especially con-sidering the numerical noise. For the purpose of computing the self-energy (which was the focus of Sec. II of Ref. 23), or computing the importance of "anharmonic effects" in the adiabatic approximation (Sec. III of Ref. 23), such numerical noise was addressed by keeping a large finite imaginary broadening (between 0.1-0.4 eV) and a fixed wavevector sampling 32x32x32. The choice of the value of the imaginary broadening was detailed in the Appendix of this paper, which was explicitly mentioned in the caption of this Table 1: "See Appendix for the values of the η parameters used". Such imaginary broadening was determined to ensure removal of the noise in the self-energy, in order to produce Fig. 1 and 2, but however did not deliver an accurate absolute value for the ZPR in Table I of Sec. II.
Indeed, in another publication of some of ours (even slightly before Ref. 23), namely Ref. 24, the proper limit of vanishing imaginary broadening with much improved wavevector sampling was performed, for C and BN (up to 125x125x125 grid). The same values for C and BN as mentioned in our current manuscript, namely -330 meV and -406 meV, were found (with one more digit, see Table VII,

C. Modification of the ZPR due to lattice parameter changes
The ZPR AHC g is computed at the fixed theoretical lattice parameters and internal coordinates minimizing the DFT energy (or fixed experimental lattice parameters in the case of Ge, GaAs and TiO 2 ). The phonon population effects on the lattice parameters, and induced change of band gap, are not taken into account in such procedure. This corresponds strictly to the harmonic approximation.
In the so-called quasi-harmonic approximation, the equilibrium lattice parameters and possibly internal coordinates are relaxed in order to minimize the free energy of the system. This makes such parameters depend on the temperature (thermal expansion), but also induces their zero-point change at T = 0 K. In turn, such changes modifies the gap, an effect to be accounted for in order to compare with the experimental ZPR g . This effect has been investigated by Garro and coworkers, Ref. 87. It is also included in several first-principles calculations of the ZPR, see e.g. Ref. 88. In the present work, we have indeed added this effect beyond the harmonic approximation: we have computed the equilibrium lattice parameters with and without zero-point motion, then evaluated the change of gap. In the first case, we have minimized the Born-Oppenheimer energy without phonon contribution, giving {R f ix i }. In the second case, the free energy, including zero-point phonon contributions computed at different lattice parameters and interpolated between them, was minimized, giving {R T =0 i }. This procedure is similar to the one described in Ref. 89, albeit extended to two dimensions (uniaxial and basal lattice parameters) in the case of wurtzite materials. The differences in band gap values for such lattice parameters form the lattice ZPR modification of the electronic structure: In such approach, one supposes that there is no crossinfluence of the modification of lattice parameters on the AHC ZPR. Also, one relies on the approximation that the internal parameters, if not determined by symmetry, are those that minimize the Born-Oppenheimer energy at the corresponding R i lattice parameters. In other words, the deviatoric thermal forces and their effect on the band gap are not included, but the deviatoric thermal stresses are taken into account, following Ref. 90. In the list of materials of Tables S2 and S3, deviatoric thermal forces only appear for wurtzite materials, that have a quite symmetric tetrahedral environment, and thus, very small expected effect.

D. A survey of zero-point renormalization computations
In complement to the references given in the introduction of the main text, we list additional ones in view of further comparisons, without pretending to be exhaustive. Refs. 15-25, 71, 88, 91-101 report calculations of the ZPR g based on the AHC approach with KS-DFT wavefunctions and eigenenergies. Refs. 18, 27-35, 92, 102-108 report calculations of the ZPR g based on the ASC approach, usually based on KS-DFT, but also sometimes based on GW, see Refs. 30, 35, and 105. As analyzed by the present generalized Fröhlich model, the ASC might be reasonably predictive for purely covalent materials (e.g. C, Si, Ge), or weakly infrared-active ones (e.g. GaAs), however, when the ZPR gFr is a significant fraction of the ZPR, the predictive capability of the ASC approaches should be questioned. Even in the covalent case, non-adiabatic corrections can be sizeable (5...15%), see Table  VII of Ref. 24. The ASC GW calculations performed until now 30,35,105 rely at most on 4x4x4 supercells, which is at the limit of convergence. Going from 4x4x4 to 5x5x5 supercells for KS-DFT indeed still brings some modification, as shown by Ref. 35. For this reason, the effect of relying on GW corrections can hardly be clarified in the present context, which focuses on the non-adiabatic effects.
In a very recent publication 109 , Engel and coworkers have proposed a new method to compute the ZPR, combining supercell calculations with an interpolation method, still in the adiabatic approximation. They applied it for 9 materials, a subset of those examined by Karsai and coworkers 35 (focusing only on covalent or weakly ionic materials). The results that they obtain do not change the overall assessment of the adiabatic approach that we present here: for some materials, the agreement with experiment is slightly improved or deteriorated.
Studies of strongly infrared-active materials using the ASC include e.g. hydrogen-containing materials HF, H 2 O, NH 3 and CH 4 in Ref. 33, where a giant electron-phonon interaction in molecular crystal and related important nonquadratic coupling had been investigated, the study of ice in Ref. 34, and the study of LiF, MgO and TiO 2 in Ref. 105, and some of the results from Ref. 35 for IR-active materials.
In addition to these references, the electron-phonon interaction has been directly incorporated in GW calculations in Refs. 110-112. These are further commented upon in the next section. Table S4 presents KS-DFT ASC data that were not reported in Table S2, as the latter only mentioned such data from Ref. 35. This table S4 still focuses on our set of materials, for which AHC data are available. For the indirect gap of Si and of SiC, as well as both the direct and indirect gaps of C, the different ASC calculations agree very well. The agreement is less good for GaAs and the direct gap of Si. Anyhow, these data confirm that AHC and ASC calculations of the fundamental band gap reasonably agree for covalent materials, but strongly disagree for strongly IR-active materials (SiC-zb, TiO 2 -t, MgO-rs).  Table S3 compares gap values from experiment and from different computations. These data are used in Fig. 3. The set of 12 materials in this table includes all those for which self-consistent GW calculations with electron-hole corrections have been performed in Ref. 8, except Ne, LiF and Ar. We did not include them in our study for the following reasons. Ne and Ar are strongly affected by weak van der Waals interactions. Much more than the other materials, they will be modified by the zero-point lattice corrections, discussed earlier. For LiF, as mentioned previously, the estimated value of α show that the present perturbative treatment is likely invalid. For all three materials, excitonic effects are strong, and should likely be included for meaningful comparison with experimental data.
Suppplementary Table III 85. The latter reference also presents E G 0 W 0 g values for many materials, including more than half of those reported in Table S3. Many other E G 0 W 0 g values for those materials have been published. The spread in such values for one material can be as large as 0.5 eV. A comparative study of published E G 0 W 0 g values is out-of-scope of the present publication, but is presented in Ref. 45. The values E G 0 W 0 g in Table S3 have the goal to show the typical size of the E G 0 W 0 g to E GWeh g correction, and point out that the ZPR correction is of the same order of magnitude for materials with light nuclei.
In addition to the AHC and ASC methodologies, there has also been attempts to incorporate directly the electron-phonon interaction inside the GW approach [110][111][112] to obtain modified gap values or ZPR g . The first publication on the subject, by Botti and Marques 110 was later shown by Lambrecht, Bhandari and van Schilfgaarde 111 to ignore the range of the Fröhlich interaction, yielding an erroneous integration over the Brillouin Zone. Even in the latter work, there is a factor of 2 difference for the Fröhlich α compared to the usual definition, used in the present work. Their ZPR g results (Table I of Ref. 111 and 112) for MgO, GaN-zb and SrTiO 3 , resp. (-219 meV, -67 meV, -404 meV) can be compared with ours (-517 meV, -171 meV, -290 meV), see Table S5. They very clearly differ, with strong underestimation with respect to ours for MgO, GaN-zb, and strong overestimation for SrTiO 3 . For SrTiO 3 , our results agree better with the experimental data (-336 meV), see Table S2. For MgO and GaN-zb, a direct measurement of ZPR g is not available. Still, experimental value for the other polymorph of GaN, namely GaN-w (-180 meV) is reported in Table S2 and matches very well our ZPR AHC g value (-189 meV) for GaN-w.

IV. FROM THE ADIABATIC AHC APPROACH TO THE GENERALIZED FRÖHLICH HAMILTONIAN
In order to help understand why the generalized Fröhlich Hamiltonian captures correctly the physics of the Fan dia-gram, at variance with the adiabatic approximation, we summarize first the discussion of the AHC case given in Ref. 24. Then, we derive the generalized Fröhlich Hamiltonian from Eq. (8) of the main text.

A. The divergence of the adiabatic AHC for IR-active materials
Under the adiabatic approximation, the phonon frequencies ω q j are neglected when compared to differences between electronic eigenenergies. In other words, one supposes that the electron (or hole) responds instantaneously to the atomic vibrations, even for those electronic transitions with a characteristic time scale (inversely proportional to the eigenenergy difference) larger than the phonon oscillation period. This simplification is of course questionable for intraband transitions with small momentum transfer, as shown in Ref. 24.
Suppressing ω q j in the denominators of Eq. (8) , and evaluating the AHC ZPR according to Eqs. (7) and (9) yields For band edges, in particular the top of the valence band or the bottom of the conduction band, and considering at present non-degenerate bands, the intraband (n = n) contribution when η = 0 has two different types of divergences for q → 0, reinforcing each others: (i) for IR-active materials, the electron-phonon interaction matrix elements of H (1) q j diverge like δ nn /q, while (ii) the denominator behaves like q 2 /2m * q . In these expressions, q is the norm of the q wavevector, and m * q is the effective mass along direction q, represented by the unit vectorq, with q = q ·q. The effective mass is the inverse of the second derivative of the eigenenergy with respect to the wavevector alongq at the band extremum.
Thus, the integrand at small q behaves like 1/q 4 . Focusing on the ZPR adiab c of the bottom of the conduction band (the handling of the top of the valence band is similar) one finds, in spherical coordinates, and introducing a spherical cut-off q c whose limit q c → 0 has to be taken: where f adiab (qq) is a smooth function of q (it does not diverge, but also does not tend to zero for q → 0). As announced 24 , ZPR adiab c (q c < q) diverges like 1/q c when q c → 0. By contrast, if phonon frequencies are not suppressed in Eq. (8), which is the non-adiabatic case, one falls back to an integral of the type (here considered for one isotropic phonon branch only, and one non-degenerate isotropic electronic band): where the O(q 2 ) contribution in the denominator comes from the curvature of ω q with respect to q, usually much smaller than the curvature of the electronic eigenenergies governed by m * q , explicitly mentioned in Eq. (19). ZPR non−adiab c (q c < q) does not diverge for q c → 0, provided ω 0 is non-zero.
Acoustic modes should not be forgotten. Their eigenfrequency ω 0 tends to zero in the q → 0 limit. However, the corresponding first-order electron-phonon matrix element does not diverge like 1/q, and also, it has been shown by Allen and Heine 9 that, thanks to translational symmetry, the contribution to the Fan term from acoustic modes is exactly cancelled by the Debye-Waller term, in the q → 0 limit. Thus no divergence arises from them.
To summarize, finite phonon frequencies at q = 0 remove the AHC divergence in the non-adiabatic case.

B. The generalized Fröhlich Hamiltonian
Instead of suppressing phonon frequencies appearing in the Fan self-energy Eq.(8), like in the adiabatic approximation, a radically different strategy is followed in the generalized Fröhlich Hamiltonian: one uses only parameters at q = 0, including crucially the phonon frequencies, and extends them to the whole q-space, using simple behavior for medium and large-q (i.e. parabolic electronic dispersion and no phonon dispersion). This corresponds to macroscopic electrostatic treatment and a continuum space hypotheses.
The model proposed by Fröhlich 37 in 1954 is the most simplified version of such hypotheses: one non-dispersive LO phonon branch, one electronic parabolic band governed by a single, isotropic effective mass, and electron-phonon interaction originating from the isotropic macroscopic dielectric interaction between the IR-active LO mode and the charged electronic carrier. No Debye-Waller contribution is included. Spin might be ignored, as only one electron is considered, without spin-phonon interaction. Explicitly, the Fröhlich Hamiltonian reads 37,38,114,115 , withĤ Fr H Fr whereâ + q andâ q are phonon creation and annihilation operators,ĉ + q andĉ q are electron creation and annihilation operators, and the electron-phonon coupling parameter is given by where V BvK is the Born-von Karman normalisation volume corresponding to the k and q samplings, q is the norm of q andq is its direction, q = q ·q, ε ∞ is the optical dielectric constant, and ε 0 is the low-frequency dielectric constant. The sums over k and q run over the whole reciprocal space. At the lowest order of perturbation theory, this yields the well-known result 37,38,116 for the polaron binding energy (that we denote ZPR Fr ): Accurate diagrammatic Monte Carlo simulations 39 have demonstrated reasonable validity of this lowest-order perturbation result in a range that extends to about α ≈ 8. Beyond this, electronic self-trapping becomes too important, and the lowest-order result deviates strongly from the exact results. Multiphonon generalisations of the Fröhlich model Hamiltonian have been considered in several occasions, but only in the cubic case, while the anisotropy of the electronic effective mass has been accounted for in an approximate way, and not altogether with the degeneracy of band edges, that has either been treated approximately, or for a spherical model, without warping, or only considering one phonon branch see e.g. Refs. 47-50, 71, and 112. Our generalized Fröhlich model includes all the features of real materials altogether: multiphonon, anisotropic, degenerate band extrema. It is derived by considering the full Eq. (8) and applying the same ansatz as the original Fröhlich model: we treat as exactly as possible the region in the Brillouin zone where the phonon frequency is large with respect to the eigenenergy differences, and, coherently, suppress the contributions from bands with a large energy difference with respect to phonon frequencies (interband contributions), keeping only intraband contributions. Like for the Fröhlich Hamiltonian, we use only parameters that describe the functions in the Brillouin zone around q = 0, and extend them to the whole reciprocal space, using the same simple behaviors for medium and large-q: parabolic electronic dispersion and no phonon dispersion. We also retain only the macroscopic Fröhlich electronphonon interaction, G = 0 component, following Vogl 51 . This is also coherent with the first-principles It has been shown also by AHC that the contribution to the Fan term from acoustic modes, for which the frequency vanishes when q → 0, is exactly cancelled by the Debye-Waller term, close to q = 0.
Having highlighted the basic ideas of the derivation of the generalized Fröhlich Hamiltonian, in what follows, we write down the mathematical expressions. We focus on the conduction band expression, ZPR c , associated to one of the possibly degenerate states, labelled c, with energy ε c and with wavevector k 0 . Similar expressions for ZPR v can be derived. For the sake of simplicity, we will suppose that k 0 is the Brillouin zone center Γ, but the equations that follow can be straightforwardly adapted to another value of k 0 . We also introduce k = k + q.
Under the above hypotheses, the matrix elements of H with Ω 0 the primitive cell volume, V BvK the Born-von-Karman normalisation volume corresponding to the k and q samplings, and the scalar product k , n |k, n P is computed from the periodic part of the Bloch functions. Sums over k and q of the type ∑ q f (q)/V BvK are to be replaced by Ω 0 /(2π) 3 d 3 q f (q) in the macroscopic limit in what follows. We have used the notation ω j0 (q) for the q → 0 limit of the j-phonon branch (ω j0 (q) = lim q→0 ω j(qq) ), while the related mode-polarity vectors p j (q), see Eq. (41) of Ref. 118, are defined from the Born effective charges and the eigendisplacements of the phonon mode, summing over the atoms labelled by κ, and the Cartesian direction of displacement α, Similarly, the optic dielectric constant alongq is obtained from the optic dielectric tensor by summing over wavevector components, ε ∞ (q) = ∑ αβqα ε ∞ αβq β , as in Eq. (56) of Ref. 58. The mode-polarity vector vanishes for all non-IRactive modes. Thus naturally only the LO modes contribute to the ZPR gFr in this generalized Fröhlich approach (for acoustic modes and TO modes, the mode-polarity vector vanishes). Still, the directional dependency of the phonon frequencies, of the mode-polarity vectors, and of the dielectric constant, is correctly taken into account.
We then restrict ourselves to the bands that are degenerate with c, and choose ak-independent complete basis for these states at k 0 = Γ, denoted |Γ, m . We also define thekdependent overlap between states m at Γ and states n alongk in the k → 0 limit: This definition supposes that the phase of the |kk, n states is chosen continuously as a function of k, for vanishing k along thek direction. The phase can always be chosen to be continuous, and this choice will have no bearing on the results obtained later. Still, it does not imply any continuity hypothesis for the phase between states n for differentk directions. The matrix s nm (k) is unitary. Note that the |Γ, m states are fixed before thek direction is known. The coefficients s nm (k) can be computed from first principles following the perturbative treatment of Refs. 117 and 59, which is the first-principles equivalent of the k.p approach, restricted to the subspace spanned by the bands that connect to the degenerate states. It delivers Luttinger-Kohn parameters D αβ j j , 119 that allows one to find s nm (k) for allk. The Luttinger-Kohn parameters can also be determined from experimental measurements. The s matrices can be thought as being the adaptation of spherical harmonics to the specific symmetry representation to which the degenerate wavefunctions belong.
With these definitions, we write the generalized Fröhlich electron-phonon interaction as This expression depends on the directionsk,k , andq, but not explicitly on their norm, except for the 1 q factor. We focus then on the denominators appearing in the selfenergy expression Eq. (8), evaluated at ω = ε c , following Eq. (9). Again we consider only the intraband terms, for bands n degenerate with c in the q → 0 limit. These bands, similarly to the c state, are unoccupied, hence f k n = 0. Thus the numerator of the second term in Eq. (8) vanishes. We find: See Refs. 117 and 59 for the detailed treatment of m * n (k ), the directional effective mass associated to one electronic band among a set of degenerate bands at an extremum.

C. Lowest-order perturbation treatment
Given the Hamiltonian, the lowest-order perturbation ZPR value can be inferred by integration over the radial q coordinate. One obtains for the conduction bands: Let us denote by n deg the degeneracy of the band edge. In the non-degenerate case (n deg = 1), only n = c has to be considered, then obviously |s n c (q)| = 1 independently ofq. Grouping terms adequately, we obtain ZPR gFr The link can then be established with Eq. (25) in case of isotropic effective mass, isotropic dielectric tensor, and single LO phonon branch. Let us further examine the degenerate case (n deg > 1). We suppose that the degeneracy is not accidental, but due to symmetry reasons. In this case, we can exploit the symmetry to show that the computation |s n c (q)| 2 can be avoided, as this factor can be replaced by 1/n deg in the integral. We will present the demonstration for the case of a three-fold degeneracy, that arises due to a cubic space group. This demonstration could be made general thanks to group theory.
At the bottom of the conduction band, located at the k 0 wavevector, but that for convenience we again choose to be Γ, we set up a basis of three orthonormalized eigenfunctions, denoted {|X , |Y , |Z }, that form an irreducible representation of the symmetry group. The eigenfunction |Γ, c for which we compute the ZPR gFr g , can be decomposed in this basis, with coefficients (u cX , u cY , u cZ ): The eigenfunctions |qq, n over which the n sum andq integral are done in Eq. (31) can be decomposed as well in this basis: |qq, n = s n X (q)|X + s n Y (q)|Y + s n Z (q)|Z Thus Eq. (31) becomes where The second-rank tensor T mm does not depend on the state c.
However the same ZPR gFr c must be obtained from Eq. (35) for any symmetrically equivalent state c . Hence, the tensor T mm must be invariant under all the symmetry operations of the group, which means where δ is the Kronecker symbol. t can be evaluated thanks to the trace of T mm , yielding ZPR gFr   Table S5 presents the ZPR for the bottom of the conduction band and the top of the valence band from the present firstprinciples calculations: either based on the AHC approach, or from feeding the few relevant first-principles parameters in the generalized Fröhlich model. These data are used in Fig. 4. Table S6 presents the parameters computed from first principles, fed in the (generalized) Fröhlich model, only for materials with one phonon branch, with cubic symmetry and edges with non-degenerate electronic state, though including both We discuss how a few-parameter Fröhlich Hamiltonian, of the type we have introduced, can be obtained from experimental data, thus constituting a few-parameter model of the ZPR of real materials.
In the Fröhlich model 37,38,114,115 four parameters are needed: the LO phonon frequency ω LO , the electronic dielectric constant ε ∞ , the low-frequency dielectric constant ε 0 , and the effective mass m * . These four parameters combine to give the α parameter mentioned in Eq. (24). The Fröhlich Hamiltonian can be expressed in the natural units of electronic effective mass and phonon frequency, to make appear α as the single parameter defining the Fröhlich Hamiltonian 114,116 . Despite its simplicity, this Fröhlich Hamiltonian has resisted analytical treatments.
One can consider the generalization of such Fröhlich Hamiltonian to arbitrary electronic dispersion, phonon dispersion or electron-phonon interaction, as mentioned in the introduction of Ref. 41. However, one has a model for a real material only if the parameters describing such dispersions and such interaction are available. Ideally, there should be as few parameters as possible if the goal is to extract the relevant physics.
In this respect, the Hamiltonian in Eqs. (1)-(5) is the proper generalization of Eqs. (19)- (23) to all band extrema (degenerate or not, isotropic or anisotropic) in all crystal types (cubic or not), because it can be parameterized with only a few numbers, and all such parameters can be found, in principle, from experiment, even if, in the present work, they were deduced from first principles.
Indeed, the electronic part, Eq. (2), is fully determined if the so-called Luttinger-Kohn parameters 119 are available. Taking the example of the top of the valence band of a cubic material, with three-fold degeneracy, only three Luttinger-Kohn parameters, called A, B, and C in Ref. 119 are needed to obtain the direction-dependent effective masses, m * n (k), of the three electronic bands, for all directions k. Similarly, the directiondependent phonon frequencies ω j0 (q) and their mode-polarity vectors p γ, j (q) can be obtained from neutron and infra-red spectra measurements and parametrized with a few numbers. For a cubic material with two atoms per unit cell, there is only one LO frequency, with mode-polarity vector directly determined from the LO-TO splitting. The optic dielectric constant along any direction q comes from the knowledge of a 3x3 symmetric tensor, easily obtained by optical measurements.
Conversely, in Ref. 52, by Verdi and Giustino, a model for the electron-phonon interaction (or vertex) is proposed (at variance with our model of the full Hamiltonian, with a simplified vertex). It is used for a very different purpose. In this very nice study, Verdi and Giustino used the Fröhlich idea and the microscopic theory developed by Vogl 51 to model the nonanalytic behaviour of the electron-phonon vertex in the q → 0 limit. This approach was instrumental to make the electronphonon interaction amenable to Fourier interpolation in an attempt to diminish the computational resources needed to capture the long-range behaviour of the scattering potentials and perform mobility calculations using Wannier functions. However, their work does not propose a ZPR model or Fröhlich Hamiltonian per se. Indeed, the electronic dispersion in Verdi and Giustino is the full first-principles dispersion in the entire Brillouin Zone and similarly for the phonon dispersion, unlike the one of the original Fröhlich approach (one parameter for each of these dispersions), and the present one, which relies on few parameters, the number of which depends on the crystal type and symmetries, and on the Luttinger-Kohn parameters.
where β = k B T (with k B the Boltzmann constant), T the temperature, Z I the canonical partition function among the quantum nuclear states m with energies E m (Z I = ∑ m exp(−β E m )), and ε c m the band edge average taken over the corresponding many-body nuclear wavefunction. See e.g. Eq. (4) of Ref. 120. Such approach allows one to perform nonperturbative calculations, and thus includes effects beyond the adiabatic AHC approach, as emphasized in Sec.XI.A.2 of Ref. 36. It also allows one to use as starting point the GW eigenenergies, instead of the DFT ones, as proposed and used by some of us 30 , and further exploited by Monserrat 105 and Karsai et al. 35 .
As detailed in Ref. 23, and illustrated in Fig. 4 of this reference, one can dissociate the harmonic approximation for the nuclear positions from the (non)quadratic behaviour of the eigenenergies as a function of the nuclear positions. While phonon anharmonicities can be safely ignored for most of the materials in the present context (except SrTiO 3 , see the METHODS section of the main text, the (non)quadratic behaviour might very important in the ASC formalism. Within the harmonic approximation, one finds Eq.
In this equation, z denotes the ensemble of all phonon coordinate amplitudes, while n denotes the ensemble of all phonon occupation numbers. As such, this expression does not diverge for IR-active materials, thanks to the sampling over coordinate amplitudes and the ε c (z) asymptotic linear behaviour, see Fig. 4 of Ref. 23 or Fig. 2 of Ref. 33. By contrast, making the quadratic approximation for the ε c (z) yields the AHC result for the Fan term, and hence, a divergence: the curvature becomes indeed infinite for IR-active longwavelength phonons. An extreme nonquadratic behavior is also observed in the degenerate case, as illustrated in Fig. 1 of Ref. 35. The connection with the adiabatic AHC result has been analyzed in detail in Ref. 21. However, obviously, Eq. (40) as well as Eq. (41) are based on the adiabatic approximation: the quantum nuclear state is not responding to the presence of an added or removed electron. Instead, Eq. (40) should have reflected such dependence to be correct. The correlation between electron and nuclei is not present in Eq. (40): the electron is affected by the phonons, but the phonons are incorrectly not affected by the addition of an electron. For IR-active materials, as explained in the main text, the lack of such correlation has important consequences for the predictive capability of the method.
As another consequence of the time-dependent atomic motion, the nonquadratic behaviour, present in the adiabatic treatment, looses its importance. The behaviour observed for degenerate electronic state, in the adiabatic approximation, exemplified by Fig. 1 of Ref. 35 might nevertheless still be relevant in the context of the Jahn-Teller effect. Anyhow, the transition between states, due to the first-order electron-phonon interaction, always involves a modification of the energy by the emission or absorption of a phonon. Hence, the corresponding initial and final energies are not identical.

B. The non-quadratic effect for IR-inactive materials
Coming back to the C-diamond case, the "anharmonic effect" in the adiabatic approximation is argued to be large in the above-mentioned Sec.III of Ref. 23, that some of us coauthored. At present, we call it more properly "non-quadratic effect" following Ref. 33. However, this study was performed within the adiabatic approximation. In the case of infrared active materials, without such "non-quadratic effects", the ZPR is infinite in the adiabatic approximation. The non-quadratic effects succeed to bring the adiabatic ZPR to a finite value, i.e. it has "infinite" impact. For C, without such "non-quadratic effects", the ZPR is finite in the adiabatic approximation, but nevertheless comes from the integration of a diverging integrand. Then the non-quadratic effects modify the adiabatic harmonic value by 40% by suppressing this diverging integrand. Still, despite yielding a finite value, the adiabatic approximation is incorrect to start with: if atomic motion effects are taken into account (and they should be taken into account from the very start, since the physics is the one of electrons not being able to follow the atomic motions), there is no diverging integrand that the non-quadratic effects would reduce afterwards. Thus, the crucial role of atom dynamics cannot be ignored even in diamond. We actually think that the estimation of non-quadratic effects in our publication Ref. 23 or in Ref. 33 cannot be transposed to the non-adiabatic treatment, even for the case of C. The non-quadratic effects in Ref. 23 or Ref. 33 remove a non-existing divergence, either global (IRactive materials), or integrable (non-IR-active materials), spuriously introduced by the adiabatic approximation.

C. GW approximation within the ASC: discussion
In such calculation of the ZPR of diamond, Ref. 30, we have obtained a 40% change of direct gap ZPR when the quantities needed in the ASC approach were obtained from GW instead of DFT. This effect appears to be comparable to the difference between AHC and experimental values reported in the present manuscript, while (as for the non-quadratic effect) C is an IR-inactive material.
In Ref. 35, a similar change was reproduced for the direct gap, but the change for the indirect (true) gap of diamond was much smaller, about 23%. For the other materials studied in Ref. 35, using the ASC approach, the ratio between GW and DFT, systematically obtained, was usually closer to unity, but not always. One can wonder whether such difference between GW and DFT would not also be obtained within the non-adiabatic AHC approach.
The estimation of such an effect within the non-adiabatic treatment is at present out of reach. Actually, as already argued in the previous sections, one can hardly justify transposing the conclusions obtained in the ASC approach to the nonadiabatic AHC approach, because the ASC approach relies on the non-quadratic coupling to eliminate the divergence of the ZPR, while for most materials, we believe such elimination to originate from the coupled dynamics of electron and phonons. In view of the agreement obtained in Fig. 2 for non-adiabatic AHC+DFT, we can simply infer that the DFT to GW effect is apparently smaller in the non-adiabatic AHC case than in the ASC case.
Also, let us emphasize that such a DFT to GW effect can hardly equal the non-adiabatic effect reported in the present manuscript. Indeed, for IR-active materials, such effect suppresses an infinity and replaces it by a finite value, which is far more important than a 40% modification. Even if the non-quadratic effect is taken into account, Fig. 2 shows that the difference between ASC and AHC is not merely 40% for many materials, but can be as large as a factor of 2 or 3.
So, we do not deny that other effects, like GW, could be important at the level of the 25% agreement that we show in Fig.  2. However, the non-adiabatic effect is not simply a 25% modification, it reduces a divergent quantity (possibly integrable) to a finite value.