Barrier crossing to the small Holstein polaron regime
aa r X i v : . [ c ond - m a t . o t h e r] A ug Barrier crossing to the small Holstein polaron regime
Peter Hamm ∗ , G. P. Tsironis + ∗ Physikalisch-Chemisches Institut, Universit¨at Z¨urich, Winterthurerstr. 190, CH-8057 Z¨urich, Switzerland + Department of Physics, University of Crete and Institute of Electronic Structure and Laser,FORTH, P.O. Box 2208, Heraklion 71003, Crete, Greece. (Dated: October 25, 2018)We investigate the dimensionality effects of the Holstein polaron from the fully quantum regime, where thecrossover between large and small polaron solutions is known to be continuous in all dimensions, into thelimit described by the semiclassical Discrete Nonlinear Schr¨odinger (DNLS) Equation, where the crossoveris continuous in 1D but discontinuous in higher dimensions. We use exact numerics on one hand and a twovariable parametrization of the Toyozawa ansatz on the other in order to probe the crossover region in allparameter regimes. We find that a barrier appears also in 1D separating the two types of solutions, seeminglyin contradiction to the common paradigm for the DNLS according to which the crossover is barrier-free. Wequantify the polaron behavior in the crossover region as a function of the exciton overlap and find that the barrierremains small in 1D and tunnelling through it is not rate-limiting.
The Holstein model has been used widely to describe trans-port properties of electrons or excitons in diverse systems,ranging from molecular crystals [1] to strongly correlatedelectron-phonon systems [2], interface charge localization inalkane layers [3], organic transistors [4, 5], proteins andprotein-like crystals [6, 7, 8] as well as DNA [9]. In thecontext of the model, the electronic degrees of freedom arecoupled to dispersionless phonons leading to the formation oflarge, i.e. much larger than the lattice spacing, or small, viz.essentially single-site, polarons. The relative magnitudes ofthree parameters of the theory, viz. electronic overlap inte-gral, phonon energy and exciton-phonon coupling determinethe specific polaron properties. Outstanding issues have beenthe nature of polarons in different regimes, their transitions aswell as the precise onset of selftrapping. The Discrete Non-linear Schr¨odinger (DNLS) equation, that is a semiclassicalapproximation to the Holstein model, gives a continuous tran-sition in 1D and a discontinuous transition in higher dimen-sions [10, 11], while most variational calculations, in partic-ular those based on the Toyozawa wavefunction, show a dis-continuity in all dimensions [12, 13, 14, 15, 16, 17, 18]. Onthe other hand, it has been shown recently with the help of anumerically exact solution of the full-quantum Holstein po-laron that the cross-over is continuous in 1D [19] as well asin higher dimensions [20, 21]. The non-existence of phasetransitions in polaron systems can also be proven on very gen-eral grounds [22]. In the present paper, we address the natureof the polaron transitions using exact numerics for the fullyquantum Holstein model on one hand, as well as a simple,physically motivated two parameter variational minimizationon the other hand. The comparison of the two approachesresolves the issue of the large to small polaron transition, pro-duces quantitative information on the nature of the selftrap-ping transition and finalizes pending issues regarding deficien-cies of standard approximations.The Holstein Hamiltonian reads: H = H ex + H ph + H ex , ph H ex = − J (cid:229) j (cid:16) B j † B j + + B j † B j − (cid:17) H ph = ¯ h w (cid:229) j (cid:16) b j † b j + / (cid:17) H ex , ph = − c (cid:229) j B j † B j (cid:16) b j † + b j (cid:17) (1)where b † j ( b j ) and B † j ( B j ) destroy (create) a phonon and anelectron (or exciton) at site j , respectively, and J , ¯ h w , c are theexcitonic overlap, phonon energy and exciton-phonon cou-pling, respectively. In what follows we use two dimension-less parameters, i.e. the exciton coupling J / ¯ h w and exciton-phonon coupling c / ¯ h w while all energies are given in unitsof one phonon-quantum ¯ h w . We also restrict our study to thewavenumber k = k = q ≡ ( b †0 + b ) , with phonon coordinates at sites different fromthe exciton site traced out. We note that the results are qualita-tively similar in all dimensions: The reduced phonon densityis centered at q ≈ c for small exciton couplings J (i.e. a smallpolaron), similar to the J = q = c [6]. For sufficientlylarge exciton couplings J , the phonon displacement shifts to q ≈
0, i.e. essentially a free exciton without any phonon-displacement. We observe that the transition between the tworegimes is smoother in 1D than it is in 2D and 3D, however,the essential point is that free exciton and small polaron solu-tions coexist in a certain parameter range, also in 1D. The twostates gradually change their relative weights, but hardly theircharacter, as the exciton coupling J is increased. If we inter-pret the reduced phonon density as one that originates froman effective potential V ( q ) , then the coexistence of two solu-tions hints to the presence of a barrier separating them. Thisconclusion is in disagreement with the semiclassical DNLScase [10, 11], which does not reveal any barrier in 1D butrather a gradual transition.The phase diagrams in Fig. 1, bottom row, summarize theseresults where we plot the expectation value of the mean dis- - - - - - - q q q
1D 2D 3D c J JJ c c a b cd fe
J J J
FIG. 1: (Color online) Results from a numerically exact diagonalization of the Holstein Hamiltonian in 1D (left), 2D (middle) and 3D (right).Top row (a-c): Phonon density of the polaron ground state as a function of the phonon-coordinate q . Shown is the reduced phonon density atthe site of the exciton. The exciton-phonon coupling was set to c = .
0. Bottom row (d-f): Phase diagram, plotting the expectation value ofthe phonon coordinate h q i in units of c for the polaron ground state. Contour lines are equidistant with spacing 0.1. In red (light gray) are thephase separation lines from Eq. 2, whereas in blue (dark gray) are those obtained from equalling the two variational minima of the Toyozawawavefunction (see Fig. 2). placement of the phonon coordinate h q i in units of c (whichis the maximum value obtained in the J = c and J , the result is qualitatively very similar in alldimensions with a continuous transition between small andlarge polaron solution. Nevertheless, it is quite evident thatthe transition becomes increasingly more abrupt in 2D and3D for large couplings (i.e. as we approach the regime of theDNLS), whereas the cross-over stays smooth in the 1D case.The phase separation lines can be fitted extremely well to ageneric relationship (Fig. 1, bottom row, red lines): c c = / a + √ a J (2)with a = a = .
34 in 2D and a = .
41 in 3D, re-spectively. The 1D value stems from Lindenberg et al. em-pirical relationship [23], whereas those for 2D and 3D co-incide with the critical couplings in the DNLS limit (when √ a J ≫ / a ) [11]. Hence, Eq. (2) does give the correct semi-classical limits expected from the DNLS analysis.Although exact, the results of Fig. 1 originate from a di-agonalization of a huge matrix and provide relatively lit-tle physical insight. We therefore employ an approximate,yet more intuitive trial function, the Toyozawa wavefunc-tion [12, 13, 14, 17], and use the numerically exact solutionas an accuracy check. The Toyozawa trial function is writtenas a Bloch state: | Y i = √ N (cid:229) l e ikl | y l i (3) with | y i ≡ (cid:229) i a i B † i | i ex (cid:213) j | q j − i i j (4)where all | y l i ≡ | y i are identical, | i ex is an exciton vacuumstate and | q i j a coherent phonon state: | q i j ≡ e q ( b † j − b j ) | i ph . (5)The Toyozawa wavefunction function can be considered aBloch-extension of symmetry breaking soliton solutions [14],where both excitons and phonons are dressing one commontrapping site l . In 1D, the norm of the wavefunction | y i is [13]: h y | y i = (cid:229) i " (cid:229) j a j a j − i (cid:213) k e − ( q k − q k − i ) / (6)and the expectation value of the Hamiltonian: h y | H | y i = − J (cid:229) i " (cid:229) j a j a j − i (cid:213) k e − ( q k − q k − i + ) / (7) + (cid:229) i " (cid:229) j a j a j − i (cid:229) l q l q l − i (cid:213) k e − ( q k − q k − i ) / − c (cid:229) i " (cid:229) j a j a j − i ( q j + q j − i ) (cid:213) k e − ( q k − q k − i ) / J =9 c =4.011D q q q h ex ba J =4.5 c =4.102D J =3 c =4.173D c FIG. 2: Variational energy of the Toyozawa wavefunction as a func-tion of q and h ex in 1D (left), 2D (middle) and 3D (right). Exci-ton coupling is J = / D and the exciton-phonon coupling c = c varc has been calculated so that both variational minima are equally deep.Contour line spacings are 0 .
05 in (a) and 0 . In more than 1D (dimension D ), the indices i and j are re-placed by { i x , i y ( , i z ) } and { j x , j y ( , j z ) } , respectively, in eithersum/product. Furthermore, the exciton term is repeated D times with the one-site shift (i.e. the q j − i + -term) in eitherdirection (for isotropic exciton coupling, as considered here,this can be reduced to one exciton term with prefactor 2 DJ ).Zhao et al. [14] have minimized the energy of the Toyozawawavefunction by effectively varying all parameters { a i } and { q i } (the calculation was done in momentum space), still leav-ing us with a multi-parameter solution. In contrast, we makean exponential ansatz, in analogy to Refs. [11, 13, 17]: a i ≡ h | i | ex q i ≡ q h | i | ph (8)that contains only three parameters. Minimizing Eq. 7 re-veals (cid:229) i q i = c as one condition [13, 14] which allows us toeliminate one of these parameters leading to h ph = ( D √ c − D √ q ) / ( D √ c + D √ q ) . Hence, we can express the variationalenergy as a function of effectively two parameters, q and h ex .Figure 2 shows the variational energy of the trial functionfor J = / D . We observe the emergence of two variationalminima in all dimensions with a barrier separating them. Thebarrier is very shallow in 1D ( ≈ .
1) while it is much morepronounced in 2D or 3D. This is the barrier which is respon-sible for the double-peak structure of the wavefunction in thenumerically exact solution, which also exists in 1D (Fig. 1a).Venzl and Fischer [13] have shown a very similar variationalenergy surface in 1D (as a function of h ex and h ph , eliminat-ing q ), however, did not address the double-minimum due toa choice of too small value of te exciton coupling, viz. J = c varc that induces the large-to-small polarontransition, we find that the numerically obtained line c varc ( J ) follows very closely the relationship of Eq. 2 (compare blueand red lines in Fig. 1d-f). This near identity of c varc and Eq. 2starts to deviate in 1D for very large exciton overlaps J > ∼ E c J / Ñ w=9 e tunnel Dc FIG. 3: (Color online) Energies of the lowest energy state for J = c varied through the cross-over re-gion in 1D. Black line: numerically exact result. Red (light gray)lines: The two minima of the variational energy. Blue (dark gray)line: Energy corrected by tunnelling between the two components.The inset focuses into the crossing region and defines the terms in-troduced in Eq. 10. and compares it with the two variational minima (where theyexist, red lines). Outside the crossing region, the variationalresult agrees reasonably well with the exact energy – hencethe two-parameter trial function Eq. 8 is quite good in thisregime. However, both deviate significantly at the crossingpoint, indicating that the Toyozawa wavefunction per se doesnot capture the physics of the crossover correctly [12, 14, 15].In particular, the variational minimum, as well as all otherproperties of the trial wavefunction, would not be analyticalat the cross-over point. This is simply due to the fact that thetwo variational minima are separated by a barrier (Fig. 2a),and the solution of an overall minimization switches abruptlywhen one of these local minima decreases below the other.Bariˇsi´c has extended the variational space by consideringa linear combination of two (or more) Toyozawa wavefunc-tions [17]: | y i ≡ c l | y ( l ) i + c s | y ( s ) i , (9)showing that this removes the discontinuity in 1D. Fig. 2amakes very clear why that is so: The two Toyozawa wave-functions will sit in either of two minima of the variationalenergy, and tunneling coupling between them will (a) lowerthe energy and thereby correct for approximately half of themismatch to the numerically exact result, and (b) render thecross-over continuous (see Fig. 3, blue line) [24]In order to understand the nature of the barrier separatingthe small and large polaron solutions we address the questionof how the crossover domain behaves in the limit J → ¥ . Mix-ing of the large and the small polaron states will occur when-ever the tunnel coupling is on the order of the energy separa-tion between them, e tunnel > ∼ | E ( l ) − E ( s ) | . Hence, the quotientof the tunnel coupling divided by the difference of the deriva-tives of the two variational energies with respect to c definesa measure of the extension Dc of the crossover region (Fig. 3,inset, shows pictorially e tunnel and Dc ): Dc ≡ e tunnel d (cid:0) E ( l ) − E ( s ) (cid:1) / d c (10) -3 -2 -1 -3 -2 -1 b
2D 1D e t unn e l a
3D 2D DJ Dc FIG. 4: (Color online) (a) Tunnel coupling e tunnel and (b) extensionof the cross-over region Dc as a function of exciton coupling J (with c = c varc ) in 1D (red), 2D (blue) and 3D (green). Figure 4a, red line, shows that the tunnel coupling e tunnel staysessentially constant in 1D for all J ’s, however, the denomina-tor in Eq. 10 decreases with J so that Dc increases roughly lin-early for large enough J (Fig. 4b, red line). Consequently, thecross-over remains continuous even for J → ¥ , leading to thecorrect 1D result in the adiabatic limit. In contrast, in higherdimensions, both the tunnel coupling e tunnel as well as Dc de-cay exponentially with J (Fig. 4, blue and green lines). Hence,the superposition of two Toyozawa wavefunctions, althoughpredicting in principle a continuous crossover, leads (in 2D,3D) to an increasingly more abrupt transition for J → ¥ , justlike the numerically exact result does (Fig. 1e,f).In conclusion, we studied the dependence of quantum po-laron regimes on lattice dimensionality employing exact nu-merics and a physically motivated variational ansatz based onthe Toyozawa wavefunction. The exact numerical solutionshows that the large-to-small polaron crossover is continuous in all dimensions although much more abrupt in two and threedimensions, especially as the parameters approach the adia-batic regime (Fig. 1e,f). In 1D, the transition is smooth andremains continuous also for large J ’s (Fig. 1d).More physical insight into the polaron features is obtainedfrom an approximate trial wavefunction, viz. the Toyozawawavefunction of Eq. 4, that leads to variational results whoseregime of validity has been analyzed extensively [14, 15].When the latter is augmented with an exponential ansatz ofEq. 8 containing effectively only two variational parameters,it leads to a very intuitive and clear picture. Despite the sim-plicity of this approach, it provides results that coincide withthe exact ones outside the crossover region, while, neverthe-less, predicting non-analytic properties at the transition pointin all dimensions. The appearence of this discontinuity waspointed out previously [17, 20] but it has been found in thiswork to be due to the presence of two, almost degenerate solu-tions coexisting in a certain parameter regime while separatedthrough a barrier. Consequently, a superposition of two Toy-ozawa wavefunctions, as in Eq.9, reveals qualitatively correctresults for the cross-over region in all dimensions and in allparameter regimes [17]. The resulting intuitive picture maybe used to asses quantitatively the properties of the cross-overregion as a function of the exciton coupling. One importantoutcome of this analysis is that a barrier separating large andsmall polaron solutions exists also in 1D, in sharp contrastto the predictions of semiclassical theories. We neverthelessfind that tunnelling through the barrier remains efficient in allparameter regimes (Fig. 4a, red line), since the barrier staysof the order of ¯ h w and hence does not hamper motion in theDNLS limit where phonons are treated classically. The Hol-stein model is a fundamental ”minimal” model for stronglyinteracting systems and thus our results may also have ram-ifications for more complex models such as for instance theHolstein-Hubbard model. [1] T. Holstein, Ann. Phys. , 325 (1959).[2] W. Z. Wang, A. R. Bishop, J. T. Gammel, and R. N. Silver,Phys. Rev. Lett. , 3284 (1998).[3] N. H. Ge, C. M. Wong, R. L. Lingle Jr., J. D. McNeill, K. J.Gaffney, and C. B. Harris, Science , 202 (1998).[4] L. Torsi, A. Dodabalapur, L. J. Rothberg, A. W. P. Fung, andH. E. Katz, Science , 1462 (1996).[5] W. A. Schoonveld, J. Wildeman, D. Fichou, P. A. Bobbert, B. J.van Wees, and T. M. Klapwijk, Nature , 977 (2000).[6] A. C. Scott, Phys. Reports , 1 (1992).[7] J. Edler, P. Hamm, and A. C. Scott, Phys. Rev. Lett. , 067403(2002).[8] P. Hamm and G. P. Tsironis, Eur. Phys. J. Sp. Topics , 303(2007).[9] S. S. Alexandre, E. Artacho, J. M. Soler, and H. Chachan, Phys.Rev. Lett. , 108105 (2003).[10] V. V. Kabanov and O. Y. Mashtakov, Phys. Rev. B , 6060(1993).[11] G. Kalosakas, S. Aubry, and G. P. Tsironis, Phys. Rev. B ,3094 (1998).[12] Y. Toyozawa, Prog. Theor. Phys. , 29 (1961). [13] G. Venzl and S. Fischer, Phys. Rev. B , 6437 (1985).[14] Y. Zhao, D. W. Brown, and K. Lindenberg, J. Chem. Phys. ,3159 (1997).[15] A. H. Romero, D. W. Brown, and K. Lindenberg, J. Chem.Phys. , 6540 (1998).[16] V. Cataudella, G. de Filippis, and G. Iadonisi, Phys. Rev. B ,15163 (1999).[17] O. S. Bariˇsi´c, Phys. Rev. B , 144301 (2002).[18] O. S. Bariˇsi´c, Europhys. Lett. , 57004 (2007).[19] J. Bonˇca, S. A. Trugman, and I. Batisti´c, Phys. Rev. B , 1633(1999).[20] L.-C. Ku, S. A. Trugman, and J. Bonˇca, Phys. Rev. B ,174306 (2002).[21] L.-C. Ku and S. A. Trugman, Phys. Rev. B , 014307 (2007).[22] B. Gerlach and H. L¨owen, Phys. Rev. B , 4291 (1987).[23] A. H. Romero, D. W. Brown, and K. Lindenberg, Phys. Rev. B , 13728 (1999).[24] The expansion coefficients c l and c s in Eq. 9 are determined byminimization through direct matrix diagonalization. The over-lap and mixing matrix elements are calculated in analogy toEqs. 6 and 7, respectively, with one coefficient a j or q k in either term stemming from | y ( l ) i , and the other, a j − i or q k − i , from | y ( s ) ii