All-coupling polaron optical response: analytic approaches beyond the adiabatic approximation
aa r X i v : . [ c ond - m a t . o t h e r] J un All-coupling polaron optical response: analytic approaches beyond the adiabaticapproximation
S. N. Klimin, ∗ J. Tempere, † and J. T. Devreese ‡ TQC, Universiteit Antwerpen, Universiteitsplein 1, B-2610 Antwerpen, Belgium (Dated: August 6, 2018)In the present work, the problem of an all-coupling analytic description for the optical conductivityof the Fr¨ohlich polaron is treated, with the goal being to bridge the gap in validity range that existsbetween two complementary methods: on the one hand the memory function formalism and onthe other hand the strong-coupling expansion based on the Franck-Condon picture for the polaronresponse. At intermediate coupling, both methods were found to fail as they do not reproduceDiagrammatic Quantum Monte Carlo results. To resolve this, we modify the memory functionformalism with respect to the Feynman-Hellwarth-Iddings-Platzman (FHIP) approach, in orderto take into account a non-quadratic interaction in a model system for the polaron. The strong-coupling expansion is extended beyond the adiabatic approximation, by including into the treatmentnon-adiabatic transitions between excited polaron states. The polaron optical conductivity that weobtain by combining the two extended methods agree well, both qualitatively and quantitatively,with the Diagrammatic Quantum Monte Carlo results in the whole available range of the electron-phonon coupling strength.
I. INTRODUCTION
The polaron, first proposed as a physical concept byL. D. Landau in the context of electrons in polar crys-tals, has become a generic notion describing a particleinteracting with a quantized bosonic field. The polaronproblem has consequently been used for a long time as atesting ground for various analytic and numerical meth-ods with applications in quantum statistical physics andquantum field theory. In condensed matter physics, thepolaron effect coming from the electron-phonon interac-tion is a necessary ingredient in the description of the DCmobility and the optical response in polar crystals (seeRef. ). Polaronic effects are manifest in many interestingsystems, such as magnetic polarons , polarons in semi-conducting polymers , and complex oxides which aredescribed in terms of the small-polaron theory . Large-polaron theory has recently been stimulated by the pos-sibility to study polaronic effects using highly tunablequantum gases: the physics of an impurity immersed inan atomic Bose-Einstein condensate can be modeled onthe basis of a Fr¨ohlich Hamiltonian. Another recent de-velopment in large-polaron physics stems from the ex-perimental advances in the determination of the bandstructure of highly polar oxides , relevant for supercon-ductivity, where the optical response of complex oxidesexplicitly shows the large-polaron features .Diagrammatic Quantum Monte Carlo (DQMC) meth-ods have been applied in recent years to numerically cal-culate the ground state energy and the optical conduc-tivity of the Fr¨ohlich polaron . Advances in computa-tional techniques such as DQMC inspired renewed studyof the key problem in polaron theory – an analytic de-scription of the polaron response. For the small-polaron optical conductivity, the all-coupling analytic theory hasbeen successfully developed showing good agreementwith the numeric results of the DQMC. However, theoptical response problem for a large polaron is not yet completely solved analytically.Asymptotically exact analytic solutions for the polaronoptical conductivity have been obtained in the limits ofweak and strong coupling . A first proposal for anall-coupling approximation for the polaron optical con-ductivity has been formulated in Ref. (below referredto as DSG), further developing the Feynman-Hellwarth-Iddings-Platzman theory (FHIP) and using the Feyn-man variational approach . However, in Ref. , it wasalready demonstrated that FHIP at large α is inconsis-tent with the Heisenberg uncertainty relations. This in-consistency is revealed in Ref. through extremely nar-row peaks of the optical conductivity at large α . Nev-ertheless, the peak positions for the polaron optical con-ductivity as obtained in Ref. have been confirmed withhigh accuracy by the DQMC calculation . This in-spired further attempts to develop analytical methodsfor the polaron optical response, especially at interme-diate and strong coupling. Among these analytic meth-ods, an extension of the DSG method has been proposedin Ref. introducing an extended memory function for-malism with a relaxation time determined from the ad-ditional sum rule for the polaron optical conductivity.Alternatively, for the strong coupling regime, the strongcoupling expansion (SCE) based on the Franck-Condonscheme for multiphonon optical conductivity has beendeveloped in Refs. .The extended memory-function formalism and thestrong coupling expansion of Ref. are complementaryto each other. The memory-function formalism is well-substantiated for small and intermediate values of α ,and the strong coupling expansion adequately describesthe opposite limit of large α . However these two meth-ods only qualitatively agree with each other and withthe DQMC data in the range of intermediate couplingstrengths. On one hand, the memory function formalismdisagrees with DQMC at large α . On the other hand,the strong-coupling expansion only qualitatively repro-duces the shape of the optical conductivity and fails atintermediate α . The main aim of the current paper isto extend both the memory function formalism and thestrong coupling expansion in order to bridge the gap thatremains between their regions of validity , such that thecombination of both methods allows to find analyticalresults in agreement with the numeric DQMC results atall coupling.We have added the following new elements in the the-ory which lead to an overlapping of the areas of appli-cability for two aforesaid analytic methods. For weakand intermediate coupling strengths, an extension of theFeynman variational principle and the memory-functionmethod for a polaron with a non-quadratic trial ac-tion has been developed. As distinct from the mem-ory function formalism of Ref. , we do not use addi-tional sum rules and relaxation times, and perform thecalculation ab initio . For intermediate and strong cou-pling strengths, the strong coupling expansion of Ref. is extended beyond the adiabatic approximation account-ing for non-adiabatic transitions between excited polaronstates. This leads to a substantial expansion of therange of validity for the strong-coupling expansion to-wards smaller α and to an overall improvement of itsagreement with DQMC.The paper is organized as follows. In Sec. II we de-scribe an all-coupling analytic description for the polaronoptical conductivity within the extended memory func-tion formalism with a non-parabolic trial action and thenon-adiabatic strong-coupling expansion. Sec. III con-tains the discussion of the obtained optical conductivityspectra and their comparison with results of other meth-ods and with the DQMC data. The discussion is followedby conclusions, Sec. IV. II. ANALYTIC METHODS FOR THEPOLARON OPTICAL CONDUCTIVITYA. Memory function formalism with anon-parabolic trial action
To generalize the memory function formalism, we startby extending Feynman’s variational approach to transla-tion invariant non-Gaussian trial actions. The electron-phonon system is described by the Fr¨ohlich Hamiltonian,using the Feynman units with ~ = 1 , the LO-phononfrequency ω LO = 1, and the band mass m b = 1,ˆ H = ˆp X q (cid:18) ˆ a + q ˆ a q + 12 (cid:19) + 1 √ V X q p √ παq (cid:0) ˆ a q + ˆ a + − q (cid:1) e i q · ˆr , (1)where ˆr is the position operator of the electron, ˆp is itsmomentum operator; ˆ a † q and ˆ a q are, respectively, the cre-ation and annihilation operators for longitudinal optical (LO) phonons of wave vector q . The electron-phononcoupling strength is described by the Fr¨ohlich couplingconstant α . As this Hamiltonian is quadratic in thephonon degrees of freedom, they can be integrated outanalytically in the path-integral approach . The remain-ing electron degree of freedom is described via an actionfunctional where the effects of electron-phonon interac-tion are contained in an influence phase: S [ r e ( τ )] = 12 β Z ˙r e ( τ ) dτ − Φ[ r e ( τ )] . (2)Here r e ( τ ) is the path of the electron, expressed in imag-inary time so as to obtain the euclidean action, and β = 1 / ( k B T ) with T the temperature. The influencephase corresponding to (1) isΦ[ r e ( τ )] = √ πα Z d q (2 π ) β Z dτ β Z dτ ′ cosh (cid:16) | τ − τ ′ | − β (cid:17) sinh( β/ × e i q · [ r e ( τ ) − r e ( τ ′ ) ] . (3)This depends on the difference in electron position at dif-ferent times, resulting in a retarded action functional. Inthe path-integral formalism, thermodynamic potentials(such as the free energy) are calculated via the partitionsum, which in turn is written as a sum over all possiblepaths r e ( τ ) of the electron that start and end in the samepoint, weighted by the exponent of the action: e − βF = Z = Z D r e e − S [ r e ( τ )] . (4)Feynman’s original variational method considers a quadratic trial action where the phonon degrees of free-dom are replaced a fictitious particle with coordinate r f ( τ ), S quad [ r e ( τ ) , r f ( τ )] = β Z (cid:20) m ˙r e m f ˙r f V ( r f − r e ) dτ, (5)interacting with the electron through a harmonic poten-tial: V ( r f − r e ) = m f ω r f − r e ) . (6)The partition sum corresponding to this model is Z quad = Z D r e Z D r f e − S quad [ r e ( τ ) , r f ( τ )] . (7)Expectation values of of functionals F [ r e ( τ )] of the elec-tron path are given by hF [ r e ( τ )] i quad = 1 Z quad Z D r e F [ r e ( τ )] × Z D r f e − S quad [ r e ( τ ) , r f ( τ )] . (8)In the above formula, it is clear that performing the pathintegral over r f exactly may simplify the result. Theresult of this integration still depends on the path r e ( τ )and is written in terms of an influence phase, Z D r f exp − β Z " m f ˙r f m f ω r f − r e ) dτ = Z f exp { Φ quad [ r e ( τ )] } , (9)where Z f is a constant independent of r e ( τ ). The influ-ence phase corresponds to a quadratic, retarded interac-tion for the electron. The model system partition sumbecomes Z quad = Z f Z D r e exp − β Z m ˙r e dτ + Φ quad [ r e ( τ )] . (10)Feynman restricted his trial action to a quadratic action,since only for case one can calculate the influence phaseanalytically, and obtainΦ quad [ r e ( τ )] = − m f ω β Z dτ β Z dτ ′ [ r e ( τ ) − r e ( τ ′ )] × cosh h ω (cid:16) | τ − τ ′ | − β (cid:17)i sinh( βω/ . (11)The essence of the Feynman variational method con-sists in writing the partition function of the true electron-phonon system (2) as Z = 1 Z f Z D r e exp { Φ[ r e ( τ )] − Φ quad [ r e ( τ )] }× Z D r f exp {− S quad [ r e ( τ ) , r f ( τ )] } . (12)Indeed, performing the path integration for the fictitiousparticle via (9) cancels Φ quad [ r e ( τ )] as well as the factor Z f , and leaves the kinetic energy contribution, restoringthe action function of the true system. The usefulness ofthe above expression lies in the fact that it can also beinterpreted as an expectation value with respect to themodel system: Z = Z quad Z f h exp { Φ[ r e ( τ )] − Φ quad [ r e ( τ )] }i quad . (13)Now using Jensen’s inequality D e F [ r e ( τ )] E > e hF [ r e ( τ )] i , (14) and taking the logarithm of the resulting expression leadsto F F + 1 β h Φ quad [ r e ( τ )] − Φ[ r e ( τ )] i quad , (15)which is the Jensen-Feynman variational inequality. Herewe introduce the notation Z quad / Z f = e − βF . Note thatwhen we write the above expression in terms of the ac-tions rather than the influence phases, one gets the morefamiliar form of the Jensen-Feynman inequality: F F + 1 β h S [ r e ( τ )] − S quad [ r e ( τ )] i quad . (16)Using the Feynman variational approach with the Gaus-sian trial action, excellent results are obtained for the po-laron ground-state energy, free energy, and the effectivemass. Moreover, this approach has been effectively usedto derive the DSG all-coupling theory for the polaron op-tical conductivity, Ref. . However, as mentioned in theintroduction, the DSG and DQMC results contradict toeach other in the range of large α . The most probablesource of this contradiction is the Gaussian form of thetrial action used in the DSG theory. Indeed, the modelsystem contains only a single frequency, leading to un-physically sharp peaks in the spectrum, subject to ther-mal broadening only . Extensions to the formalism have tried to overcome this problem by including an ad-hoc broadening of the energy level, chosen in such as wayas to comply with the sum rules. A remarkable successin the problem of the polaron optical response has beenachieved in the recent work , where the all-coupling po-laron optical conductivity is calculated using the generalquadratic trial action instead of the Feynman model witha single fictitious particle. The resulting optical conduc-tivity is in good agreement with DQMC results in theweak- and intermediate-coupling regimes and is qualita-tively in line with DQMC even at extremely strong cou-pling, resolving the issue of the linewidth in the FHIPapproach.In the literature, there are attempts to re-formulatethe Feynman variational approach avoiding retarded trialactions. For example, Cataudella et al. introduce anextended action which contains the coordinates of theelectron, the fictitious particle, and the phonons. Thisaction, however, is not exactly equivalent to the actionof the electron-phonon system, and hence the results ob-tained in need verification. In Ref. , we introduced anextended action/Hamiltonian for an electron-phonon sys-tem and reformulated the Feynman variational methodin the Hamiltonian representation. This method leadsto the same result as the Feynman variational approach.However the method of Ref. reproduces the strong cou-pling limit for the polaron energy only when using aGaussian trial action.In the current work, we propose to extend the Feynmanvariational approach to trial systems with non-parabolicinteractions between an electron and a fictitious parti-cle. The difficulty with using non-Gaussian trial actionsis that the influence phase (here Φ quad ) can only be com-puted analytically for quadratic action functionals. How-ever, quantum-statistical expectation values (such as theone in the Jensen-Feynman inequality) can be calculatedfor non-quadratic model systems by other means, in par-ticular if the spectrum of eigenvalues and eigenfunctionscan be found. So, what we propose is to focus on keep-ing the influence phase for a quadratic model system inthe expressions, while at the same time allowing for non-Gaussian potentials for the expectation values .Consider a (non-quadratic) variational trial action S var [ r e ( τ ) , r f ( τ )] = β Z (cid:20) m ˙r e m f ˙r f U ( r f − r e ) dτ (17)with a general potential U . Since S quad [ r e ( τ ) , r f ( τ )] = S var [ r e ( τ ) , r f ( τ )]+ β Z [ V ( r f − r e ) − U ( r f − r e )] dτ, (18)we can rewrite (12) to: Z = 1 Z f Z D r e Z D r f exp { Φ[ r e ( τ )] − Φ quad [ r e ( τ )] }× exp − β Z [ V ( r f − r e ) − U ( r f − r e )] dτ × exp {− S var [ r e ( τ ) , r f ( τ )] } . (19)Expectation values of functionals of r e ( τ ) and r f ( τ ) withrespect to the non-quadratic variational model systemare given by hF [ r e , r f ] i var = 1 Z var Z D r e Z D r f F [ r e , r f ] × exp {− S var [ r e ( τ ) , r f ( τ )] } . (20)This allows to interpret expression (19) as Z = Z var Z f h exp { Φ[ r e ( τ )] − Φ quad [ r e ( τ )] − R β [ V ( r f − r e ) − U ( r f − r e )] dτ oE var (21)With Z var / Z f = e − βF var and using again Jensen’s in-equality, we arrive at: F F var + 1 β h Φ quad [ r e ( τ )] − Φ[ r e ( τ )] i var + h V ( r f − r e ) − U ( r f − r e ) i var (22) We have used that for time-independent potentials, β R β dτ = 1. When U = V , and we choose a quadraticinteraction potential, this restores the original Jensen-Feynman variational principle. However, deviations froma quadratic potential result in two changes. Firstly, thereis an additional term corresponding to the expectationvalue of the difference between the chosen variational po-tential and the quadratic one. Secondly, the expectationvalues are to be calculated with respect to the chosenvariational potential U rather than with respect to thequadratic potential. We again emphasize that the cal-culation of such expectation values does not require aquadratic potential. Thus the new variational inequality (22) is an extension of the Feynman – Jensen inequality.It is important for the calculations that S var is trans-lation invariant but non-retarded action, so that all ex-pressions in the variational functional (22) have the sameform in both representations – path integral and stan-dard quantum mechanics. Apart from the parametersappearing in the trial action S var , the inequality (22) stillcontains as variational parameters m f and ω , inheritedfrom the “auxiliary” quadratic action S quad and appear-ing in Φ quad and V ( r f − r e ). These do not depend onthe parameters of the electron-phonon interaction at all,and therefore the minimization of the free energy withrespect to m f and ω can be performed before the mini-mization of the whole variational functional with respectto the remaining parameters of the non-Gaussian trialaction S var .The extended Jensen-Feynman inequality (22), despitehaving more variational parameters, does not lead in gen-eral to a substantially lower polaron free energy than theoriginal Feynman result, except in the extremely strongcoupling regime, where the present variational functionalanalytically tends (for T = 0) to the exact strong cou-pling limit obtained by Miyake . However, its advantagewith respect to the original Feynman treatment is in cal-culating the optical conductivity. A physically reasonablechoice of the trial interaction potential U ( r f − r ) is nolonger restricting to a single frequency oscillator. Accord-ing to Refs. , the self-consistent potential for an elec-tron induced by the lattice polarization is parabolic nearthe bottom and Coulomb-like at large distances. There-fore, for the calculation of the optical conductivity, wechoose a trial potential in the piecewise form, stitchingtogether a parabolic and a Coulomb-like potential in acontinuous way. The spectrum of internal states of themodel system with this potential necessarily consists ofan infinite number non-equidistant energy levels with theenergies E n < E >
0. Accounting for transitions between allthese levels, one must expect a significant broadening ofthe peak absorption.The polaron optical conductivity is calculated usingthe memory-function formalism which differs from thatapplied in Refs. by using the non-quadratic trialaction described above. The optical conductivity of agas of interacting polarons within the memory-functiontechnique, is given by the formula which is structurallysimilar to the polaron optical conductivity , σ (Ω) = e n m b i Ω − χ (Ω) / Ω , (23)where n = N/V is the carrier density. The memoryfunction in the non-quadratic setting is given by χ (Ω) = 23 ~ m b Z d q (2 π ) q | V q | Z ∞ dt (cid:0) e i Ω t − (cid:1) × Im cos h ω (cid:16) t + i ~ β (cid:17)i sinh (cid:16) β ~ ω (cid:17) D e i q · r ( t ) e − i q · r E var , (24)where ω is the LO phonon frequency, and the corre-lation function (cid:10) e i q · r ( t ) e − i q · r (cid:11) var is calculated with thequantum states of the trial Hamiltonian correspondingto S var . In the particular case T = 0 we apply the for-mula following from (24) χ (Ω) = 13 π ~ m b lim δ → + Z ∞ dq | V q | q ∞ Z dt e − δt × (cid:0) e i Ω t − (cid:1) Im (cid:16) e − iω t D e i q · r ( t ) e − i q · r E var (cid:17) . (25)Rather than computing the correlation function (cid:10) e i q · r ( t ) e − i q · r (cid:11) var as a path integral through (20), wechoose to evaluate it in the equivalent Hamiltonian for-malism. In this Hamiltonian framework, (25) is writtenas a sum over the eigenstates of the trial Hamiltonianfor the electron and the fictitious particle interactingthrough the potential U ,ˆ H var = ˆp ˆp f m f + U ( ˆr f − ˆr ) . (26)Thus the correlation function is given by: D e i q · r ( t ) e − i q · r E var = X k ′ ; l ′ ,n ′ ,m ′ e i t ~ (cid:20) ( ε , − ε l ′ ,n ′ ) − ~ ( k ′ ) M (cid:21) × (cid:12)(cid:12)(cid:10) ψ ;0 , , (cid:12)(cid:12) e i q · r (cid:12)(cid:12) ψ k ′ ; l ′ ,n ′ ,m ′ (cid:11)(cid:12)(cid:12) , (27)where M = 1 + m f is the total mass of the trial systemand | ψ k ; l,n,m i are the wave functions of the trial system.The wave function | ψ k ; l,n,m i is factorized as a product ofa plane wave for the center-of-mass motion (with center-of-mass coordinate R ) and a wave function for the rela-tive motion | ϕ l,n,m i (with the coordinate vector ρ of therelative motion), | ψ k ; l,n,m i = 1 √ V e i k · R | ϕ l,n,m i , (28) | ϕ l,n,m i = R l,n ( ρ ) Y l,m ( θ, ϕ ) . (29) The quantum numbers for the Hamiltonian ˆ H var are themomentum k , the quanta l, m related to to angular mo-mentum, and a nodal quantum number n for the relativemotion wavefunction. The quantum numbers l, n deter-mine the energy ε l,n associated with the relative motionbetween electron and fictitious particle (including boththe discrete and continuous parts of the energy spec-trum). When substituting (27) into the memory functionwe arrive at the result χ (Ω) = 13 π ~ m b Z ∞ dq | V q | q ∞ Z dte − δt × X k ′ ; l ′ ,n ′ ,m ′ (cid:12)(cid:12)(cid:10) ψ ;0 , , (cid:12)(cid:12) e i q · r (cid:12)(cid:12) ψ k ′ ; l ′ ,n ′ ,m ′ (cid:11)(cid:12)(cid:12) × (cid:0) e i Ω t − (cid:1) Im e − i t ~ (cid:20) ( ε l ′ ,n ′ − ε , ) + ~ ( k ′ ) M + ~ ω (cid:21) ! , (30)with δ → +0. Further, the Feynman units are used,where ~ = 1, ω = 1, and the band mass m b = 1. Inthese units, the squared modulus | V q | is | V q | = 2 √ παq The memory function is then given by the formula χ (Ω) = 2 √ α π Z ∞ dq q × X k ′ ; l ′ ,n ′ ,m ′ (cid:12)(cid:12)(cid:10) ψ ;0 , , (cid:12)(cid:12) e i q · r (cid:12)(cid:12) ψ k ′ ; l ′ ,n ′ ,m ′ (cid:11)(cid:12)(cid:12) × ∞ Z dte − δt (cid:0) e i Ω t − (cid:1) × Im (cid:18) e − it h ( ε l ′ ,n ′ − ε , ) + M ( k ′ ) +1 i (cid:19) . (31)The matrix element in (31) is a product of two matrixelements: (cid:10) ψ k ; l,n,m (cid:12)(cid:12) e i q · r (cid:12)(cid:12) ψ k ′ ; l ′ ,n ′ ,m ′ (cid:11) = 1 V D e − i kR (cid:12)(cid:12) e i q · R (cid:12)(cid:12) e i k ′ R E (cid:10) ϕ l,n,m (cid:12)(cid:12) e iµ q · ρ (cid:12)(cid:12) ϕ l ′ ,n ′ ,m ′ (cid:11) , (32)where µ is the reduced mass of the trial system. The firstmatrix element is1 V D e − i kR (cid:12)(cid:12) e i q · R (cid:12)(cid:12) e i k ′ R E = δ k ′ , k − q . (33)This eliminates the integration over the final electron mo-mentum k ′ and reduces the memory function to the ex-pression χ (Ω) = 2 √ α π Z ∞ dq q ∞ Z dte − δt × X l ′ ,n ′ ,m ′ (cid:12)(cid:12)(cid:10) ϕ , , (cid:12)(cid:12) e iµ q · ρ (cid:12)(cid:12) ϕ l ′ ,n ′ ,m ′ (cid:11)(cid:12)(cid:12) × (cid:0) e i Ω t − (cid:1) Im (cid:18) e − it (cid:16) q M + ε l ′ ,n ′ − ε , +1 (cid:17) (cid:19) . (34)The summation over m ′ is performed explicitly: X m,m ′ (cid:12)(cid:12)(cid:10) ϕ l,n,m (cid:12)(cid:12) e iµ q · ρ (cid:12)(cid:12) ϕ l ′ ,n ′ ,m ′ (cid:11)(cid:12)(cid:12) = (2 l + 1) (2 l ′ + 1)2 Z ∞ ρ dρ Z ∞ dρ ′ ( ρ ′ ) × R l,n ( ρ ) R l ′ ,n ′ ( ρ ) R l,n ( ρ ′ ) R l ′ ,n ′ ( ρ ′ ) × Z π sin ( µq | ρ − ρ ′ | ) µq | ρ − ρ ′ |× P l (cos θ ) P l ′ (cos θ ) sin θdθ. (35)The modulus | ρ − ρ ′ | is expressed as | ρ − ρ ′ | = q ρ + ( ρ ′ ) − ρρ ′ cos θ. (36)Hence we can use the expansion of sin ( µq | ρ − ρ ′ | ) µq | ρ − ρ ′ | throughthe Legendre polynomials:sin ( µq | ρ − ρ ′ | ) µq | ρ − ρ ′ | = ∞ X l ′′ =0 (2 l ′′ + 1) × j l ′′ ( µqρ ) j l ′′ ( µqρ ′ ) P l ′′ (cos θ ) . (37)The integral of the product of three Legendre polynomi-als is expressed through the 3 j -symbol: Z π P l ′′ (cos θ ) P l (cos θ ) P l ′ (cos θ ) sin θdθ = 2 (cid:18) l l ′ l ′′ (cid:19) . (38)Therefore we find that X m,m ′ (cid:12)(cid:12)(cid:10) ϕ l,n,m (cid:12)(cid:12) e iµ q · ρ (cid:12)(cid:12) ϕ l ′ ,n ′ ,m ′ (cid:11)(cid:12)(cid:12) = ∞ X l ′′ =0 (2 l + 1) (2 l ′ + 1) (2 l ′′ + 1) × (cid:18) l l ′ l ′′ (cid:19) S q ( l, n | l ′′ | l ′ , n ′ ) , (39)where S q ( l, n | l ′′ | l ′ , n ′ ) is the matrix element with radialwave functions for the trial system, S q ( l, n | l ′′ | l ′ , n ′ ) ≡ Z ∞ R l,n ( ρ ) R l ′ ,n ′ ( ρ ) j l ′′ ( µqρ ) ρ dρ. (40) For l = 0 the result of the summation over intermediatestates is reduced to the formula X m ′ (cid:12)(cid:12)(cid:10) ϕ ,n, (cid:12)(cid:12) e iµ q · ρ (cid:12)(cid:12) ϕ l ′ ,n ′ ,m ′ (cid:11)(cid:12)(cid:12) = (2 l ′ + 1) S q (0 , | l ′ | l ′ , n ′ ) . (41)After this summation and the integration over time, thememory function takes the form χ (Ω) = √ α π Z ∞ dq q X l,n (2 l + 1) S q (0 , | l | l, n ) × (cid:18) − Ω q,l,n + iδ −
1Ω + Ω q,l,n + iδ + 2Ω q,l,n (cid:19) . (42)( δ → +0)with the transition frequency for transitions between theground and excited states of the trial system accompa-nied by an emission of a phonon:Ω q,l,n ≡ q M + ε l,n − ε , + 1 . (43)Expression (42) is used for the numerical calculationof the polaron optical conductivity within the extendedmemory function formalism. B. Non-adiabatic strong coupling expansion
Next, we describe the strong coupling approach and itsextension beyond the adiabatic approximation, denotedbelow as the non-adiabatic SCE . Here, the goal is to takenon-adiabatic transitions between different excited levelsof a polaron into account in the formalism. The notationsin this subsection are the same as in Ref. . The polaronoptical conductivity in the strong coupling regime is rep-resented by the Kubo formula,Re σ (Ω) = Ω2 Z ∞−∞ e i Ω t f zz ( t ) dt, (44)with the dipole-dipole correlation function f zz ( t ) = X n,l,m, X n ′ ,l ′ ,m ′ , X n ′′ ,l ′′ ,m ′′ h ψ n,l,m | ˆ z | ψ n ′′ ,l ′′ ,m ′′ i× h ψ n ′ ,l ′ ,m ′ | ˆ z | ψ i× D ph (cid:12)(cid:12)(cid:12)D ψ (cid:12)(cid:12)(cid:12) e it ˆ H ′ (cid:12)(cid:12)(cid:12) ψ n,l,m E × D ψ n ′′ ,l ′′ ,m ′′ (cid:12)(cid:12)(cid:12) e − it ˆ H ′ (cid:12)(cid:12)(cid:12) ψ n ′ ,l ′ ,m ′ E(cid:12)(cid:12)(cid:12) ph E . (45)where | ψ n,l,m i are the polaron states as obtained withinthe strong coupling ansatz in Ref. . The transformedHamiltonian ˆ H ′ of the electron-phonon system after thestrong coupling unitary transformation takes the formˆ H ′ = ˆ H ′ + ˆ W (46)with the termsˆ H ′ = ˆp X q | f q | + V a ( ˆr ) + X q (cid:18) ˆ b + q ˆ b q + 12 (cid:19) , (47)ˆ W = X q (cid:16) ˆ w q ˆ b q + ˆ w ∗ q ˆ b + q (cid:17) . (48)Here, w q are the amplitudes of the renormalized electron-phonon interactionˆ w q = p √ παq √ V (cid:0) e i q · ˆr − ρ q , (cid:1) , (49)where ρ q , is the expectation value of the operator e i q · ˆr with the trial electron wave function | ψ i : ρ q , = (cid:10) ψ (cid:12)(cid:12) e i q · ˆr (cid:12)(cid:12) ψ (cid:11) , (50)and V a ( ˆr ) is the self-consistent potential energy for theelectron, V a ( ˆr ) = − X q √ παq V ρ − q , e i q · ˆr . (51)The eigenstates of the Hamiltonian ˆ H ′ are the prod-ucts of the electron wave functions and those of thephonon vacuum | ψ n,l,m i | ph i . The dipole-dipole correla-tion function f zz ( t ) given by (66) is simplified within theadiabatic approximation for the ground state and usingthe selection rules for the dipole transition matrix ele-ments and the symmetry properties of the polaron Hamil-tonian, as in Ref. . The correlation function, using theinteraction representation, takes the form, f zz ( t ) = X n ′ ,n h ψ | ˆ z | ψ n, , i h ψ n ′ , , | ˆ z | ψ i e − i Ω n, t × h ψ n, , |h ph |U ( t ) | ph i| ψ n ′ , , i , (52)with the Franck-Condon transition frequency,Ω n, ≡ ε n, − ε , , (53)and the evolution operator U ( t ) = T exp (cid:20) − i Z t ds ˆ W ( s ) (cid:21) , (54)where T is the time-ordering symbol, and ˆ W ( s ) is the in-teraction Hamiltonian in the interaction representation,ˆ W ( s ) = e i ˆ H ′ s ˆ W e − i ˆ H ′ s . (55)As found in early works on the strong-coupling Fr¨ohlichpolaron (see, for review, Refs. ), the energy differ-ences between different excited FC states for a strongcoupling polaron are much smaller than the energy dif-ference between the ground and lowest excited FC state.Therefore we keep here the adiabatic approximation for the ground state and, consequently, for the transition be-tween the ground and excited states. On the contrary,the adiabatic approximation for the transitions between different excited states is not applied in (52), as distinctfrom the calculation in Ref. .The matrix elements for the dipole transitions from theground state to other excited states than | ψ , , i (i. e., h ψ | z | ψ n, , i with n = 1) have small relative oscillatorstrengths with respect to h ψ | z | ψ , , i (of order ∼ − ).Therefore further on we consider the next-to-leading or-der nonadiabatic corrections for the contribution to (52)with n = n ′ = 1 and the adiabatic expression for thecontribution with the other terms. In this approxima-tion, (52) takes the form f zz ( t ) = X n |h ψ | ˆ z | ψ n, , i| e − i Ω n, t × h ψ n, , |h ph |U ( t ) | ph i| ψ n, , i . (56)The exact averaging over the phonon variables is per-formed by the disentangling of the evolution operator (inanalogy with ). As a result, we obtain the formula f zz ( t ) = X n |h ψ | z | ψ n, , i| e − i Ω n, t × D ψ n, , (cid:12)(cid:12)(cid:12) T e exp (cid:16) ˆΦ (cid:17)(cid:12)(cid:12)(cid:12) ψ n, , E (57)with the “influence phase” (assuming ~ = 1 and ω = 1)ˆΦ = − Z t ds Z s ds ′ e − i ( s − s ′ ) X q ˆ w q ( s ) ˆ w + q ( s ′ ) , (58)and T e the time-ordering symbol with respect to the elec-tron degrees of freedom. The correlation function (57) isthe basis expression for the further treatment.The next approximation is the restriction to theleading-order semi-invariant expansion: D ψ n, , (cid:12)(cid:12)(cid:12) T e exp (cid:16) ˆΦ (cid:17)(cid:12)(cid:12)(cid:12) ψ n, , E ≈ exp D ψ n, , (cid:12)(cid:12)(cid:12) T e (cid:16) ˆΦ (cid:17)(cid:12)(cid:12)(cid:12) ψ n, , E . (59)As shown in Ref. , this approximation accounts of thestatic Jahn-Teller effect, and it works well, because thedynamic Jahn-Teller effect appears to be very small. Theinfluence phase is invariant under spatial rotations sothat D ψ n, , (cid:12)(cid:12)(cid:12) T e (cid:16) ˆΦ (cid:17)(cid:12)(cid:12)(cid:12) ψ n, , E = D ψ n, , (cid:12)(cid:12)(cid:12) T e (cid:16) ˆΦ (cid:17)(cid:12)(cid:12)(cid:12) ψ n, , E = D ψ n, , − (cid:12)(cid:12)(cid:12) T e (cid:16) ˆΦ (cid:17)(cid:12)(cid:12)(cid:12) ψ n, , − E . (60)Hence the correlation function (57) can be simplified to f zz ( t ) = X n |h ψ | ˆ z | ψ n, , i| × exp (cid:16) − i Ω n, t − X q X n ′ ,l ′ ,m ′ ,m |h ψ n, ,m | ˆ w q | ψ n ′ ,l ′ ,m ′ i| × − iω n ′ ,l ′ ; n, t − e − iω n ′ ,l ′ ; n, t ω n ′ ,l ′ ; n, ! . (61)with the notation ω n ′ ,l ′ ; n, ≡ ε n ′ ,l ′ − ε n, . (62)In our previous treatments of the strong coupling po-laron optical conductivity, we neglected the matrix el-ements for ˆ w q between the electron energy levels withdifferent energies. Physically, this approximation corre-sponds to the adiabatic approximation. Here, we go be-yond this approximation, taking into account the transi-tions between different excited states but still assuming(as before) that the adiabatic approximation holds forthe transitions between the ground and excited states.In other words, we keep in the sum P n ′ ,l ′ ,m ′ ,m in (61)only the excited states.Introducing parameters related to the extension of theHuang-Rhys parameter used in Ref. : S n ′ ,l ; n, ≡ ω n ′ ,l ; n, X q X m ′ ,m |h ψ n, ,m | ˆ w q | ψ n ′ ,l,m ′ i| , (63)the correlation function is rewritten as follows: f zz ( t )= X n |h ψ | z | ψ n, , i| exp − i Ω n, t − X n ′ ,l S n ′ ,l ; n, × (cid:0) − iω n ′ ,l ; n, t − e − iω n ′ ,l ; n, t (cid:1)i . (64)The states | ψ n ′ ,l,m ′ i can be subdivided to two groups:(1) the states | ψ , ,m ′ i with the energy level ε , , (2) thehigher energy states with ( n ′ , l ) = (1 , . Taking into account thesecond group of states provides the step beyond the adi-abatic approximation – this is the focus of the presenttreatment. We denote the parameters corresponding tothe adiabatic approximation by S n ≡ S n, n, ≡ X q X m ′ ,m |h ψ n, ,m | ˆ w q | ψ n, ,m ′ i| . (65)Correspondingly, the correlation function (64) takes the form f zz ( t ) = X n |h ψ | z | ψ n, , i| × exp h − i Ω n, t − S n (cid:0) − it − e − it (cid:1) − X ( n ′ ,l ) =( n, S n ′ ,l ; n, × (cid:0) − iω n ′ ,l ; n, t − e − iω n ′ ,l ; n, t (cid:1)(cid:3) . (66)Following Refs. , the factor (cid:0) − it − e − it (cid:1) is ex-panded in powers of t :1 − it − e − it = 12 t + O (cid:0) t (cid:1) . (67)This approximation is appropriate in the strong couplingregime, where the phonon frequency is small with re-spect to other frequencies, such as the Franck-Condonfrequency Ω , , which increases as Ω , ∝ α at large α . In the case when non-adiabatic terms are non takeninto account, the expansion (67) provides an envelopeof the optical conductivity spectrum obtained in .The other factor, (cid:0) − iω n ′ ,l ; n, t − e − iω n ′ ,l ; n, t (cid:1) , shouldnot be expanded in the same way, because the frequen-cies ω n ′ ,l ; n, ( n ′ , l ) = (1 ,
1) also increase in the strongcoupling limit as α . Therefore we keep this factor as is,without expansion. As a result, in the strong couplingregime the correlation function f zz ( t ) is: f zz ( t ) = X n |h ψ | z | ψ n, , i| × exp (cid:18) − δS n − i ˜Ω n, t − S n t + X ( n ′ ,l ) =( n, S n ′ ,l ; n, e − iω n ′ ,l ; n, t . (68)with the parameters: δS n ≡ X ( n ′ ,l ) =(1 , S n ′ ,l ; n, , (69) δ Ω n ≡ X ( n ′ ,l ) =(1 , S n ′ ,l ; n, ω n ′ ,l ; n, , (70)˜Ω n, ≡ Ω n, − δ Ω n . (71)The parameter δS n plays a role of the Debye-Waller fac-tor and ensures the fulfilment of the f -sum rule for theoptical conductivity. The parameter δ Ω n is the shiftof the Franck-Condon frequency to a lower value dueto phonon-assisted transitions to higher energy states.The exponent can be expanded, yielding a description interms of multiphonon processes:exp X ( n ′ ,l ) =( n, S n ′ ,l ; n, e − iω n ′ ,l ; n, t = X { p n ′ ,l ≥ } Y ( n ′ ,l ) =( n, S p n ′ ,l ; n, n ′ ,l ; n, p n ′ ,l ; n, ! × e − i P n ′ ,l p n ′ ,l ; n, ω n ′ ,l ; n, t , (72)where the sum P { p n ′ ,l } is performed over all combina-tions { p n ′ ,l ≥ } .With the expansion (72), the polaron optical conduc-tivity takes the form:Re σ (Ω)= Ω X n |h ψ | z | ψ n, , i| e − δS n r π S n × X { p n ′ ,l ; n, ≥ } Y ( n ′ ,l ) =( n, S p n ′ ,l ; n, n ′ ,l ; n, p n ′ ,l ; n, ! × exp − (cid:16) ˜Ω n, + P n ′ ,l p n ′ ,l ; n, ω n ′ ,l ; n, − Ω (cid:17) S n . (73)In formula (73), the term where all p n ′ ,l ; n, = 0 corre-sponds to the adiabatic approximation and exactly re-produces the result of Ref. . The other terms representthe non-adiabatic contributions to Re σ (Ω), and are cor-rection terms to the previously found results. III. RESULTS AND DISCUSSIONS
The polaron optical conductivity derived in the abovesection is in line with the physical understanding of theunderlying processes for the polaron optical response,achieved in early works and summarized in Ref. .It is based on the concept of the polaron excitations ofthree types: • Relaxed Excited States (RES) for which the lat-tice polarization is adapted to the electronic distri-bution; • Franck-Condon states (FC) where the lattice polar-ization is “frozen”, adapted to the polaron groundstate; • Scattering states characterized by the presence ofreal phonons along with the polaron.The polaron RES can be formed when the electron-phonon coupling is strong enough, for α ' .
5. At weakcoupling, the polaron optical response at zero tempera-ture is due to transitions from the polaron ground state to scattering states. In other words, the optical absorptionspectrum of a weak-coupling polaron is determined bythe absorption of radiation energy, which is re-emitted inthe form of LO phonons. At stronger couplings, the con-cept of the polaron relaxed excited states first introducedin Ref. becomes of key importance. In the range of suf-ficiently large α when the polaron RES are formed, theabsorption of light by a polaron occurs through transi-tions from the ground state to RES which can be accom-panied by the emission of different numbers n ≥ . Within the framework of formalismsbased on the memory function (MF), we compare thefollowing theories: • The original DSG method of Ref. , where the ex-pectation value in 31 is calculated with respect toa Gaussian trial action. This will be denoted by MF-1 in the figures. • The extended MF formalism of , where an ad-hocbroadening with a strength determined from sumrules is included in (23). This will be denoted by MF-2 . • The current non-quadratic MF formalism, basedon the extension of the Jensen-Feynman inequalityintroduced in this paper, denoted by
MF-new .Among the strong-coupling expansions (SCE), we distin-guish: • The strong-coupling result in the adiabatic approx-imation, as obtained in Ref. . This will be denotedhere by SCE-1 . • The adiabatic approximation of Ref. , which usesmore accurate trial polaron states. This will be de-noted by SCE-2 . • The current non-adiabatic strong coupling expan-sion, denoted by
SCE-new .The subsequent figures show the results for increas-ing α . In Figure 1, the optical conductivity is shown forsmall coupling, α = 1 , and for α = 3 , .
25 which corre-spond to the dynamic regime where the RES starts toplay a role. In this regime, analytic solutions are pro-vided by the various memory function formalisms listedabove, and we compare them to DQMC numeric data .At weak coupling ( α = 1 , panel ( a )), all the approachesbased on the memory function give results in agreementwith DQMC. For α = 3 (panel ( b )), the current methodgives a better fit to the DQMC result that the other two0 c a = 5.25 MF- new MF-2 MF-1 DQMC R e s ( w ) [ i n un it s n e / ( m b w L O )] W / w b a = 3 MF- new MF-2 MF-1 DQMC R e s ( w ) [ i n un it s n e / ( m b w L O )] W / w a = 1 MF- new MF-2 MF-1 DQMC R e s ( w ) [ i n un it s n e / ( m b w L O )] W / w a FIG. 1: Polaron optical conductivity calculated for α = 1 ( a ), α = 3 ( b ) and α = 5 .
25 ( c ) within the present non-quadraticMF formalism (denoted in the figure as MF- new ), comparedwith the polaron optical conductivity calculated within theextended memory-function formalism (MF-2) of Ref. , theresults of the memory-function approach using the Feynmanparabolic trial action (MF-1), and the diagrammatic quan-tum Monte Carlo (DQMC) . methods. For a stronger coupling, α = 5 .
25 (panel ( c ))the MF-2 approach substantially improves the originalresult MF-1, but the optical conductivity spectrum cal-culated within the new non-quadratic MF formalism liescloser to the DQMC data than either of the other two.Fig. 2 demonstrates the behavior of the polaron op-tical conductivity spectra in the intermediate couplingregime, for α = 6 . α = 7. In this regime, the ex-isting memory function approaches ( MF-1 , MF-2 ) as wellas the existing strong coupling expansions (
SCE-1 , SCE- 2 ) do not provide satisfactory results. The new memoryfunction approach and the new strong coupling expansionare in much better agreement with the DQMC data. W / w a = 7 SCE- new SCE-2 SCE-1 MF-2 MF-1 DQMC R e s ( w ) [ i n un it s n e / ( m b w L O )] W RES W FC b a W / w a = 6.5 MF- new SCE- new
MF-2 MF-1 DQMC R e s ( w ) [ i n un it s n e / ( m b w L O )] FIG. 2: Polaron optical conductivity calculated for α = 6 . a ) and α = 7 ( b ) using different analytic approaches:the non-quadratic MF formalism (MF- new ), the extendedmemory-function formalism of Ref. (MF-2), the memory-function approach with the Feynman parabolic trial action (MF-1), the non-adiabatic strong-coupling expansion (de-noted at the figure as SCE- new ), the adiabatic strong-coupling expansions of Refs. (SCE-1 and SCE-2). Theresults are compared to DQMC data of Refs. . This range of coupling parameters is where one wouldwant to cross over from using a memory function basedapproach to a strong coupling expansion. Whereas theexisting methods do not allow to bridge this gap at in-termediate coupling, the extensions that we have pro-posed here are suited to implement such a cross-over.The present memory-function approach with the non-parabolic trial action leads to a relatively small exten-sion of the range of α where the polaron optical conduc-tivity compares well with the DQMC data, namely from α ≈ . α ≈ .
5. For α / .
5, the memory-functionapproach with the non-parabolic trial action providesa better agreement with DQMC than all other knownapproximations. Remarkably, the optical conductivityspectra as given by the non-quadratic MF formalism andthe non-adiabatic SCE are both in better agreement withthe Monte Carlo data than any of the preceding analyti-cal methods. For α = 6 .
5, the polaron optical conductiv-ity calculated within non-quadratic MF formalism and1the non-adiabatic SCE lie rather close to each other. Wecan conclude therefore that the ranges of validity of thosetwo approximations overlap, despite the fact that theseapproximations are based on different assumptions.The maximum of the optical conductivity spectrumprovided by the non-quadratic MF formalism for α = 6 . , where it was suggested thatthese two features in the DQMC spectra can correspondphysically to the dynamic (RES) and the Franck-Condoncontributions. The present results are in line with thatphysical picture.In Fig. 2 ( b ), the arrows indicate the FC transitionfrequency for the transition to the first excited FC stateΩ , ≡ Ω FC and the RES transition frequency Ω RES fora strong coupling polaron as calculated in Ref. . Wecan see that both the shape and the position of the max-imum of the optical conductivity band obtained withinthe adiabatic approximation in Refs. are rather farfrom those for the DQMC data. Taking into accountnon-adiabatic transitions drastically improves the agree-ment of the strong coupling approximation with DQMC,even for α = 7, which, strictly speaking, is not yetthe strong coupling regime. The value α = 7 can berather estimated as an intermediate coupling. However,even at this intermediate coupling strength, the resultsof present approach lie much closer to the DQMC datathan those obtained within all other aforesaid analyticmethods. Also a substantial improvement of the agree-ment between the strong coupling expansion and DQMCis clearly expressed in Fig. 3, where the polaron opticalconductivity spectra are shown for the strong couplingregime for α = 8 to α = 9. For strong couplings, thenon-adiabatic SCE accurately reproduces both the peakposition and the overall shape of the DQMC spectra. Fi-nally, we see that the results of the non-adiabatic SCEremain accurate also in the extremely strong couplingregime, as shown in Fig. 4. IV. CONCLUSIONS
In the present work, we have modified two basic an-alytic methods for the polaron optical conductivity inorder to extend their ranges of applicability for theelectron-phonon coupling constant in such a way thatthese ranges overlap. The memory function formalismusing a trial action for a model two-particle system hasbeen extended to work with non-quadratic interactionpotentials in the model system. This method combinesthe translation invariance of the trial system, which is W / w W FC W RES a = 9 SCE- new SCE-2 SCE-1 DQMC c R e s ( w ) [ i n un it s n e / ( m b w L O )] W RES W FC a = 8.5 SCE- new SCE-2 SCE-1 DQMC b R e s ( w ) [ i n un it s n e / ( m b w L O )] W / w W / w W RES a = 8 SCE- new SCE-2 SCE-1 DQMC a R e s ( w ) [ i n un it s n e / ( m b w L O )] W FC FIG. 3: Polaron optical conductivity calculated for α = 8 ( a ), α = 8 . b ) and α = 9 ( c ) within several analytic strong cou-pling approaches and compared to DQMC data of Refs. .The notations are the same as in Fig. 2. one of the main advantages of the Feynman variationalapproach, with a more realistic interaction between theelectron and the fictitious particle. This extension leadsto a substantial improvement of the polaron optical con-ductivity for small and intermediate coupling strengthswith respect to the preceding known versions of the mem-ory function approach.The other method is the strong-coupling expansion,and we have extended it beyond the Franck-Condonadiabatic approximation by taking into account non-adiabatic transitions between different excited polaronstates. As a result, the modified non-adiabatic strong-coupling expansion appears now to be in good agreementwith the numerical DQMC data in a wide range of α
15 20 25 30 35 400.000.050.100.150.20 a = 15 SCE- new SCE-2 SCE-1 DQMC b R e s ( w ) [ i n un it s n e / ( m b w L O )] W / w W RES W FC
10 15 20 25 300.000.050.100.150.20 W / w a = 13 SCE- new SCE-2 SCE-1 DQMC a R e s ( w ) [ i n un it s n e / ( m b w L O )] W FC W RES
FIG. 4: Polaron optical conductivity in the extremely strongcoupling regime, for α = 13 ( a ) and α = 15 ( b ). The notationsare the same as in Fig. 2. from intermediate coupling strength to the strong cou-pling limit. For the intermediate coupling value α = 6 . α the agreement between the results of thenon-adiabatic SCE and DQMC becomes gradually bet- ter. At very strong coupling, even the preceding adi-abatic SCE is already sufficiently good, so that theimprovement due to the non-adiabatic transitions, e. g.,for α = 15, is relatively small. However, for a slightlyweaker coupling, e. g., for α = 9, we can observea drastically improved agreement with DQMC for thepresent non-adiabatic SCE as compared to the adiabaticapproximation. We can conclude that at present, thestrong coupling approximation taking into account non-adiabatic contributions provides the best agreement withthe DQMC results for α ' . α ≤
15, available for DQMC.In summary, extending the MF and SCE formalismsleads to an overlapping of the areas of α where these twoanalytic methods are applicable. These analytic methodshave been verified, appearing to be in good agreementwith numeric DQMC data at all α available for DQMC.We therefore possess the analytic description of the po-laron optical response which embraces the whole rangeof the coupling strength. Acknowledgments
We thank A. S. Mishchenko for valuable discussionsand the DQMC data for the polaron optical conduc-tivity, and V. Cataudella for the details of the EMFFmethod. Discussions with F. Brosens and D. Sels aregratefully acknowledged. This research has been sup-ported by the Flemish Research Foundation (FWO-Vl), project nrs. G.0115.12N, G.0119.12N, G.0122.12N,G.0429.15N, by the Scientific Research Network of theResearch Foundation-Flanders, WO.033.09N, and by theResearch Fund of the University of Antwerp. ∗ On leave from Department of Theoretical Physics, StateUniversity of Moldova, MD-2009 Kishinev, Republic ofMoldova. † Lyman Laboratory of Physics, Harvard University ‡ Technische Universiteit Eindhoven, P.O. Box 513, 5600 MBEindhoven, The Netherlands. L. D. Landau,
Phys. Z. Sowjetunion , 664 (1933) [Englishtranslation in Collected Papers , Gordon and Breach, NewYork, 1965, pp. 67-68]. A. S. Alexandrov and J. T. Devreese,
Advances in PolaronPhysics (Springer, 2009). R. von Helmolt, J. Wecker, B. Holzapfel, L. Schultz, andK. Samwer, Phys. Rev. Lett. , 2331 (1993). H. Sirringhaus et al ., Nature (London) , 685 (1999). T. Holstein, Ann. Phys. (N.Y.) , 325 (1959). M. Setvin, C. Franchini, X. Hao, M. Schmid, A. Janotti,M. Kaltak, C. G. Van de Walle, G. Kresse, and U. Diebold,Phys. Rev. Lett. , 086402 (2014). X. Hao, Z. Wang, M. Schmid, U. Diebold, and C. Fran-chini, Phys. Rev. B , 085204 (2015). J. Vlietinck, W. Casteels, K. Van Houcke, J. Tempere, J.Ryckebusch, and J. T. Devreese, New J. Phys. , 033023(2015). W. Meevasana, X. J. Zhou, B. Moritz, C.-C. Chen, R. H.He, S.-I. Fujimori, D. H. Lu, S.-K. Mo, R. G. Moore, F.Baumberger, T. P. Devereaux, D. van der Marel, N. Na-gaosa, J. Zaanen and Z.-X. Shen, New Journal of Physics , 023004 (2010). J. L. M. van Mechelen, D. van der Marel, C. Grimaldi, A.B. Kuzmenko, N. P. Armitage, N. Reyren, H. Hagemann, and I. I. Mazin, Phys. Rev. Lett. , 226403 (2008). J. T. Devreese, S. N. Klimin, J. L. M. van Mechelen, andD. van der Marel, Phys. Rev. B , 125119 (2010). A. S. Mishchenko, N. V. Prokof’ev, A. Sakamoto, andB. V. Svistunov, Phys. Rev. B , 6317 (2000). A. S. Mishchenko, N. Nagaosa, N. V. Prokof’ev, A.Sakamoto, and B. V. Svistunov, Phys. Rev. Lett. ,236401 (2003). G. L. Goodvin, A. S. Mishchenko, and M. Berciu, Phys.Rev. Lett. , 076403 (2011). E. Kartheuser, R. Evrard, and J. Devreese Phys. Rev. Lett. , 94-97 (1969). J. Devreese, J. De Sitter, and M. Goovaerts, Phys. Rev. B , 2367 (1972). J. T. Devreese, in
Polarons in Ionic Crystals and PolarSemiconductors (North-Holland, Amsterdam, 1972), pp.83 – 159. G. De Filippis, V. Cataudella, A. S. Mishchenko, C. A.Perroni, and J. T. Devreese, Phys. Rev. Lett. , 136405(2006). V. L. Gurevich, I. G. Lang, and Yu. A. Firsov, Sov. Phys.Solid State , 918 (1962). J. Devreese, W. Huybrechts, and L. Lemmens, Phys. Sta-tus Solidi B , 77 (1971). S. N. Klimin and J. T. Devreese, Phys. Rev. B , 035201(2014). R. P. Feynman, R. W. Hellwarth, C. K. Iddings, and P. M.Platzman, Phys. Rev. , 1004 (1962). R. P. Feynman, Phys. Rev. , 660 (1955). D. Sels and F. Brosens, Phys. Rev. E , 012124 (2014). D. Sels and F. Brosens, Phys. Rev. E , 042110 (2014). D. Sels, arXiv:1605.04998 (2016). V. Cataudella, G. De Filippis, and C.A. Perroni, “SinglePolaron Properties in Different Electron-Phonon Models”,in:
Polarons in Advanced Materials , ed. by A. S. Alexan-drov, Springer Series in Materials Science, Volume 103,2007, pp. 149-189. S. N. Klimin and J. T. Devreese, Solid State Communica-tions , 144 (2011). S. J. Miyake, J. Phys. Soc. Japan , 181 (1975). S. I. Pekar,
Untersuchungen ¨uber die Elektronentheorie derKristalle (Akademie-Verlag, Berlin, 1954). G. R. Allcock, in
Polarons and Excitons , edited by C. G.Kuper and G. D. Whitfield (Oliver and Boyd, Edinburgh,1963), pp. 45 – 70. A. S. Mishchenko, N. Nagaosa, and N. Prokof’ev, Phys.Rev. Lett. , 166402 (2014). F. M. Peeters and J. T. Devreese, Phys. Rev. B , 6051(1983). R. P. Feynman, Phys. Rev.84