Polaron-induced band renormalization due to linear and quadratic electron-phonon coupling

TL;DR

This paper introduces a nonperturbative analytical method to study polaron-induced band renormalization in electron-phonon systems, accounting for both linear and quadratic couplings, and applies it to organic crystals like naphthalene.

Contribution

It develops a novel, explicit analytical approach to analyze band narrowing due to electron-phonon interactions beyond linear coupling, applicable in strong-coupling regimes.

Findings

Derived explicit analytical expressions for eigenvalue spectra.

Demonstrated band narrowing as a function of temperature.

Applied the model to naphthalene with ab initio coupling parameters.

Abstract

We present a novel approach to electron-lattice interaction beyond the linear-coupling regime. Based on the solution of a Holstein-Peierls-type model, we derive explicit analytical expressions for the eigenvalue spectrum of the Hamiltonian, resulting in a narrowing of bands as a function of temperature. Our approach enables the intuitive interpretation in terms of quasiparticles, i.e. polaron bands and dressed-phonon frequencies. Being nonperturbative, the formalism also applies in the strong-coupling case. We apply it to the organic crystal naphthalene, with the coupling strengths obtained by \textit{ab initio} calculations.

Equations

$$\displaystyle H\$$

$$=\sum _{mn}{\epsilon }_{mn}{a}_m^{†}{a}_n+\sum _{\mathbf{Q}}\hslash {\omega }_{\mathbf{Q}}(b_{\mathbf{Q}}^{†}{b}_{\mathbf{Q}}+\text{textonehalf})$$

$$+\sum _{mn\mathbf{Q}}\hslash {\omega }_{\mathbf{Q}}{g}_{\mathbf{Q}mn}{a}_m^{†}{a}_n(b_{\mathbf{Q}}^{†}+b_{-\mathbf{Q}})$$

$$+\sum _{mn\mathbf{Q}}\hslash {\omega }_{\mathbf{Q}}{f}_{\mathbf{Q}mn}{a}_m^{†}{a}_n(b_{\mathbf{Q}}^{†}+b_{-\mathbf{Q}})(b_{-\mathbf{Q}}^{†}+b_{\mathbf{Q}}).$$

$$\displaystyle\tilde{H}\,=\,e^{S}\,H\,e^{S^{{\dagger}}{}}\$$

$$S\equiv \sum _{mn}{C}_{mn}{a}_m^{†}{a}_n$$

$$C_{mn}\equiv \sum _{\mathbf{Q}}\left[{\alpha }_{\mathbf{Q}mn}\left(b_{\mathbf{Q}}^{†}-b_{-\mathbf{Q}}\right)+{\beta }_{\mathbf{Q}mn}\left(b_{\mathbf{Q}}^{†}{b}_{-\mathbf{Q}}^{†}-b_{\mathbf{Q}}{b}_{-\mathbf{Q}}\right)\right]$$

$${\alpha }_{\mathbf{Q}mn}\equiv -\frac{1}{4}{\left[{\left((\mathbb{I}+4f_{\mathbf{Q}}{)}^{-1/4}-\mathbb{I}\right)}^{-1}(\mathbb{I}+4f_{\mathbf{Q}}{)}^{-1}\text{ln}(\mathbb{I}+4f_{\mathbf{Q}})g_{\mathbf{Q}}\right]}_{mn}$$

$$\displaystyle\beta_{{\bf{Q}}mn}\equiv\frac{1}{8}\big{[}\textrm{ln}(\mathbb{I}+4f_{\bf{Q}})\big{]}_{mn}$$

$$\tilde{H}=\sum _{\mathbf{k}}{a}_{\mathbf{k}}^{†}{a}_{\mathbf{k}}\left({\epsilon }_0+{\tilde{\epsilon }}_{\mathbf{k}}^{\prime }+{\eta }_{\mathbf{k}}+{\chi }_{\mathbf{k}}\right)+\sum _{\mathbf{Q}}\hslash {\omega }_{\mathbf{Q}}\left(b_{\mathbf{Q}}^{†}{b}_{\mathbf{Q}}+\text{textonehalf}\right)$$

$${\tilde{\epsilon }}_{mn}^{\prime }$$

$${\eta }_{mn}$$

$${\chi }_{mn}$$

$${\lambda }_{\mathbf{Q}mn}$$

$${\Lambda }_{\mathbf{k}}$$

$$a_{\mathbf{k}}$$

$$\displaystyle\langle\ \tilde{\varepsilon}^{\prime}_{mn}\ \rangle\ =\ \varepsilon^{\prime}_{mn}\ \ \textrm{exp}\Big{[}-{\sum_{\bf{Q}}}^{\prime}\mathbb{G}_{{\bf{Q}}mn}\ (\nicefrac{{1}}{{2}}+N_{\bf{Q}})$$

$$\displaystyle\qquad\ \ \ \qquad\quad\qquad\ -{\sum_{\bf{Q}}}^{\prime}\mathbb{F}_{{\bf{Q}}mn}\ ((1+N_{\bf{Q}})^{2}+N^{2}_{\bf{Q}})\ \ \Big{]}\ ,$$

$$\displaystyle\mathbb{G}_{{\bf{Q}}mn}\ \equiv\ |\ \alpha_{{\bf{Q}}mm}-\alpha_{{\bf{Q}}nn}\ |^{2}\ +\ \sum_{l\neq m,n}\ \big{(}\ |\,\alpha_{{\bf{Q}}ml}\,|^{2}+|\,\alpha_{{\bf{Q}}nl}\,|^{2}\ \big{)}\,,$$

$$\displaystyle\mathbb{F}_{{\bf{Q}}mn}\ \equiv\ |\ \beta_{{\bf{Q}}mm}-\beta_{{\bf{Q}}nn}\ |^{2}\ +\ \sum_{l\neq m,n}\ \big{(}\ |\,\beta_{{\bf{Q}}ml}\,|^{2}+|\,\beta_{{\bf{Q}}nl}\,|^{2}\ \big{)}\,,$$

$${\epsilon }_{\pm }(\mathbf{k})=$$

$$+2{\epsilon }_{a{b}^{\prime }}\left[\cos \left(\mathbf{k}\cdot (\mathbf{a}+\mathbf{b})\right)+\cos \left(\mathbf{k}\cdot (\mathbf{a}-\mathbf{b})\right)\right]$$

$$+2{\epsilon }_{ac}\cos \left(\mathbf{k}\cdot (\mathbf{a}+\mathbf{c})\right)+2{\epsilon }_{a{c}^{\prime }}\cos \left(\mathbf{k}\cdot (\mathbf{a}-\mathbf{c})\right)$$

$$+2{\epsilon }_{b{c}^{\prime }}\left[\cos \left(\mathbf{k}\cdot (\mathbf{b}+\mathbf{c})\right)+\cos \left(\mathbf{k}\cdot (\mathbf{b}-\mathbf{c})\right)\right]$$

$$\pm 2{\epsilon }_{ab}\left[\cos \left(\mathbf{k}\cdot \frac{\mathbf{a}+\mathbf{b}}{2}\right)+\cos \left(\mathbf{k}\cdot \frac{\mathbf{a}-\mathbf{b}}{2}\right)\right]$$

$$\pm 2{\epsilon }_{abc}\left[\cos \left(\mathbf{k}\cdot \frac{\mathbf{a}+\mathbf{b}+2\mathbf{c}}{2}\right)+\cos \left(\mathbf{k}\cdot \frac{\mathbf{a}-\mathbf{b}+2\mathbf{c}}{2}\right)\right]$$

$$\pm 2{\epsilon }_{ab{c}^{\prime }}\left[\cos \left(\mathbf{k}\cdot \frac{\mathbf{a}+\mathbf{b}-2\mathbf{c}}{2}\right)+\cos \left(\mathbf{k}\cdot \frac{\mathbf{a}-\mathbf{b}-2\mathbf{c}}{2}\right)\right].$$

$$g_{\nu ,mn}=\frac{1}{\hslash {\omega }_{\nu }}\frac{\delta {\epsilon }_{mn}}{\delta U_{\nu }};f_{\nu ,mn}=\frac{1}{2\hslash {\omega }_{\nu }}\frac{{\delta }^2{\epsilon }_{mn}}{\delta U_{\nu }^2}.$$

Full text

The coupling between electronic and vibrational degrees of freedom is one of the most prominent examples of fermion–boson interaction in condensed-matter systems. Its theoretical description and analysis is central to the understanding of many physical phenomena like heat and charge transport, or superconductivity Mahan (1990); Crespi and Cohen (1993). The complex and possibly strong interaction between electrons and phonons makes perturbative approaches often problematic. A way to obtain non-perturbative solutions of the system Hamiltonian is to perform transformations which allow for a description of the system in terms of independent dressed particles instead of coupled bare particles. This is the concept behind the polaron transformation, which leads to analytical solutions in the linear-coupling regime Mahan (1990). Analogous methods for the nonlinear nonlocal electron-phonon coupling (EPC) have not been available so far.

It is known from model Hamiltonians Adolphs and Berciu (2013a) and from applications of Green-function theory Cardona and Thewalt (2005); Giustino et al. (2010); Cannuccia and Marini (2011); Kawai et al. (2014); Gonze et al. (2011); Antonius et al. (2014a); Poncé et al. (2014) as initiated in early work by Allen, Heine, and Cardona Allen and Heine (1976); Allen and Cardona (1981, 1983) that the renormalization of electronic bands due to quadratic electron-phonon coupling is often of the same order as the linear one. Since these two contributions are often opposite in sign Gonze et al. (2011), the neglect of quadratic coupling can lead to significant errors in electronic-structure calculations considering vibrational effects. The relevance of quadratic EPC terms has also been pointed out for properly describing a number of phenomena, like band renormalization and spin crossover D’Avino et al. (2011), charge transfer Pecile et al. (1989) or the Jahn-Teller effect Koizumi et al. (2000), and many materials like graphene Lozovik and Sokolik (2010); Ziegler et al. (2011); Basko (2008), SrTiO3 Liu et al. (2001), molecular systems Gonze et al. (2011) or systems with torsional motions Coropceanu et al. (2007). Recent research indicates dramatic effects of even small nonlinear EPC in model Hamiltonians Li et al. (2015); Li and Johnston (2015). Despite the prominent character of the quadratic electron-phonon interaction, just a few non-perturbative model approaches exist in the literature Munn and Silbey (1978); Riseborough (1984); Entin-Wohlman et al. (1985); Olsen (2009); Adolphs and Berciu (2013b); Adolphs and Berciu (2014).

In this work, we present a novel formalism for the non-perturbative treatment of linear and quadratic electron-phonon interaction. We derive an analytical solution that facilitates the understanding of the essential underlying physics in an intuitive way. It also allows for obtaining accurate solutions in the strong coupling regime, where a perturbative treatment fails, and for tackling large systems that would not be numerically affordable not (a). Moreover, the nonlocality of the Hamiltonian makes our procedure quite generally applicable, and, thus for instance, suited to describe phenomena like charge transport Coropceanu et al. (2009); Tu et al. (2018); Ortmann et al. (2009).

Our theory is based upon a quadratic Holstein-Peierls-type model for the interaction between electrons and phonons. An analogous approach, including linear coupling only, has been successfully applied previously Hannewald et al. (2004a). Using a tight-binding description for crystalline materials, the corresponding Hamiltonian takes the form

[TABLE]

The electronic subsystem is described by $H_{el}$ where the fermionic operators $a_{m}$ ($a^{{{\dagger}}}_{m}$) correspond to the annihilation (creation) of electrons at lattice site ${\bf R}_{m}$ with on-site energies $\varepsilon_{mm}$ and the strength of their interatomic coupling being governed solely by the transfer integrals $\varepsilon_{mn}$. The phononic subsystem is given by $H_{ph}$ where the bosonic operators $b_{\bf{Q}}$ ($b^{{{\dagger}}}_{\bf{Q}}$) describe the annihilation (creation) of a phonon mode $\nu$ with wave vector $\bf{q}$ at frequency $\omega_{\bf{Q}}$ [$\bf{Q}\equiv(\bf{q},\nu)$]. The interaction between the electronic and vibrational degrees of freedom is governed by $H_{el-ph}$. The electron-phonon matrix elements $g_{{\bf{Q}}mn}$ and $f_{{\bf{Q}}mn}$ (note that $f_{{\bf{Q}}}$, $g_{{\bf{Q}}}$ are matrices) can be determined using the methods presented in Refs. Giustino et al. (2010); Gonze et al. (2011); Poncé et al. (2014) or Hannewald et al. (2004a), or by fitting to experimental data Zhugayevych and Tretiak (2015). These coefficients correspond to the linear and quadratic electron-phonon coupling, respectively, because $(b^{\dagger}_{\bf{Q}}+b_{-\bf{Q}})$ is proportional to the Fourier transform of the nuclear displacements SM . We explicitly take into account both the local ($m=n$) Holstein-like terms as well as the nonlocal ($m\neq n$) Peierls-like terms, which may lead to remarkably different results than considering just the local ones SM ; Hannewald et al. (2004a).

Our approach relies on a unitary transformation that makes the Hamiltonian (Eq. Polaron-induced band renormalization due to linear and quadratic electron-phonon coupling) diagonal and thus facilitates an analytical solution. To achieve this goal in a non-perturbative manner, all the terms of the Hamiltonian that do not conserve the number of phonons are removed. By this procedure, we control the complexity that originates from the nonlocal Peierls-type coupling in $H_{el-ph}$. To this end, we adopt approximations that have been applied in a variant of the polaron transformation Mahan (1990); Stefanucci and van Leeuwen (2013) that has turned out successful in the linear-coupling case Hannewald et al. (2004a, b). This approach provides an analytical solution and allows for an intuitive interpretation within the quasiparticle picture.

We propose the following canonical transformation of the Hamiltonian (Eq. Polaron-induced band renormalization due to linear and quadratic electron-phonon coupling):

[TABLE]

In Eq. 2, the functions of matrices are represented by the corresponding Taylor expansion. By performing this transformation, the Hamiltonian (Eq. Polaron-induced band renormalization due to linear and quadratic electron-phonon coupling) becomes SM :

[TABLE]

where

[TABLE]

$N$ is the number of unit cells, and $\varepsilon_{0}$ is the value of the diagonal entries of the matrix of transfer integrals ($\varepsilon_{0}\equiv\varepsilon_{mm}$) Hannewald et al. (2004a). The different terms of $\tilde{H}$ (Eq. 3) correspond to dressed electrons and phonons. $a^{{\dagger}}_{\bf k}a_{\bf k}$ and $b^{{\dagger}}_{\bf Q}b_{\bf Q}$ represent the number operators of electrons and phonons, respectively. Since $H$ (Eq. Polaron-induced band renormalization due to linear and quadratic electron-phonon coupling) contains terms which are not proportional to $b^{{\dagger}}_{\bf Q}b_{\bf Q}$, it is not straightforward to assign particle numbers to the eigenvalues of $H$. Conversely, our approach characterizes every eigenvalue of $\tilde{H}$ (Eq. 3) after the evaluation of $\langle\tilde{\varepsilon}^{\prime}_{\bf{k}}\rangle$ by the number of dressed electrons and phonons, respectively. (A more detailed explanation of the meaning of the terms of (Eq. 3) is presented in the Supplemental Material.) With this evaluation, the transformed Hamiltonian $\tilde{H}$ (Eq. 3) becomes diagonal, which is explained below (see Eq. 5). The advantage of our transformation (Eq. 2a) becomes clear by noticing what follows: $H$ contains terms proportional to $(b^{{\dagger}}_{\bf{Q}}b^{{\dagger}}_{-\bf{Q}}{-}b_{\bf{Q}}b_{-\bf{Q}})$ and to $(b^{{\dagger}}_{\bf{Q}}{+}b_{-\bf{Q}})$, which do not conserve the number of phonons. The accurate calculation of the eigenvalues of the Holstein-Peierls Hamiltonian by numerical means requires very large basis sets, which makes such calculation numerically very demanding Zhugayevych and Tretiak (2015). Indeed, as shown in the Supplemental Material SM for a 2D square lattice Li et al. (2015); Li and Johnston (2015), the numerical solution of this model system, including local and nonlocal electron-vibrational couplings, requires high phonon populations in the basis sets. In contrast, our analytical solver presented here provides accurate solutions very efficiently.

We proceed now by evaluating Eq. 4a, which will allow us to quantify the polaron transfer integrals and to make $\tilde{H}$ diagonal. If the matrices $f_{{\bf{Q}}}$, $g_{{\bf{Q}}}$ commute with $\varepsilon^{\prime}$, this leads to $\tilde{\varepsilon}^{\prime}_{mn}=\varepsilon^{\prime}_{mn}$ making $\tilde{H}$ diagonal without any approximation. If they do not commute, we can make $\tilde{H}$ diagonal by replacing $\tilde{\varepsilon}^{\prime}_{mn}$ by its thermal average $\langle\tilde{\varepsilon}^{\prime}_{mn}\rangle$. Such an approach was employed to solve analogous problems for the linear Hannewald et al. (2004a, b); Zhugayevych and Tretiak (2015) and quadratic cases Munn and Silbey (1978); Entin-Wohlman et al. (1985). By proceeding in this way, an analytical solution becomes feasible, while the behavior of the system is still reliably reproduced. We note that the approximation of thermal averaging is commonly used in literature Allen and Heine (1976); Dahm and Ueda (2007); Tamura et al. (2012); Li et al. (2012); De Filippis et al. (2015). In our case it is particularly harmless Hannewald et al. (2004a, b), because the equations depending on phonon operators are proportional to commutators of $\varepsilon^{\prime}$ the matrices $\alpha_{\bf{Q}}$ and $\beta_{\bf{Q}}$ SM , which are expected to be very small or zero. The thermal averaging can also be avoided if applying perturbation theory on top of $\tilde{H}$, as done in Refs. Munn and Silbey (1978); Entin-Wohlman et al. (1985).

For the explicit evaluation of the thermal averages $\langle\tilde{\varepsilon}^{\prime}_{mn}\rangle$, we apply the Baker-Campbell-Hausdorff (BCH) theorem Hannewald et al. (2004a); Dyinkin (1947) to Eq. 4a. This yields an infinite series that can be evaluated by truncating it at a finite order SM . Nonetheless, an analytical expression for the entire series is also possible, leading to an expression for the band narrowing as a function of temperature

[TABLE]

where $N_{\bf{Q}}$ is the Bose distribution and ′ in the summations means that one must omit the contributions of the ${\bf{Q}}$’s which satisfy $[f_{\bf{Q}},\varepsilon]=0$ (for $\mathbb{F}_{\bf{Q}}$), and $[f_{\bf{Q}},\varepsilon]=[g_{\bf{Q}},\varepsilon]=0$ (for $\mathbb{G}_{\bf{Q}}$) (this happens, e.g., if ${\bf{q}}=0$). In the derivation of (Eq. 5) SM , we have used the approximation $[g_{\bf{Q}},g_{\bf{Q}^{\prime}}]=[g_{\bf{Q}},f_{\bf{Q}^{\prime}}]=[f_{\bf{Q}},f_{\bf{Q}^{\prime}}]=0$ for all $g_{\bf{Q}},g_{\bf{Q}^{\prime}}$, and we only consider the most important contributions ($\varepsilon_{mn}$, $\varepsilon_{nm}$ and $\varepsilon_{mm}$) in the evaluation of the commutators of $f_{\bf Q}$ and $g_{\bf Q}$ with $\varepsilon$. Proceeding this way, corresponds to a generalization of approximations that where shown to have a very small effect Hannewald et al. (2004a, b). The terms $(1+N_{{\bf{Q}}})$ and ($\nicefrac{{1}}{{2}}+N_{\bf{Q}})$ in Eq. 5 can be identified to correspond to the emission and absorption of one phonon, respectively. They arise from the linear EPC. Analogously, the terms in $(1{+}N_{{\bf{Q}}})^{2}$ and $N_{{\bf{Q}}}^{2}$ correspond to emission and absorption processes of two phonons, which arise from the quadratic EPC term of $H$. The evaluation of $\langle\tilde{\varepsilon}^{\prime}_{mn}\rangle$, either using Eq. 5 or by truncating the BCH expansion of $\varepsilon^{\prime}_{mn}$ SM , makes the transformed Hamiltonian $\tilde{H}$ diagonal and provides analytical (symbolic) expressions for its eigenvalue spectrum.

We illustrate the strength of the formalism presented above by applying it to the molecular organic semiconductor naphthalene, which crystallizes in a monoclinic structure (space group $P2_{1/a}$) with two non-equivalent molecules in the unit cell. It has been subject of previous investigations to demonstrate band-narrowing based on a tight-binding model, but including linear coupling only Hannewald et al. (2004a). All the ab initio calculations are performed with the VASP code Kresse and Furthmüller (1996a, b), that implements the projector augmented wave (PAW) method. Hard pseudopotentials Blöchl (1994); Kresse and Joubert (1999) are used for C and H atoms; exchange-correlation effects are described with the local-density approximation Perdew and Zunger (1981). A planewave cutoff energy of 1200 eV and a $6\times 8\times 6$ ${{\bf k}}$-mesh have been employed in the ground-state calculation. We adopt the experimental lattice parameters (at 5K) Capelli et al. (2006) of a=8.08 Å, b=5.93 Å, c=8.63 Å, and $\beta=124.7^{\circ}$, and relax all atomic coordinates until the maximum force is smaller than 0.0001 eV/Å. The vibrational properties at the $\Gamma$ point are obtained through density-functional perturbation theory Baroni et al. (2001), providing the phonon frequencies, $\omega_{\bf Q}$, and eigenvectors (mass-weighted normal modes) $U_{\bf Q}$ Poncé et al. (2014). Our tight-binding Hamiltonian includes the on-site energy and the ten transfer integrals between nearest and next-nearest neighbors, corresponding to ${\mathbf{R}_{m}-\mathbf{R}_{n}}=0,$ ${\pm\mathbf{a},\,}$ ${\pm\mathbf{b},\,}$ ${\pm\mathbf{c},\,}$ ${\pm(\mathbf{a}\pm\mathbf{b}),\,}$ ${\pm(\mathbf{a+c}),\,}$ ${\pm(\mathbf{a-c}),\,}$ ${\pm(\mathbf{b}\pm\mathbf{c}),\,}$ ${\pm(\mathbf{a}\pm\mathbf{b})/2,\,}$ ${\pm({{\bf a}/2}\pm{{\bf b}/2}+\mathbf{c}),\,}$ and ${\pm({{\bf a}/2}\pm{{\bf b}/2}-\mathbf{c})}$, where a, b, and c are the lattice vectors. The electronic eigenvalues provided by this model Katz et al. (1963) are given by:

[TABLE]

The $\pm$ signs in Eq. Polaron-induced band renormalization due to linear and quadratic electron-phonon coupling arise due to the fact that there are two naphthalene molecules per unit cell Katz et al. (1963). To obtain the parameters of Eq. Polaron-induced band renormalization due to linear and quadratic electron-phonon coupling, we proceed in a manner analogous to the one presented in Hannewald et al. (2004a), though more general. The 11 $\varepsilon_{mn}$ parameters used in the tight-binding Hamiltonian are obtained by a least-square fit to the ab initio band structures (including 264 indepedent $\bf k$-points). The coupling coefficients $g_{\mathbf{Q}mn}$ and $f_{\mathbf{Q}mn}$ correspond to the first and second derivatives of $\varepsilon_{mn}$, respectively, with respect to structural deformations according to each phonon eigenvector, given by

[TABLE]

We restrict the calculation of coefficients, to $\Gamma$-point phonons considering the nine intermolecular phonon branches ($\nu$) with lowest phonon frequencies not (b) as also done in Ref. Hannewald et al. (2004a). For the calculation of the coupling coefficients (Eq. 7), we perform finite-difference displacements of the nuclear positions along the phonon eigenvectors. This results in approximately parabolic curves for each of the tight-binding coefficients as a function of displacement. Each of these curves is fitted to a cubic function, whose slope and curvature are used to evaluate $g_{\nu,mn}$ and $f_{\nu,mn}$, respectively. We then calculate the coupling coefficients at finite $\bf q$ by assuming the eigenmodes to be dispersionless Hannewald et al. (2004a), i.e., $\omega_{{\bf{Q}}}\equiv\omega_{{\bf{q}}\nu}=\omega_{\nu}$, $g_{{\bf{Q}}mn}=g_{\nu,mn}(e^{-i{\bf{q}}{\bf{R}}_{m}}+e^{-i{\bf{q}}{\bf{R}}_{n}})/2$, $f_{{\bf{Q}}mn}=f_{\nu,mn}(1+\textrm{cos}({i{\bf{q}}({\bf{R}}_{m}-{\bf{R}}_{n})})/2$ SM .

In Fig. 1 we display the band widths of the top valence band and the lowest conduction band (for brevity termed HOMO and LUMO, respectively) as predicted by our new formalism. For comparison, we show the results obtained by purely linear electron-vibrational coupling Hannewald et al. (2004a). The effect of the quadratic coupling gives rise to significant changes, which are a substantial reduction of the LUMO (24% at T=0) and an increase of the HOMO band width. While the absolute amounts only vary very little along the whole temperature range, the relative changes are significant. For example, the variation of the HOMO band width by second-order coupling is about 45% at 300K. This is in line with previous works on electron-phonon interaction, which showed that the contribution of the quadratic coupling to band renormalization must not be ignored Antonius et al. (2014b); Giustino et al. (2010); Gonze et al. (2011). Here we demonstrate that this also holds for polaron-induced band-narrowing effects. Interestingly, we learn that the quadratic coupling can either increase or decrease the band width, also in line with the effect of the quadratic coupling on band renormalization Gonze et al. (2011). For this system, we find a low contribution of the quadratic nonlocal terms of the Hamiltonian ($f_{{\bf Q}mn}$ with $m\neq n$) not (b).

To summarize this work, we have presented a nonperturbative theory of electron-phonon interaction based on an analytical solution of the quadratic Holstein-Peierls-type Hamiltonian. Despite the inherent complexity of the problem, our results allow for an intuitive interpretation in terms of quasiparticles, i.e., independent polarons and dressed phonons. As all the related material parameters can be obtained from ab initio calculations, our theory represents an important step towards the quantitative analysis of temperature-dependent band narrowing in real materials beyond perturbational approaches. Moreover, our approach could be extended towards polaron-mobility theories Ortmann et al. (2009, 2011) including quadratic electron-phonon scattering processes.