A generalised Davydov-Scott model for polarons in linear peptide chains
AA generalised Davydov-Scott model for polarons in linear peptide chains
April 2017
J. Luo * B. M. A. G. Piette † *[email protected] † [email protected] of Mathematical Sciences, Durham University a r X i v : . [ c ond - m a t . o t h e r] A p r bstract We present a one-parameter family of mathematical models describing the dynamics of polaronsin periodic structures, such as linear polypeptides, which, by tuning the model parameter,can be reduced to the Davydov or the Scott model. We describe the physical significanceof this parameter and, in the continuum limit, we derive analytical solutions which representstationary polarons. On a discrete lattice, we compute stationary polaron solutions numerically.We investigate polaron propagation induced by several external forcing mechanisms. We showthat an electric field consisting of a constant and a periodic component can induce polaronmotion with minimal energy loss. We also show that thermal fluctuations can facilitate theonset of polaron motion. Finally, we discuss the bio-physical implications of our results.PACS number(s): 71.38.-k, 63.20.kd, 05.60.-k, 05.40.-a
Introduction
The polaron, a quasi-particle formed by the coupling of an electron to a vibrating lattice,was first theorised by L. D. Landau in 1933 [1]. In essence, polaron formation is a process ofelectron self-trapping. The presence of the electron causes localised distortions in the naturalvibrational mode of the lattice, a.k.a. the lattice phonon. In return, if the electromagneticinteraction between the electron and lattice is appropriate, then the phonon distortions canlower the potential well for the electron, thus trapping the electron.Some twenty years after the inception of the polaron concept, a mathematical description ofit was formalised by H. Fr¨ohlich [2] and subsequently T. Holstein [3, 4]. Since then, propertiesof the Fr¨ohlich-Holstein polaron have been well studied, with some authors hypothesising anapplication of dynamical polarons as electron transporters in conductive material [5–7]. In the1970s, A. S. Davydov used the basis of polaron theory to explain some biological processes [8].Specifically he proposed that, in an α -helical protein, a certain intramolecular oscillator caninteract with the peptide chain in a way similar to an electron interacting with a crystallattice. Davydov suggested that this interaction could lead to the localisation and propagationof vibrational energy in the α -helix. Later, A. C. Scott modified Davydov’s theory, takinginto account the internal geometry of peptide units [9]. Some authors argued that, given thepolarisability of peptide units, electron self-trapping is also possible in proteins, and it can bedescribed by the same mathematical model that Davydov and Scott used [10, 11]. It shouldtherefore be possible that polaronic transport of electrons may take place in proteins, too.Recently, L. S. Brizhik et al. reported on the properties of static and dynamical polarons insimple molecular chains, and adverted to the applicability of their results to electron transportin biomolecules such as proteins [12–14]. Their studies were based on the Davydov-Scott model.In the current study, we propose a generalisation to the Davydov-Scott model, and use itto explore the properties of polarons in a linear peptide chain. In the generalised model,there is an extra parameter which represents the extent to which the electron-polypeptideinteraction is spatially symmetric. In section 2, we describe our model and explain why theextra parameter is necessary. We also give physical interpretations of all other parameters inthe model, justifying the choices of their values where possible. Then, we derive a set of coupleddynamical equations which govern the electron and phonon parts of the polaron, as well ashow they interact. In section 3 we look at solutions to our equations which are stationary,and thus deduce properties of static polarons admissible by our model, such as the polaron’sbinding energy. The process of solving the equations is carried out analytically as well asnumerically. By the former approach, a closed-form expression for the solution is found, butits use is limited, because the solution process involves a few approximations and simplifyingassumptions. By the numerical approach, no convenient expression for the solution is possible,but the method solves the equations directly without simplifications. We compare the resultsproduced by the two different methods.Section 4 concerns dynamical polarons. We discover that it is possible to use a suitable externalforcing to displace the stationary polaron, and to sustain its motion in such a way that itsenergy remains highly stable. We investigate how the polaron’s motion depends upon ourforcing parameters. We use only numerical methods to obtain our results in section 4, as wellas those in section 5, where we consider how the polaron’s motion is affected by temperatureof the environment. For this part, the external forcing from section 4 remains in place, but wealso utilise a parameter which controls the magnitude of the thermal effect. To account for therandom nature of thermal fluctuations, we repeat each numerical simulation many times over,taking the average of the results. Finally, we conclude by discussing the physical realisabililty1f our mathematical model, particularly how the external forcing which we study in section 4may be realised. We also briefly discuss the generalisability of our model to studying electrontransport by polarons in α -helices. In both Davydov’s and Scott’s models, the Hamiltonian for a system of excitons interacting withone-dimensional lattice phonons is written in Fr¨ohlich-Holstein form ˆ H = ˆ H e + ˆ H p + ˆ H int , whereˆ H e , ˆ H p and ˆ H int represent energy contributions from the exciton, phonon and interaction parts,respectively [2, 3, 9, 15, 16]. We adopt this Hamiltonian for our model, and following [12–14] weconsider an additional external Hamiltonian , ˆ H ext , so that our Hamiltonian takes the formˆ H = ˆ H e + ˆ H p + ˆ H int + ˆ H ext , (1)where ˆ H e describes a tight-binding electron, the stretching and compressing of hydrogen bondsin the peptide chain are phonon oscillations described by ˆ H p , ˆ H int accounts for the electron-phonon interaction, and ˆ H ext represents the effect of an external electric field. We assume thatthe peptide chain consists of N + 1 identical units and N identical hydrogen bonds. In thetight-binding approximation, we haveˆ H e = N (cid:88) n =0 J ˆ A † n ˆ A n − N − (cid:88) n =0 J (cid:16) ˆ A † n +1 ˆ A n + ˆ A † n ˆ A n +1 (cid:17) . (2)The subscript n in eq. (2) labels peptide units, which are the unit cells of our lattice. ˆ A † n and ˆ A n are local electron creation and annihilation operators, respectively. J is the potential energyof a localised electron. Modelling each unit as a point-dipole, we assume the nearest-neighbourdipole interaction energy is a constant and write it as − J [17–19]. The external Hamiltonian,ˆ H ext = − N (cid:88) n =0 qE ( t ) R ( n − n ) ˆ A † n ˆ A n , (3)models the effect of an electric field with strength E ( t ) on the potential energy of a localisedelectron with charge − q . The potential energy due to E ( t ) is set to zero at some arbitrary n , and R = 4 . H p , is a classical one. In the harmonic approximation, the hydrogen bonds aremodelled as Hookean springs with force constant K , and therefore ˆ H p takes the formˆ H p = N (cid:88) n =0 P n M + N − (cid:88) n =0 M Ω ( U n +1 − U n ) , (4)where M = 1 . × − kg is the average mass of a peptide unit in a membrane α -helix [20],and we have defined Ω := (cid:112) K/M . U n and P n are, respectively, the displacement and conjugatemomentum of the n th unit. Thus, the first and second sums in the expression for ˆ H p represent,respectively, the kinetic and potential energies of the lattice. We take the value of Ω to be thenatural angular frequency of slow phonons in an α -helix, Ω = 5 . × s − [21–23]. To derivethe interaction Hamiltonian, ˆ H int , Davydov and Scott assumed that the energy of an on-siteexcitation depends on lattice deformations in its vicinity. For us, the local deformations are S n := U n +1 − U n , (5)2amely the amount by which the lengths of hydrogen bonds deviate from equilibrium. ByDavydov and Scott’s assumption, if we write the electron energy at site n in a Taylor expansion,the first two terms are J + χG n ( S n , S n − ), where χ is a constant, G n is a linear function, and (cid:12)(cid:12) χG n /J (cid:12)(cid:12) (cid:28)
1. Then the interaction Hamiltonian is ˆ H int = (cid:80) Nn =0 χG n ˆ A † n ˆ A n . Davydov assumedthat S n and S n − have equal influence on local excitation energies [15], so G n = ( S n + S n − ) / H int isˆ H Davint = ( χ/ U − U ) ˆ A † ˆ A + N − (cid:88) n =1 ( U n +1 − U n − ) ˆ A † n ˆ A n + ( U N − U N − ) ˆ A † N ˆ A N ] . (6)Davydov’s model is therefore spatially symmetric, since ˆ A † n ˆ A n is coupled equally to U n +1 andto U n − . Scott modified Davydov’s model by opting for the antisymmetric G n ( S n , S n − ) = S n instead [9]. The reason is, while both authors let an intra-peptide C=O oscillator take therole of the exciton, Davydov neglected the internal geometry of the peptide units, but Scottpointed out that every unit has its C=O pair immediately adjacent to the next hydrogen bondin the chain. This leads Scott to assume that ˆ A † n ˆ A n is coupled to, without loss of generality, S n , and not S n − . Scott therefore hadˆ H Scoint = N − (cid:88) n =0 χ ( U n +1 − U n ) ˆ A † n ˆ A n . (7)Since we are modelling an electron as opposed to an intra-peptide oscillator, we cannot useScott’s argument to justify assuming that ˆ A † n ˆ A n is coupled to U n +1 and not to U n − . Norshould we assume that on-site energies are affected equally by deformations on both sides,as Davydov did. We therefore propose G ( S n , S n − ) = χ r S n + χ l S n − , taking without loss ofgenerality χ r > ≤ χ l ≤ χ r . Then,ˆ H int = χ r ( U − U ) ˆ A † ˆ A + χ l ( U N − U N − ) ˆ A † N ˆ A N + N − (cid:88) n =1 [ χ r ( U n +1 − U n ) + χ l ( U n − U n − )] ˆ A † n ˆ A n . (8)By defining χ := χ r + χ l , β = χ r − χ l χ r + χ l ∈ [0 , , (9)we can write ˆ H int = χ β )( U − U ) ˆ A † ˆ A + χ − β )( U N − U N − ) ˆ A † N ˆ A N + N − (cid:88) n =1 χ (cid:2) ( U n +1 − U n − ) + β ( U n +1 + U n − − U n ) (cid:3) ˆ A † n ˆ A n . (10)We treat χ as an adjustable parameter. Setting β = 0 ( χ l = χ r ) gives us the symmetric modelof Davydov as per eq. (6), whilst setting β = 1 ( χ l = 0) produces the antisymmetric model ofScott as per eq. (7). The larger β is, the less spatial symmetry our model possesses. Indeed,for β ∈ [0 , n -( n + 1) coupling strength to n -( n −
1) coupling strength is givenby χ r /χ l = (1 + β ) / (1 − β ), and this ratio is strictly increasing with β .We write the electronic state of the system as a linear superposition of local excitations [12], | Ψ( t ) (cid:105) = N (cid:88) n =0 α n ( t ) ˆ A † n | vac (cid:105) , (11)3here | vac (cid:105) is the vacuum state, and α n ∈ C is the probability amplitude for an electronlocalised at the n th site, subject to the normalisation condition, N (cid:88) n =0 | α n | = 1 . (12)We proceed to derive dynamical equations for α n and U n . By equating coefficients of ˆ A † n | vac (cid:105) on both sides of the Schr¨odinger equation, i (cid:126) d | Ψ (cid:105) / d t = ( ˆ H e + ˆ H int + ˆ H ext ) | Ψ (cid:105) , we obtain i (cid:126) d α n d t = (cid:20) J + χ S n + S n − ) + χ β ( S n − S n − ) (cid:21) α n − J ( α n +1 + α n − ) − eE ( t ) R ( n − n ) α n . (13)We have defined S − = S N = 0 , α − = α N +1 = 0 , (14)so that eq. (13) holds for all n including the boundary terms ( n = 0 , N ). Equations for U n arederived from classical Hamilton equations, d U n / d t = ∂H cla /∂P n and d P n / d t = − ∂H cla /∂U n ,where H cla := (cid:104) Ψ | ( ˆ H p + ˆ H int ) | Ψ (cid:105) . These equations are M d U n d t = ( S n − S n − ) + χ (cid:20)(cid:16) | α n +1 | + | α n | (cid:17) − (cid:16) | α n | + | α n − | (cid:17)(cid:21) − χ β (cid:20)(cid:16) | α n +1 | −| α n | (cid:17) − (cid:16) | α n | −| α n − | (cid:17)(cid:21) . (15)In order that eq. (15) holds at the boundaries, we have set α = α N = 0 . (16)This boundary condition is justified because we expect the probability distribution | α n | to behighly localised with half-width of O (1), and because we will be working with large latticeswith N (cid:29)
1. We also impose the following boundary condition on U n , representing a peptidechain which is fixed at one end. U = d U d t = 0 . (17)Next, we introduce the gauge transformation, α n ( t ) = ψ n ( t ) exp (cid:20) − it (cid:126) ( J − J ) (cid:21) , (18)which sets J = 2 in eq. (13). Physically this represents a shift in the arbitrary referencevalue from which energy is measured. Combining the 2 α n term with the J term in eq. (13),we obtain the discrete Laplacian, − J ( α n +1 + α n − − α n ). Meanwhile, to account for theinteraction between the peptide chain and its environment, we need to add Langevin terms tothe r.h.s. of eq. (15) [13, 14, 24, 25]. They are, a damping term describing energy dissipationdue to friction, − Γ d U n / d t , where Γ is the viscous damping coefficient; and a stochastic term F n ( t ), describing random forces due to thermal fluctuations. Specifically, F n ( t ) is normally-distributed with zero mean and correlation function (cid:104) F m ( t ) F n ( t (cid:48) ) (cid:105) = 2Γ k B Θ δ m,n δ ( t − t (cid:48) ), where k B is the Boltzmann constant and Θ is the temperature of the environment.4caling time by Ω − and length by R gives us the following non-dimensionalised dynamicalequations for ψ n and u n := U n /R , for n = 0 , , . . . , N . i ˙ ψ n = σ (cid:2) ( s n + s n − ) + β ( s n − s n − ) (cid:3) ψ n − ρ ( ψ n +1 + ψ n − − ψ n ) − (cid:15) ( τ )( n − n ) ψ n , (19a)¨ u n = ( s n − s n − ) + δ (cid:2) ( c n − c n − ) − β ( g n − g n − ) (cid:3) − γ ˙ u n + f n ( τ ) , (19b)where we have defined s n := u n +1 − u n , g n := | ψ n +1 | −| ψ n | , c n := | ψ n +1 | + | ψ n | , (20)and where the overdot denotes differentiation with respect to dimensionless time τ , and ρ = J (cid:126) Ω , σ = Rχ (cid:126) Ω , (cid:15) = qER (cid:126) Ω , δ = χ M R Ω , γ = Γ M Ω , f n = F n M R Ω . (21)Equation (19) holds subject to the boundary conditions (14), (16) and (17), as well as thenormalisation condition, N (cid:88) n =0 | ψ n | = 1 . (22)It is easily verifiable that setting β = 0 and β = 1 in eq. (19) produces Davydov’s andScott’s dynamical equations, respectively [9, 15]. ρ is known as an adiabaticity parameter ,as it is the ratio of the characteristic time scale of phonon vibrations to that of electronicphase variations [26]. We fix J , following [12–14] (which used a different scaling), at ρ = 2 . M, R and Ω are fixed, the ratio σ/δ = 1880 is constant. The range of δ whichwe consider throughout this study correspond to χ ∼ O (10 − ) Newtons, agreeing with [9].Finally, we take γ = 0 .
05, agreeing with [12–14] up to different scaling factors.
We derive stationary polaron solutions to eq. (19), subject to zero electric field ( (cid:15) = 0) and zerotemperature ( f n = 0). We consider analytical and numerical methods separately, and comparethe results. When f n = 0 and ˙ u n = ¨ u n = 0, eq. (19b) becomes s n − s n − = δ (cid:2) β ( g n − g n − ) − ( c n − c n − ) (cid:3) , (23)which holds if s n ≡ u n +1 − u n = δ ( βg n − c n ) = δ (cid:104) ( β − | ψ n +1 | − ( β + 1) | ψ n | (cid:105) . (24)Putting eq. (24) into eq. (19a) and requiring (cid:15) = 0 gives us i ˙ ψ n = − σδ (cid:104)(cid:0) − β (cid:1) | ψ n +1 | + (cid:0) − β (cid:1) | ψ n − | + 2 (cid:0) β (cid:1) | ψ n | (cid:105) ψ n − ρ ( ψ n +1 + ψ n − − ψ n ) . (25)5efining ∆ ψ n := ψ n +1 + ψ n − − ψ n , ∆ | ψ n | := | ψ n +1 | + | ψ n − | − | ψ n | , (26) λ := 4 σδρ ≡ χ M Ω J , η := σδρ (cid:0) − β (cid:1) ≡ λ (cid:0) − β (cid:1) , (27)we can rewrite eq. (25) as iρ − ˙ ψ n + ∆ ψ n + λ | ψ n | ψ n + η ∆ | ψ n | ψ n = 0 . (28)We note that, since M, Ω and J are all fixed, the parameter λ inherits the adjustability of χ .Now, in a stationary state, the time dependence of ψ n can be at most a variation of the phasefactor. Following [15], we consider the ansatz ψ n ( τ ) = exp ( iρH τ + ikξ ) φ ( ξ ) (cid:12)(cid:12) ξ =( n − N/ R , (29)where ξ is a real, continuous variable in the domain [ − N R/ , N R/ φ is a real, twice-differentiable function, and H and k are constants. In particular, H is an energy eigenvalue,in the sense that iρ − ˙ ψ n = − H ψ n . (30)In the limit N (cid:29) R becomes small compared to the domain size, which enables us to invokethe continuum approximation, ψ n ± = exp (cid:0) iρH τ + ik ( ξ ± R ) (cid:1) (cid:32) φ ( ξ ) ± Rφ (cid:48) ( ξ ) + R φ (cid:48)(cid:48) ( ξ ) + O ( R ) (cid:33)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ξ =( n − N/ R , (31)implying | ψ n ± | = φ ( ξ ) ± Rφ ( ξ ) φ (cid:48) ( ξ ) + R φ ( ξ ) φ (cid:48)(cid:48) ( ξ ) + R (cid:0) φ (cid:48) ( ξ ) (cid:1) + O ( R ) (cid:12)(cid:12)(cid:12) ξ =( n − N/ R . (32)Putting eqs. (29) to (32) into eq. (28), then dividing the result by exp ( iρH τ + ikξ ) andretaining terms up to O ( R ), we obtain0 = − H φ ( ξ ) + (cid:104) cos( kR ) (cid:0) φ ( ξ ) + R φ (cid:48)(cid:48) ( ξ ) (cid:1) − φ ( ξ ) + i sin( kR ) (cid:0) Rφ (cid:48) ( ξ ) (cid:1)(cid:105) + λφ ( ξ ) + ηR (cid:104) φ ( ξ ) φ (cid:48)(cid:48) ( ξ ) + 2 (cid:0) φ (cid:48) ( ξ ) (cid:1) (cid:105)(cid:12)(cid:12)(cid:12)(cid:12) ξ =( n − N/ R . (33)The last term in eq. (33) is equivalent to ηR d ( φ ( ξ )) / d ξ . Equating imaginary parts ofeq. (33) gives us k = 0. After the scaling x := ξ/R , the real part of eq. (33) becomes − H φ + φ xx + λφ + η ( φ ) xx φ = 0 when x = n − N/ . (34)The subscript x denotes differentiation with respect to x . We seek φ ( x ) which satisfies eq. (34) for all x , not just when x = n − N/
2. Then, from such a φ ( x ) we will be able to recover thediscrete solution ψ n ( τ ) via ξ = xR and eq. (29). Further to being globally defined, we requirethat φ ( x ) has vanishing derivatives at infinity, and satisfies the normalisation condition, (cid:90) ∞−∞ φ ( x ) d x = 1 . (35)6f η = 0 (i.e. β = 1), then eq. (34) reduces to the nonlinear Schr¨odinger equation with a cubicnonlinearity, which has a well-known solution satisfying all the above constraints [9], H = λ / , φ ( x ) = ± (cid:114) λ λx x. (36)Consider η > β < φ ) xx ≡ φφ xx + 2( φ x ) , we rewrite eq. (34) as − H φ + φ xx (cid:0) ηφ (cid:1) + λφ + 2 η ( φ x ) φ = 0 . (37)We dedicate the remainder of this section to analysing eq. (37). It is an autonomous equationfor φ ( x ), which allows us to define h ( φ ) := φ x , and write φ xx ≡ d( φ x )d φ φ x = hh φ . (38)We define y ( φ ) := h , so that y φ = 2 hh φ , and multiply eq. (37) by 2 to obtain (cid:0) ηφ (cid:1) y φ + 4 ηφy = 2 H φ − λφ . (39)The l.h.s. of eq. (39) is the total derivative of (1 + 2 ηφ ) y with respect to φ . We therefore have y ( φ ) = (cid:82) (cid:0) H φ − λφ (cid:1) d φ ηφ = H φ − λφ / C ηφ . (40)The integration constant C is determined by considering the limit x → ∞ , in which φ → y ≡ ( φ x ) →
0. We therefore have C = 0. Now we note that, if H ≤
0, then the r.h.s.of eq. (40) is negative whenever φ (cid:54) = 0, so it cannot equal the l.h.s. which is ( φ x ) . Thus, if H ≤ φ ( x ) satisfying eq. (40) is identically zero. We therefore require H > φ , we obtain(2 φφ x ) = 4 H φ − λφ ηφ . (41)We then define Φ := φ , and eq. (41) becomes(Φ x ) = 4 H Φ − λ Φ η Φ . (42)If eq. (42) has a solution which is globally non-negative and twice-differentiable, has vanishingderivatives at infinity, and satisfies (cid:90) ∞−∞ Φ( x ) d x = 1 , (43)then we claim that Φ( x ) must attain its global upper bound of 2 H /λ at some finite x , andthat every local maximum of Φ must also be a global maximum. The proof of this claim isas follows. Since Φ cannot be identically zero, and since lim x →±∞ Φ( x ) = 0, Φ( x ) must haveat least one turning point, at some finite x and non-zero Φ. But we observe from eq. (42)that Φ x vanishes if and only if Φ = 0 or 2 H /λ . Therefore, Φ( x ) must attain its global upperbound of 2 H /λ at least once, and no other local maximum value is possible. This concludesthe proof. We further propose that wherever Φ( x ) attains its maximum value, say at x = x max ,the second derivative Φ xx does not vanish there. The proof is as follows. On the one hand,we have Φ xx ≡ φφ xx + 2( φ x ) ; when x = x max , we also have φ x = Φ x / (2 φ ) = 0, thereforeΦ xx = 2 φφ xx . On the other hand, eq. (34) is equivalent not only to eq. (37) but also to7 φ Φ xx = H φ − φ xx − λφ . It follows that, at x = x max , we have (1 / (2 φ ) + ηφ )Φ xx = H φ − λφ ,and therefore Φ xx = (2 H Φ − λ Φ ) / (1 + 2 η Φ). Since Φ( x max ) = 2 H /λ , it follows thatΦ xx ( x max ) = − H / ( λ + 4 ηH ) <
0, as required. A corollary of this proposition is thatthere must exist some neighbourhood of x max containing no maxima of Φ( x ) other than x max itself. Without loss of generality, let x max = 0. Suppose the corollary is false, so that everyneighbourhood of 0 contains some non-zero x which is a maximum of Φ( x ). Then, theremust exist some sequence x n approaching 0 such that Φ( x ) has a maximum at every x n , withΦ( x n ) = Φ(0). But this leads to a contradiction. Indeed, for every x n we have the Taylorexpansion Φ( x n ) = Φ(0) + Φ xx (0) x n / O ( x n ), where the first derivative is absent becauseΦ( x ) has a maximum at 0. It then follows that Φ xx (0) = lim n →∞ x n ) − Φ(0)) /x n = 0,which contradicts the previous proposition. Therefore the corollary is proven. Since 0 has amaxima-free neighbourhood, we say that Φ( x ) has an isolated maximum at 0.We note that we can indeed require that x = 0 is a maximum of Φ( x ), because eq. (42) istranslationally invariant: if Φ( x ) is a solution then so is Φ( x − c ) for any constant c . We exploitthis invariance, requiring that Φ( x ) satisfiesΦ := Φ(0) = 2 H /λ. (44)Now we claim that there exists b >
0, which may be infinite, such that lim x → b Φ( x ) = 0, andΦ( x ) (cid:54) = 0 for all x ∈ (0 , b ). The proof of this claim is as follows. If Φ( x ) (cid:54) = 0 for all x ∈ (0 , ∞ ),then we are done. If Φ( x ) = 0 for some x ∈ (0 , ∞ ), then the set { x ∈ (0 , ∞ ) : Φ( x ) = 0 } musthave a minimum. If it does not, then there would be a sequence x n > n → ∞ , x n → x n ) →
0; but this would contradict the continuity of Φ( x ) at x = 0. Thus,letting b equal the least positive zero of Φ( x ), then we are done.Next, we propose that no other solution on [0 , b ) exists, and the proof is as follows. If Φ( x )has any maxima in (0 , ∞ ), then the set M := { x ∈ (0 , ∞ ) : x is a maximum of Φ( x ) } musthave a minimum, because otherwise we would have a contradiction to the fact that x = 0 isan isolated maximum of Φ( x ). Let x = min M , and suppose x < b . Since Φ( x ) = Φ(0), andsince Φ( x ) has no maximum in the interval (0 , x ), and since Φ( x ) is continuous, it must attainits minimum value at some point b (cid:48) ∈ (0 , x ). But that implies Φ( b (cid:48) ) = 0, where b (cid:48) < x < b ,contradicting the fact that Φ( x ) (cid:54) = 0 for all x ∈ (0 , b ). Therefore, we must have x > b .Since Φ x vanishes only at maxima and minima, it follows that Φ x is non-vanishing on (0 , b ),and therefore Φ( x ) is strictly decreasing on [0 , b ). That is, any solution to eq. (42) on [0 , b )satisfying all the aforementioned contraints must also satisfyΦ x = − g (Φ) := − (cid:115) H Φ − λ Φ η Φ = − (cid:112) H Φ (cid:115) − Φ / Φ η Φ , ≤ x < b, Φ ≥ Φ > . (45)On any closed interval [Φ , Φ ] ⊂ (0 , Φ ), the function g (Φ) is continuous and non-zero, sothe reciprocal function 1 /g (Φ) is continuous and bounded, and therefore Riemann integrable.But g (Φ) approaches 0 as Φ → Φ , meaning 1 /g (Φ) becomes unbounded. Thus, integration of1 /g (Φ) on the interval [Φ , Φ ] is not trivial. Luckily Φ = Φ is an integrable singularity of thefunction 1 /g (Φ), because the Puiseux series of 1 /g (Φ) about Φ is O ((Φ − Φ ) − / ). Therefore,for any Φ ∈ (0 , Φ ], eq. (45) is equivalent to (cid:90) Φ Φ g (Φ) dΦ = x (Φ ) − x (Φ ) = x (Φ ) , (46)which determines a unique x (Φ ) ∈ [0 , b ). The l.h.s. of eq. (46) is a strictly decreasing functionof Φ , meaning x (Φ ) has a unique inverse function which is also strictly decreasing, Φ ( x ), on8he domain x ∈ [0 , b ). But Φ( x ) is an existing function satisfying eq. (45) for all x ∈ [0 , b ).Therefore, we must have Φ ( x ) = Φ( x ) for all x ∈ [0 , b ), and the uniquess of Φ( x ) follows.In summary, we have so far established the following. If eq. (42) has a solution which is globallynon-negative and twice-differentiable, has vanishing derivatives at infinity, and has the propertythat its integral over R is 1, then eq. (42) has a solution, say Φ( x ), which satisfies all the aboveconstraints as well as the condition (44), and there exists some b > x ) is strictly decreasing on [0 , b ), and lim x → b − Φ( x ) = 0, and Φ( x ) is the uniquesolution on [0 , b ). Moreover, using exactly the same arguments as above, it can also be shownthat there exists some a < x ) is strictly increasing on( a, x → a + Φ( x ) = 0, and Φ( x ) is the unique solution on ( a, interval ofuniqueness , ( a, b ), Φ x is given byΦ x = G ( x, Φ) := − sgn( x ) g (Φ) , (47)where g is defined by eq. (45), and sgn is the sign function.Now we describe a method which, given λ and η , determines the unique Φ( x ) on ( a, b ), andalso determines a, b, H in the process. Indeed we will show that for any λ and η , the intervalof uniqueness for Φ( x ) must be ( a, b ) = R . The fact that ( a, b ) = R shall have the followingsubtle consequence. Note that the derivation of eq. (42) involved a multiplication by Φ ≡ φ .Thus, the deduction from eq. (42) back to eq. (37) holds on the condition that Φ (cid:54) = 0. Since a and b are the smallest (in magnitude) zeros of Φ( x ), we see that the equivalence betweeneqs. (37) and (42) breaks down outside the interval ( a, b ). That is to say, eqs. (37) and (42)are equivalent globally if and only if ( a, b ) = R .The method is as follows. For x ∈ [0 , b ), consider the coordinate transformation, Z (Φ) := arsech (cid:0) Y (Φ) (cid:1) , where Y (Φ) := (cid:114) ΦΦ = (cid:114) λ Φ2 H . (48)Φ( x ) is a bijection from [0 , b ) to (0 , Φ ], Y (Φ) is a bijection from (0 , Φ ] to (0 , ,
1] to [0 , ∞ ). Therefore, all the coordinatetransformations are invertible. For x ∈ ( a, Z with respect to x we find, for all x ∈ ( a, b ), Z x = Z Y · Y Φ · Φ x = − Y √ − Y · √ Φ Φ · (cid:0) − sgn( x ) g (Φ) (cid:1) = sgn( x )2Φ (cid:112) − Φ / Φ · g (Φ)= sgn( x ) √ H √ η Φ , (49)where we have used definition (45) of g (Φ). Moreover, by definition we have Y = sech Z , andit follows that 2 η Φ = 2 η Φ Y = 2 η (2 H /λ )sech Z . Defining ν := 4 ηH /λ, (50)we rewrite eq. (49) as Z x = sgn( x ) √ H (cid:112) ν sech Z . (51)9ue to eq. (48), we have Z ( x = 0) = 0. We can therefore solve eq. (51) as follows.sgn( x ) (cid:112) H (cid:90) x d˜ x = (cid:90) Z (cid:112) ν sech ˜ Z d ˜ Z, (52)implying sgn( x ) (cid:112) H x = arsinh (cid:18) sinh Z √ ν (cid:19) + √ ν arctan (cid:32) √ ν sinh Z (cid:112) ν + cosh Z (cid:33) . (53)Now we can determine the values of a and b . In the limit Z → + ∞ , the definition of thecoordinate transformations, as per eq. (48), dictates that we must have either x → a or x → b .At the same time, eq. (53) dictates that we must have x → ±∞ , because the arctan functionon the r.h.s. of eq. (53) is bounded, whilst the arsinh function diverges to + ∞ . It thereforefollows that ( a, b ) = R .The next step is to rewrite eq. (53) as an expression for x in terms of Φ, so that we can invertthe expression to find Φ( x ) for x ∈ R . By definition (48) we have cosh Z = 1 /Y = Φ / Φ,and it follows that sinh Z = cosh Z − / Φ) −
1. Since Z is by definition non-negative,we must take the positive square root, sinh Z = (cid:112) (Φ / Φ) −
1. Then eq. (53) becomessgn( x ) (cid:112) H x = arsinh (cid:115) − (Φ / Φ )(1 + ν ) (Φ / Φ ) + √ ν arctan (cid:115) ν (cid:0) − (Φ / Φ ) (cid:1) ν (Φ / Φ ) . (54)We claim that, given Φ > x ∈ R , eq. (54) uniquely determines a value of Φ >
0. Theproof is as follows. If x = 0, then immediately from eq. (54) we have Φ = Φ , and we are done.If x (cid:54) = 0, consider the function G (Φ) := arsinh (cid:115) − (Φ / Φ )(1 + ν ) (Φ / Φ ) + √ ν arctan (cid:115) ν (cid:0) − (Φ / Φ ) (cid:1) ν (Φ / Φ ) − sgn( x ) (cid:112) H x, (55)where x and Φ are parameters. Differentiating eq. (55) with respect to Φ, we findd G dΦ = − (cid:115) ν (Φ / Φ )1 − (Φ / Φ ) < > . (56)This means G is strictly decreasing for Φ >
0. Since G (Φ) → ∞ in the limit Φ →
0, and G (Φ ) = − sgn( x ) √ H x <
0, and G is continuous, we must have G vanishing at exactly onevalue of Φ ∈ (0 , Φ ). This concludes the proof. Moreover, we observe that in eq. (54) the l.h.s.is invariant under x (cid:55)→ − x . Thus, on R we have Φ( − x ) = Φ( x ).In practice, given any Φ > x ∈ R , we can compute Φ( x ) by locating the zero of G (Φ).However, the value of Φ cannot be freely chosen. Instead, it is determined by the normalisationcondition (43) which, since Φ( x ) is an even function on R , now reads1 = 2 (cid:90) ∞ Φ( x ) d x = 2 (cid:90) ∞ Z =0 Φ Z + x d Z, (57)where the Z + x is the positive- x branch of Z x , as per eq. (51). It then follows that1 = 2 (cid:90) ∞ Φ (cid:112) ν sech Z √ H d Z. (58)10sing Φ = Φ sech Z , we deduce √ H = (cid:90) ∞ sech Z (cid:112) ν sech Z d Z (59a)= 12 + (1 + ν ) arctan √ ν √ ν . (59b)Multiplying eq. (59b) by 2 √ ν , replacing H by λ Φ /
2, and replacing ν by 4 ηH /λ = 2 η Φ , weobtain the following transcendental equation for Φ . (cid:112) λη = (cid:112) η Φ + (1 + 2 η Φ ) arctan (cid:112) η Φ . (60)To show that exactly one solution to eq. (60) exists, we consider the function F (Φ ) := (cid:112) η Φ + (1 + 2 η Φ ) arctan (cid:112) η Φ − (cid:112) λη. (61)Differentiating F (Φ ) with respect to Φ , we findd F dΦ = 2 η (cid:18) √ η Φ + arctan (cid:112) η Φ (cid:19) > > . (62)This means F (Φ ) is strictly increasing for Φ >
0. Since lim Φ → F (Φ ) = −√ λη <
0, and F (Φ ) → ∞ in the limit Φ → ∞ , and F is continuous, we must have F vanishing at exactlyone value of Φ >
0. In practice, given parameters λ and η , we can compute Φ by locating thezero of F (Φ ), and Φ uniquely determines the energy eigenvalue, H = λ Φ /
2. We can thenfeed the value Φ into eq. (54), and then for every x ∈ R we can find Φ( x ) by means we havealready described. In summary, given λ and η , eqs. (44), (54) and (60) together constitute ananalytical solution to eq. (42); and as we have already proven, it must be the unique globalsolution to eq. (42) which satisfies all the constraints we have imposed.We note that if the parameter η →
0, we should recover the solution to the nonlinearSchr¨odinger equation, given by eq. (36); and indeed we do. Firstly, in the limit η → ν →
0, which means we cannot use eq. (60) to determine Φ , because the derivationof eq. (60) involved a multipication by √ ν . Instead, we must extract Φ from eq. (59). In thelimit ν →
0, eq. (59a) is simply (cid:112) λ/ (8Φ ) = (cid:82) ∞ sech Z d Z = 1. It follows that Φ = λ/ H = λ Φ / λ /
16, againagreeing with eq. (36). Finally, when ν →
0, eq. (54) is simplysgn( x ) (cid:112) H x = arsinh (cid:112) (Φ / Φ) − , (63)which is equivalent to Φ / Φ = 1 + sinh ( √ H x ) = cosh ( λx/ sech ( λx/ H provides a link between Φ( x ) and the binding energy of the stationarypolaron. By eqs. (11), (18) and (29), where k = 0 and ξ = xR , the polaron state is written interms of local excitations as | Ψ (cid:105) = (cid:80) Nn =0 α n ˆ A † n | vac (cid:105) , and in the limit N (cid:29)
1, we have α n = φ ( n − N/
2) exp (cid:20) − it (cid:126) ( J − J − H J ) (cid:21) , φ ( x ) = Φ( x ) . (64)Thus, the stationary | Ψ (cid:105) solves i (cid:126) d | Ψ (cid:105) / d t = ( ˆ H e + ˆ H int ) | Ψ (cid:105) as well as satisfying i (cid:126) d | Ψ (cid:105) / d t =( J − J − H J ) | Ψ (cid:105) . By definition, the polaron’s binding energy, E b , is its total internal energymeasured with respect to J . In units of J , we have E b := (cid:104) Ψ | ˆ H e + ˆ H p + ˆ H int | Ψ (cid:105) − J J . (65)11n expression for (cid:104) ˆ H p (cid:105) in terms of Φ( x ) can be found by using eqs. (4) and (24). Since thepolaron is stationary, the kinetic part of (cid:104) ˆ H p (cid:105) is zero, so we have (cid:104) ˆ H p (cid:105) J = M Ω R J N − (cid:88) n =0 ( u n +1 − u n ) = σδ ρ N − (cid:88) n =0 (cid:104) ( β − | ψ n +1 | − ( β + 1) | ψ n | (cid:105) = λ N − (cid:88) n =0 (cid:104) ( β − | ψ n +1 | + ( β + 1) | ψ n | (cid:105) + η N − (cid:88) n =0 | ψ n +1 | | ψ n | . (66)It then follows that E b = − − H + (cid:104) ˆ H p (cid:105) J = − − H + λ N − (cid:88) n =0 (cid:104) ( β − Φ( n + 1 − N/ + ( β + 1) Φ( n − N/ (cid:105) + η N − (cid:88) n =0 Φ( n + 1 − N/
2) Φ( n − N/ . (67)We have made use of definition (27) of λ, η , the fact that | ψ n | = | φ n | for all n , as well asthe fact that | φ n | is approximated by Φ( n − N/ symmetry parameter β , and the effective couplingparameter λ . Recall that the former is a measure of the spatial symmetry of the electron-phonon interaction, and the latter measures the strength of this interaction. These are theonly two parameters that affect the stationary polaron’s physical properties (as η is merely aconvenient combination of β and λ ).Figure 1a shows how Φ and the half-width of the polaron varies with β and λ . We define thehalf-width as the distance between the two x -values at which Φ( x ) = Φ /
2, and it is a measureof how localised the polaron is. As one would expect, the half-width is negatively correlatedwith Φ , which is the maximum height of Φ( x ). The figure shows Φ increasing with λ , andhalf-width decreasing with λ , and the rate of change of each quantity is greater given largervalues of β . That is to say, the more spatially asymmetric the electron-lattice interaction is,the more influential λ is. The figure also has the following implication on the accuracy of Φ( x )as an approximation to the discrete stationary solution to eq. (19). In a discrete solution, ψ n = exp( iρH τ ) φ n , the physical interpretation of | ψ n | is the probability of the electron beinglocalised at the n th lattice site. Therefore, the normalisation condition is defined in terms of asum, (cid:80) Nn =0 | ψ n | = 1, and consequently we must have | ψ n | ≤ n . When a continuumsolution Φ( x ) is used to approximate the discrete one, we have the relation Φ ≡ max | ψ n | .Thus, any continuum solution with Φ > β = 1,Φ exceeds 1 if λ is greater than 8, since Φ = λ/
8. On the other hand, when β = 0, wecomputed Φ for λ up to 100, and Φ remains less than 0.6.In fig. 1b we see that H increases with λ whilst the polaron’s binding energy gains magnitude,meaning the larger λ is the more energy is required to break up the polaron. Once again, thelarger β is, the more rapidly these quantities vary with λ . We note that the thick (black) curvefor H , corresponding to β = 1, is exactly the graph of H = λ /
16, as per eq. (36). Comparingfigs. 1a and 1b, we see that a polaron which is more strongly bound has a larger Φ and asmaller half-width, i.e., it is more localised. 12 a) (b)Figure 1: (Colour online.) The height Φ of the analytical Φ( x ) solution ((a), right axis), thepolaron half-width ((a), left axis), the energy eigenvalue H ((b), right axis), and the polaronbinding energy ((b), left axis), according to analytical solutions Φ( x ). The dependenceof each quantity upon β and λ is represented by a family of curves. The thick (black)curve always corresponds to β = 1, and as β decreases towards 0, the thin (blue) curves,corresponding to β = 0 . , . , . . . , . ,
0, become either steeper or shallower.
The gradient of curves in fig. 1b vary with β , and the variation is more pronounced when β is close to 1. This suggests that the system is highly sensitive to variations in β when β islarge, but not so when β is small. Moreover, as λ decreases, curves corresponding to differentvalues of β begin to converge; specifically this happens when λ ≈
1. This suggests that when λ is small, the extent to which the electron-phonon interaction is spatially symmetric has littlebearing on the physical properties of stationary polarons. In this section we solve eqs. (28) and (30) directly, using a numerical scheme, but not withoutthe help of analytical results from section 3.1. We then compare the resulting stationarypolaron states with the ones we obtained via continuum approximation.Expanding eq. (28) using the definitions of ∆ ψ n and ∆ | ψ n | , we have − H ψ n + ( ψ n +1 + ψ n − − ψ n ) + λ | ψ n | ψ n + η (cid:16) | ψ n +1 | + | ψ n − | − | ψ n | (cid:17) = 0 . (68)Any solution ψ n to eq. (68) is an attractor of the following map [27]. ψ n (cid:55)→ H ( ψ n ) (cid:107)H ( ψ n ) (cid:107) , (69)where H ( ψ n ) := ( ψ n +1 + ψ n − − ψ n ) + λ | ψ n | ψ n + η (cid:16) | ψ n +1 | + | ψ n − | − | ψ n | (cid:17) , (70a) (cid:107)H ( ψ n ) (cid:107) = (cid:118)(cid:117)(cid:117)(cid:116) N (cid:88) n =0 H ( ψ n ) . (70b)13 a) (b)Figure 2: (Colour online.) (a) The maximum localisation probability max | ψ n | (rightaxis), and polaron binding energy (left axis), as functions of β and λ . The curves for β = 0 largely overlap those for β = 0 .
2, so for practical reasons they are plotted onseparate scales. For each value of β , we are interested only in those λ for which max | ψ n | is not too close to the extremes, i.e. 0 and 1.(b) Thick lines: numerical stationary solutions, | ψ n | (right axis), and the associated u n (left axis). Solutions are shifted along the n -axis, to avoid overlap. From left toright: β = 0 , . , . , . ,
1. Dotted (red) lines: λ = 1 .
0. Solid (black) lines: λ = 2 . λ = 4 .
4. Thin lines: approximate solutions, by analytical methodsof section 3.1. From left to right: β = 0 , . , . , . ,
1; all with λ = 2 . We take the approximate solution from section 3.1 as initial guess, and repeatedly apply eq. (69)until we reach convergence. When converged, ψ n is the stationary solution to eq. (68), and (cid:107)H ( ψ n ) (cid:107) is equal to H . In practice, on a grid with N = 200, convergence is typically reachedwithin O (10 ) iterations, which amounts to O (10 ) seconds of computing time. We havecomputed stationary solutions for various β and λ , and some results are presented in fig. 2.Figure 2a contains information about two key aspects of the stationary polaron state: theelectron probability distribution, and the polaron binding energy. Qualitatively speaking, it isin agreement with predictions of the continuum approximation, as per fig. 1: as λ increases, thepolaron becomes more localised, and more strongly bound. Moreover, the effect of increasing λ is more profound given larger values of β . However, further comparison between figs. 1aand 2a reveals a noteworthy difference. When β = 1, Φ as a linear function of λ , whereasfig. 2a suggests that max | ψ n | , which is approximated by Φ , is not linearly dependent on λ .In fact, given any β , max | ψ n | grows significantly faster with λ than fig. 1a predicts. Despitethat, the growth of max | ψ n | in fig. 2a eventually stalls, when λ becomes sufficiently large. Thisis a manifestation of a fundamental difference between the continuum and discrete equations,which we discussed in section 3.1: the continuum equations place no limit on how large Φ canbe, whereas the discrete system limits max | ψ n | to 1.Figure 2b shows a selection of | ψ n | solutions. Comparing all the dotted (red) lines, whichcorrespond to λ = 1 . β , we see that they are essentially identical. Thisconfirms the belief that when λ is close to 1, systems with different β -values unify. Thefigure also shows some stationary solutions to the other half of eq. (19), namely u n , which isexpressed in terms of the stationary | ψ n | solution as per eq. (24). Recall that physically u n represents the displacement of the n th molecule from its equilibrium position. In order that thepoint-dipole model for lattice units is valid, the lattice distortion must satisfy the condition14 u n +1 − u n | (cid:28) | u n +1 − u n | ∼ O (10 − ).Comparing all the dashed (blue) lines in fig. 2b, which correspond to λ = 4 . β , enables us to make the following observation. When β = 0, the u n solution is centred atthe location of max | ψ n | , in the sense that its graph is rotationally symmetric about n = 80.This agrees with our intuition that when β = 0, i.e. when the electron-phonon interactionis spatially symmetric, the electron in the stationary state causes equal lattice distortion toits left and right. As β increases, the maximum magnitude of u n remains the same, but thecentre of u n shifts away from the location of max | ψ n | , in response to the decrease in spatialsymmetry. When β = 1, the molecule at the location of max | ψ n | ( n = 120 in this case) isbarely displaced, whereas molecules to the right of this point are displaced considerably. Now,the potential energy in the lattice is a sum over terms of the form ( u n +1 − u n ) , which is thesquare of the gradient of the u n graph at site n . In the steady state, this gradient is zero exceptat a few sites around the location of max | ψ n | , and it is clear that solutions corresponding tolarger values of β have steeper gradients there. We therefore conclude that, in the stationarystate, systems with greater spatial asymmetry store more potential energy in the lattice.In fig. 2b we also see a comparison between some | ψ n | solutions and their counterpart con-tinuum approximations, Φ( x ). In particular, we look at the thick solid (black) lines and theiraccompanying thin solid (black) lines. The comparison reveals that, fixing λ , in this case λ = 2 .
6, Φ( x ) is a more accurate approximant to | ψ n | when β is smaller. As β approaches 1, itbecomes apparent that Φ( x ) under-estimates the height of the | ψ n | profile. If λ is sufficientlylarge, however, Φ( x ) becomes an over-estimate of the profile height. This all comes down onceagain to the fact that the continuum equations do not limit the height of the Φ( x ) solution. In this part of the study we explore properties of polarons which propagate along the peptidechain, under an external forcing (cid:15) ( τ ), and zero temperature ( f n ( τ ) = 0). Physically, (cid:15) ( τ ) mayrepresent the strength of a time-dependent electric field. We solve eq. (19) as an initial valueproblem, using the stationary ψ n and u n solutions which we computed in section 3.2 as theinitial configuration of the system. We prescribe a suitable (cid:15) ( τ ), setting n to the locationwhere the stationary | ψ n | attains its maximum. Then we integrate the system forward intime using the 4 th -order Runge-Kutta method. To ensure numerical stability, the integrationtime-step is set at ∆ τ = 0 . | ψ n | attains its maximum at lattice site ¯ n ,then polaron position is the vertex location of the parabola extrapolated from three points:(¯ n, | ψ ¯ n | ) , (¯ n − , | ψ ¯ n − | ) , (¯ n + 1 , | ψ ¯ n +1 | ). We note that, if the polaron is dynamical, then thestationary solution given by eq. (64) is no longer valid, and therefore we cannot take eq. (67)as the expression for the binding energy. Instead, the binding energy E b as per eq. (65) willbe computed directly from the numerical solutions.15 .1 Constant or periodic electric fields The most obvious choice of (cid:15) ( τ ) is a constant, (cid:15) ( τ ) = ¯ (cid:15) > τ ≥ . (71)Using moderately-localised stationary states (with max | ψ n | ≈ .
6) as initial conditions, wecomputed polaron trajectories under various values of ¯ (cid:15) . Our results show that, given β = 1and λ = 3 .
0, a constant forcing of any ¯ (cid:15) ∼ O (10 − ) induces nothing but small oscillationsof the polaron around its initial position. An example of trajectory is presented in fig. 3. Figure 3: (Colour online.) Some polaron trajec-tories, given β = 1 and λ = 3 .
0, under either aconstant or a periodic forcing (cid:15) . Solid (black) line: (cid:15) = 0 .
1. Dashed (blue) line: (cid:15) = 0 . πτ /T ), T = 500. 1000 units of τ equals 1.8 nanoseconds. As ¯ (cid:15) is increased beyond 0.1, we findthat eventually the forcing does be-come strong enough to dislodge theelectron from its potential well, andpropel the polaron along the peptidechain. However, as the polaron prop-agates, the magnitude of its bindingenergy decreases rapidly, and the po-laron “delocalises”, i.e. breaks up intounbound components, within severalhundred time units. A direct mani-festation of the polaron’s energy lossand eventual delocalisation is that the | ψ n | profile loses height and gains localpeaks at lattice sites far away from theglobal maximum. For the sake of con-sistency, throughout the remainder ofthis study we shall say that a polaronhas delocalised if its maximum height drops to below 0.1, as it must then be the case thatother local peaks have magnitudes comparable to the global maximum. Figure 4 shows anexample of a constant forcing large enough to cause polaron displacement, and it illustratesthe resultant rapid delocalisation of the polaron. (a) (b)Figure 4: Polaron propagation under β = 1 , λ = 3 . (cid:15) = 0 . τ .(b) The | ψ n | profile of the polaron upon delocalisation. Figure 4a shows the trajectory of a polaron which, within roughly 600 time units, is displacedby just over 300 lattice sites. Its binding energy steadily decreases in magnitude, until thepolaron delocalises at τ ≈ ψ n | , at the time of delocalisation. This | ψ n | profile has evolved from an initial configurationpossessing a maximum height of 0.64, and no local peaks apart from the global maximum.If the polaron’s binding energy decreases in magnitude, then the polaron’s ability to transportenergy is diminished. Beyond the example shown in fig. 4, all our results are consistent withthe hypothesis that, regardless of β and λ , a constant (cid:15) causes the polaron to undergo eithersmall periodic oscillations, or rapid losses in energy. We would like to find ways to displacethe polaron without significant energy loss. Therefore, we must look for forms of (cid:15) other thanconstants. The next most natural choice of (cid:15) is periodic, (cid:15) ( τ ) = A sin 2 πτT for τ ≥ , (72)where A is the amplitude and T is the period. Physically this may represent an electromagneticplane wave which is monochromatic, i.e. coherent. Under periodic (cid:15) ( τ ) with A up to 0.2,regardless of β and λ we find that the polaron simply oscillates about its initial position. Thepolaron’s oscillatory motion has a period which coincides with T , and an amplitude whichis positively correlated with A . An example of such trajectories is shown in fig. 3. Whilethe polaron remains highly stable over time, its position averaged over its periodic remainsconstant. Thus, if we want polarons which transport energy from one lattice site to another,we must again look for an alternative form of (cid:15) ( τ ). Having studied the effects of constant forcing and periodic forcing in section 4.1, and discoveredthat neither serves to displace the polaron with minimal energy loss, in this section we considerforcing of the form (cid:15) ( τ ) = ¯ (cid:15) + A sin 2 πτT . (73)Equation (73) represents the combination of the two types of forcing considered previously,with a constant component and a sinusoidal one. One may also think of (cid:15) ( t ) as a mean-shiftedperiodic forcing (MSPF). In particular, the mean ¯ (cid:15) is chosen to be lower than the constantforcing (cid:15) which is required to displace the polaron, in the manner of fig. 4. Therefore, ¯ (cid:15) onits own would not give the electron enough energy to escape its potential well. But we hopethat the component A can periodically push the electron energy over the threshold, resultingin polaron motion. Another possible advantage of this setup is that A may periodically lowerthe electron energy, slowing it down and giving the lattice time to “catch up”, thus makingthe polaron motion more sustainable than it would be under a constant forcing.Mathematically, (cid:15) ( τ ) depends on three independent parameters, ¯ (cid:15), A and T . Before investigat-ing the effect of each of these parameters, we present fig. 5, which is a direct comparison withfig. 4.We have replaced the constant forcing (cid:15) = 0 .
15, which resulted in fig. 4, with an MSPFwhich has the same maximum amplitude as before. The difference is that now this maximumamplitude is reached once every period T . Figure 5a shows that, within roughly 10 periods,the polaron is displaced by nearly 400 lattice sites. Contrary to the uniform manner in whichthe polaron moves in fig. 4a, now the polaron moves towards one end of the peptide chain andthen the other, within each period of (cid:15) ( τ ). The overall displacement of the polaron is due tothe fact that each movement to one end of the chain is larger than the subsequent swing backthe other way. We note that while the polaron moves slightly further compared to fig. 4a, its17 a) (b)Figure 5: Polaron propagation under β = 1 , λ = 3 . (cid:15) ( τ ) = 0 .
025 + 0 .
125 sin(2 πτ / τ .(b) The height (right axis) and half-width (left axis) of the | ψ n | profile as functions of τ . lifetime , i.e. the amount of time elapsed before delocalisation, is much longer. Overall, thepolaron in fig. 5 propagates with a lower (average) velocity , V , defined by V = average position over final complete period of motion − initial positionnumber of complete periods × T , (74)where the numerator is the displacement of the polaron, which we denote by D . Our resultsshow that, of the parameters ¯ (cid:15), A and T , the dominant factor which determines the polaron’svelocity is the constant component ¯ (cid:15) . We will discuss this in more depth in relation to fig. 6.Figure 5a also shows how the polaron’s binding energy, E b , varies in time. Following an initialdrop in magnitude, E b mostly oscillates between − .
75 and − .
5, until another sharp decreasein magnitude leading up to delocalisation at τ ≈ n = 300, E b is about − .
6. In fig. 5a, the polaron’s average position over the 6 th period is roughly 300, and theaverage binding energy over this period is − .
2. That is to say, under the MSPF, the polaronis carrying twice as much energy when it reaches n = 300, compared to when it reaches n = 300under the constant forcing. Even though the constant forcing gets the polaron to n = 300 inless time, we consider the MSPF a better mechanism for polaron transport, because the polaronbinding energy is more stable. Indeed, the same can be said when the destination n is anythinglarger than 30. If the destination is n <
30, then the constant forcing takes the polaron to n in such a small amount of time that it causes no more variation in binding energy than theMSPF does. In general, all our results are consistent with the hypothesis that, by splitting aconstant forcing into constant and sinusoidal components, we lower the polaron’s velocity butincrease its stability and lifetime. We say that, compared to the constant forcing, the MSPF isa better long-distance tranport mechanism, where speed can be sacrificed for energy efficiency.In fig. 5b we see another aspect of the polaron’s motion, namely, how the height and half-widthof the | ψ n | profile vary with time. Following an initial decrease, the profile height, max | ψ n | ,mostly oscillates between 0.2 and 0.4, until a sudden drop to 0.1, leading to delocalisation.Meanwhile, the half-width mostly oscillates between 1.5 and 4, following an initial growth.The peaks in the half-width, as well as the troughs in max | ψ n | , occur precisely when thepolaron turns from moving in one direction to moving in the other. This suggests that whenthe polaron accelerates, it “spreads out”, and so the half-width widens and max | ψ n | drops.We observe this phenomenon in all our results.Based on our observations, we theorise that a polaron’s directed motion may be explainedphysically as follows. Since the forcing (cid:15) ( τ ) is the effect of an electric field, it first-and-foremost18 a) ¯ (cid:15) = 0 . , T = 500. (b) ¯ (cid:15) = 0 . , T = 500.(c) ¯ (cid:15) = 0 . , T = 2000.Figure 6: (Colour online.) Some polaron trajectories under the MSPF, (cid:15) ( τ ) = ¯ (cid:15) + A sin(2 πτ /T ). Dotted (red) lines: A = 0 .
10; solid (black) lines: A = 0 .
15; dashed (blue) lines: A = 0 .
20. Each figure contains 9 trajectories,whose initial positions have been shifted to avoid overlap. Every trajectorystarting from position 200 correspond to a polaron with symmetry param-eter β = 0 and effective coupling parameter λ = 7 .
6. Trajectories startingfrom position 600 correspond to β = 0 . λ = 4 .
9, and those startingfrom position 1000 correspond to β = 1 and λ = 3 . λ has to be variedwith β , in order to keep the initial | ψ n | profiles unchanged. In this case, allinitial conditions have max | ψ n | = 0 . provides the electron with extra energy. This is evident in the dramatic energy variation duringthe first period of (cid:15) ( τ ) (see fig. 5a). Following this, it becomes much easier for the electronto overcome the significantly diminished polaron binding energy, E b . This is why the onset ofpolaron motion always follows a drastic drop in magnitude of E b . Whenever | (cid:15) ( τ ) | becomeslarge enough to give the electron sufficient energy to overcome E b , the electron is dislodgedfrom its potential well and propelled along the lattice. If the electron-lattice coupling is strongenough, then the lattice distortion can keep up with the electron, and so the polaron can remainintact. Whenever | (cid:15) ( τ ) | drops below the binding threshold, the electron-lattice interaction slowsdown the electron and causes its probability distribution to spread out. This is why the half-width of | ψ n | always peaks at times when the polaron’s instantaneous velocity is zero. If | (cid:15) ( τ ) | remains below the binding threshold for long enough, then the polaron’s position can plateau,as is seen in fig. 6, particularly in the lowermost solid (black) lines in figs. 6b and 6c. If (cid:15) ( τ ) hasa large enough periodic component A , then it is possible for | (cid:15) ( τ ) | to overcome the thresholdtwice per period: once with (cid:15) >
0, once with (cid:15) <
0. If the electron moves towards large n in the (cid:15) > n in the (cid:15) < (cid:15) (cid:54) = 019nsures that the electron always spends more time moving one way than the other, hence theoverall directedness of the polaron trajectories. This point is most clearly demonstrated by thetrajectories in fig. 6c, where the period T = 2000.Within each of figs. 6a to 6c, we can compare the trajectories represented by the same line type,but have different starting positions. This reveals the effect of varying the spatial symmetryof electron-lattice interaction. A greater spatial asymmetry, i.e. a larger β , causes the polaronto be more susceptible to displacement. Also within each of figs. 6a to 6c, we can comparetrajectories starting from the same position, but are represented by different line types. Thisreveals the effect of varying the forcing amplitude A . The larger A is, the more the polaronoscillates back and forth during each period of motion. We can also compare trajectories infigs. 6a and 6b which have identical line types and starting positions. This suggests that theoverall velocity of the polaron is determined by ¯ (cid:15) , in the sense that the larger ¯ (cid:15) is, the more thepolaron moves per unit time. Finally, by comparing trajectories in figs. 6b and 6c which haveidentical line types and starting positions, we hope to see the effect of varying T . However,this comparison is not particularly enlightening. We therefore present fig. 7, which not onlyprovides more insight into the effect of T , but also helps to quantify our observations, andreinforce our theories, about the effects of ¯ (cid:15) and A .For several combinations of ¯ (cid:15) and T , we examine how the polaron’s lifetime, displacement andvelocity vary with A . The results are displayed together, in fig. 7, so that the effects of ¯ (cid:15) and T can be gauged also. Firtly we consider the lifetime, τ d . We computed all our numericalsolutions up to τ = 50000, which is several times larger than the typical lifetime of a polaronthat moves under the MSPF. If the polaron is not displaced by the MSPF, then it is effectivelypermanent, in the sense that its energy oscillates instead of dissipating over time, and it wouldhave a lifetime far exceeding 50000. Thus, in fig. 7, the lifetime of a permanent (undisplaced)polaron is represented as τ d = 50000. For each combination of ¯ (cid:15) and T , there exists some critical amplitude , A = A c , below which the polaron is undisplaced by the MSPF. At A = A c ,the combined magnitude of the forcing, (cid:15) comb := ¯ (cid:15) + A , becomes large enough to displace thepolaron, and τ d drops sharply. This drop can sometimes result in a lifetime of only severalthousand time units - see for instance the bottom-left subfigure in fig. 7, corresponding to¯ (cid:15) = 0 .
1. When ¯ (cid:15) is smaller, say ¯ (cid:15) = 0 .
005 (top-left subfigure), the drop in lifetime is lessdramatic. As A increases beyond A c , the polaron’s lifetime drops further, if only slightly.Next, we look at the polaron’s displacement, D . When A is small, the polaron does not movebarring small oscillations, the types of which we saw in fig. 3. As A reaches critical value A c , the polaron turns from being quasi-stationary to moving by several hundred lattice sitesduring its lifetime. Evidently, the value of A c is independent of T . Note that we only considerthe displacement of polarons whose lifetimes are at least 2 T , and we set the displacement ofpolarons with shorter lifetimes to zero - see for instance the dashed (blue) lines in the centreand bottom-middle subfigures.Whilst the value of A c does not depend on T , the amount of displacement caused by A c does.However, it is unclear from our results what their correlation is. As A increases beyond A c ,the qualitative behaviour of D is that it decreases. This is due to the fact that increasing A causes the polaron to delocalise more quickly, and therefore the polaron has less time tomove. To understand how A affects the amount of polaron displacement per unit time , weexamine the polaron’s (average) velocity, V , as per definition (74). When A is small, V is zero.As A reaches critical value A c , the velocity becomes typically O (10 − ). Exactly what valuethis critical velocity V c takes depends on ¯ (cid:15) - the larger ¯ (cid:15) is, the larger V c is. As A increasesbeyond A c , sometimes V simply decays away - see for instance the top-right subfigure, where¯ (cid:15) = 0 . V grows before its decay - see for instance the middle-right20 igure 7: (Colour online.) Polaron lifetime τ d , displacement D , and velocity V ,under the MSPF, (cid:15) = ¯ (cid:15) + A sin(2 πτ /T ), with symmetry parameter fixed at β = 0 . A . Top row: ¯ (cid:15) = 0 . (cid:15) = 0 .
03. Bottom row: ¯ (cid:15) = 0 .
1. Dotted (red) lines: T = 100. Solid(black) lines: T = 500. Dashed (blue) lines: T = 2000. and bottom-right subfigures, where ¯ (cid:15) = 0 .
03 and 0.1 respectively. Such behaviour is possiblewhen the polaron lifetime decays with A more quickly than the displacement does. Whenthis happens, there may exist some optimal amplitude , A = A m , at which the polaron attainsmaximum velocity, V m . A m may coincide with A c - see for instance the top-right subfigure.Meanwhile, the middle-right and bottom-right illustrate clearly that, for different values of T ,the critical A c remains the same, whereas the optimal A m changes. Qualitatively speaking, thelarger T is, the smaller A m is.Whilst the value of A c does not depend on T , it does depend on ¯ (cid:15) - we see this by comparingany row of subfigures in fig. 7 to any other row. But how do A c and ¯ (cid:15) correlate? Our resultsshow that, as ¯ (cid:15) grows, A c drops, but crucially the combined magnitude (cid:15) comb = ¯ (cid:15) + A c remainsroughly constant. Specifically, in the top row we see ¯ (cid:15) = 0 .
005 and A c = 0 . (cid:15) = 0 .
030 and A c = 0 . (cid:15) = 0 .
100 and A c = 0 . (cid:15) comb = 0 .
172 when A reaches critical. Recall that, when using astraightforward constant forcing (cid:15) = ¯ (cid:15) , there is also a threshold value for ¯ (cid:15) , below which the21olaron simply exhibits small oscillations, and above which the polaron moves at high speedbut delocalises very quickly. It is noteworthy that this threshold is (cid:15) = 0 .
154 (given β = 0 . (cid:15) = 0 .
154 causes polaron displacement, (cid:15) = 0 .
153 does not; and if one wishes to add on aperiodic component A sin(2 πτ /T ) in order to move the polaron, one needs A ≥ . (cid:15) + A far exceed what ¯ (cid:15) is required on its own to move the polaron. This phenomenon isobserved across all values of β .In practice, then, what would make a good combination of forcing parameters, which propelthe polaron with decent speed but does not cause large energy dissipation too quickly? First ofall, a large ¯ (cid:15) results in a large velocity but an energetically unstable polaron which delocalisesrapidly - so rapidly that it may move the polaron less far in its lifetime than a small ¯ (cid:15) does.The middle column of fig. 7 precisely illustrates this point. Meanwhile, a small ¯ (cid:15) results in long-living polarons which can move very far, because of how stable they are, but they would takemore time to reach the same destination, compared to polarons under a large ¯ (cid:15) . On balance,a moderate value of ¯ (cid:15) such as 0.03 is preferable. Secondly, once a ¯ (cid:15) is chosen, it remains tochoose A and T , and it is obvious that the ideal choice of A is the optimal amplitude, A = A m .Meanwhile, if T is small, such as T = 100, A m would be large. On the other hand, if T is large,such as T = 2000, the value of the maximum velocity would be small. We observe both of theseextremes very clearly in the middle-right subfigure of fig. 7. Once again, these observationsare not specific to β = 0 .
6, but universal for all values of β . Overall, we believe that the bestMSPF parameters which we have tested are such combinations where ¯ (cid:15) ≈ . T ≈ A ≈ A m which, given β = 0 .
6, is A m = 0 . A m and β in section 4.3.We note one anomaly which we observe in fig. 7 but did not expect. When ¯ (cid:15) = 0 .
005 (top row),if T = 2000 and A is large enough, then the displacement (and therefore velocity) can take largenegative values, meaning the polaron moves in the opposite direction to what we expected, andwith large speeds. Whilst we are uncertain as to what causes this counter displacement, it iscertainly another reason to reject small ¯ (cid:15) and large T when choosing forcing parameters. β To produce fig. 7, we fixed β = 0 .
6. How would the figure have looked if β had been different?Our results show that qualitatively it would exhibit the same behaviour, characterised bycritical amplitudes A c , and optimal amplitudes A m . Quantitatively, the values of A c and A m would change. It is therefore natural to investigate how they change with β . After all, ourgeneralisation to the Davydov-Scott model is manifest in the extra parameter β . Firstly weestablish the following preliminary result.Recall that the stationary polaron, upon which we impose the MSPF, is characterised bytwo quantities: its probability distribution, specifically its maximum localisation probability,max | ψ n | , and its binding energy. These are in turn determined by the symmetry parameter β and effective coupling parameter λ . As β varies, so does the value of λ required to keepmax | ψ n | constant. This correlation is shown in fig. 8. It is clear that λ ( β ) is a decreasingfunction. We have made sure that whenever we altered β we also took λ = λ ( β ), so that all ofour moving polarons begin as stationary states which share the same probability distribution.An alternative would have been to take whatever λ is required to keep the binding energyconstant. Our results show that if we had decided to keep the binding energy constant at, say, − .
5, then max | ψ n | would have been 0.53 at β = 0, or 0.81 at β = 1. It is a central feature ofour model that two stationary polarons with the same maximum localisation probability need22ot have the same binding energy, and vice versa. Figure 8: For β in [0 , λ ( β ), the value of λ required to keepmax | ψ n | = 0 .
64; left axis: bindingenergy of stationary polaron resultingfrom β and λ ( β ). For example, when β = 0 . λ = λ (0 .
6) = 5 .
0, thebinding energy is − .
47. Figure 9: Critical amplitude A c (rightaxis, solid line), optimal amplitude A m (right axis, dashed line), criticalvelocity V c (left axis, solid line), andoptimal velocity V m (left axis, dashedline), as functions of β . Parameters:¯ (cid:15) = 0 . , T = 500. Having established the relation λ ( β ), we can study how A c and A m depend on β . We do soby fixing ¯ (cid:15) and T , and working out what A c and A m are for various { β, λ ( β ) } . For instance,in fig. 7 we saw that if β = 0 . λ = λ (0 .
6) = 5 .
0, ¯ (cid:15) = 0 .
03 and T = 500, then A c = 0 .
142 and A m = 0 . (cid:15) and T , and vary β ? The result is displayed in fig. 9. When β = 0,we have A c = 0 .
12, and when β = 1 we have A c = 0 .
1. In fact, A c is minimal when β = 1, whichsuggests that a system with antisymmetric electron-phonon interaction is most conducive topolaron displacement by MSPF. One might expect that a system with symmetric interactionwould be least conducive to polaron displacement, and therefore A c should be maximal when β = 0. This is not the case. We observe that A c is maximal when β = 0 .
6, which models asystem with moderately asymmetric electron-phonon interaction. Meanwhile, polaron velocityproduced by critical forcing, V c , is maximal when β = 0 . A m , varies little when β is less than 0.7, but decays sharply when β increases beyond 0.7, to such an extent that it almost equals the critical amplitude A c . Theoptimal velocity, V m , is typically of the same order of magnitude as V c . When β = 0 . V m and V c almost coincide.Earlier, based on fig. 7, we asserted that the period T of the MSPF has little effect on the valueof A c . Indeed, our results show that, if we had produced fig. 9 with T fixed at either 100 or2000, the A c curve would have been virtually unaffected. We also conjectured that the criticalamplitude A c is negatively and linearly correlated with ¯ (cid:15) , so that the combined magnitude (cid:15) comb = ¯ (cid:15) + A c remains constant as ¯ (cid:15) varies. This is supported by our results. Indeed, fig. 9is produced with ¯ (cid:15) fixed at 0.03; but if we had produced fig. 9 with ¯ (cid:15) fixed at either 0.005 or0.1, the A c curve would simply have been shifted along the vertical axis, by an amount equalto the difference between the new ¯ (cid:15) and 0.03. 23 Dynamical polarons in non-zero temperature
We study the effect of random fluctuations which result from non-zero temperatures in theenvironment surrounding the lattice. The randomised forcing on the lattice is represented bythe normally-distributed f n ( τ ), which by definition (21) must have, for τ ≥
0, the followingfirst and second moments. (cid:104) f n ( τ ) (cid:105) = 0 , (75a) (cid:104) f m ( τ ) f n ( τ + ∆ τ ) (cid:105) = 2 γθδ mn ∆ τ , (75b)where θ is the dimensionless temperature, θ = k B Θ M R Ω , (76)and Θ is the temperature. Compared to the dimensional F n ( t ) that we introduced in section 2,the τ which now appears in f n ( τ ) is a discrete index. We have followed the standard procedureof replacing δ ( t − t (cid:48) ) by 1 / ∆ τ , up to non-dimensionalisation constants. Beginning with astationary polaron, which we computed numerically in section 3.2, we integrate the systemof eq. (19) forward in time from τ = 0. Using a random number generation algorithm, wegenerate a new vector f n before each integration step. If θ is large, we find that it can causelarge distortions in the lattice and rapid delocalisation of the polaron, due to excessive energyinput to the system. For appropriate values of θ , we see that the polaron’s binding energyundergoes small fluctuations, but on the average it tends to shift towards zero. After sometime, the binding energy stabilises. For example, given β = 0 . λ = 5 .
0, the stationarypolaron has binding energy − .
47. Integrating from τ = 0 with θ = 0 . τ ∼ O (10 ) the binding energy settles, on average, around − .
15. The period of time requiredfor a polaron to reach such a thermal equilibrium is the thermalisation phase of the polarondynamics. During this phase, the forcing (cid:15) ( τ ) on the electron is kept at zero, but the electronnevertheless undergoes small fluctuations around its initial position, due to its coupling to thethermalised lattice. Our results show that, irrespective of β , we are unable to raise θ > . θ induces excessive lattice distortions which cause the polaron to delocalisebefore reaching thermal equilibrium.In each simulation, we integrate the system, with (cid:15) ( τ ) = 0, until the polaron reaches thermalequilibrium. Then we reset τ = 0, and “turn on” the forcing (cid:15) ( τ ) for τ ≥
0. We examine how thepolaron subsequently moves, under combinations of (cid:15) ( τ ) and f n ( τ ). Since the thermalisationphase raises the polaron energy, we expect that a thermalised polaron would be easier todisplace, in the sense that it would require a smaller (cid:15) ( τ ) to displace it, compared to the zerotemperature case. Indeed, our results confirm this. To obtain our results in this section, everydynamical simulation, with a set of chosen parameters { β, λ ( β ) , ¯ (cid:15), A, T, θ } , is run 100 times,and averages of quantities such as polaron lifetime and displacement are then taken.Figure 10 is to be compared directly with fig. 7, which contained results for β = 0 . θ = 0.Specifically, fig. 10 is to be compared with the solid (black) lines in the middle row of subfiguresin fig. 7, for which two of the parameters in (cid:15) = ¯ (cid:15) + A sin(2 πτ /T ) were fixed: ¯ (cid:15) = 0 .
03, and T = 500. We saw that, given said parameter values, the critical amplitude was A c = 0 . θ in the system, we define A c to be the smallest A for which theaverage polaron displacement (over 100 simulations) exceeds 10 lattice sites. According to thisdefinition, when β = 0 . , ¯ (cid:15) = 0 . , T = 500 and θ = 0 . A c = 0 . θ = 0. Fixing said values of ¯ (cid:15), T and θ , we find thatthe value of A c depends on β in a manner shown in fig. 11. That is, A c is minimal when β = 1,24 igure 10: (Colour online.) From left to right: polaron lifetime τ d , displacement D , and velocity V , under the MSPF, (cid:15) = ¯ (cid:15) + A sin(2 πτ /T ), and the thermalforcing, f n ( τ ) with temperature θ . The horizontal axis is A . The symmetryparameter fixed at β = 0 .
6. ¯ (cid:15) = 0 .
03 and T = 500 are fixed. Black lines: θ = 0 . θ = 0 . A c (right axis) and critical velocity V c (left axis), as functions of β . Param-eters: ¯ (cid:15) = 0 . , T = 500 , θ = 0 . θ c asa function of β and ¯ (cid:15) , A = 0. suggesting that an antisymmetric electron-phonon interaction makes it easiest to displace thepolaron. Meanwhile, A c is maximal when β ≈ .
5, suggesting that, counter-intuitively, whatmakes displacing the polaron most difficult is not a symmetric electron-phonon interaction,but a moderately asymmetric one. Indeed, this A c ( β ) function is very similar to the one infig. 9, where we also had ¯ (cid:15) = 0 . , T = 500 fixed, but θ = 0. Now with θ = 0 . A c ( β )curve in fig. 11 is significantly lower. This means that, regardless of β , a non-zero temperaturemakes it easier to displace the polaron by the forcing (cid:15) = ¯ (cid:15) + A sin(2 πτ /T ), in the sense thata smaller combined magnitude ¯ (cid:15) + A is required. It is also noteworthy that, under a non-zerotemperature, the onset of polaron motion is more gradual, in the sense that a critical amplituderesults in a very small velocity, V c . Indeed, comparing V c in fig. 9 with V c in fig. 11, we seethat the latter is 2 orders of magnitude smaller.25ore can be said about fig. 10. When we raise the temperature to θ = 0 . A = 0. This suggests that,given β = 0 . (cid:15) = 0 .
03, there exists some critical temperature θ = θ c between 0.0001and 0.0005, for which just the combination of (cid:15) ( τ ) = ¯ (cid:15) and f n ( τ ) is sufficient to displace thepolaron, and no periodic component in (cid:15) ( τ ) is needed. θ c is critical in the sense that, if θ isany lower than θ c , then the combination of (cid:15) ( τ ) = ¯ (cid:15) and f n ( τ ) does not energise the polaronenough to move it, and a non-zero A is required. Indeed our results show that, given β = 0 . (cid:15) = 0 .
03, the critical temperature is θ c = 0 . θ c changes as we vary β and ¯ (cid:15) , and the results are shown in fig. 12. We observe thequalitative trend that, the larger ¯ (cid:15) is, the less thermal energy is required to make up for theextra energy that the polaron needs in order to move. We also observe that, in general, thelarger β is, the less thermal energy is required to displace the polaron. This fits in nicely withour understanding that, when β is close to 1, we have an electron-phonon interaction which isbiased towards one end of the lattice, making the electron more susceptible to displacement. Figure 13: A polaron trajectory (right axis), and thecorresponding time-evolution of the polaron’s bind-ing energy (left axis), given β = 0 . , ¯ (cid:15) = 0 . , A = 0,and θ = θ c = 0 . In fig. 13 we present a typical po-laron trajectory when θ = θ c . Underthis critical temperature, some simu-lations would produce no polaron dis-placement at all, but most trajecto-ries would be similar to that in fig. 13,clearly showing a directed movement.We propose to explain the shape ofthese trajectories as follows. First ofall, the combined magnitude of theMSPF would be much lower than whatis required to move the polaron underzero temperature. In fig. 13 for exam-ple, we have ¯ (cid:15) + A = 0 .
03, whereasthe critical value under zero temper-ature, as we discovered in section 4.2,is ¯ (cid:15) + A = 0 . (cid:15) . But before the electron has time to movefar, the random forces may have further distorted the local lattice sites in such a way that ahigh potential barrier is restored. This then traps the electron again, giving the polaron timeto recover its integrity, before the next random time at which the electron jumps out of itspotential well. This explains why a trajectory under the critical temperature appears jagged,showing the polaron “hopping” one or two sites at a time, in stark contrast with the smoothand regular polaron motion exhibited in fig. 6.26 Discussions and conclusions
In this study we have presented a new mathematical model describing polaron dynamics inlinear peptide chains. The model is dependent on a symmetry parameter, β , which measuresthe extent to which the interaction between the polaron’s electron and phonon componentsis spatially symmetric. We have shown that when β takes its extreme values, 0 and 1, themodel reduces to existing ones for which it was assumed that the electron-phonon interactionwas, respectively, symmetric and antisymmetric. We have justified the physical neccessity ofincluding β in the model, in that one should not simply assume the electron to be coupledequally strongly to lattice points on either side, or to be coupled only to the lattice point onone side. Instead, the spatial symmetry should be determined by the adjustable parameter β .Apart from β , we have also identified two composite parameters which are most vital to theintrinsic properties of the polaron. Firstly there is the adiabaticity parameter, ρ , measuringthe characteristic time scale separation between the electron and phonon, which we justifiablyfixed throughout the study. Then there is the effective coupling parameter, λ , measuring thestrength of the electron-phonon interaction. The combination of β and λ determines the twoaspects of the stationary polaron: its maximum localisation probability, and its binding energy.We have computed both of these quantities as functions of β and λ . Moreover, in the infinitelattice limit, we have obtained stationary polaron solutions by analytically integrating thesystem, and the results are in good agreement with our numerical solutions on a finite lattice.Our main results relate to using an external forcing to displace the polaron, in a manner whichcauses minimal energy loss and which, crucially, is directed. Such polaron dynamics could beachieved only if the electron is dislodged from its self-trapping potential well, and the locallattice distortions propagate coherently with the electron, and some mechanism exists whichensures the electron always moves towards one end of the lattice. If the second condition is notmet, then over time the electron probability density function would become broader, leading todelocalisation of the polaron. We have found that a constant external force, ¯ (cid:15) , on the electronis insufficient for displacing the polaron, unless ¯ (cid:15) is larger than some threshold value, but thenthe forcing causes rapid energy loss and delocalisation. We have also found that a sinusoidalforce, A sin(2 πτ /T ), on the electron is never sufficient for displacing the polaron, throughoutthe range of A that we tested. We then combined the constant and sinusoidal forces, resultingin the mean-shifted periodic forcing (MSPF), (cid:15) = ¯ (cid:15) + A sin(2 πτ /T ). We have discovered that,for each ¯ (cid:15) which is insufficient on its own to displace the polaron, there is some critical value A c ,such that the polaron is displaced if and only if A ≥ A c . There is also an optimal value A m , suchthat the polaron attains maximum velocity at A = A m . The value of A c is irrespective of theperiod T , whilst A m is negatively correlated with T . As ¯ (cid:15) is decreased, A c increases, in such away that the combined magnitude ¯ (cid:15) + A c remains constant. This suggests that there is a certainamount of extra energy that the electron needs in order to overcome the polaron binding, andhow much of it comes from the constant or sinusoidal part is inconsequential, as long as thetwo parts combine to a large enough overall amplitude. Nevertheless, the split between ¯ (cid:15) and A does determine the manner in which the polaron propagates, specifically its velocity andstability. The velocity is predominantly determined by ¯ (cid:15) , and positively correlated with it; butthe stability of the polaron is negatively correlated with ¯ (cid:15) . By comparing three sets of { ¯ (cid:15), A } with the same combined ¯ (cid:15) + A , namely { . , . } , { . , . } , { . , . } (while keepingall other parameters fixed), we found that { . , . } produces optimal balance betweenpolaron velocity and stability.We have examined how the aforementioned phenomena depends upon β . To do so, we neededa way of isolating the effect of varying β . This posed a difficulty, because if we were to fix λ β then both the maximum localisation probability and binding energy of the station-ary polaron would change. We would then be comparing dynamical behaviours of dissimilarpolarons. We therefore decided to vary λ with β , in a way that allowed us to generate a setof stationary polarons, one for each combination of { β, λ ( β ) } , such that they all had the samemaximum localisation probability. Then we launched these polarons using the same externalforcing and compared the results. We have found that β = 1, representing a spatially antisym-metric electron-phonon interaction, produces a polaron which is easiest to move, in the sensethat the least amount of forcing is required. We have also found that the symmetric model, β = 0, does not make the polaron most difficult to move ( β ≈ . β = 0 model which pushes the electron towardsone end of the lattice, despite it being coupled to the other end equally strongly.We have also studied the MSPF under non-zero temperatures, θ >
0. The manifestationof thermal effects is random forces on the lattice points. We have found that a non-zero θ facilitates polaron propagation, in the sense that it lowers the critical amplitude A c , for anygiven ¯ (cid:15) . Moreover, a non-zero θ results in a gradual onset of polaron motion, meaning therate of change of polaron velocity with respect to A near A = A c is small, compared to theonset under θ = 0. Our results have also shown that, whenever there is polaron propagation,whether θ = 0 or θ >
0, the relative displacements between neighbouring lattice points remainunder O (10 − ). This is a necessary condition which allows us to model the lattice points aspoint dipoles.Some of the choices of parameters in the MSPF may be justified physically as follows. It iswell known that across the plasma membrane of a living cell, a resting membrane potential ismaintained by intercellular chemical processes [28]. It is also well known that within the plasmamembrane there exist highly stable transmembrane regions of proteins, for instance the humanprolactin receptor 2N7I [29], and the rat monoamine oxidase A 1O5W [30], both of which are α -helical structures spanning the entire membrane width. Given a constant potential differenceof ∆ V across a linear, homogeneous, isotropic dielectric medium with constant width d , theeffective electric field inside the medium is given by E = ∆ V / ( κd ), where κ is the dielectricconstant (a.k.a. relative permittivity) of the medium [31]. Assuming the plasma membrane issuch a medium, we can then attribute the physical origin of ¯ (cid:15) to the resting membrane potential,and calculate ¯ (cid:15) using eq. (21), namely ¯ (cid:15) = qE R/ ( (cid:126) Ω). For many cells the values of the E and d are well established. As an example, one may look at human red blood cells (erythrocytes),one of the most widely studied cells in nature, and point to [32] for the value E = − . d = 78˚A. However, the value of κ for a membrane is highly contentious,due to the fact that it depends sensitively upon a large variety of biophysical attributes ofthe membrane, such as hydration [35], pH value [36], and structural stability [37]. In a recentreview, it was reported that the value of κ in literature ranges from 1 to 40 [38]. Feeding thesevalues of E , d and κ into the equation for ¯ (cid:15) , we find that ¯ (cid:15) ranges from 0.0033 to 0.13. Wehave taken care to ensure that in this study the values of ¯ (cid:15) falls strictly within this range. For aphysical origin of the periodic term, A sin(2 πτ /T ), one could look to common electromagneticradiations which fill the environment around us in the modern age, such as the radiation fromtelecommunication transmitters. In particular, the values of T which we have considered, 100,500 and 2000, respectively match the frequencies of the IEEE 802.11ad protocal Wi-Fi band,the K u band frequencies for satellite communications and broadcasting, and the UHF bandfrequencies for cellular communications [39,40]. However, the amplitudes of the aforementionedradiations are much smaller than the values of A for which we have observed noteworthy results.For instance, treating the mobile telephone transmitter as an omni-directional dipole with peakpower P , we can estimate the amplitude ˜ A of its output waves at operational distance d , byusing the well-known formula P/ (4 πd ) = (cid:15) c ˜ A /κ , where (cid:15) , c, κ are the vacuum permittivity,28peed of light, and relative permittivity of the medium, respectively. Feeding P = 1W [41] and κ ≤
40 [38] into the formula, we obtain dimensionless A ≤ × − / ( d/ metres). This meansthat in order to obtain A = 0 .
1, one needs the operational distance d to be O (10 − ) metres,which is unrealistic. It is therefore clear that, in a real cell environment, the effect on a polarondue to a combination of resting membrane potential and random thermal forces are dominantover any external electromagnetic radiation that may commonly be present. This highlights theimportance of our observation that, given a constant electric field ¯ (cid:15) , there exists some criticaltemperature θ c , such that the polaron undergoes directed drift if and only if θ ≥ θ c . In otherwords, a combination of ¯ (cid:15) and θ can be sufficient for displacing the polaron, in a manner whichcauses minimal energy loss and which is directed, just as a combination of ¯ (cid:15) and A sin(2 πτ /T )can. It has been reported that stationary polarons formed on 3-dimensional lattices, such asan α -helix, are more strongly bound and therefore can be stabler during propagation [42]. Itis our hope that our model can be adapted to study such 3-D systems, and that the stabilisingeffect of the helical geometry will enable us to raise θ , from our current values of O (10 )K tophysiological temperatures. Acknowledgement
We would both like to thank Larissa Brizhik for kindly answering some questions.
References [1] L. D. Landau.
Phys. Z. Sowjet. , 3:664, 1933.[2] H. Fr¨ohlich.
Proc. R. Soc. Lond. A , 215:291, 1952.[3] T. Holstein.
Ann. Phys. , 8:325, 1959.[4] T. Holstein.
Ann. Phys. , 8:343, 1959.[5] A. J. Heeger, S. Kivelson, J. R. Schrieffer, and W.-P. Su.
Rev. Mod. Phys. , 60:781, 1988.[6] N. K. Voulgarakis and G. P. Tsironis.
Phys. Rev. B , 63:014302, 2000.[7] L. S. Brizhik, A. A. Eremko, B. M. A. G. Piette, and W. J. Zakrzewski.
Phys. Rev. B ,68:104301, 2003.[8] A. S. Davydov.
Phys. Scripta. , 20:387, 1979.[9] A. C. Scott.
Phys. Rep. , 217:1, 1992.[10] G. N. Chuev and V. D. Lakhno.
J. Theor. Biol. , 163:51, 1993.[11] E. M. Conwell and S. V. Rakhmanova.
Proc. Natl. Acad. Sci. USA , 97:4556, 2000.[12] L. S. Brizhik, A. A. Eremko, B. M. A. G. Piette, and W. J. Zakrzewski.
J. Phys. Condens.Matter , 20:255242, 2008.[13] L. S. Brizhik, A. A. Eremko, B. M. A. G. Piette, and W. J. Zakrzewski.
J. Phys. Condens.Matter , 22:155105, 2010.[14] L. S. Brizhik, A. A. Eremko, B. M. A. G. Piette, and W. J. Zakrzewski.
Phys. Rev. E ,89:062905, 2014. 2915] A. S. Davydov.
Sov. Phys. Usp. , 25:898, 1982.[16] A. S. Davydov.
Solitons in Molecular Systems . Kluwer Academic Publishers, 2nd edition,1991.[17] J. C. Eilbeck, P. S. Lomdahl, and A. C. Scott.
Phys. Rev. B , 30:4703, 1984.[18] A. L. Lehninger, D. L. Nelson, and M. M. Cox.
Principles of Biochemistry . WorthPublishers, New York, 2nd edition, 1993.[19] W. Barford.
J. Chem. Phys. , 126:134905, 2007.[20] D. N. Langelaan, M. Wieczorek, C. Blouin, and J. K. Rainey.
J. Chem. Inf. Model. ,50:2213, 2010.[21] K. G. Brown, S. C. Erfurth, E. W. Small, and W. L. Peticolas.
Proc. Natl. Acad. Sci.USA , 69:1467, 1972.[22] K.-C. Chou.
Biochem. J. , 209:573, 1983.[23] K.-C. Chou.
Biophys. J. , 45:881, 1984.[24] D. S. Lemons and A. Gythiel.
Am. J. Phys. , 65:1079, 1997.[25] T. Schlick.
Molecular Modeling and Simulation . Springer, 2nd edition, 2010.[26] D Hennig.
Phys. Rev. E , 64:041908, 2001.[27] G. Kalosakas, S. Aubry, and G. P. Tsironis.
Phys. Rev. B , 58:3094, 1998.[28] M. Luckey.
Membrane Structrural Biology . CUP, 2nd edition, 2014.[29] K. Bugge, E. Papaleo, G. W. Haxholm, J. T. S. Hopper, C. V. Robinson, J. G. Olsen,K. Lindorff-Larsen, and B. B. Kragelund.
Nat. Commun. , 7:11578, 2016.[30] J. Ma, M. Yoshimura, E. Yamashita, A. Nakagawa, A. Ito, and T. Tsukihara.
J. Mol.Biol. , 338:103, 2004.[31] Jackson J. D.
Classical Electrodynamics . Wiley, New York, 3rd edition, 1999.[32] K. Cheng, H. C. Haspel, M. L. Vallano, B. Osotimehin, and M. Sonenberg.
J. MembraneBiol. , 56:191, 1980.[33] L. McCaughan and S. Krimm.
Science , 207:1481, 1980.[34] R. M. Hochmuth, C. A. Evans, H. C. Wiles, and J. T. McCown.
Science , 220:101, 1983.[35] D. L. Mobley, K. A. Dill, and J. D. Chodera.
J. Phys. Chem. B , 112:938, 2008.[36] V. Z. Spassov and L. Yan.
Protein Science , 17:1955, 2008.[37] S. Vicatos, M. Roca, and A. Warshel.
Proteins , 77:670, 2009.[38] L. Li, C. Li, Z. Zhang, and E. Alexov.
J. Chem. Theory Comput. , 9:2126, 2013.[39] IEEE Std 802.11ad-2012 (Amendment to IEEE Std 802.11-2012, as amended by IEEE Std802.11ae-2012 and IEEE Std 802.11aa-2012), 2012. doi: 10.1109/IEEESTD.2012.6392842.[40] IEEE Std 521-2002 (Revision of IEEE Std 521-1984), 2003. doi:10.1109/IEEESTD.2003.94224.[41] S. L¨onn, U. Forss´en, P. Vecchia, A. Ahlbom, and M. Feychting.
Occup. Environ. Med. ,61:769, 2004. 3042] L. S. Brizhik, A. A. Eremko, B. M. A. G. Piette, and W. J. Zakrzewski.