Fully Anharmonic, Non-Perturbative Theory of Vibronically Renormalized Electronic Band Structures

We develop a first-principles approach for the treatment of vibronic interactions in solids that overcomes the main limitations of state-of-the-art electron-phonon coupling formalisms. In particular, anharmonic effects in the nuclear dynamics are accounted to all orders via ab initio molecular dynamics simulations. This non-perturbative, self-consistent approach evaluates the response of the wave functions along the computed anharmonic trajectory; thus it fully considers the coupling between nuclear and electronic degrees of freedom. We validate and demonstrate the merits of the concept by calculating temperature-dependent spectral functions and band gaps for silicon and the cubic perovskite SrTiO3, a strongly anharmonic material featuring soft modes. In the latter case, our approach reveals that anharmonicity and higher-order vibronic couplings can contribute substantially to the electronic-structure at finite-temperatures, noticeably affecting macroscopic properties, such as absorption coefficients as well as thermal and electrical conductivities.

Electronic band structures are a fundamental concept in material science used to qualitatively understand and quantitatively assess optical and electronic properties of materials, e.g., charge carrier mobilities and absorption spectra of semiconductors. Over the last decade, two pivotal advancements have paved the way towards predictive, quantitative ab initio calculations of electronic band structures: improvements in the treatment of electronic exchange and correlation [1,2] and the inclusion of electron-phonon interactions via perturbative manybody formalisms based on the Allen-Heine theory [3]. The latter approach has been widely used to calculate temperature-dependent effects on the electronic structure stemming from the nuclear motion [4][5][6][7][8][9][10][11][12][13][14][15][16][17][18]. However, such perturbative calculations rely on two approximations. a) The nuclear motion is approximated in a harmonic model which is equivalent to the concept of phonons and b) the vibronic interaction between electronic and nuclear degrees of freedom is treated by perturbation theory in terms of electron-phonon coupling. In both approximations, interactions at finite temperatures T are thus described via truncated Taylor expansions, using derivatives computed at the static equilibrium geometry, i.e., for the total energy minimum corresponding to the atomic geometry R eq obtained in the classical T = 0 K limit. Clearly, both these approximations are problematic whenever large deviations from R eq occur, e.g., at elevated temperatures and for soft bonded atoms. Despite several efforts to mitigate either one of these approximations [19][20][21][22][23][24][25][26][27], a consistent computational approach that accounts on equal footing for both anharmonic effects in the nuclear motion and the full vibronic coupling is still lacking.
In this work, we fill this gap by deriving a fully anharmonic, non-perturbative first-principles theory of vibronic coupling. We demonstrate the validity of the approach and of its implementation in the all-electron, numeric atomic orbitals code FHI-aims [28] by comput-ing temperature-dependent spectral functions and band gaps for silicon and cubic SrTiO 3 . In the former, largely harmonic case, our approach reproduces perturbationtheory data, thus confirming that anharmonic effects and higher-order vibronic couplings are negligible in silicon. However, this is not the case for the perovskite SrTiO 3 . Besides confirming the role of strong anharmonic effects in SrTiO 3 [29,30], we also show that the anharmonic motion results in a breakdown of the perturbative electron-phonon coupling model. Accounting for anharmonic, higher-order vibronic coupling (AViC) is thus essential to obtain the correct temperature dependence of the electronic-structure in this material.
In the following, the energy R l of the electronic state |ψ R l is obtained by solving the Schrödinger equation where H R el is the electronic Hamiltonian of the system at the atomic geometry R. This may be a Kohn-Sham Hamiltonian with a certain exchangecorrelation functional. For readability, we use the generalized index l to indicate both the band index n and the wave vector k. The temperature dependence of R l is evaluated within the Born-Oppenheimer approximation via the canonical ensemble average at temperature T : Here, k B is the Boltzmann constant, Z the canonical partition function, P the momenta of the nuclei, and E(R, P ) the total energy of the combined electronic and nuclear system. For the evaluation of Eq. (1), the stateof-the-art formalism [31] resorts to the two perturbative approximations mentioned above. When the harmonic approximation to the potential energy surface ( which allows for a straightforward evaluation of the phase-space integral [23-25, 32, 33]. When the dependence of the electronic states on the nuclear motion is truncated up to second order in the atomic displacements R l ≈ pt,R l , then the ensemble average in Eq. (2) yields the perturbative Allen-Heine energies pt,R l ha T . In this work, we rely on neither of the two approximations. First, ab initio molecular dynamics (aiMD) trajectories with length t 0 are used to evaluate the canonical ensemble average in Eq. (1) as time (t) average This accounts for the full anharmonicity of the PES. Second, the dependence of the electronic eigenenergies on the nuclear positions is explicitly evaluated by solving H at each aiMD step R(t). All orders of coupling between electronic and nuclear degrees of freedom are included by these means. This involves reexpanding in terms of the wave functions at equilibrium |ψ eq m . With that, one obtains: In this form, it is evident that Eq. (5) not only incorporates the first non-vanishing derivatives of H

R(t) el
− H eq el as perturbative formalisms, but all orders. Similarly, all orders of couplings with the nuclear motion -not just quadratic terms -are captured via the coefficients p ml , which describe the intricate R(t)-dependence of the wave functions along the aiMD. Accordingly, all orders of AViC are statistically captured by these means. Our approach, named stAViC in the following, is thus valid even when the (harmonic) phonon ansatz is inappropriate.
Note that the practical evaluation of Eq. (3) with aiMD corresponds to the classical, high-temperature limit. For the purpose of this work, this is a justified approximation, as substantiated below. Furthermore, aiMD simulations have to be performed in rather large supercells in order to capture the full vibrational spectrum. The electronic energies l = N K and wave functions ψ l = ψ N K obtained in the aiMD are associated to band indices N and wave vectors K spanning a reduced Brillouin zone (BZ), where capital letters are chosen to indicate supercell quantities.
of silicon as function of temperature obtained via the stAViC ∆ g MD T (orange) and via the non-perturbative harmonic approach: ∆ g ha-cl T (grey) and ∆ g ha-qm T (red). All calculations were performed using DFT-LDA and 6×6×6 supercells containing 432 atoms. Perturbative harmonic calculations (blue, [12]) and experimental data (black, [34]) are shown as well.
This "BZ folding" hinders a direct comparison with the energy levels computed in the static limit at R eq , for which the wave vectors k of eq l = eq nk and ψ eq l = ψ eq nk span the full fundamental BZ. To obtain a meaningful band structure in the fundamental BZ also for the aiMD case, the expansion coefficients introduced in Eq. (4) are used to "unfold" the states ψ N K . To this aim, we consider the spectral function expressed in the Lehman representation [38]: Compared to Eqs. (4)-(5), in which the perturbed eigenvalue N K is obtained by expressing its wave function as superposition of equilibrium states ψ eq nk , Eq. (6) reflects the inverse relationship: Each perturbed eigenvalue N K contributes to all states nk in the fundamental BZ, whereby p determines the strength of this contribution.
For a given aiMD configuration R(t), we obtain the momentum-resolved spectral function A nk (E) by summing over n in Eq. (6): The spectral weight P  [36]; the respective band gap in the static limit (3.568 eV) was determined via linear regression [37] from the high T > 800 K data. equilibrium states with wave vector k [39,40]. The momentum resolved spectral function in thermodynamic equilibrium A R(t) k (E) T is eventually obtained by computing the thermodynamic average of the A R(t) k (E) along R(t) using Eq. (3). Quasi-particle peaks and thus band gaps g T are extracted from A R(t) k (E) T by scanning over the energy axis. Details of the implementation can be found in the Suppl. Material.
As a validation of our approach, Fig. 1 shows the temperature-dependence of the band gap renormalization ∆ g T = g T − g ha-cl 0K of bulk silicon (DFT-LDA). Our stAViC calculations ∆ g MD T are in excellent agreement with reference data ∆ pt g ha−qm T obtained with the perturbative, harmonic formalism [12] for T 400 K. The discrepancies between the two approaches at lower temperatures are exclusively caused by quantum-nuclear effects not captured in aiMD. In Fig. 1, this is demonstrated by comparing non-perturbative, harmonic data obtained by evaluating Eq. (2) with Monte Carlo sampling [23,41] using classical ∆ g ha-cl T and quantummechanical ∆ g ha-qm T statistics. In both cases, anharmonic effects are thus neglected, while higher-order vibronic couplings are included via Eq. (5). The fact that the anharmonic ∆ g MD T and the harmonic approach ∆ g ha-cl T almost coincide in the classical limit proves that anharmonic effects are indeed negligible for silicon. Similarly, higher-order vibronic couplings are negligible here, given that the non-perturbative ∆ g ha-qm T and the perturbative data ∆ pt g ha-qm T follow closely each other.
For other materials, however, AViCs are often not negligible, as we demonstrate for the perovskite SrTiO 3 . At T = 0 K, this material exhibits a tetragonal structure (c/a = 0.998, space group I4/mcm), in which the individual tetrahedra are slightly tilted with respect to each other [43]. Above 105 K [44] and up to its melting point at 2300 K [36], SrTiO 3 exhibits a cubic high-symmetry structure (space group P m3m), in which all tetrahedra appear to be aligned as shown in Fig. 2(a). This cubic structure does not correspond to a minimum, but to a saddle point of the PES and thus features imaginary phonon frequencies. Even in the cubic lattice (c/a = 1), the tetrahedra favor a tilted arrangement in the static limit, corresponding to the minima in Fig. 2(a). Thermodynamic hopping between these wells results, on average, in an apparent alignment of the tetrahedra, in close analogy to other vibrationally-stabilized materials [45][46][47][48]. Perturbative approaches cannot capture this complex dynamics: If the saddle point with aligned tetrahedra is chosen as the static equilibrium R eq , phonon modes with imaginary frequencies have to be "frozen in" [49] and their coupling to the electronic-structure is neglected. If one of the minima with tilted tetrahedra is chosen as R eq , both the harmonic approximation for the PES and the parabolic electron-phonon model be- (E) T of cubic SrTiO3 along Γ-R (128 k-points) for T = 300 and 1200 K calculated non-perturbatively on the harmonic (left) and anharmonic PES (right) using DFT-PBE and 30 uncorrelated geometries in a 5 × 5 × 5 supercell containing 625 atoms. The band structure in static equilibrium is shown as white lines. The respective densities of states (DOS) at high symmetry point R is shown on the right as magenta (harmonic) and black (aiMD).
come not only inaccurate, but even qualitatively wrong at elevated temperatures, at which multiple minima are explored, as exemplified by the parabolic fits shown in Fig. 2(a). This breakdown of the harmonic model and of the perturbative electron-phonon coupling model has a direct impact on the calculated thermodynamic properties of SrTiO 3 .
The temperature-dependence of the band gap renormalization of SrTiO 3 is shown in Fig. 2(b). All calculations were performed at the PBE level and van-der-Waals interactions were included using the Tkatchenko-Scheffler method [50]. Corrections [51] for long-range polar effects [52] that are not fully captured within the finite aiMD supercells are included. Thermal lattice expansion and the associated, non-negligible band-gap opening of, e.g., 154 meV at 1000 K, are also accounted for non-perturbatively. All these aspects are detailed in the Suppl. Material. As discussed for Fig. 1, the fact that ∆ g ha-cl T and ∆ g ha-qm T become comparable for T > 500 K, implies that the use of classical aiMD is justified in this regime. In contrast to silicon, distinct deviations between harmonic ∆ g ha-cl T and anhar-monic ∆ g MD T data are observed for SrTiO 3 , leading to an additional renormalization in stAViC as large as 147 meV at 600 K and 260 meV at 1200 K. With respect to the perturbative, harmonic data, this corresponds to a remarkable increase of 18 % and 27 %, respectively. With respect to experiment [36], stAViC improves the agreement significantly and quantitatively reproduces the measured high-temperature slope.
More insights can be obtained by comparing the spectral functions calculated by sampling the harmonic and anharmonic PES, as done in Fig. 3 for two different temperatures. At 300 K, the two approaches yield very similar results. At 1200 K, however, the spectral functions exhibit a strikingly different behaviour, in particular at the R-point. Qualitatively, this can rationalized from the fact that the most anharmonic, imaginary phonon frequencies also occur at R. Note that the valence band maximum (VBM) lies at the R-point as well; its different upshift in the two approaches is thus the root for the different temperature-dependence of the band gaps observed in Fig. 2 (b). Furthermore, the spectral functions in Fig. 3 reveal that not only the quasi-particle peaks, but also the curvatures (effective masses) and linewidths are substantially altered by AViCs. Clearly, this calls for more systematic investigations along these lines.
In this work, we have derived and validated a fully anharmonic, non-perturbative theory of the vibronic interactions in solids that overcomes the two main approximations (harmonic and electron-phonon coupling model) used in perturbative state-of-the-art formalisms. The presented stAViC methodology lays the founding for further developments, e.g., for elucidating the influence of many-body effects [53] on AViCs and for studying the role of AViCs on absorption coefficients, Auger recombination rates, and transport coefficients [31]. For the latter case, our work suggest that AViCs can play a substantial role, given that they significantly influence effective masses, linewidths (lifetimes), and band-gaps, hence also altering free carrier densities n c (T ). For instance, the simple model n c (T ) ∝ exp [− g T /(2k B T )] reveals that the intrinsic n c in SrTiO 3 increases by two orders of magnitude at 1000 K due to the band-gap's temperature dependence. Since a large fraction of this effect is driven by AViCs, correctly capturing this physics, as done in the stAViC formalism, is pivotal to obtain a correct description of refractory materials.
CC thanks Hagen-Henrik Kowalski and Florian Knoop for fruitful discussions. This project was supported by TEC1p (the European Research Council (ERC) Horizon 2020 research and innovation programme, grant agreement No. 740233), BigMax (the Max Planck Societys Research Network on Big-Data-Driven Materials-Science), and the NOMAD pillar of the FAIR-DI e.V. association. All the electronic-structure theory calculations produced in this project are available on the NOMAD repository: http://dx.doi.org/10.17172/NOMAD/2020.03.18-1.