# Polaron Mass and Electron-Phonon Correlations in the Holstein Model

Authors: Marco Zoli
Journal: Advances in Condensed Matter Physics
 (2010)

---

## Abstract

The Holstein
Molecular Crystal Model is investigated by a
strong coupling perturbative method which,
unlike the standard Lang-Firsov approach,
accounts for retardation effects due to the
spreading of the polaron size. The effective
mass is calculated to the second perturbative
order in any lattice dimensionality for a broad
range of (anti)adiabatic regimes and
electron-phonon couplings. The crossover from a
large to a small polaron state is found in all
dimensionalities for adiabatic and intermediate
adiabatic regimes. The phonon dispersion largely
smoothes such crossover which
is signalled by polaron mass enhancement and on-site localization
of the correlation function. The notion of self-trapping together
with the conditions for the existence of light polarons, mainly in
two- and three-dimensions, is discussed. By the imaginary time
path integral formalism I show how nonlocal electron-phonon
correlations, due to dispersive phonons, renormalize downwards thee-ph coupling justifying the possibility for
light and essentially small 2D Holstein polarons.


---

## Body

## 1. Introduction

The interest for phonons and lattice distortions in High-Temperature Superconductors (HTSc) is today more than alive [1, 2]. While the microscopic origin of the pairing mechanism in cuprate HTSc has not yet been unravelled, evidence has been provided [3–8] that electron-lattice interactions and local lattice fluctuations are correlated with the onset of the superconducting transition [9–11]. Recent investigations [12] would suggest a similar role for the lattice also in the newly discovered layered pnictide-oxyde quaternary superconducting compounds [13–15]. In fact, the discovery of HTSc in copper oxides [16] was motivated by the knowledge that the Jahn-Teller effect [17–19] is strong in Cu2+. In nonlinear molecules a lattice distortion lifts the degeneracy of the electronic states and lowers the overall ground state energy: the energy gain is the Jahn-Teller stabilization energy which can be measured optically. As vibrational and electronic energies are of the same order, the nuclear motion cannot be separated from the electrons and the combined electron-lattice object becomes a mobile Jahn-Teller polaron [20] which can be described in terms of the nonlinear Holstein Hamiltonian [21]. Over the last twenty years, the focus on the HTSc has largely contributed to trigger the study of the polaron properties [22–25]. Earlier and more generally, polarons had become a significant branch in condensed matter physics [26–29] as Landau introduced the concept of an electron which can be trapped by digging its own hole in an ionic crystal [30]. A strong coupling of the electron to its own lattice deformation implies violation of the Migdal theorem [31] and polaron collapse of the electron band [32] with the appearance of time retarded interactions in the system. A great advance in the field came after Feynman used the path integrals [33] to calculate energy and effective mass of the Fröhlich polaron [34, 35] by a variational method [36]. In the statistical path integral, the quantum mechanical electron motion is represented by a real space path x(τ) being function of an imaginary time τ whose range has an upper bound given by the inverse temperature. As a key feature of the Feynman treatment, the phonon degrees of freedom can be integrated out exactly (being the action quadratic and linear in the oscillator coordinates) but leave a substantial trace as a time retarded potential naturally emerges in the exact integral action. The electron, at any imaginary time, interacts with its position at a past time. This self interaction mirrors the fact that the lattice needs some time to readjust to the deformation induced by the electron motion and follow it. For decades, the Fröhlich polaron problem has been extensively treated by path integral techniques [37–39] which, remarkably, can be applied for any value of the coupling constant [40]. A review work on the Fröhlich polaron is in [41].With regard to small polarons, path integral investigations started with the groundbreaking numerical work by De Raedt and Lagendijk [42–44] who derived fundamental properties for the Holstein polaron. While a sizable electron-phonon coupling is a requisite for polaron formation [45–47], also the dimensionality and degree of adiabaticity of the physical system could essentially determine the stability and behavior of the polaron states [48–62]. When the characteristic phonon energy is smaller than the electronic energy the system is set in the adiabatic regime. In materials such as the HTSc, colossal magnetoresistance oxides [63–65], organic molecular crystals [66], DNA molecules [67, 68] and finite quantum structures [69, 70], intermediate adiabatic polarons are relevant as carriers are strongly coupled to high energy optical phonons. For the HTSc systems, a path integral description had been proposed for polaron scattering by anharmonic potentials, due to lattice structure instabilities, as a possible mechanism for the nonmetallic behavior of the c-axis resistivity [71, 72].It has been questioned whether small (bi)polarons could indeed account for highTc due to their large effective mass [73] but, later on, it was recognized that dispersive phonons renormalize the effective coupling in the Holstein model yielding much lighter masses [74]. Also long range Fröhlich e-ph interactions provide a pairing mechanism and bipolaron mass consistent with high Tc [75–78] while accounting, among others, for angle-resolved photoemission spectroscopy data in cuprates HTSc [79, 80].In view of the relevance of the polaron mass issue, I review in this article some work [81] on the Holstein polaron with dispersive phonons which examines the notion of self-trapping and estimates some polaron properties in the parameter space versus lattice dimensionality. The latter is introduced in the formalism by modelling the phonon spectrum through a force constant approach which weighs the first neighbors intermolecular shell. This is done in the spirit of the work by Holstein [82] on the small polaron motion for a one-dimensional narrow band system. Calculating the polaron hopping between localized states due to multiphonon scattering processes in perturbation theory, Holstein first pointed out that the phonon dispersion was in fact a vital ingredient of the theory whereas a dispersionless spectrum would yield a meaningless divergent hopping probability. Time dependent perturbative theory, the hopping integral being the perturbation, applies as long as the interaction time is shorter than the polaron lifetime on a state. Computing the time integral for the hopping probability, I had earlier shown [53] that such interaction time is shorter in higher dimensionality thus making the time dependent perturbative method more appropriate in the latter conditions. For a given dimensionality, the interaction time is also reduced by increasing the strength of the intermolecular forces hence, the phonon dispersion. Consistently also the crossover temperature between band-like motion (at low T) and hopping motion (at high T) [83] shifts upwards, along the T axis, in higher dimensional systems with larger coordination numbers. Hopping polaron motion is the relevant transport mechanism in several systems including DNA chains [84, 85].It is worth emphasizing that our previous investigations on the Holstein polaron and the present paper assume dimensional effects as driven by the lattice while the electron transfer integral is a scalar quantity. In a different picture, first proposed by Emin [86] for quasi 1D solids in the adiabatic limit and recently generalized by variational exact diagonalization techniques [87], the electron transfer integral is instead a vector whose components are switched on in order to treat anisotropic polaron properties. The system dimensionality is then increased by tuning the electronic subsystem parameters, the 3D case corresponding to an isotropic hopping integral, whereas the phonon spectrum is dispersionless. While the method to introduce the system dimensionality in [87] differs from ours, some trends regarding the polaron mass in 1D and 3D are apparently at variance with ours as it will be pointed out in the following discussion.Being central to polaron studies, the effective mass behavior ultimately reflects the abovementioned property of the polaron problem: when the electron moves through the lattice, it induces and drags a lattice deformation which however does not follow instantaneously, the retardation causing a spread in the size of the quasiparticle. The electron-phonon correlation function provides a measure of the deformation around the instantaneous position of the electron and may be used to quantify the polaron size. Several refined theoretical tools, including quantum Monte Carlo calculations [42], density matrix renormalization group [88], variational methods [89, 90] and exact diagonalization techniques [91, 92] have been applied to the Holstein Hamiltonian to compute such polaron characteristics covering various regimes of electron-phonon coupling and adiabatic ratio. While a qualitative and often quantitative [88, 89] agreement among several numerical methods has emerged especially with regard to ground state polaron properties, we feel that analytical investigations may provide a useful insight also in the specific regime of intermediate adiabaticity and e-ph coupling for which the leading orders of perturbation theory traditionally fail in the weak coupling [93] or become less accurate in the strong coupling [94, 95].Extremely interesting is that range ofe-ph couplings which are sufficiently strong to allow for polaron formation but not too strong to prevent the polaron from moving through the lattice. When the relevant energy scales for electrons and phonons are not well separated, such range of couplings produces a composite particle which is not trapped by its lattice deformation. It is this intermediate regime which forms the focus of the present work.The retardation effect and its consequences are here treated by means of a variational analytical method based on the Modified Lang Firsov (MLF) transformation [96] which considerably improves the standard Lang Firsov (LF) [97] approach on which strong coupling perturbation theory (SCPT) is based. Section 2 provides the generalities of the MLF method applied to the dispersive Holstein Hamiltonian. In Section 3, the polaron mass is calculated in the (anti)adiabatic regimes covering a broad range of cases in parameter space. The spreading of the adiabatic polaron size over a few lattice sites is shown in Section 4 through computation of the electron-phonon correlations. In Section 5, I give further physical motivations for the existence of light Holstein polarons: using a path-integral method, I show how the electron-phonon coupling is renormalized downwards in momentum space by nonlocal correlations arising from the dispersive nature of the phonon spectrum. Some final remarks are in Section 6.


## 2. Modified Lang-Firsov Method

The Holstein diatomic molecular model was originally cast [21] in the form of a discrete nonlinear Schrödinger equation for electrons whose probability amplitude at a molecular site depends on the interatomic vibration coordinates [98]. The nonlinearities are tuned by the electron-phonon coupling g, whose strength drives the crossover between a large and a small polaron for a given value of the adiabatic parameter. The polaron radius is measured with respect to the lattice constant [99, 100]. When the size of the lattice distortion is of the same order of (or less than) the lattice constant the polaron has a small radius and the discreteness of the lattice must be taken into account.In second quantization the dimension dependent Holstein Hamiltonian with dispersive harmonic optical phonons reads(1)H=-t∑〈ij〉ci†cj+g∑ini(bi†+bi)+∑qωqbq†bq,
wheret is the hopping integral and the first sum is over z nearest neighbors. ci† and ci are the real space electron creation and annihilation operators at the i-site, ni(=ci†ci) is the number operator, bi† and bi are the phonons creation and annihilation operators. bq† is the Fourier transform of bi† and ωq is the frequency of the phonon with vector momentum q. Unlike the Su-Schrieffer-Heeger Hamiltonian [101, 102] a paradigmatic model in polymer physics, in  (1) the electron hopping does not depend on the relative displacement between adjacent molecular sites hence, the phonon created by bi† is locally coupled to the electronic density.The LF transformation uses a phonon basis of fixed displacements (at the electron residing site) which diagonalizes the Hamiltonian in  (1) in absence of hopping. The hopping term is then treated as a perturbation [74, 94]. However the standard LF approach does not account for the retardation between the electron and the lattice deformations which induces a spread in the size of the polaron. Precisely this effect becomes important for the intermediate e-ph coupling values which may be appropriate for some HTSc.The idea underlying the MLF transformation [103] is that to consider the displacements of the oscillators at different sites around an electron in order to describe the retardation effect. For the present case of dispersive phonon the MLF transformation, applied to the Hamiltonian in  (1), reads
(2)H̃=eSHe-S,S=∑qλqnq(b-q†-bq),nq=1N∑inie-iq·Ri=1N∑kck+q†ck,
where Ri are the lattice vectors, and λq are the variational parameters which represent the shifts of the equilibrium positions of the oscillators (quantized ion vibrations) with momentum q. The conventional Lang-Firsov transformation is recovered by setting λq=g/ωq.Explicitly, the MLF transformed Holstein Hamiltonian in  (2) is
(3)H̃=-ϵp∑ini-tp∑ijci†cj×exp[1N∑qλqbq†(eiq·Ri-eiq·Ri)]×exp[-1N∑qλqbq(e-iq·Ri-e-iq·Ri)]+∑qωqbq†bq+∑q(g-λqωq)nq(b-q†+bq),
where the polaron self-energy ϵp is
(4)ϵp=1N∑q(2g-λqωq)λq,
and the polaronic hopping is
(5)tp=texp[-1N∑qλq2(1-γqz)],γq=2∑i=x,y,zcosqi.The coordination numberz is twice the system dimensionality.Looking at (3), it is clear that the unperturbed Hamiltonian H0 can be taken as
(6)H0=-ϵp∑ini+∑qωqbq†bq,
while the remaining part of the Hamiltonian (H̃-H0) in the MLF basis is considered as the perturbation part.The energy eigenstates ofH0 are given by
(7)|ϕi,{nq}〉=ci†|0〉e|nq1,nq2,nq3,…〉ph,
where i is the electron site, and nq1, nq2, and nq3 are the phonon occupation numbers in the phonon momentum states q1, q2, and q3, respectively. The lowest energy eigenstate of the unperturbed Hamiltonian has no phonon excitations, that is, nq=0 for all q. The ground state has an energy E00=-ϵp and is N-fold degenerate, where N is the number of sites in the system. The perturbation lifts the degeneracy, and to first order in t the ground state energy of the 3D-polaron with momentum k is given by
(8)E0(k)=-ϵp-tpγk.The second-order correction to the ground-state energy of the polaron with momentumk is given by
(9)E0(2)(k)=∑k′‍∑{nq}1∑qnqωq×|〈{nq},k′|H̃-H0|k,{0}〉|2,
where, in principle, intermediate states having all possible phonon numbers contribute to  (9).By minimizing the zone center ground state energy, the variational parametersλq are obtained as
(10)λq=gωq+ztp(1-γq/z),
and, by  (10), the one phonon matrix element between the ground state |k=0,{nq=0}〉 and the first excited state 〈1q,k′| vanishes. Then, the one phonon excitation process yields no contribution to  (9). By  (10) one gets the λq's for the 1D, 2D, and 3D systems and, via  (5), the narrowing of the polaron band [104].


## 3. Polaron Mass

Given the formal background, the polaron mass is calculated both for the Lang-Firsov and for the Modified Lang-Firsov method to the second order in SCPT. Generally the second order correction (i) makes a relevant contribution to the ground state energy, (ii) does not affect the bandwidth, (iii) introduces the mass dependence on the adiabatic parameter which would be absent in the first order. Altogether the polaron landscape introduced by the second order SCPT is much more articulated than the simple picture suggested by the first order of SCPT which nonetheless retains its validity towards the antiadiabatic limit. In first order SCPT, ground state energy, bandwidth, and effective mass appear as equivalent, interchangeable properties to describe the polaron state while band narrowing and abrupt mass enhancement arecoincident signatures of small polaron formation in parameter space. However, there is no compelling physical reason for such coincidence to occur as the bandwidth is in fact a zone edge quantity whereas the effective mass is a zone center quantity. Thus, by decoupling onset of the band narrowing and self trapping of the polaron mass, the second order of SPCT provides the framework for a richer polaron structure in momentum space.As emphasized in Section1, dispersive phonons have a relevant role in the Holstein model: a dispersionless spectrum would in fact predict larger polaron bandwidths in lower dimensionality and yield a divergent site jump probability for the small polaron in time dependent perturbation theory [82, 105]. Here, a lattice model is assumed in which first neighbors molecular sites interact through a force constants pair potential. Then, the optical phonon spectrum is given in 1D, 2D (square lattice), and 3D (simple cubic lattice), respectively, by
(11)ω1D2(q)=α+γM+1Mα2+2αγcosq+γ2,ω2D2(q)=α+2γM+1Mα2+2αγg(q)+γ2(2+h(q)),ω3D2(q)=α+3γM+1Mα2+2αγj(q)+γ2(3+l(q)),g(q)=cosqx+cosqy,h(q)=2cos(qx-qy),j(q)=cosqx+cosqy+cosqz,l(q)=2cos(qx-qy)+2cos(qx-qz)+2cos(qy-qz),
where α is the intramolecular force constant and γ is the intermolecular first neighbors force constant. M is the reduced molecular mass.  Thus, the intra- and intermolecular energies are ω0=2α/M and ω1=γ/M, respectively. Some care should be taken in setting the phonon energies as the second order perturbative term grows faster than the first order by increasing ω1 [74]. Hence too large dispersions may cause a breakdown of the SCPT. In terms of ω0, the dimensionless parameter zt/ω0 defines the adiabatic (zt/ω0>1) and the antiadiabatic (zt/ω0<1) regime. Some other choices are found in literature with t/ω0>1 or t/ω0>1/4 defining the adiabatic regime. As shown below, such discrepancies may lead to significantly different interpretations of the polaron behavior in parameter space mainly for higher dimensions.Hereafter I takeω0=100 meV and tune t to select a set of zt/ω0 which sample both antiadiabatic and adiabatic regime without reaching the adiabatic limit. The dynamics of the lattice is in fact central to our investigation. A moderate to strong  range of e-ph  couplings is assumed (g/ω0≽1) so that the general conditions for polaron formation are fulfilled throughout the range of adiabatic ratios [22, 95].Figure1 plots the ratio of the one-dimensional polaron mass to the bare band mass, against the e-ph  coupling, calculated both in the Lang-Firsov scheme and in the Modified Lang-Firsov expression. The intermediate regime 2t=ω0 is assumed in Figure 1(a). While at very strong couplings the MLF plots converge towards the LF predictions, a remarkably different behavior between the LF and the MLF mass shows up for moderate g. The LF method overestimates the polaron mass for g∈[~1-2], and mostly, it does not capture the rapid mass increase found instead in the MLF description. Note that, around the crossover, the MLF polaron mass is of order ten times the bare band mass choosing ω1=60 meV. Large intermolecular energies enhance the phonon spectrum thus reducing the effective masses in both methods. In the MLF method, large ω1 tend also to smooth the mass behavior in the crossover region.Ratio of the one-dimensional polaron mass to the bare band mass versuse-ph coupling according to the Lang-Firsov and the Modified Lang-Firsov methods. The adiabatic parameter is set at (a) the intermediate value, 2t/ω0=1; (b) a fully adiabatic regime, 2t/ω0=2; and (c) an antiadiabatic regime, 2t/ω0=0.25. ω0=100 meV and ω1 (in units meV) are the intramolecular and intermolecular energies of the diatomic molecular chain, respectively.
(a)(b)(c)Figure1(b) presents the case of an adiabatic regime: the discrepancies between LF and MLF plots are relevant for a broad range of e-ph couplings. The latter is more extended for larger intermolecular energies as also seen in Figure 1(a). Mass renormalization is poor in the MLF curves up to the crossover which is clearly signalled by a sudden although continuous mass enhancement whose abruptness is significantly smoothed for the largest values of intermolecular energies.Figure1(c) shows a fully antiadiabatic case in which the LF and MLF plots practically overlap throughout the whole range of couplings. This confirms that the LF method is appropriate in the antiadiabatic limit which is essentially free from retardation effects. Unlike the previous cases the MLF plots are always smooth versus g indicating that no self-trapping event occurs in the antiadiabatic regime. Here the polaron does not trap as it is always small for the whole range of couplings.It should be reminded that the concept ofself-trapping traditionally indicates an abrupt transition between an infinite size state at weak e-ph couplings and a finite (small) size polaron at strong e-ph couplings. In one-dimension, according to the traditional adiabatic polaron theory [106–108], the polaron solution is always the ground state of the system and no self-trapping occurs. Instead, in higher dimensionality a minimum coupling strength is required to form finite size polarons, hence self-trapped polarons can exist at couplings larger than that minimum.However, these conclusions have been critically reexamined in the recent polaron literature and the same notion ofpolaron size has been questioned due to the complexity of the polaron quasiparticle itself [109, 110]. Certainly, as a shrinking of the polaron size yields a weight increase, the polaron mass appears as the most reliable indicator of the self-trapping transition. The latter, as the results in Figure 1 suggest, is rather a crossover essentially dependent on the degree of adiabaticity of the system and crucially shaped by the internal structure of the phonon cloud which I have modelled by tuning the intermolecular forces. In this view, self trapping events can be found also in the parameter space of 1D systems. This amounts to say that also finite size polarons can self-trap if a sudden, although continuous, change in their effective mass occurs for some values of the e-ph couplings in some portions of the intermediate/adiabatic regime. The continuity of such event follows from the fact that the ground state energy is analytic function of g for optical phonon dispersions [111]; hence the possibility of phase transitions is ruled out in the Holstein model. At finite temperature phase transitions cannot occur as the free energy is smooth for any phonon dispersion [112–114].As fluctuations in the lattice displacements around the electron site are included in the MLF variational wavefunction, the calculated polaron mass should not display discontinuities by varying the Hamiltonian parameters through the crossover [115]. Mathematically the crossover points are selected through the simultaneous occurrence of a maximum in the first logarithmic derivative and a zero in the second logarithmic derivative of the MLF polaron mass with respect to the coupling parameter. Such inflection points, corresponding to the points of the most rapid increase for m*, are marked in Figure 2 on the plots of the mass ratios computed for a broad range of (anti)adiabatic parameters both in one, two, and three dimensions.Ratio of the Modified Lang-Firsov polaron mass to the bare band mass versuse-ph coupling in (a) 1D, (b) 2D, and (c) 3D. A set of twelve zt/ω0 values ranging from the antiadiabatic to the adiabatic regime is considered. From left to right: zt/ω0=0.25,0.5,0.75,1.0,1.25,1.5,1.75,2.0,2.25,2.5,2.75,3.0. ω0=100 meV. The diamonds mark the occurrence of the self-trapping event.
(a)(b)(c)In 1D, see Figure2(a), the crossover occurs for g values between ~1.8 and 2.3 and the corresponding self-trapped polaron masses are of order ~5–50 times the bare band mass thus suggesting that relatively light small polarons can exist in 1D molecular solids with high phonon spectrum. The onset of the self-trapping line is set at the intermediate value 2t=ω0 and the self-trapped mass values grow versus g by increasing the degree of adiabaticity. There is no self-trapping in the fully antiadiabatic regime as the electron and the dragged phonon cloud form a compact unit also at moderate e-ph  couplings. Then, the mass increase is always smooth in the antiadiabatic regime.Some significant results are found in 2D as shown in Figure2(b): (i) for a given g and adiabaticity ratio, the 2D mass is lighter than the 1D mass and the 2D LF limit is attained at a value which is roughly one order of magnitude smaller than in 1D; (ii) the crossover region is shifted upwards along the g axis with the self trapping events taking place in the range, g~2.2–2.6; the masses are of order ~5–10 times the bare band mass; (iii) the line of self-trapping events changes considerably with respect to the 1D plot; the marked curve is parabolic in 2D with an extended descending branch starting at the intermediate value 4t/ω0=1; (iv) in the deep adiabatic regime, the lattice dimensionality smoothens the mass increase versus g.This effect is even more evident in 3D, see Figure2(c), as there are no signs of abrupt mass increase even for the largest values of the adiabatic parameter. At the crossover, 3D masses are of order ~5–10 times the bare band mass with the self trapping points lying in the range, g~2.5–2.9. At very large couplings the mass ratio becomes independent of t and converges towards the LF value. In this region (and for the choice ω1=60 meV) the 3D Lang-Firsov mass is one order of magnitude smaller than the 2D mass. As the coordination number grows versus dimensionality, large intermolecular forces are more effective in hardening the 3D phonon spectrum thus leading to lighter 3D polaron masses than 2D ones.The self trapping transition appears in Figure2 as smoother in higher dimensions thus contradicting the trend found by previous investigations, markedly by exact diagonalization techniques [60, 87] and variational calculations [116]. Two reasons may be invoked to resolve the discrepancy, the first being more physical and the second more technical. (1) From  (11), it can be easily seen that ωd2[q=(0,0,0)]-ωd2[q=(π,π,π)]=2dω12. Then, for a given ω1, the phonon band is more dispersive in higher dimensionality d. Hence increasing the system dimension corresponds to attribute larger weight to the intermolecular forces which ultimately smoothen the crossover as made evident in Figure 1. (2) The works in [60, 87, 116] take t/ω0 as adiabatic ratio in any dimension whereas our calculation spans the same range of zt/ω0 values in any dimension. Thus, for instance, the fully adiabatic zt/ω0=3 plots in Figure 2 would correspond to an antiadiabatic case in 3D (t/ω0=1/2) and a slightly adiabatic case in 1D (t/ω0=3/2) following the definition in [60, 87, 116]. Consistently, the 3D antiadiabatic polaron mass is smoother than the 1D adiabatic polaron mass. Vice versa, assuming the same ratio, that is, t/ω0=1/2, one should compare the rightmost plot in Figure 2(c) with the fourth from left in Figure 2(a): the latter is smoother than the former as the 1D case is now antiadiabatic while the 3D case is fully adiabatic.Nonetheless, the findings displayed in Figure2 agree with density matrix renormalization group [88], variational Hilbert space [60] methods and weak coupling perturbation theory [57] in predicting a lighter mass in higher dimensionality.


## 4. Electron-Phonon Correlations

Within the MLF formalism one may also compute the electron-phonon correlation functions in the polaron ground state. This offers a measure of the polaron size as electron and phonons displacements can be taken at different neighbors sites. The on-siteχ0, the first neighbor site χ1 and the second neighbor site χ2 correlation functions are plotted in Figure 3 against the e-ph  coupling both in 1D and 2D. For g~1, there is some residual quasiparticle weight (not appreciable on this scale) associated with χ3 and χ4 which however tend to vanish by increasing the coupling. An adiabatic regime is assumed to point out the self-trapping event which occurs when χ0~1 and χn~0 for n≥1. This happens for g~2 in 1D and g~2.5 in 2D for the lowest case of ω1 here considered. By studying the correlation functions for two values of ω1 it is seen how the intermolecular forces smooth the crossover and extend the polaron size over a larger range of g. The case ω1=60 meV, which allows a comparison with Figure 2, is intermediate between the two plots presented in Figure 3. In such case the self trapping appears at g~2.3 in 1D in fair agreement with the corresponding inflection point in the 1D effective mass (see fifth curve from right in Figure 2(a)). In 2D the polaron size shrinks at slightly larger g than in 1D and the on-site localization occurs at g somewhat larger than the mass inflection point.Electron-phonon correlation functions: on-siteχ0, first neighbor site χ1 and second neighbor site χ2 versus e-ph coupling in (a) 1D and (b) 2D. An adiabatic regime is considered and two values of intermolecular energy are assumed.
(a)(b)Altogether the diamond marked loci displayed in Figure2 indicate that strong mass renormalization is accompanied by on-site (or two sites) polaron localization thus identifying the self trapping line with the formation of small polarons. Before the crossover takes place, for instance, at 1≤g≤2 in Figure 3(a), the 1D polaron spreads over a few (two to four) lattice sites. This is the also the range of couplings in which the polaron band shows the largest deviations [91] from the cosine-like band of standard LF perturbation theory [94] due to the importance of longer (than nearest neighbors) hopping processes. The notion of intermediate polaron seems to us as the most appropriate for such mobile and relatively light objects.


## 5. Path Integral Method

The results obtained so far can be put on sound physical bases by applying the space-time path integral method to the dispersive Holstein Hamiltonian. The method permits to incorporate the effect of the electron-phonon correlations in a momentum dependent effectivee-ph coupling.The phonons operators in (1) can be generally written in terms of the isotropic displacement field un as
(12)bi†+bi=1N∑q2Mω(q)∑nexp(iq·(Ri-Rn))un.Then, thee-ph  term in (1) transforms as follows:
(13)He-ph=gN∑q2Mω(q)∑i,nci†ciunexp(iq·(Ri-Rn)).The sum overn spans all nth neighbors of the Ri lattice site in any dimensionality. It is clear that the phonon dispersion introduces e-ph real space correlations which would be absent in a dispersionless model. Note that  (13) contains the same physics as the e-ph term in (1). Fourier transforming the atomic displacement field and taking the lattice constant |a|=1, from (13), I get for a linear chain and a square lattice, respectively:
(14)Hde-ph=gN3/2∑ici†ci∑q,q′2Mω(q)exp(iq′·Ri)×uq′Sd(q′-q),S1D(q′-q)≡1+2∑n=1n*cos(n(q′-q)),S2D(q′-q)≡1+2∑n=1n*[cos(nΔqx)+cos(nΔqy)]+2∑m,n=1n*[cos(mΔqx+nΔqy)+cos(mΔqx-nΔqy)],Δqx≡qx′-qx,Δqy≡qy′-qy.While, in principle, the sum overn should cover all the N sites in the lattice, the cutoff n* permits to monitor the behavior of the coupling term as a function of the range of the e-ph  correlations. In particular, in 1D the integer n numbers of the neighbors shells up to n* while in the square lattice, the n*=1 term includes the second neighbors shell, the sum up to n*=2 includes the fifth neighbors shells, n*=3 covers the ninth shell, and so on. Switching off the interatomic forces, ω(q)=ω0, one would recover from (14) a local e-ph  coupling model with Sd≡1. As no approximation has been done at this stage (14) is general.


### 5.1. Semiclassical Holstein Model

I apply to the Holstein Hamiltonian space-time mapping techniques [117–119] to write the general path integral for an electron particle in a bath of dispersive phonons. The method has been used to treat also e-ph polymer models [120, 121] in which the electron hopping causes a shift in the atomic displacements and, as a consequence, the vertex function depends on both the electronic and the phononic wave vector. Such shift is however absent in the Holstein model as (1) makes clear [122, 123].I introducex(τ) and y(τ′) as the electron coordinates at the i and j lattice sites, respectively, and the electronic Hamiltonian (first term in (1)) transforms into
(15)He(τ,τ′)=-t(c†(x(τ))c(y(τ′))+c†(y(τ′))c(x(τ))),
whereτ and τ′ are continuous variables (∈[0,β]) in the Matsubara Green's functions formalism with β being the inverse temperature hence the electron hops are not constrained to first neighbors sites. After setting τ′=0,  y(0)≡0, the electron operators are thermally averaged over the ground state of the Hamiltonian. As a result, the average energy per lattice site, he(τ)≡〈He(τ)〉/N, associated to electron hopping reads (in d dimensions):
(16)he(τ)=-t(G[-x(τ),-τ]+G[x(τ),τ]),G[x(τ),τ]=1β∫dkπdexp[ik·x(τ)]∑nexp(-iℏνnτ)iℏνn-ϵk,
whereνn are the fermionic Matsubara frequencies and ϵk is the electron dispersion relation.Consider now thee-ph term. The spatial e-ph correlations contained in (14) are mapped onto the time axis introducing the τ dependence in the displacement field: uq→uq(τ). Assuming periodic atomic particle paths: uq(τ+β)=uq(τ), the displacement is expanded in NF Fourier components:
(17)uq(τ)=uo+∑n=1NF2((ℜun)qcos(ωnτ)-(ℑun)qsin(ωnτ))
with ωn=2nπ/β. I take a semiclassical version of the Holstein Hamiltonian [124], assuming that the phonon coordinates in  (14) as classical variables interacting with quantum mechanical fermion operators. Such approximation may affect the thermodynamics of the system as the quantum lattice fluctuations in fact play a role mainly for intermediate values of the e-ph coupling [89, 125].Averaging (14) on the electronic ground state, the e-ph energy per lattice site is defined as
(18)〈Hde-ph〉N=∑q〈Hde-ph〉q,〈Hde-ph〉q=gN3/22Mω(q)∑q′ρq′uq′Sd(q′-q),ρq′=1N∑i〈ci†ci〉exp(iq′·Ri).Then, on the base of (17) and (18), I identify the perturbing source current [126] for the Holstein model with the τ dependent averaged e-ph  Hamiltonian term:
(19)j(τ)=∑qjq(τ),jq(τ)≡〈Hde-ph〉q.Note that the Holstein source current does not depend on the electron path coordinates. The time dependence is incorporated only in the atomic displacements. This property will allow us to disentangle phonon and electron degrees of freedom in the path integral and in the total partition function.After these premises, one can proceed to write the general path integral for an Holstein electron in a bath of dispersive phonons. Assuming a mixed representation, the electron paths are taken in real space while the phonon paths are in momentum space. Thus, the electron path integral reads(20)〈x(β)∣x(0)〉=∏q〈x(β)∣x(0)〉q,〈x(β)∣x(0)〉q=∫Duq(τ)exp[-∫0βdτM2(u̇q2(τ)+ω2(q)uq2(τ))]×∫Dx(τ)exp[-∫0βdτ(m2ẋ2(τ)+he(τ)-jq(τ))],
where m is the electron mass. The perturbing current in (20) is integrated over τ using (17) and (19). The result is
(21)∫0βdτjq(τ)=βuogd(q),gd(q)=gN22Mω(q)∑q′ρq′Sd(q′-q),
where gd(q) is thus a time averaged e-ph potential.The total partition function can be derived from  (20) by imposing the closure condition both on the phonons (see (17)) and on the electron paths, x(β)=x(0). Using (21), I obtain
(22)ZT=∏q∮Duqexp[βuogd(q)-∫0βdτM2(u̇q2(τ)+ωq2uq2(τ))]×∮Dxexp[-∫0βdτ(m2ẋ2(τ)+he(τ))],
where∮Duq and ∮Dx are the measures of integration which normalize the kinetic terms in the phonon field and electron actions, respectively [126].The phonon degrees of freedom in (22) can be integrated out analytically yielding:
(23)ZT=∏qP(q)×∮Dx(τ)×exp[-∫0βdτ(m2ẋ2(τ)+he(τ))],P(q)=1βω(q)exp[(gd(q)λM)22πω(q)2]×∏n=1NF(2nπ)2(2nπ)2+(βω(q))2
being λM=πℏ2β/M.Equation (23) is the final analytical result from which the thermodynamics of the model can be computed [123]. The exponential function in P(q) embodies the effect of the non local correlations due to the dispersive nature of the phonon spectrum. Phonon and electron contributions to the partition function are decoupled although the effective potential gd(q) carries a dependence on the electron density profile in momentum space through the function ρq.


### 5.2. Electron-Phonon Coupling

The time (temperature) averagede-ph potential in (21) is computed in the case of a linear chain and of a square lattice. As an example, I take the electron density profile as ρq=ρocos(q) in 1D and ρq=ρocos(qx)cos(qy) for the 2D system. Since the momentum integration runs over qi∈[0,π/2],  ρo represents in both cases the total electron density. More structured density functions containing longer range oscillations are not expected to change the trend of the results hereafter presented. Low-energy phonon spectra parameters are here assumed, ω0=20 meV and ω1=10 meV. Setting g=3, I take a strong Holstein coupling although the general trend of the results holds for any value of g.In Figure4(a), g1D(q)/ρo is plotted for three choices of the cutoff n* to emphasize how the potential depends on the e-ph interaction range. The constant value of the potential obtained for ω1=0 is also reported on. For short correlations (n*=4) there is a range of wave vectors in which the effective coupling becomes even larger than in the dispersionless case: this is due to the fact that the approximation related to a short cutoff is not sufficiently accurate. In fact, extending the range of the e-ph correlations, the effective potential is progressively and substantially reduced for all momenta with respect to the dispersionless Holstein model. Numerical convergence for g1D(q) is found at the value n*=24 which corresponds to 48 lattice sites along the chain. Then, the 1D coupling renormalization is q-dependent and generally larger for q~π/2 where the phonon dispersion generates stronger e-ph correlations. Exact diagonalization techniques [91] had found that multiphonon states with longer range hopping processes substantially soften the polaron band narrowing with respect to the predictions of standard SCPT with the largest deviations occurring at the polaron momentum ~π/2. In fact strong, nonlocal e-ph correlations should favour hopping on further than neighbors sites occupied by phonons belonging to the lattice deformation.(a) Time averagede-ph coupling (in units meVÅ-1  versus wave vector for a linear chain. n* represents the cutoff on the e-ph correlations. ρ0 is the electron density. ω0=20 meV and ω1=10 meV.  g=3. The dispersionless e-ph coupling is obtained for ω1=0. (b) Time averaged e-ph coupling in two dimensions versus the y component of the momentum. The cases n*=1,  n*=2, and n*=3 imply that the correlation range includes the second, the fifth and the ninth neighbors shell, respectively. The input parameters are as in (a).
(a)(b)The projections of the two-dimensionale-ph potential along the y component of the wave vector is shown in Figure 4(b) for three values of the cutoff. For n*=1 the correlation range is extended to the second neighbor shell thus including 8 lattice sites. For n*=2 and n*=3 the normalization is over 24 and 48 lattice sites, respectively. Then, the three values of n* in 2D span as many lattice sites as the three values of n* in 1D, respectively. I have made this choice to normalize consistently the potential for the linear chain and the square lattice. There is a strong renormalization for the 2D effective potential with respect to the dispersionless case for any value of the wave vector. Remarkably, the reduction is already evident for n*=1. By extending the correlation range this tendency becomes more pronounced for qy close to the center and to the edge of the reduced Brillouin zone. An analogous behavior is found by projecting the 2D potential along the qx axis. The 2D potential stabilizes by including the 9th neighbor shell (n*=3) in the correlation range.Then, an increased range for thee-ph correlations leads to a very strong reduction of the effective coupling, an effect which is qualitatively analogous to that one would get by hardening the phonon spectrum. This is the physical reason underlying the fact that the Holstein polaron mass can be light.


## 5.1. Semiclassical Holstein Model

I apply to the Holstein Hamiltonian space-time mapping techniques [117–119] to write the general path integral for an electron particle in a bath of dispersive phonons. The method has been used to treat also e-ph polymer models [120, 121] in which the electron hopping causes a shift in the atomic displacements and, as a consequence, the vertex function depends on both the electronic and the phononic wave vector. Such shift is however absent in the Holstein model as (1) makes clear [122, 123].I introducex(τ) and y(τ′) as the electron coordinates at the i and j lattice sites, respectively, and the electronic Hamiltonian (first term in (1)) transforms into
(15)He(τ,τ′)=-t(c†(x(τ))c(y(τ′))+c†(y(τ′))c(x(τ))),
whereτ and τ′ are continuous variables (∈[0,β]) in the Matsubara Green's functions formalism with β being the inverse temperature hence the electron hops are not constrained to first neighbors sites. After setting τ′=0,  y(0)≡0, the electron operators are thermally averaged over the ground state of the Hamiltonian. As a result, the average energy per lattice site, he(τ)≡〈He(τ)〉/N, associated to electron hopping reads (in d dimensions):
(16)he(τ)=-t(G[-x(τ),-τ]+G[x(τ),τ]),G[x(τ),τ]=1β∫dkπdexp[ik·x(τ)]∑nexp(-iℏνnτ)iℏνn-ϵk,
whereνn are the fermionic Matsubara frequencies and ϵk is the electron dispersion relation.Consider now thee-ph term. The spatial e-ph correlations contained in (14) are mapped onto the time axis introducing the τ dependence in the displacement field: uq→uq(τ). Assuming periodic atomic particle paths: uq(τ+β)=uq(τ), the displacement is expanded in NF Fourier components:
(17)uq(τ)=uo+∑n=1NF2((ℜun)qcos(ωnτ)-(ℑun)qsin(ωnτ))
with ωn=2nπ/β. I take a semiclassical version of the Holstein Hamiltonian [124], assuming that the phonon coordinates in  (14) as classical variables interacting with quantum mechanical fermion operators. Such approximation may affect the thermodynamics of the system as the quantum lattice fluctuations in fact play a role mainly for intermediate values of the e-ph coupling [89, 125].Averaging (14) on the electronic ground state, the e-ph energy per lattice site is defined as
(18)〈Hde-ph〉N=∑q〈Hde-ph〉q,〈Hde-ph〉q=gN3/22Mω(q)∑q′ρq′uq′Sd(q′-q),ρq′=1N∑i〈ci†ci〉exp(iq′·Ri).Then, on the base of (17) and (18), I identify the perturbing source current [126] for the Holstein model with the τ dependent averaged e-ph  Hamiltonian term:
(19)j(τ)=∑qjq(τ),jq(τ)≡〈Hde-ph〉q.Note that the Holstein source current does not depend on the electron path coordinates. The time dependence is incorporated only in the atomic displacements. This property will allow us to disentangle phonon and electron degrees of freedom in the path integral and in the total partition function.After these premises, one can proceed to write the general path integral for an Holstein electron in a bath of dispersive phonons. Assuming a mixed representation, the electron paths are taken in real space while the phonon paths are in momentum space. Thus, the electron path integral reads(20)〈x(β)∣x(0)〉=∏q〈x(β)∣x(0)〉q,〈x(β)∣x(0)〉q=∫Duq(τ)exp[-∫0βdτM2(u̇q2(τ)+ω2(q)uq2(τ))]×∫Dx(τ)exp[-∫0βdτ(m2ẋ2(τ)+he(τ)-jq(τ))],
where m is the electron mass. The perturbing current in (20) is integrated over τ using (17) and (19). The result is
(21)∫0βdτjq(τ)=βuogd(q),gd(q)=gN22Mω(q)∑q′ρq′Sd(q′-q),
where gd(q) is thus a time averaged e-ph potential.The total partition function can be derived from  (20) by imposing the closure condition both on the phonons (see (17)) and on the electron paths, x(β)=x(0). Using (21), I obtain
(22)ZT=∏q∮Duqexp[βuogd(q)-∫0βdτM2(u̇q2(τ)+ωq2uq2(τ))]×∮Dxexp[-∫0βdτ(m2ẋ2(τ)+he(τ))],
where∮Duq and ∮Dx are the measures of integration which normalize the kinetic terms in the phonon field and electron actions, respectively [126].The phonon degrees of freedom in (22) can be integrated out analytically yielding:
(23)ZT=∏qP(q)×∮Dx(τ)×exp[-∫0βdτ(m2ẋ2(τ)+he(τ))],P(q)=1βω(q)exp[(gd(q)λM)22πω(q)2]×∏n=1NF(2nπ)2(2nπ)2+(βω(q))2
being λM=πℏ2β/M.Equation (23) is the final analytical result from which the thermodynamics of the model can be computed [123]. The exponential function in P(q) embodies the effect of the non local correlations due to the dispersive nature of the phonon spectrum. Phonon and electron contributions to the partition function are decoupled although the effective potential gd(q) carries a dependence on the electron density profile in momentum space through the function ρq.


## 5.2. Electron-Phonon Coupling

The time (temperature) averagede-ph potential in (21) is computed in the case of a linear chain and of a square lattice. As an example, I take the electron density profile as ρq=ρocos(q) in 1D and ρq=ρocos(qx)cos(qy) for the 2D system. Since the momentum integration runs over qi∈[0,π/2],  ρo represents in both cases the total electron density. More structured density functions containing longer range oscillations are not expected to change the trend of the results hereafter presented. Low-energy phonon spectra parameters are here assumed, ω0=20 meV and ω1=10 meV. Setting g=3, I take a strong Holstein coupling although the general trend of the results holds for any value of g.In Figure4(a), g1D(q)/ρo is plotted for three choices of the cutoff n* to emphasize how the potential depends on the e-ph interaction range. The constant value of the potential obtained for ω1=0 is also reported on. For short correlations (n*=4) there is a range of wave vectors in which the effective coupling becomes even larger than in the dispersionless case: this is due to the fact that the approximation related to a short cutoff is not sufficiently accurate. In fact, extending the range of the e-ph correlations, the effective potential is progressively and substantially reduced for all momenta with respect to the dispersionless Holstein model. Numerical convergence for g1D(q) is found at the value n*=24 which corresponds to 48 lattice sites along the chain. Then, the 1D coupling renormalization is q-dependent and generally larger for q~π/2 where the phonon dispersion generates stronger e-ph correlations. Exact diagonalization techniques [91] had found that multiphonon states with longer range hopping processes substantially soften the polaron band narrowing with respect to the predictions of standard SCPT with the largest deviations occurring at the polaron momentum ~π/2. In fact strong, nonlocal e-ph correlations should favour hopping on further than neighbors sites occupied by phonons belonging to the lattice deformation.(a) Time averagede-ph coupling (in units meVÅ-1  versus wave vector for a linear chain. n* represents the cutoff on the e-ph correlations. ρ0 is the electron density. ω0=20 meV and ω1=10 meV.  g=3. The dispersionless e-ph coupling is obtained for ω1=0. (b) Time averaged e-ph coupling in two dimensions versus the y component of the momentum. The cases n*=1,  n*=2, and n*=3 imply that the correlation range includes the second, the fifth and the ninth neighbors shell, respectively. The input parameters are as in (a).
(a)(b)The projections of the two-dimensionale-ph potential along the y component of the wave vector is shown in Figure 4(b) for three values of the cutoff. For n*=1 the correlation range is extended to the second neighbor shell thus including 8 lattice sites. For n*=2 and n*=3 the normalization is over 24 and 48 lattice sites, respectively. Then, the three values of n* in 2D span as many lattice sites as the three values of n* in 1D, respectively. I have made this choice to normalize consistently the potential for the linear chain and the square lattice. There is a strong renormalization for the 2D effective potential with respect to the dispersionless case for any value of the wave vector. Remarkably, the reduction is already evident for n*=1. By extending the correlation range this tendency becomes more pronounced for qy close to the center and to the edge of the reduced Brillouin zone. An analogous behavior is found by projecting the 2D potential along the qx axis. The 2D potential stabilizes by including the 9th neighbor shell (n*=3) in the correlation range.Then, an increased range for thee-ph correlations leads to a very strong reduction of the effective coupling, an effect which is qualitatively analogous to that one would get by hardening the phonon spectrum. This is the physical reason underlying the fact that the Holstein polaron mass can be light.


## 6. Conclusions

According to the traditional notion of small polaron, the strong electron coupling to the lattice deformation implies a polaron collapse of the electron bandwidth (zt) and remarkable mass renormalization.  Thus, bipolaronic theories have been thought of being inconsistent with superconductivity at high Tc as the latter is inversely proportional to the bipolaron effective mass. Motivated by these issues, I have investigated the conditions for the existence of light polarons in the Holstein model and examined the concept of self-trapping versus dimensionality for a broad range of (anti)adiabatic regimes. The self trapping events correspond to the points of the most rapid increase for the polaron effective mass versus the e-ph coupling. A modified version of the Lang-Firsov transformation accounts for the spreading of the polaron size due to retardation effects which are relevant mainly for moderate e-ph couplings. Assuming large optical phonon frequencies (ω0), it is found that polaron masses display an abrupt although continuous increase driven by the e-ph coupling in any dimension, included the 1D case. Such self trapping events occur in the adiabatic and intermediate (zt/ω0~1) regimes provided that the e-ph coupling is very strong (g/ω0≥2). The second-order of strong coupling perturbation theory permits to decouple two fundamental properties of the polaron landscape, namely, band narrowing and effective mass enhancement. I have considered a range of e-ph couplings (g/ω0≥1) in which the conditions for polaron formation are fulfilled. For a given value of the adiabatic ratio, the onset of the band narrowing anticipates, along the g axis, the crossover to the heavy polaron state which pins the quasiparticle essentially on one lattice site. In such buffer zone, between band narrowing onset and self trapping, the charge carriers can be appropriately defined as intermediate polarons.In all dimensions, there is room for non trapped intermediate polarons in the region of the moderate couplings. Such polarons spreads over a few lattice sites and their real space extension grows with dimensionality. Thus, analytical methods can suitably describe some polaron properties also in that interesting intermediate window ofe-ph couplings (i.e., 1≤g/ω0≤2 for the 1D system) in which adiabatic polarons are not large but not yet on-site localized. In the square lattice, at the onset of the self trapping state, the polaron over bare electron mass ratios are ~5–10 and even smaller for lower couplings.The physical origins for the possibility of light polarons have been analysed treating the dispersive Holstein Hamiltonian model by the imaginary time path integral method. The perturbing source current is the time averagede-ph Hamiltonian term. The partition function of the model has been derived integrating out the phonon degrees of freedom in semiclassical approximation. Focusing on the role of the intermolecular forces, I have quantitatively determined the effects of the dispersive spectrum on the time-averaged e-ph coupling which incorporates the nonlocality effects. It is found that, in the square lattice, the e-ph coupling is strongly renormalized downwards and this phenomenon is already pronounced by the (minimal) inclusion of the second neighbors intermolecular shell. This may explain why 2D Holstein polarons can be light while the real space range of the e-ph correlations is relatively short.



---

*Source: 815917-2009-10-15.xml*