A generalised Davydov-Scott model for polarons in linear peptide chains

TL;DR

This paper introduces a flexible mathematical model for polarons in linear peptide chains, unifying previous models and exploring their dynamics under various external influences with potential biological implications.

Contribution

It develops a one-parameter family of models that generalize Davydov and Scott models, providing analytical and numerical solutions for polaron behavior.

Findings

Electric fields can induce low-loss polaron motion.

Thermal fluctuations facilitate polaron mobility.

Analytical stationary solutions derived in the continuum limit.

Abstract

We present a one-parameter family of mathematical models describing the dynamics of polarons in linear periodic structures such as polypeptides. By tuning the parameter, we are able to recover the Davydov and the Scott models. We describe the physical significance of this parameter. In the continuum limit, we derive analytical solutions which represent stationary polarons. On a discrete lattice, we compute stationary polaron solutions numerically. We investigate polaron propagation induced by several external forcing mechanisms. We show that an electric field consisting of a constant and a periodic component can induce polaron motion with minimal energy loss. We also show that thermal fluctuations can facilitate the onset of polaron motion. Finally, we discuss the bio-physical implications of our results.

Equations

$$\hat{H}={\hat{H}}_{\text{e}}+{\hat{H}}_{\text{p}}+{\hat{H}}_{\text{int}}+{\hat{H}}_{\text{ext}},$$

$${\hat{H}}_{\text{e}}=\sum _{n=0}^N{J}_0{\hat{A}}_n^{†}{\hat{A}}_n-\sum _{n=0}^{N-1}{J}_1\left({\hat{A}}_{n+1}^{†}{\hat{A}}_n+{\hat{A}}_n^{†}{\hat{A}}_{n+1}\right).$$

$${\hat{H}}_{\text{ext}}=-\sum _{n=0}^{N}qE(t)R\left(n-n_0\right){\hat{A}}_n^{†}{\hat{A}}_n,$$

$${\hat{H}}_{\text{p}}=\sum _{n=0}^N\frac{P_n^2}{2M}+\sum _{n=0}^{N-1}M{\Omega }^2\text{leavevmode}\frac{{\left(U_{n+1}-U_n\right)}^2}{2},$$

$$S_n\text{mathchar}58=U_{n+1}-U_n,$$

$${\hat{H}}_{\text{int}}^{\text{Dav}}=(\chi /2)[(U_1-U_0){\hat{A}}_0^{†}{\hat{A}}_0+\sum _{n=1}^{N-1}(U_{n+1}-U_{n-1}){\hat{A}}_n^{†}{\hat{A}}_n+(U_N-U_{N-1}){\hat{A}}_N^{†}{\hat{A}}_N].$$

$${\hat{H}}_{\text{int}}^{\text{Sco}}=\sum _{n=0}^{N-1}\chi (U_{n+1}-U_n){\hat{A}}_n^{†}{\hat{A}}_n.$$

$$\displaystyle\hat{H}_{\textnormal{int}}=\leavevmode\nobreak\$$

$$+\sum _{n=1}^{N-1}[{\chi }_r(U_{n+1}-U_n)+{\chi }_l(U_n-U_{n-1})]{\hat{A}}_n^{†}{\hat{A}}_n.$$

$$\chi \text{mathchar}58={\chi }_r+{\chi }_l,\beta =\frac{{\chi }_r-{\chi }_l}{{\chi }_r+{\chi }_l}\in \left[0,1\right],$$

$$\displaystyle\hat{H}_{\textnormal{int}}=\leavevmode\nobreak\$$

$$+\sum _{n=1}^{N-1}\frac{\chi }{2}\left[\left(U_{n+1}-U_{n-1}\right)+\beta \left(U_{n+1}+U_{n-1}-2U_n\right)\right]{\hat{A}}_n^{†}{\hat{A}}_n.$$

$$|\Psi (t)⟩=\sum _{n=0}^N{\alpha }_n(t){\hat{A}}_n^{†}|\text{vac}⟩,$$

$$\sum _{n=0}^N{\left|{\alpha }_n\right|}^2=1.$$

$$i\hslash \frac{\text{d}{\alpha }_n}{\text{d}t}=$$

$$-eE(t)R\left(n-n_0\right){\alpha }_n.$$

$$S_{-1}=S_N=0,{\alpha }_{-1}={\alpha }_{N+1}=0,$$

$$M\frac{{\text{d}}^2{U}_n}{\text{d}{t}^2}=$$

$$-\frac{\chi }{2}\beta \left[\left({\left|{\alpha }_{n+1}\right|}^2-{\left|{\alpha }_n\right|}^2\right)-\left({\left|{\alpha }_n\right|}^2-{\left|{\alpha }_{n-1}\right|}^2\right)\right].$$

$${\alpha }_0={\alpha }_N=0.$$

$$U_0=\frac{\text{d}{U}_0}{\text{d}t}=0.$$

$${\alpha }_n(t)={\psi }_n(t)\exp \left[-\frac{it}{\hslash }\left(J_0-2J_1\right)\right],$$

$$i{\dot{\psi }}_n$$

$${\ddot{u}}_n$$

$$s_n\text{mathchar}58=u_{n+1}-u_n,g_n\text{mathchar}58={\left|{\psi }_{n+1}\right|}^2-{\left|{\psi }_n\right|}^2,c_n\text{mathchar}58={\left|{\psi }_{n+1}\right|}^2+{\left|{\psi }_n\right|}^2,$$

$$\rho =\frac{J_1}{\hslash \Omega },\sigma =\frac{R\chi }{2\hslash \Omega },ϵ=\frac{qER}{\hslash \Omega },\delta =\frac{\chi }{2MR{\Omega }^2},\gamma =\frac{\Gamma }{M\Omega },f_n=\frac{F_n}{MR{\Omega }^2}.$$

$$\sum _{n=0}^N{\left|{\psi }_n\right|}^2=1.$$

$$s_n-s_{n-1}=\delta \left[\beta \left(g_n-g_{n-1}\right)-\left(c_n-c_{n-1}\right)\right],$$

$$s_n\equiv u_{n+1}-u_n=\delta \left(\beta g_n-c_n\right)=\delta \left[\left(\beta -1\right){\left|{\psi }_{n+1}\right|}^2-\left(\beta +1\right){\left|{\psi }_n\right|}^2\right].$$

$$i{\dot{\psi }}_n=$$

$$-\rho \left({\psi }_{n+1}+{\psi }_{n-1}-2{\psi }_n\right).$$

$$\Delta {\psi }_n$$

$$\lambda$$

$$i{\rho }^{-1}{\dot{\psi }}_n+\Delta {\psi }_n+\lambda {\left|{\psi }_n\right|}^2{\psi }_n+\eta \Delta {\left|{\psi }_n\right|}^2{\psi }_n=0.$$

Full text

A generalised Davydov-Scott model for polarons in linear peptide chains

April 2017

J. Luo*

B. M. A. G. Piette$\dagger$

*jingxi.luo@durham.ac.uk

$\dagger$b.m.a.g.piette@durham.ac.uk

Department of Mathematical Sciences, Durham University

Abstract

We present a one-parameter family of mathematical models describing the dynamics of polarons in 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 significance of this parameter and, in the continuum limit, we derive analytical solutions which represent stationary polarons. On a discrete lattice, we compute stationary polaron solutions numerically. We investigate polaron propagation induced by several external forcing mechanisms. We show that an electric field consisting of a constant and a periodic component can induce polaron motion with minimal energy loss. We also show that thermal fluctuations can facilitate the onset 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

1 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 of electron self-trapping. The presence of the electron causes localised distortions in the natural vibrational mode of the lattice, a.k.a. the lattice phonon. In return, if the electromagnetic interaction between the electron and lattice is appropriate, then the phonon distortions can lower the potential well for the electron, thus trapping the electron.

Some twenty years after the inception of the polaron concept, a mathematical description of it was formalised by H. Fröhlich [2] and subsequently T. Holstein [3, 4]. Since then, properties of the Fröhlich-Holstein polaron have been well studied, with some authors hypothesising an application of dynamical polarons as electron transporters in conductive material [5, 6, 7]. In the 1970s, A. S. Davydov used the basis of polaron theory to explain some biological processes [8]. Specifically he proposed that, in an $\alpha$-helical protein, a certain intramolecular oscillator can interact with the peptide chain in a way similar to an electron interacting with a crystal lattice. Davydov suggested that this interaction could lead to the localisation and propagation of vibrational energy in the $\alpha$-helix. Later, A. C. Scott modified Davydov’s theory, taking into account the internal geometry of peptide units [9]. Some authors argued that, given the polarisability of peptide units, electron self-trapping is also possible in proteins, and it can be described by the same mathematical model that Davydov and Scott used [10, 11]. It should therefore 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 in simple molecular chains, and adverted to the applicability of their results to electron transport in biomolecules such as proteins [12, 13, 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 it to 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-polypeptide interaction is spatially symmetric. In section 2, we describe our model and explain why the extra parameter is necessary. We also give physical interpretations of all other parameters in the model, justifying the choices of their values where possible. Then, we derive a set of coupled dynamical equations which govern the electron and phonon parts of the polaron, as well as how 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’s binding energy. The process of solving the equations is carried out analytically as well as numerically. By the former approach, a closed-form expression for the solution is found, but its use is limited, because the solution process involves a few approximations and simplifying assumptions. By the numerical approach, no convenient expression for the solution is possible, but the method solves the equations directly without simplifications. We compare the results produced by the two different methods.

Section 4 concerns dynamical polarons. We discover that it is possible to use a suitable external forcing to displace the stationary polaron, and to sustain its motion in such a way that its energy remains highly stable. We investigate how the polaron’s motion depends upon our forcing parameters. We use only numerical methods to obtain our results in section 4, as well as those in section 5, where we consider how the polaron’s motion is affected by temperature of the environment. For this part, the external forcing from section 4 remains in place, but we also utilise a parameter which controls the magnitude of the thermal effect. To account for the random 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 realisabililty of our mathematical model, particularly how the external forcing which we study in section 4 may be realised. We also briefly discuss the generalisability of our model to studying electron transport by polarons in $\alpha$-helices.

2 The model and dynamical equations

In both Davydov’s and Scott’s models, the Hamiltonian for a system of excitons interacting with one-dimensional lattice phonons is written in Fröhlich-Holstein form $\hat{H}=\hat{H}_{\textnormal{e}}+\hat{H}_{\textnormal{p}}+\hat{H}_{\textnormal{int}}$, where $\hat{H}_{\textnormal{e}},\hat{H}_{\textnormal{p}}$ and $\hat{H}_{\textnormal{int}}$ represent energy contributions from the exciton, phonon and interaction parts, respectively [2, 3, 15, 16, 9]. We adopt this Hamiltonian for our model, and following [12, 13, 14] we consider an additional external Hamiltonian, $\hat{H}_{\textnormal{ext}}$, so that our Hamiltonian takes the form

[TABLE]

where $\hat{H}_{\textnormal{e}}$ describes a tight-binding electron, the stretching and compressing of hydrogen bonds in the peptide chain are phonon oscillations described by $\hat{H}_{\textnormal{p}}$, $\hat{H}_{\textnormal{int}}$ accounts for the electron-phonon interaction, and $\hat{H}_{\textnormal{ext}}$ represents the effect of an external electric field. We assume that the peptide chain consists of $N+1$ identical units and $N$ identical hydrogen bonds. In the tight-binding approximation, we have

[TABLE]

The subscript $n$ in eq. 2 labels peptide units, which are the unit cells of our lattice. $\hat{A}^{\dagger}_{n}$ and $\hat{A}_{n}$ are local electron creation and annihilation operators, respectively. $J_{0}$ is the potential energy of a localised electron. Modelling each unit as a point-dipole, we assume the nearest-neighbour dipole interaction energy is a constant and write it as $-J_{1}$ [17, 18, 19]. The external Hamiltonian,

[TABLE]

models the effect of an electric field with strength $E(t)$ on the potential energy of a localised electron with charge $-q$. The potential energy due to $E(t)$ is set to zero at some arbitrary $n_{0}$, and $R=4.5\textnormal{\r{A}}$ is the equilbrium lattice spacing. Since the electron mass is several orders smaller than the mass of a peptide unit, we take a semi-classical approach where the phonon Hamiltonian, $\hat{H}_{\textnormal{p}}$, is a classical one. In the harmonic approximation, the hydrogen bonds are modelled as Hookean springs with force constant $K$, and therefore $\hat{H}_{\textnormal{p}}$ takes the form

[TABLE]

where $M=1.774\times 10^{-25}\textnormal{kg}$ is the average mass of a peptide unit in a membrane $\alpha$-helix [20], and we have defined $\Omega\mathrel{\mathop{\mathchar 58\relax}}=\sqrt{K/M}$. $U_{n}$ and $P_{n}$ are, respectively, the displacement and conjugate momentum of the $n\textsuperscript{th}$ unit. Thus, the first and second sums in the expression for $\hat{H}_{\textnormal{p}}$ represent, respectively, the kinetic and potential energies of the lattice. We take the value of $\Omega$ to be the natural angular frequency of slow phonons in an $\alpha$-helix, $\Omega=5.5\times 10^{12}\textnormal{s}^{-1}$ [21, 22, 23]. To derive the interaction Hamiltonian, $\hat{H}_{\textnormal{int}}$, Davydov and Scott assumed that the energy of an on-site excitation depends on lattice deformations in its vicinity. For us, the local deformations are

[TABLE]

namely the amount by which the lengths of hydrogen bonds deviate from equilibrium. By Davydov and Scott’s assumption, if we write the electron energy at site $n$ in a Taylor expansion, the first two terms are $J_{0}+\chi G_{n}(S_{n},S_{n-1})$, where $\chi$ is a constant, $G_{n}$ is a linear function, and $\mathinner{\!\left\lvert\chi G_{n}/J_{0}\right\rvert}\ll 1$. Then the interaction Hamiltonian is $\hat{H}_{\textnormal{int}}=\sum_{n=0}^{N}\chi G_{n}\hat{A}^{\dagger}_{n}\hat{A}_{n}$. Davydov assumed that $S_{n}$ and $S_{n-1}$ have equal influence on local excitation energies [15], so $G_{n}=\left(S_{n}+S_{n-1}\right)/2$, and $\hat{H}_{\textnormal{int}}$ is

[TABLE]

Davydov’s model is therefore spatially symmetric, since $\hat{A}^{\dagger}_{n}\hat{A}_{n}$ is coupled equally to $U_{n+1}$ and to $U_{n-1}$. Scott modified Davydov’s model by opting for the antisymmetric $G_{n}(S_{n},S_{n-1})=S_{n}$ instead [9]. The reason is, while both authors let an intra-peptide C=O oscillator take the role of the exciton, Davydov neglected the internal geometry of the peptide units, but Scott pointed out that every unit has its C=O pair immediately adjacent to the next hydrogen bond in the chain. This leads Scott to assume that $\hat{A}^{\dagger}_{n}\hat{A}_{n}$ is coupled to, without loss of generality, $S_{n}$, and not $S_{n-1}$. Scott therefore had

[TABLE]

Since we are modelling an electron as opposed to an intra-peptide oscillator, we cannot use Scott’s argument to justify assuming that $\hat{A}^{\dagger}_{n}\hat{A}_{n}$ is coupled to $U_{n+1}$ and not to $U_{n-1}$. Nor should 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-1})=\chi_{r}S_{n}+\chi_{l}S_{n-1}$, taking without loss of generality $\chi_{r}>0$ and $0\leq\chi_{l}\leq\chi_{r}$. Then,

[TABLE]

By defining

[TABLE]

we can write

[TABLE]

We treat $\chi$ as an adjustable parameter. Setting $\beta=0$ ($\chi_{l}=\chi_{r}$) gives us the symmetric model of Davydov as per eq. 6, whilst setting $\beta=1$ ($\chi_{l}=0$) produces the antisymmetric model of Scott as per eq. 7. The larger $\beta$ is, the less spatial symmetry our model possesses. Indeed, for $\beta\in[0,1)$, the ratio of $n$-$(n+1)$ coupling strength to $n$-$(n-1)$ coupling strength is given by $\chi_{r}/\chi_{l}=(1+\beta)/(1-\beta)$, and this ratio is strictly increasing with $\beta$.

We write the electronic state of the system as a linear superposition of local excitations [12],

[TABLE]

where $\ket{\textnormal{vac}}$ is the vacuum state, and $\alpha_{n}\in\mathbb{C}$ is the probability amplitude for an electron localised at the $n$th site, subject to the normalisation condition,

[TABLE]

We proceed to derive dynamical equations for $\alpha_{n}$ and $U_{n}$. By equating coefficients of $\hat{A}^{\dagger}_{n}\ket{\textnormal{vac}}$ on both sides of the Schrödinger equation, $i\hbar\leavevmode\nobreak\ \textnormal{d}\ket{\Psi}/\textnormal{d}t=(\hat{H}_{\textnormal{e}}+\hat{H}_{\textnormal{int}}+\hat{H}_{\textnormal{ext}})\ket{\Psi}$, we obtain

[TABLE]

We have defined

[TABLE]

so that eq. 13 holds for all $n$ including the boundary terms ($n=0,N$). Equations for $U_{n}$ are derived from classical Hamilton equations, $\textnormal{d}U_{n}/\textnormal{d}t=\partial H_{\textnormal{cla}}/\partial P_{n}$ and $\textnormal{d}P_{n}/\textnormal{d}t=-\partial H_{\textnormal{cla}}/\partial U_{n}$, where $H_{\textnormal{cla}}\mathrel{\mathop{\mathchar 58\relax}}=\braket{\Psi}{(\hat{H}_{\textnormal{p}}+\hat{H}_{\textnormal{int}})}{\Psi}$. These equations are

[TABLE]

In order that eq. 15 holds at the boundaries, we have set

[TABLE]

This boundary condition is justified because we expect the probability distribution $\mathinner{\!\left\lvert\alpha_{n}\right\rvert}^{2}$ to be highly localised with half-width of $\mathcal{O}(1)$, and because we will be working with large lattices with $N\gg 1$. We also impose the following boundary condition on $U_{n}$, representing a peptide chain which is fixed at one end.

[TABLE]

Next, we introduce the gauge transformation,

[TABLE]

which sets $J_{0}=2$ in eq. 13. Physically this represents a shift in the arbitrary reference value from which energy is measured. Combining the $2\alpha_{n}$ term with the $J_{1}$ term in eq. 13, we obtain the discrete Laplacian, $-J_{1}\left(\alpha_{n+1}+\alpha_{n-1}-2\alpha_{n}\right)$. Meanwhile, to account for the interaction between the peptide chain and its environment, we need to add Langevin terms to the r.h.s. of eq. 15 [13, 14, 24, 25]. They are, a damping term describing energy dissipation due to friction, $-\Gamma\leavevmode\nobreak\ \textnormal{d}U_{n}/\textnormal{d}t$, where $\Gamma$ 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 $\braket{F_{m}(t)F_{n}(t^{\prime})}=2\Gamma k_{B}\Theta\delta_{m,n}\delta(t-t^{\prime})$, where $k_{B}$ is the Boltzmann constant and $\Theta$ is the temperature of the environment.

Scaling time by $\Omega^{-1}$ and length by $R$ gives us the following non-dimensionalised dynamical equations for $\psi_{n}$ and $u_{n}\mathrel{\mathop{\mathchar 58\relax}}=U_{n}/R$, for $n=0,1,\dots,N$.

[TABLE]

where we have defined

[TABLE]

and where the overdot denotes differentiation with respect to dimensionless time $\tau$, and

[TABLE]

Equation 19 holds subject to the boundary conditions (14), (16) and (17), as well as the normalisation condition,

[TABLE]

It is easily verifiable that setting $\beta=0$ and $\beta=1$ in eq. 19 produces Davydov’s and Scott’s dynamical equations, respectively [15, 9]. $\rho$ is known as an adiabaticity parameter, as it is the ratio of the characteristic time scale of phonon vibrations to that of electronic phase variations [26]. We fix $J_{1}$, following [12, 13, 14] (which used a different scaling), at $\rho=2.1$. Moreover, since $M,R$ and $\Omega$ are fixed, the ratio $\sigma/\delta=1880$ is constant. The range of $\delta$ which we consider throughout this study correspond to $\chi\sim\mathcal{O}(10^{-11})$ Newtons, agreeing with [9]. Finally, we take $\gamma=0.05$, agreeing with [12, 13, 14] up to different scaling factors.

3 Stationary polaron solutions

We derive stationary polaron solutions to eq. 19, subject to zero electric field ($\epsilon=0$) and zero temperature ($f_{n}=0$). We consider analytical and numerical methods separately, and compare the results.

3.1 Analytical results

When $f_{n}=0$ and $\dot{u}_{n}=\ddot{u}_{n}=0$, eq. 19b becomes

[TABLE]

which holds if

[TABLE]

Putting eq. 24 into eq. 19a and requiring $\epsilon=0$ gives us

[TABLE]

Defining

[TABLE]

we can rewrite eq. 25 as

[TABLE]

We note that, since $M,\Omega$ and $J_{1}$ are all fixed, the parameter $\lambda$ inherits the adjustability of $\chi$. Now, in a stationary state, the time dependence of $\psi_{n}$ can be at most a variation of the phase factor. Following [15], we consider the ansatz

[TABLE]

where $\xi$ is a real, continuous variable in the domain $[-NR/2,NR/2]$, $\phi$ is a real, twice-differentiable function, and $H_{0}$ and $k$ are constants. In particular, $H_{0}$ is an energy eigenvalue, in the sense that

[TABLE]

In the limit $N\gg 1$, $R$ becomes small compared to the domain size, which enables us to invoke the continuum approximation,

[TABLE]

implying

[TABLE]

Putting eqs. 29, 30, 31 and 32 into eq. 28, then dividing the result by $\exp\left(i\rho H_{0}\tau+ik\xi\right)$ and retaining terms up to $\mathcal{O}(R^{2})$, we obtain

[TABLE]

The last term in eq. 33 is equivalent to $\eta R^{2}\leavevmode\nobreak\ \textnormal{d}^{2}(\phi(\xi))^{2}/\textnormal{d}\xi^{2}$. Equating imaginary parts of eq. 33 gives us $k=0$. After the scaling $x\mathrel{\mathop{\mathchar 58\relax}}=\xi/R$, the real part of eq. 33 becomes

[TABLE]

The subscript $x$ denotes differentiation with respect to $x$. We seek $\phi(x)$ which satisfies eq. 34 for all $x$, not just when $x=n-N/2$. Then, from such a $\phi(x)$ we will be able to recover the discrete solution $\psi_{n}(\tau)$ via $\xi=xR$ and eq. 29. Further to being globally defined, we require that $\phi(x)$ has vanishing derivatives at infinity, and satisfies the normalisation condition,

[TABLE]

If $\eta=0$ (i.e. $\beta=1$), then eq. 34 reduces to the nonlinear Schrödinger equation with a cubic nonlinearity, which has a well-known solution satisfying all the above constraints [9],

[TABLE]

Consider $\eta>0$ (i.e. $\beta<1$). Since $(\phi^{2})_{xx}\equiv 2\phi\phi_{xx}+2(\phi_{x})^{2}$, we rewrite eq. 34 as

[TABLE]

We dedicate the remainder of this section to analysing eq. 37. It is an autonomous equation for $\phi(x)$, which allows us to define $h(\phi)\mathrel{\mathop{\mathchar 58\relax}}=\phi_{x}$, and write

[TABLE]

We define $y(\phi)\mathrel{\mathop{\mathchar 58\relax}}=h^{2}$, so that $y_{\phi}=2hh_{\phi}$, and multiply eq. 37 by 2 to obtain

[TABLE]

The l.h.s. of eq. 39 is the total derivative of $(1+2\eta\phi^{2})y$ with respect to $\phi$. We therefore have

[TABLE]

The integration constant $C$ is determined by considering the limit $x\rightarrow\infty$, in which $\phi^{2}\rightarrow 0$ and $y\equiv(\phi_{x})^{2}\rightarrow 0$. We therefore have $C=0$. Now we note that, if $H_{0}\leq 0$, then the r.h.s. of eq. 40 is negative whenever $\phi\neq 0$, so it cannot equal the l.h.s. which is $(\phi_{x})^{2}$. Thus, if $H_{0}\leq 0$ then the only $\phi(x)$ satisfying eq. 40 is identically zero. We therefore require $H_{0}>0$. Multiplying eq. 40 by $4\phi^{2}$, we obtain

[TABLE]

We then define $\Phi\mathrel{\mathop{\mathchar 58\relax}}=\phi^{2}$, and eq. 41 becomes

[TABLE]

If eq. 42 has a solution which is globally non-negative and twice-differentiable, has vanishing derivatives at infinity, and satisfies

[TABLE]

then we claim that $\Phi(x)$ must attain its global upper bound of $2H_{0}/\lambda$ at some finite $x$, and that every local maximum of $\Phi$ must also be a global maximum. The proof of this claim is as follows. Since $\Phi$ cannot be identically zero, and since $\lim_{x\rightarrow\pm\infty}\Phi(x)=0$, $\Phi(x)$ must have at least one turning point, at some finite $x$ and non-zero $\Phi$. But we observe from eq. 42 that $\Phi_{x}$ vanishes if and only if $\Phi=0$ or $2H_{0}/\lambda$. Therefore, $\Phi(x)$ must attain its global upper bound of $2H_{0}/\lambda$ at least once, and no other local maximum value is possible. This concludes the proof. We further propose that wherever $\Phi(x)$ attains its maximum value, say at $x=x_{\textnormal{max}}$, the second derivative $\Phi_{xx}$ does not vanish there. The proof is as follows. On the one hand, we have $\Phi_{xx}\equiv 2\phi\phi_{xx}+2(\phi_{x})^{2}$; when $x=x_{\textnormal{max}}$, we also have $\phi_{x}=\Phi_{x}/(2\phi)=0$, therefore $\Phi_{xx}=2\phi\phi_{xx}$. On the other hand, eq. 34 is equivalent not only to eq. 37 but also to $\eta\phi\Phi_{xx}=H_{0}\phi-\phi_{xx}-\lambda\phi^{3}$. It follows that, at $x=x_{\textnormal{max}}$, we have $(1/(2\phi)+\eta\phi)\Phi_{xx}=H_{0}\phi-\lambda\phi^{3}$, and therefore $\Phi_{xx}=(2H_{0}\Phi-2\lambda\Phi^{2})/(1+2\eta\Phi)$. Since $\Phi(x_{\textnormal{max}})=2H_{0}/\lambda$, it follows that $\Phi_{xx}(x_{\textnormal{max}})=-4H_{0}^{2}/(\lambda+4\eta H_{0})<0$, as required. A corollary of this proposition is that there must exist some neighbourhood of $x_{\textnormal{max}}$ containing no maxima of $\Phi(x)$ other than $x_{\textnormal{max}}$ itself. Without loss of generality, let $x_{\textnormal{max}}=0$. Suppose the corollary is false, so that every neighbourhood of 0 contains some non-zero $x$ which is a maximum of $\Phi(x)$. Then, there must exist some sequence $x_{n}$ approaching 0 such that $\Phi(x)$ has a maximum at every $x_{n}$, with $\Phi(x_{n})=\Phi(0)$. But this leads to a contradiction. Indeed, for every $x_{n}$ we have the Taylor expansion $\Phi(x_{n})=\Phi(0)+\Phi_{xx}(0)x_{n}^{2}/2+\mathcal{O}(x_{n}^{3})$, where the first derivative is absent because $\Phi(x)$ has a maximum at 0. It then follows that $\Phi_{xx}(0)=\lim_{n\rightarrow\infty}2(\Phi(x_{n})-\Phi(0))/x_{n}^{2}=0$, which contradicts the previous proposition. Therefore the corollary is proven. Since 0 has a maxima-free neighbourhood, we say that $\Phi(x)$ has an isolated maximum at 0.

We note that we can indeed require that $x=0$ is a maximum of $\Phi(x)$, because eq. 42 is translationally invariant: if $\Phi(x)$ is a solution then so is $\Phi(x-c)$ for any constant $c$. We exploit this invariance, requiring that $\Phi(x)$ satisfies

[TABLE]

Now we claim that there exists $b>0$, which may be infinite, such that $\lim_{x\rightarrow b}\Phi(x)=0$, and $\Phi(x)\neq 0$ for all $x\in(0,b)$. The proof of this claim is as follows. If $\Phi(x)\neq 0$ for all $x\in(0,\infty)$, then we are done. If $\Phi(x)=0$ for some $x\in(0,\infty)$, then the set $\{x\in(0,\infty)\mathrel{\mathop{\mathchar 58\relax}}\Phi(x)=0\}$ must have a minimum. If it does not, then there would be a sequence $x_{n}>0$ such that, as $n\rightarrow\infty$, $x_{n}\rightarrow 0$ and $\Phi(x_{n})\rightarrow 0$; but this would contradict the continuity of $\Phi(x)$ at $x=0$. Thus, letting $b$ equal the least positive zero of $\Phi(x)$, then we are done.

Next, we propose that no other solution on $[0,b)$ exists, and the proof is as follows. If $\Phi(x)$ has any maxima in $(0,\infty)$, then the set $\mathcal{M}\mathrel{\mathop{\mathchar 58\relax}}=\{x\in(0,\infty)\mathrel{\mathop{\mathchar 58\relax}}x\leavevmode\nobreak\ \textnormal{is a maximum of}\leavevmode\nobreak\ \Phi(x)\}$ must have a minimum, because otherwise we would have a contradiction to the fact that $x=0$ is an isolated maximum of $\Phi(x)$. Let $x_{1}=\min\mathcal{M}$, and suppose $x_{1}<b$. Since $\Phi(x_{1})=\Phi(0)$, and since $\Phi(x)$ has no maximum in the interval $(0,x_{1})$, and since $\Phi(x)$ is continuous, it must attain its minimum value at some point $b^{\prime}\in(0,x_{1})$. But that implies $\Phi(b^{\prime})=0$, where $b^{\prime}<x_{1}<b$, contradicting the fact that $\Phi(x)\neq 0$ for all $x\in(0,b)$. Therefore, we must have $x_{1}>b$. Since $\Phi_{x}$ vanishes only at maxima and minima, it follows that $\Phi_{x}$ is non-vanishing on $(0,b)$, and therefore $\Phi(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

[TABLE]

On any closed interval $[\Phi_{1},\Phi_{2}]\subset(0,\Phi_{0})$, the function $g(\Phi)$ is continuous and non-zero, so the reciprocal function $1/g(\Phi)$ is continuous and bounded, and therefore Riemann integrable. But $g(\Phi)$ approaches 0 as $\Phi\rightarrow\Phi_{0}$, meaning $1/g(\Phi)$ becomes unbounded. Thus, integration of $1/g(\Phi)$ on the interval $[\Phi_{1},\Phi_{0}]$ is not trivial. Luckily $\Phi=\Phi_{0}$ is an integrable singularity of the function $1/g(\Phi)$, because the Puiseux series of $1/g(\Phi)$ about $\Phi_{0}$ is $\mathcal{O}((\Phi-\Phi_{0})^{-1/2})$. Therefore, for any $\Phi_{1}\in(0,\Phi_{0}]$, eq. 45 is equivalent to

[TABLE]

which determines a unique $x(\Phi_{1})\in[0,b)$. The l.h.s. of eq. 46 is a strictly decreasing function of $\Phi_{1}$, meaning $x(\Phi_{1})$ has a unique inverse function which is also strictly decreasing, $\Phi_{1}(x)$, on the domain $x\in[0,b)$. But $\Phi(x)$ is an existing function satisfying eq. 45 for all $x\in[0,b)$. Therefore, we must have $\Phi_{1}(x)=\Phi(x)$ for all $x\in[0,b)$, and the uniquess of $\Phi(x)$ follows.

In summary, we have so far established the following. If eq. 42 has a solution which is globally non-negative and twice-differentiable, has vanishing derivatives at infinity, and has the property that its integral over $\mathbb{R}$ is 1, then eq. 42 has a solution, say $\Phi(x)$, which satisfies all the above constraints as well as the condition (44), and there exists some $b>0$ which may be infinite such that $\Phi(x)$ is strictly decreasing on $[0,b)$, and $\lim_{x\rightarrow b^{-}}\Phi(x)=0$, and $\Phi(x)$ is the unique solution on $[0,b)$. Moreover, using exactly the same arguments as above, it can also be shown that there exists some $a<0$ which may be infinite such that $\Phi(x)$ is strictly increasing on $(a,0]$, and $\lim_{x\rightarrow a^{+}}\Phi(x)=0$, and $\Phi(x)$ is the unique solution on $(a,0]$. On the interval of uniqueness, $(a,b)$, $\Phi_{x}$ is given by

[TABLE]

where $g$ is defined by eq. 45, and sgn is the sign function.

Now we describe a method which, given $\lambda$ and $\eta$, determines the unique $\Phi(x)$ on $(a,b)$, and also determines $a,b,H_{0}$ in the process. Indeed we will show that for any $\lambda$ and $\eta$, the interval of uniqueness for $\Phi(x)$ must be $(a,b)=\mathbb{R}$. The fact that $(a,b)=\mathbb{R}$ shall have the following subtle consequence. Note that the derivation of eq. 42 involved a multiplication by $\Phi\equiv\phi^{2}$. Thus, the deduction from eq. 42 back to eq. 37 holds on the condition that $\Phi\neq 0$. Since $a$ and $b$ are the smallest (in magnitude) zeros of $\Phi(x)$, we see that the equivalence between eqs. 42 and 37 breaks down outside the interval $(a,b)$. That is to say, eqs. 42 and 37 are equivalent globally if and only if $(a,b)=\mathbb{R}$.

The method is as follows. For $x\in[0,b)$, consider the coordinate transformation,

[TABLE]

$\Phi(x)$ is a bijection from $[0,b)$ to $(0,\Phi_{0}]$, $Y(\Phi)$ is a bijection from $(0,\Phi_{0}]$ to $(0,1]$, and the inverse sech function, arsech, is a bijection from $(0,1]$ to $[0,\infty)$. Therefore, all the coordinate transformations are invertible. For $x\in(a,0]$, we consider exactly the same transformations as eq. 48. Differentiating $Z$ with respect to $x$ we find, for all $x\in(a,b)$,

[TABLE]

where we have used definition (45) of $g(\Phi)$. Moreover, by definition we have $Y=\textnormal{sech}Z$, and it follows that $2\eta\Phi=2\eta\Phi_{0}Y^{2}=2\eta(2H_{0}/\lambda)\textnormal{sech}^{2}Z$. Defining

[TABLE]

we rewrite eq. 49 as

[TABLE]

Due to eq. 48, we have $Z(x=0)=0$. We can therefore solve eq. 51 as follows.

[TABLE]

implying

[TABLE]

Now we can determine the values of $a$ and $b$. In the limit $Z\rightarrow+\infty$, the definition of the coordinate transformations, as per eq. 48, dictates that we must have either $x\rightarrow a$ or $x\rightarrow b$. At the same time, eq. 53 dictates that we must have $x\rightarrow\pm\infty$, because the $\arctan$ function on the r.h.s. of eq. 53 is bounded, whilst the arsinh function diverges to $+\infty$. It therefore follows that $(a,b)=\mathbb{R}$.

The next step is to rewrite eq. 53 as an expression for $x$ in terms of $\Phi$, so that we can invert the expression to find $\Phi(x)$ for $x\in\mathbb{R}$. By definition (48) we have $\cosh^{2}Z=1/Y^{2}=\Phi_{0}/\Phi$, and it follows that $\sinh^{2}Z=\cosh^{2}Z-1=(\Phi_{0}/\Phi)-1$. Since $Z$ is by definition non-negative, we must take the positive square root, $\sinh Z=\sqrt{(\Phi_{0}/\Phi)-1}$. Then eq. 53 becomes

[TABLE]

We claim that, given $\Phi_{0}>0$ and $x\in\mathbb{R}$, eq. 54 uniquely determines a value of $\Phi>0$. The proof is as follows. If $x=0$, then immediately from eq. 54 we have $\Phi=\Phi_{0}$, and we are done. If $x\neq 0$, consider the function

[TABLE]

where $x$ and $\Phi_{0}$ are parameters. Differentiating eq. 55 with respect to $\Phi$, we find

[TABLE]

This means $\mathcal{G}$ is strictly decreasing for $\Phi>0$. Since $\mathcal{G}(\Phi)\rightarrow\infty$ in the limit $\Phi\rightarrow 0$, and $\mathcal{G}(\Phi_{0})=-\textnormal{sgn}(x)\sqrt{H_{0}}\leavevmode\nobreak\ x<0$, and $\mathcal{G}$ is continuous, we must have $\mathcal{G}$ vanishing at exactly one value of $\Phi\in(0,\Phi_{0})$. This concludes the proof. Moreover, we observe that in eq. 54 the l.h.s. is invariant under $x\mapsto-x$. Thus, on $\mathbb{R}$ we have $\Phi(-x)=\Phi(x)$.

In practice, given any $\Phi_{0}>0$ and $x\in\mathbb{R}$, we can compute $\Phi(x)$ by locating the zero of $\mathcal{G}(\Phi)$. However, the value of $\Phi_{0}$ cannot be freely chosen. Instead, it is determined by the normalisation condition (43) which, since $\Phi(x)$ is an even function on $\mathbb{R}$, now reads

[TABLE]

where the $Z_{x}^{+}$ is the positive-$x$ branch of $Z_{x}$, as per eq. 51. It then follows that

[TABLE]

Using $\Phi=\Phi_{0}\leavevmode\nobreak\ \textnormal{sech}^{2}Z$, we deduce

[TABLE]

Multiplying eq. 59b by $2\sqrt{\nu}$, replacing $H_{0}$ by $\lambda\Phi_{0}/2$, and replacing $\nu$ by $4\eta H_{0}/\lambda=2\eta\Phi_{0}$, we obtain the following transcendental equation for $\Phi_{0}$.

[TABLE]

To show that exactly one solution to eq. 60 exists, we consider the function

[TABLE]

Differentiating $\mathcal{F}(\Phi_{0})$ with respect to $\Phi_{0}$, we find

[TABLE]

This means $\mathcal{F}(\Phi_{0})$ is strictly increasing for $\Phi_{0}>0$. Since $\lim_{\Phi_{0}\rightarrow 0}\mathcal{F}(\Phi_{0})=-\sqrt{\lambda\eta}<0$, and $\mathcal{F}(\Phi_{0})\rightarrow\infty$ in the limit $\Phi_{0}\rightarrow\infty$, and $\mathcal{F}$ is continuous, we must have $\mathcal{F}$ vanishing at exactly one value of $\Phi_{0}>0$. In practice, given parameters $\lambda$ and $\eta$, we can compute $\Phi_{0}$ by locating the zero of $\mathcal{F}(\Phi_{0})$, and $\Phi_{0}$ uniquely determines the energy eigenvalue, $H_{0}=\lambda\Phi_{0}/2$. We can then feed the value $\Phi_{0}$ into eq. 54, and then for every $x\in\mathbb{R}$ we can find $\Phi(x)$ by means we have already described. In summary, given $\lambda$ and $\eta$, eqs. 54, 60 and 44 together constitute an analytical solution to eq. 42; and as we have already proven, it must be the unique global solution to eq. 42 which satisfies all the constraints we have imposed.

We note that if the parameter $\eta\rightarrow 0$, we should recover the solution to the nonlinear Schrödinger equation, given by eq. 36; and indeed we do. Firstly, in the limit $\eta\rightarrow 0$, we have $\nu\rightarrow 0$, which means we cannot use eq. 60 to determine $\Phi_{0}$, because the derivation of eq. 60 involved a multipication by $\sqrt{\nu}$. Instead, we must extract $\Phi_{0}$ from eq. 59. In the limit $\nu\rightarrow 0$, eq. 59a is simply $\sqrt{\lambda/(8\Phi_{0})}=\int_{0}^{\infty}\textnormal{sech}^{2}Z\leavevmode\nobreak\ \textnormal{d}Z=1$. It follows that $\Phi_{0}=\lambda/8$, agreeing with eq. 36. Then eq. 44 determines the eigenvalue $H_{0}=\lambda\Phi_{0}/2=\lambda^{2}/16$, again agreeing with eq. 36. Finally, when $\nu\rightarrow 0$, eq. 54 is simply

[TABLE]

which is equivalent to $\Phi_{0}/\Phi=1+\sinh^{2}(\sqrt{H_{0}}x)=\cosh^{2}(\lambda x/4)$, so $\Phi=\Phi_{0}\leavevmode\nobreak\ \textnormal{sech}^{2}(\lambda x/4)$, agreeing with eq. 36 once more.

The eigenvalue $H_{0}$ provides a link between $\Phi(x)$ and the binding energy of the stationary polaron. By eqs. 11, 18 and 29, where $k=0$ and $\xi=xR$, the polaron state is written in terms of local excitations as $\ket{\Psi}=\sum_{n=0}^{N}\alpha_{n}\hat{A}^{\dagger}_{n}\ket{\textnormal{vac}}$, and in the limit $N\gg 1$, we have

[TABLE]

Thus, the stationary $\ket{\Psi}$ solves $i\hbar\leavevmode\nobreak\ \textnormal{d}\ket{\Psi}/\textnormal{d}t=(\hat{H}_{\textnormal{e}}+\hat{H}_{\textnormal{int}})\ket{\Psi}$ as well as satisfying $i\hbar\leavevmode\nobreak\ \textnormal{d}\ket{\Psi}/\textnormal{d}t=(J_{0}-2J_{1}-H_{0}J_{1})\ket{\Psi}$. By definition, the polaron’s binding energy, $E_{\textnormal{b}}$, is its total internal energy measured with respect to $J_{0}$. In units of $J_{1}$, we have

[TABLE]

An expression for $\braket{\hat{H}_{\textnormal{p}}}$ in terms of $\Phi(x)$ can be found by using eqs. 4 and 24. Since the polaron is stationary, the kinetic part of $\braket{\hat{H}_{\textnormal{p}}}$ is zero, so we have

[TABLE]

It then follows that

[TABLE]

We have made use of definition (27) of $\lambda,\eta$, the fact that $\mathinner{\!\left\lvert\psi_{n}\right\rvert}=\mathinner{\!\left\lvert\phi_{n}\right\rvert}$ for all $n$, as well as the fact that $\mathinner{\!\left\lvert\phi_{n}\right\rvert}^{2}$ is approximated by $\Phi(n-N/2)$. Figure 1 shows how various aspects of the stationary polaron depends upon the symmetry parameter $\beta$, and the effective coupling parameter $\lambda$. 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 the only two parameters that affect the stationary polaron’s physical properties (as $\eta$ is merely a convenient combination of $\beta$ and $\lambda$).

Figure 1(a) shows how $\Phi_{0}$ and the half-width of the polaron varies with $\beta$ and $\lambda$. We define the half-width as the distance between the two $x$-values at which $\Phi(x)=\Phi_{0}/2$, and it is a measure of how localised the polaron is. As one would expect, the half-width is negatively correlated with $\Phi_{0}$, which is the maximum height of $\Phi(x)$. The figure shows $\Phi_{0}$ increasing with $\lambda$, and half-width decreasing with $\lambda$, and the rate of change of each quantity is greater given larger values of $\beta$. That is to say, the more spatially asymmetric the electron-lattice interaction is, the more influential $\lambda$ is. The figure also has the following implication on the accuracy of $\Phi(x)$ as an approximation to the discrete stationary solution to eq. 19. In a discrete solution, $\psi_{n}=\exp(i\rho H_{0}\tau)\phi_{n}$, the physical interpretation of $\mathinner{\!\left\lvert\psi_{n}\right\rvert}^{2}$ is the probability of the electron being localised at the $n$th lattice site. Therefore, the normalisation condition is defined in terms of a sum, $\sum_{n=0}^{N}\mathinner{\!\left\lvert\psi_{n}\right\rvert}^{2}=1$, and consequently we must have $\mathinner{\!\left\lvert\psi_{n}\right\rvert}^{2}\leq 1$ for all $n$. When a continuum solution $\Phi(x)$ is used to approximate the discrete one, we have the relation $\Phi_{0}\equiv\max{\mathinner{\!\left\lvert\psi_{n}\right\rvert}^{2}}$. Thus, any continuum solution with $\Phi_{0}>1$ cannot be reliable as an approximant. When $\beta=1$, $\Phi_{0}$ exceeds 1 if $\lambda$ is greater than 8, since $\Phi_{0}=\lambda/8$. On the other hand, when $\beta=0$, we computed $\Phi_{0}$ for $\lambda$ up to 100, and $\Phi_{0}$ remains less than 0.6.

In fig. 1(b) we see that $H_{0}$ increases with $\lambda$ whilst the polaron’s binding energy gains magnitude, meaning the larger $\lambda$ is the more energy is required to break up the polaron. Once again, the larger $\beta$ is, the more rapidly these quantities vary with $\lambda$. We note that the thick (black) curve for $H_{0}$, corresponding to $\beta=1$, is exactly the graph of $H_{0}=\lambda^{2}/16$, as per eq. 36. Comparing figs. 1(a) and 1(b), we see that a polaron which is more strongly bound has a larger $\Phi_{0}$ and a smaller half-width, i.e., it is more localised.

The gradient of curves in fig. 1(b) vary with $\beta$, and the variation is more pronounced when $\beta$ is close to 1. This suggests that the system is highly sensitive to variations in $\beta$ when $\beta$ is large, but not so when $\beta$ is small. Moreover, as $\lambda$ decreases, curves corresponding to different values of $\beta$ begin to converge; specifically this happens when $\lambda\approx 1$. This suggests that when $\lambda$ is small, the extent to which the electron-phonon interaction is spatially symmetric has little bearing on the physical properties of stationary polarons.

3.2 Numerical solutions

In this section we solve eqs. 28 and 30 directly, using a numerical scheme, but not without the help of analytical results from section 3.1. We then compare the resulting stationary polaron states with the ones we obtained via continuum approximation.

Expanding eq. 28 using the definitions of $\Delta\psi_{n}$ and $\Delta\mathinner{\!\left\lvert\psi_{n}\right\rvert}^{2}$, we have

[TABLE]

Any solution $\psi_{n}$ to eq. 68 is an attractor of the following map [27].

[TABLE]

where

[TABLE]

We take the approximate solution from section 3.1 as initial guess, and repeatedly apply eq. 69 until we reach convergence. When converged, $\psi_{n}$ is the stationary solution to eq. 68, and $\|\mathcal{H}(\psi_{n})\|$ is equal to $H_{0}$. In practice, on a grid with $N=200$, convergence is typically reached within $\mathcal{O}(10^{5})$ iterations, which amounts to $\mathcal{O}(10^{1})$ seconds of computing time. We have computed stationary solutions for various $\beta$ and $\lambda$, and some results are presented in fig. 2.

Figure 2(a) contains information about two key aspects of the stationary polaron state: the electron probability distribution, and the polaron binding energy. Qualitatively speaking, it is in agreement with predictions of the continuum approximation, as per fig. 1: as $\lambda$ increases, the polaron becomes more localised, and more strongly bound. Moreover, the effect of increasing $\lambda$ is more profound given larger values of $\beta$. However, further comparison between figs. 1(a) and 2(a) reveals a noteworthy difference. When $\beta=1$, $\Phi_{0}$ as a linear function of $\lambda$, whereas fig. 2(a) suggests that $\max\mathinner{\!\left\lvert\psi_{n}\right\rvert}^{2}$, which is approximated by $\Phi_{0}$, is not linearly dependent on $\lambda$. In fact, given any $\beta$, $\max\mathinner{\!\left\lvert\psi_{n}\right\rvert}^{2}$ grows significantly faster with $\lambda$ than fig. 1(a) predicts. Despite that, the growth of $\max\mathinner{\!\left\lvert\psi_{n}\right\rvert}^{2}$ in fig. 2(a) eventually stalls, when $\lambda$ becomes sufficiently large. This is 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 $\Phi_{0}$ can be, whereas the discrete system limits $\max\mathinner{\!\left\lvert\psi_{n}\right\rvert}^{2}$ to 1.

Figure 2(b) shows a selection of $\mathinner{\!\left\lvert\psi_{n}\right\rvert}^{2}$ solutions. Comparing all the dotted (red) lines, which correspond to $\lambda=1.0$ at various values of $\beta$, we see that they are essentially identical. This confirms the belief that when $\lambda$ is close to 1, systems with different $\beta$-values unify. The figure also shows some stationary solutions to the other half of eq. 19, namely $u_{n}$, which is expressed in terms of the stationary $\mathinner{\!\left\lvert\psi_{n}\right\rvert}^{2}$ 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 the point-dipole model for lattice units is valid, the lattice distortion must satisfy the condition $\mathinner{\!\left\lvert u_{n+1}-u_{n}\right\rvert}\ll 1$ [19]. This condition is indeed fulfilled in the stationary polaron state, since according to fig. 2(b) we have $\max\mathinner{\!\left\lvert u_{n+1}-u_{n}\right\rvert}\sim\mathcal{O}(10^{-2})$.

Comparing all the dashed (blue) lines in fig. 2(b), which correspond to $\lambda=4.4$ at various values of $\beta$, enables us to make the following observation. When $\beta=0$, the $u_{n}$ solution is centred at the location of $\max\mathinner{\!\left\lvert\psi_{n}\right\rvert}^{2}$, in the sense that its graph is rotationally symmetric about $n=80$. This agrees with our intuition that when $\beta=0$, i.e. when the electron-phonon interaction is spatially symmetric, the electron in the stationary state causes equal lattice distortion to its left and right. As $\beta$ increases, the maximum magnitude of $u_{n}$ remains the same, but the centre of $u_{n}$ shifts away from the location of $\max\mathinner{\!\left\lvert\psi_{n}\right\rvert}^{2}$, in response to the decrease in spatial symmetry. When $\beta=1$, the molecule at the location of $\max\mathinner{\!\left\lvert\psi_{n}\right\rvert}^{2}$ ($n=120$ in this case) is barely 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})^{2}$, which is the square of the gradient of the $u_{n}$ graph at site $n$. In the steady state, this gradient is zero except at a few sites around the location of $\max\mathinner{\!\left\lvert\psi_{n}\right\rvert}^{2}$, and it is clear that solutions corresponding to larger values of $\beta$ have steeper gradients there. We therefore conclude that, in the stationary state, systems with greater spatial asymmetry store more potential energy in the lattice.

In fig. 2(b) we also see a comparison between some $\mathinner{\!\left\lvert\psi_{n}\right\rvert}^{2}$ solutions and their counterpart continuum approximations, $\Phi(x)$. In particular, we look at the thick solid (black) lines and their accompanying thin solid (black) lines. The comparison reveals that, fixing $\lambda$, in this case $\lambda=2.6$, $\Phi(x)$ is a more accurate approximant to $\mathinner{\!\left\lvert\psi_{n}\right\rvert}^{2}$ when $\beta$ is smaller. As $\beta$ approaches 1, it becomes apparent that $\Phi(x)$ under-estimates the height of the $\mathinner{\!\left\lvert\psi_{n}\right\rvert}^{2}$ profile. If $\lambda$ is sufficiently large, however, $\Phi(x)$ becomes an over-estimate of the profile height. This all comes down once again to the fact that the continuum equations do not limit the height of the $\Phi(x)$ solution.

4 Dynamical polarons in zero temperature

In this part of the study we explore properties of polarons which propagate along the peptide chain, under an external forcing $\epsilon(\tau)$, and zero temperature ($f_{n}(\tau)=0$). Physically, $\epsilon(\tau)$ may represent the strength of a time-dependent electric field. We solve eq. 19 as an initial value problem, using the stationary $\psi_{n}$ and $u_{n}$ solutions which we computed in section 3.2 as the initial configuration of the system. We prescribe a suitable $\epsilon(\tau)$, setting $n_{0}$ to the location where the stationary $\mathinner{\!\left\lvert\psi_{n}\right\rvert}^{2}$ attains its maximum. Then we integrate the system forward in time using the 4th-order Runge-Kutta method. To ensure numerical stability, the integration time-step is set at $\Delta\tau=0.01$.

As we forward-integrate the system, we keep track of several scalar quantities associated with the polaron, such as its half-width and binding energy. Most importantly, we keep track of the polaron’s position, defined as follows. If $|\psi_{n}|^{2}$ attains its maximum at lattice site $\bar{n}$, then polaron position is the vertex location of the parabola extrapolated from three points: $(\bar{n},|\psi_{\bar{n}}|^{2}),(\bar{n}-1,|\psi_{\bar{n}-1}|^{2}),(\bar{n}+1,|\psi_{\bar{n}+1}|^{2})$. We note that, if the polaron is dynamical, then the stationary 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_{\textnormal{b}}$ as per eq. 65 will be computed directly from the numerical solutions.

4.1 Constant or periodic electric fields

The most obvious choice of $\epsilon(\tau)$ is a constant,

[TABLE]

Using moderately-localised stationary states (with $\max\mathinner{\!\left\lvert\psi_{n}\right\rvert}^{2}\approx 0.6$) as initial conditions, we computed polaron trajectories under various values of $\bar{\epsilon}$. Our results show that, given $\beta=1$ and $\lambda=3.0$, a constant forcing of any $\bar{\epsilon}\sim\mathcal{O}(10^{-2})$ induces nothing but small oscillations of the polaron around its initial position. An example of trajectory is presented in fig. 3.

As $\bar{\epsilon}$ is increased beyond 0.1, we find that eventually the forcing does become strong enough to dislodge the electron from its potential well, and propel the polaron along the peptide chain. However, as the polaron propagates, the magnitude of its binding energy decreases rapidly, and the polaron “delocalises”, i.e. breaks up into unbound components, within several hundred time units. A direct manifestation of the polaron’s energy loss and eventual delocalisation is that the $\mathinner{\!\left\lvert\psi_{n}\right\rvert}^{2}$ profile loses height and gains local peaks at lattice sites far away from the global maximum. For the sake of consistency, throughout the remainder of this study we shall say that a polaron has delocalised if its maximum height drops to below 0.1, as it must then be the case that other local peaks have magnitudes comparable to the global maximum. Figure 4 shows an example of a constant forcing large enough to cause polaron displacement, and it illustrates the resultant rapid delocalisation of the polaron.

Figure 4(a) shows the trajectory of a polaron which, within roughly 600 time units, is displaced by just over 300 lattice sites. Its binding energy steadily decreases in magnitude, until the polaron delocalises at $\tau\approx 600$. Meanwhile, fig. 4(b) shows the electron probability distribution, $\mathinner{\!\left\lvert\psi_{n}\right\rvert}^{2}$, at the time of delocalisation. This $\mathinner{\!\left\lvert\psi_{n}\right\rvert}^{2}$ profile has evolved from an initial configuration possessing 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 transport energy is diminished. Beyond the example shown in fig. 4, all our results are consistent with the hypothesis that, regardless of $\beta$ and $\lambda$, a constant $\epsilon$ causes the polaron to undergo either small periodic oscillations, or rapid losses in energy. We would like to find ways to displace the polaron without significant energy loss. Therefore, we must look for forms of $\epsilon$ other than constants. The next most natural choice of $\epsilon$ is periodic,

[TABLE]

where $A$ is the amplitude and $T$ is the period. Physically this may represent an electromagnetic plane wave which is monochromatic, i.e. coherent. Under periodic $\epsilon(\tau)$ with $A$ up to 0.2, regardless of $\beta$ and $\lambda$ we find that the polaron simply oscillates about its initial position. The polaron’s oscillatory motion has a period which coincides with $T$, and an amplitude which is positively correlated with $A$. An example of such trajectories is shown in fig. 3. While the polaron remains highly stable over time, its position averaged over its periodic remains constant. Thus, if we want polarons which transport energy from one lattice site to another, we must again look for an alternative form of $\epsilon(\tau)$.

4.2 Periodic electric fields with non-zero mean

Having studied the effects of constant forcing and periodic forcing in section 4.1, and discovered that neither serves to displace the polaron with minimal energy loss, in this section we consider forcing of the form

[TABLE]

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 $\epsilon(t)$ as a mean-shifted periodic forcing (MSPF). In particular, the mean $\bar{\epsilon}$ is chosen to be lower than the constant forcing $\epsilon$ which is required to displace the polaron, in the manner of fig. 4. Therefore, $\bar{\epsilon}$ on its own would not give the electron enough energy to escape its potential well. But we hope that the component $A$ can periodically push the electron energy over the threshold, resulting in polaron motion. Another possible advantage of this setup is that $A$ may periodically lower the electron energy, slowing it down and giving the lattice time to “catch up”, thus making the polaron motion more sustainable than it would be under a constant forcing.

Mathematically, $\epsilon(\tau)$ depends on three independent parameters, $\bar{\epsilon},A$ and $T$. Before investigating the effect of each of these parameters, we present fig. 5, which is a direct comparison with fig. 4.

We have replaced the constant forcing $\epsilon=0.15$, which resulted in fig. 4, with an MSPF which has the same maximum amplitude as before. The difference is that now this maximum amplitude is reached once every period $T$. Figure 5(a) shows that, within roughly 10 periods, the polaron is displaced by nearly 400 lattice sites. Contrary to the uniform manner in which the polaron moves in fig. 4(a), now the polaron moves towards one end of the peptide chain and then the other, within each period of $\epsilon(\tau)$. The overall displacement of the polaron is due to the fact that each movement to one end of the chain is larger than the subsequent swing back the other way. We note that while the polaron moves slightly further compared to fig. 4(a), its lifetime, i.e. the amount of time elapsed before delocalisation, is much longer. Overall, the polaron in fig. 5 propagates with a lower (average) velocity, $V$, defined by

[TABLE]

where the numerator is the displacement of the polaron, which we denote by $D$. Our results show that, of the parameters $\bar{\epsilon},A$ and $T$, the dominant factor which determines the polaron’s velocity is the constant component $\bar{\epsilon}$. We will discuss this in more depth in relation to fig. 6.

Figure 5(a) also shows how the polaron’s binding energy, $E_{\textnormal{b}}$, varies in time. Following an initial drop in magnitude, $E_{\textnormal{b}}$ mostly oscillates between $-0.75$ and $-1.5$, until another sharp decrease in magnitude leading up to delocalisation at $\tau\approx 4900$. There is an important observation to be made here. In fig. 4(a), we see that when the polaron reaches lattice site $n=300$, $E_{\textnormal{b}}$ is about $-0.6$. In fig. 5(a), the polaron’s average position over the 6th period is roughly 300, and the average binding energy over this period is $-1.2$. That is to say, under the MSPF, the polaron is carrying twice as much energy when it reaches $n=300$, compared to when it reaches $n=300$ under the constant forcing. Even though the constant forcing gets the polaron to $n=300$ in less time, we consider the MSPF a better mechanism for polaron transport, because the polaron binding energy is more stable. Indeed, the same can be said when the destination $n$ is anything larger 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 the MSPF does. In general, all our results are consistent with the hypothesis that, by splitting a constant forcing into constant and sinusoidal components, we lower the polaron’s velocity but increase its stability and lifetime. We say that, compared to the constant forcing, the MSPF is a better long-distance tranport mechanism, where speed can be sacrificed for energy efficiency.

In fig. 5(b) we see another aspect of the polaron’s motion, namely, how the height and half-width of the $\mathinner{\!\left\lvert\psi_{n}\right\rvert}^{2}$ profile vary with time. Following an initial decrease, the profile height, $\max\mathinner{\!\left\lvert\psi_{n}\right\rvert}^{2}$, 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\mathinner{\!\left\lvert\psi_{n}\right\rvert}^{2}$, occur precisely when the polaron turns from moving in one direction to moving in the other. This suggests that when the polaron accelerates, it “spreads out”, and so the half-width widens and $\max\mathinner{\!\left\lvert\psi_{n}\right\rvert}^{2}$ drops. We observe this phenomenon in all our results.

Based on our observations, we theorise that a polaron’s directed motion may be explained physically as follows. Since the forcing $\epsilon(\tau)$ is the effect of an electric field, it first-and-foremost provides the electron with extra energy. This is evident in the dramatic energy variation during the first period of $\epsilon(\tau)$ (see fig. 5(a)). Following this, it becomes much easier for the electron to overcome the significantly diminished polaron binding energy, $E_{\textnormal{b}}$. This is why the onset of polaron motion always follows a drastic drop in magnitude of $E_{\textnormal{b}}$. Whenever $|\epsilon(\tau)|$ becomes large enough to give the electron sufficient energy to overcome $E_{\textnormal{b}}$, the electron is dislodged from its potential well and propelled along the lattice. If the electron-lattice coupling is strong enough, then the lattice distortion can keep up with the electron, and so the polaron can remain intact. Whenever $|\epsilon(\tau)|$ drops below the binding threshold, the electron-lattice interaction slows down the electron and causes its probability distribution to spread out. This is why the half-width of $\mathinner{\!\left\lvert\psi_{n}\right\rvert}^{2}$ always peaks at times when the polaron’s instantaneous velocity is zero. If $|\epsilon(\tau)|$ 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. 6(b) and 6(c). If $\epsilon(\tau)$ has a large enough periodic component $A$, then it is possible for $|\epsilon(\tau)|$ to overcome the threshold twice per period: once with $\epsilon>0$, once with $\epsilon<0$. If the electron moves towards large $n$ in the $\epsilon>0$ case, then it will move towards small $n$ in the $\epsilon<0$ case. This explains the backwards swing exhibited by some polaron trajectories during each period of motion. The fact that $\bar{\epsilon}\neq 0$ ensures that the electron always spends more time moving one way than the other, hence the overall directedness of the polaron trajectories. This point is most clearly demonstrated by the trajectories in fig. 6(c), where the period $T=2000$.

Within each of figs. 6(a), 6(b) and 6(c), we can compare the trajectories represented by the same line type, but have different starting positions. This reveals the effect of varying the spatial symmetry of electron-lattice interaction. A greater spatial asymmetry, i.e. a larger $\beta$, causes the polaron to be more susceptible to displacement. Also within each of figs. 6(a), 6(b) and 6(c), we can compare trajectories starting from the same position, but are represented by different line types. This reveals the effect of varying the forcing amplitude $A$. The larger $A$ is, the more the polaron oscillates back and forth during each period of motion. We can also compare trajectories in figs. 6(a) and 6(b) which have identical line types and starting positions. This suggests that the overall velocity of the polaron is determined by $\bar{\epsilon}$, in the sense that the larger $\bar{\epsilon}$ is, the more the polaron moves per unit time. Finally, by comparing trajectories in figs. 6(b) and 6(c) which have identical 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 only provides more insight into the effect of $T$, but also helps to quantify our observations, and reinforce our theories, about the effects of $\bar{\epsilon}$ and $A$ .

For several combinations of $\bar{\epsilon}$ and $T$, we examine how the polaron’s lifetime, displacement and velocity vary with $A$. The results are displayed together, in fig. 7, so that the effects of $\bar{\epsilon}$ and $T$ can be gauged also. Firtly we consider the lifetime, $\tau_{\textnormal{d}}$. We computed all our numerical solutions up to $\tau=50000$, which is several times larger than the typical lifetime of a polaron that moves under the MSPF. If the polaron is not displaced by the MSPF, then it is effectively permanent, in the sense that its energy oscillates instead of dissipating over time, and it would have a lifetime far exceeding 50000. Thus, in fig. 7, the lifetime of a permanent (undisplaced) polaron is represented as $\tau_{\textnormal{d}}=50000$. For each combination of $\bar{\epsilon}$ and $T$, there exists some critical amplitude, $A=A_{\textnormal{c}}$, below which the polaron is undisplaced by the MSPF. At $A=A_{\textnormal{c}}$, the combined magnitude of the forcing, $\epsilon_{\textnormal{comb}}\mathrel{\mathop{\mathchar 58\relax}}=\bar{\epsilon}+A$, becomes large enough to displace the polaron, and $\tau_{\textnormal{d}}$ drops sharply. This drop can sometimes result in a lifetime of only several thousand time units - see for instance the bottom-left subfigure in fig. 7, corresponding to $\bar{\epsilon}=0.1$. When $\bar{\epsilon}$ is smaller, say $\bar{\epsilon}=0.005$ (top-left subfigure), the drop in lifetime is less dramatic. As $A$ increases beyond $A_{\textnormal{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 move barring small oscillations, the types of which we saw in fig. 3. As $A$ reaches critical value $A_{\textnormal{c}}$, the polaron turns from being quasi-stationary to moving by several hundred lattice sites during its lifetime. Evidently, the value of $A_{\textnormal{c}}$ is independent of $T$. Note that we only consider the displacement of polarons whose lifetimes are at least $2T$, and we set the displacement of polarons with shorter lifetimes to zero - see for instance the dashed (blue) lines in the centre and bottom-middle subfigures.

Whilst the value of $A_{\textnormal{c}}$ does not depend on $T$, the amount of displacement caused by $A_{\textnormal{c}}$ does. However, it is unclear from our results what their correlation is. As $A$ increases beyond $A_{\textnormal{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 to move. To understand how $A$ affects the amount of polaron displacement per unit time, we examine the polaron’s (average) velocity, $V$, as per definition (74). When $A$ is small, $V$ is zero. As $A$ reaches critical value $A_{\textnormal{c}}$, the velocity becomes typically $\mathcal{O}(10^{-2})$. Exactly what value this critical velocity $V_{\textnormal{c}}$ takes depends on $\bar{\epsilon}$ - the larger $\bar{\epsilon}$ is, the larger $V_{\textnormal{c}}$ is. As $A$ increases beyond $A_{\textnormal{c}}$, sometimes $V$ simply decays away - see for instance the top-right subfigure, where $\bar{\epsilon}=0.005$. Sometimes, however, $V$ grows before its decay - see for instance the middle-right and bottom-right subfigures, where $\bar{\epsilon}=0.03$ and 0.1 respectively. Such behaviour is possible when the polaron lifetime decays with $A$ more quickly than the displacement does. When this happens, there may exist some optimal amplitude, $A=A_{\textnormal{m}}$, at which the polaron attains maximum velocity, $V_{\textnormal{m}}$. $A_{\textnormal{m}}$ may coincide with $A_{\textnormal{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_{\textnormal{c}}$ remains the same, whereas the optimal $A_{\textnormal{m}}$ changes. Qualitatively speaking, the larger $T$ is, the smaller $A_{\textnormal{m}}$ is.

Whilst the value of $A_{\textnormal{c}}$ does not depend on $T$, it does depend on $\bar{\epsilon}$ - we see this by comparing any row of subfigures in fig. 7 to any other row. But how do $A_{\textnormal{c}}$ and $\bar{\epsilon}$ correlate? Our results show that, as $\bar{\epsilon}$ grows, $A_{\textnormal{c}}$ drops, but crucially the combined magnitude $\epsilon_{\textnormal{comb}}=\bar{\epsilon}+A_{\textnormal{c}}$ remains roughly constant. Specifically, in the top row we see $\bar{\epsilon}=0.005$ and $A_{\textnormal{c}}=0.167$, in the middle row we have $\bar{\epsilon}=0.030$ and $A_{\textnormal{c}}=0.142$, and the bottom row shows $\bar{\epsilon}=0.100$ and $A_{\textnormal{c}}=0.072$, each case giving $\epsilon_{\textnormal{comb}}=0.172$ when $A$ reaches critical. Recall that, when using a straightforward constant forcing $\epsilon=\bar{\epsilon}$, there is also a threshold value for $\bar{\epsilon}$, below which the polaron simply exhibits small oscillations, and above which the polaron moves at high speed but delocalises very quickly. It is noteworthy that this threshold is $\epsilon=0.154$ (given $\beta=0.6$), which is significantly lower than the critical combined amplitude of 0.172. In other words, $\epsilon=0.154$ causes polaron displacement, $\epsilon=0.153$ does not; and if one wishes to add on a periodic component $A\sin(2\pi\tau/T)$ in order to move the polaron, one needs $A\geq 0.019$, making $\bar{\epsilon}+A$ far exceed what $\bar{\epsilon}$ is required on its own to move the polaron. This phenomenon is observed across all values of $\beta$.

In practice, then, what would make a good combination of forcing parameters, which propel the polaron with decent speed but does not cause large energy dissipation too quickly? First of all, a large $\bar{\epsilon}$ results in a large velocity but an energetically unstable polaron which delocalises rapidly - so rapidly that it may move the polaron less far in its lifetime than a small $\bar{\epsilon}$ does. The middle column of fig. 7 precisely illustrates this point. Meanwhile, a small $\bar{\epsilon}$ results in long-living polarons which can move very far, because of how stable they are, but they would take more time to reach the same destination, compared to polarons under a large $\bar{\epsilon}$. On balance, a moderate value of $\bar{\epsilon}$ such as 0.03 is preferable. Secondly, once a $\bar{\epsilon}$ is chosen, it remains to choose $A$ and $T$, and it is obvious that the ideal choice of $A$ is the optimal amplitude, $A=A_{\textnormal{m}}$. Meanwhile, if $T$ is small, such as $T=100$, $A_{\textnormal{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 these extremes very clearly in the middle-right subfigure of fig. 7. Once again, these observations are not specific to $\beta=0.6$, but universal for all values of $\beta$. Overall, we believe that the best MSPF parameters which we have tested are such combinations where $\bar{\epsilon}\approx 0.030$, $T\approx 500$, and $A\approx A_{\textnormal{m}}$ which, given $\beta=0.6$, is $A_{\textnormal{m}}=0.157$. We will discuss the relationship between $A_{\textnormal{m}}$ and $\beta$ in section 4.3.

We note one anomaly which we observe in fig. 7 but did not expect. When $\bar{\epsilon}=0.005$ (top row), if $T=2000$ and $A$ is large enough, then the displacement (and therefore velocity) can take large negative values, meaning the polaron moves in the opposite direction to what we expected, and with large speeds. Whilst we are uncertain as to what causes this counter displacement, it is certainly another reason to reject small $\bar{\epsilon}$ and large $T$ when choosing forcing parameters.

4.3 The relevance of $\beta$

To produce fig. 7, we fixed $\beta=0.6$. How would the figure have looked if $\beta$ had been different? Our results show that qualitatively it would exhibit the same behaviour, characterised by critical amplitudes $A_{\textnormal{c}}$, and optimal amplitudes $A_{\textnormal{m}}$. Quantitatively, the values of $A_{\textnormal{c}}$ and $A_{\textnormal{m}}$ would change. It is therefore natural to investigate how they change with $\beta$. After all, our generalisation to the Davydov-Scott model is manifest in the extra parameter $\beta$. Firstly we establish the following preliminary result.

Recall that the stationary polaron, upon which we impose the MSPF, is characterised by two quantities: its probability distribution, specifically its maximum localisation probability, $\max\mathinner{\!\left\lvert\psi_{n}\right\rvert}^{2}$, and its binding energy. These are in turn determined by the symmetry parameter $\beta$ and effective coupling parameter $\lambda$. As $\beta$ varies, so does the value of $\lambda$ required to keep $\max\mathinner{\!\left\lvert\psi_{n}\right\rvert}^{2}$ constant. This correlation is shown in fig. 9. It is clear that $\lambda(\beta)$ is a decreasing function. We have made sure that whenever we altered $\beta$ we also took $\lambda=\lambda(\beta)$, so that all of our moving polarons begin as stationary states which share the same probability distribution. An alternative would have been to take whatever $\lambda$ is required to keep the binding energy constant. Our results show that if we had decided to keep the binding energy constant at, say, $-2.5$, then $\max\mathinner{\!\left\lvert\psi_{n}\right\rvert}^{2}$ would have been 0.53 at $\beta=0$, or 0.81 at $\beta=1$. It is a central feature of our model that two stationary polarons with the same maximum localisation probability need not have the same binding energy, and vice versa.

Having established the relation $\lambda(\beta)$, we can study how $A_{\textnormal{c}}$ and $A_{\textnormal{m}}$ depend on $\beta$. We do so by fixing $\bar{\epsilon}$ and $T$, and working out what $A_{\textnormal{c}}$ and $A_{\textnormal{m}}$ are for various $\{\beta,\lambda(\beta)\}$. For instance, in fig. 7 we saw that if $\beta=0.6$, $\lambda=\lambda(0.6)=5.0$, $\bar{\epsilon}=0.03$ and $T=500$, then $A_{\textnormal{c}}=0.142$ and $A_{\textnormal{m}}=0.157$. What if we fix $\bar{\epsilon}$ and $T$, and vary $\beta$? The result is displayed in fig. 9. When $\beta=0$, we have $A_{\textnormal{c}}=0.12$, and when $\beta=1$ we have $A_{\textnormal{c}}=0.1$. In fact, $A_{\textnormal{c}}$ is minimal when $\beta=1$, which suggests that a system with antisymmetric electron-phonon interaction is most conducive to polaron displacement by MSPF. One might expect that a system with symmetric interaction would be least conducive to polaron displacement, and therefore $A_{\textnormal{c}}$ should be maximal when $\beta=0$. This is not the case. We observe that $A_{\textnormal{c}}$ is maximal when $\beta=0.6$, which models a system with moderately asymmetric electron-phonon interaction. Meanwhile, polaron velocity produced by critical forcing, $V_{\textnormal{c}}$, is maximal when $\beta=0.87$.

The optimal amplitude, $A_{\textnormal{m}}$, varies little when $\beta$ is less than 0.7, but decays sharply when $\beta$ increases beyond 0.7, to such an extent that it almost equals the critical amplitude $A_{\textnormal{c}}$. The optimal velocity, $V_{\textnormal{m}}$, is typically of the same order of magnitude as $V_{\textnormal{c}}$. When $\beta=0.87$, $V_{\textnormal{m}}$ and $V_{\textnormal{c}}$ almost coincide.

Earlier, based on fig. 7, we asserted that the period $T$ of the MSPF has little effect on the value of $A_{\textnormal{c}}$. Indeed, our results show that, if we had produced fig. 9 with $T$ fixed at either 100 or 2000, the $A_{\textnormal{c}}$ curve would have been virtually unaffected. We also conjectured that the critical amplitude $A_{\textnormal{c}}$ is negatively and linearly correlated with $\bar{\epsilon}$, so that the combined magnitude $\epsilon_{\textnormal{comb}}=\bar{\epsilon}+A_{\textnormal{c}}$ remains constant as $\bar{\epsilon}$ varies. This is supported by our results. Indeed, fig. 9 is produced with $\bar{\epsilon}$ fixed at 0.03; but if we had produced fig. 9 with $\bar{\epsilon}$ fixed at either 0.005 or 0.1, the $A_{\textnormal{c}}$ curve would simply have been shifted along the vertical axis, by an amount equal to the difference between the new $\bar{\epsilon}$ and 0.03.

5 Dynamical polarons in non-zero temperature

We study the effect of random fluctuations which result from non-zero temperatures in the environment surrounding the lattice. The randomised forcing on the lattice is represented by the normally-distributed $f_{n}(\tau)$, which by definition (21) must have, for $\tau\geq 0$, the following first and second moments.

[TABLE]

where $\theta$ is the dimensionless temperature,

[TABLE]

and $\Theta$ is the temperature. Compared to the dimensional $F_{n}(t)$ that we introduced in section 2, the $\tau$ which now appears in $f_{n}(\tau)$ is a discrete index. We have followed the standard procedure of replacing $\delta(t-t^{\prime})$ by $1/\Delta\tau$, up to non-dimensionalisation constants. Beginning with a stationary polaron, which we computed numerically in section 3.2, we integrate the system of eq. 19 forward in time from $\tau=0$. Using a random number generation algorithm, we generate a new vector $f_{n}$ before each integration step. If $\theta$ is large, we find that it can cause large distortions in the lattice and rapid delocalisation of the polaron, due to excessive energy input to the system. For appropriate values of $\theta$, we see that the polaron’s binding energy undergoes small fluctuations, but on the average it tends to shift towards zero. After some time, the binding energy stabilises. For example, given $\beta=0.6$ and $\lambda=5.0$, the stationary polaron has binding energy $-2.47$. Integrating from $\tau=0$ with $\theta=0.0003$, we find that after $\tau\sim\mathcal{O}(10^{4})$ the binding energy settles, on average, around $-2.15$. The period of time required for a polaron to reach such a thermal equilibrium is the thermalisation phase of the polaron dynamics. During this phase, the forcing $\epsilon(\tau)$ on the electron is kept at zero, but the electron nevertheless undergoes small fluctuations around its initial position, due to its coupling to the thermalised lattice. Our results show that, irrespective of $\beta$, we are unable to raise $\theta>0.001$, because such a large $\theta$ induces excessive lattice distortions which cause the polaron to delocalise before reaching thermal equilibrium.

In each simulation, we integrate the system, with $\epsilon(\tau)=0$, until the polaron reaches thermal equilibrium. Then we reset $\tau=0$, and “turn on” the forcing $\epsilon(\tau)$ for $\tau\geq 0$. We examine how the polaron subsequently moves, under combinations of $\epsilon(\tau)$ and $f_{n}(\tau)$. Since the thermalisation phase raises the polaron energy, we expect that a thermalised polaron would be easier to displace, in the sense that it would require a smaller $\epsilon(\tau)$ to displace it, compared to the zero temperature case. Indeed, our results confirm this. To obtain our results in this section, every dynamical simulation, with a set of chosen parameters $\{\beta,\lambda(\beta),\bar{\epsilon},A,T,\theta\}$, 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 $\beta=0.6$ and $\theta=0$. Specifically, fig. 10 is to be compared with the solid (black) lines in the middle row of subfigures in fig. 7, for which two of the parameters in $\epsilon=\bar{\epsilon}+A\sin(2\pi\tau/T)$ were fixed: $\bar{\epsilon}=0.03$, and $T=500$. We saw that, given said parameter values, the critical amplitude was $A_{\textnormal{c}}=0.142$. When we have a non-zero $\theta$ in the system, we define $A_{\textnormal{c}}$ to be the smallest $A$ for which the average polaron displacement (over 100 simulations) exceeds 10 lattice sites. According to this definition, when $\beta=0.6,\bar{\epsilon}=0.03,T=500$ and $\theta=0.0001$, we see in fig. 10 that $A_{\textnormal{c}}=0.121$, which is significantly lower than the case of $\theta=0$. Fixing said values of $\bar{\epsilon},T$ and $\theta$, we find that the value of $A_{\textnormal{c}}$ depends on $\beta$ in a manner shown in fig. 12. That is, $A_{\textnormal{c}}$ is minimal when $\beta=1$, suggesting that an antisymmetric electron-phonon interaction makes it easiest to displace the polaron. Meanwhile, $A_{\textnormal{c}}$ is maximal when $\beta\approx 0.5$, suggesting that, counter-intuitively, what makes displacing the polaron most difficult is not a symmetric electron-phonon interaction, but a moderately asymmetric one. Indeed, this $A_{\textnormal{c}}(\beta)$ function is very similar to the one in fig. 9, where we also had $\bar{\epsilon}=0.03,T=500$ fixed, but $\theta=0$. Now with $\theta=0.0001$, the $A_{\textnormal{c}}(\beta)$ curve in fig. 12 is significantly lower. This means that, regardless of $\beta$, a non-zero temperature makes it easier to displace the polaron by the forcing $\epsilon=\bar{\epsilon}+A\sin(2\pi\tau/T)$, in the sense that a smaller combined magnitude $\bar{\epsilon}+A$ is required. It is also noteworthy that, under a non-zero temperature, the onset of polaron motion is more gradual, in the sense that a critical amplitude results in a very small velocity, $V_{\textnormal{c}}$. Indeed, comparing $V_{\textnormal{c}}$ in fig. 9 with $V_{\textnormal{c}}$ in fig. 12, we see that the latter is 2 orders of magnitude smaller.

More can be said about fig. 10. When we raise the temperature to $\theta=0.0005$, we find that the polaron is displaced (on average) by hundreds of sites even if $A=0$. This suggests that, given $\beta=0.6$ and $\bar{\epsilon}=0.03$, there exists some critical temperature $\theta=\theta_{\textnormal{c}}$ between 0.0001 and 0.0005, for which just the combination of $\epsilon(\tau)=\bar{\epsilon}$ and $f_{n}(\tau)$ is sufficient to displace the polaron, and no periodic component in $\epsilon(\tau)$ is needed. $\theta_{\textnormal{c}}$ is critical in the sense that, if $\theta$ is any lower than $\theta_{\textnormal{c}}$, then the combination of $\epsilon(\tau)=\bar{\epsilon}$ and $f_{n}(\tau)$ does not energise the polaron enough to move it, and a non-zero $A$ is required. Indeed our results show that, given $\beta=0.6$ and $\bar{\epsilon}=0.03$, the critical temperature is $\theta_{\textnormal{c}}=0.00032$. Furthermore, we have investigated how $\theta_{\textnormal{c}}$ changes as we vary $\beta$ and $\bar{\epsilon}$, and the results are shown in fig. 12. We observe the qualitative trend that, the larger $\bar{\epsilon}$ is, the less thermal energy is required to make up for the extra energy that the polaron needs in order to move. We also observe that, in general, the larger $\beta$ is, the less thermal energy is required to displace the polaron. This fits in nicely with our understanding that, when $\beta$ is close to 1, we have an electron-phonon interaction which is biased towards one end of the lattice, making the electron more susceptible to displacement.

In fig. 13 we present a typical polaron trajectory when $\theta=\theta_{\textnormal{c}}$. Under this critical temperature, some simulations would produce no polaron displacement at all, but most trajectories would be similar to that in fig. 13, clearly showing a directed movement. We propose to explain the shape of these trajectories as follows. First of all, the combined magnitude of the MSPF would be much lower than what is required to move the polaron under zero temperature. In fig. 13 for example, we have $\bar{\epsilon}+A=0.03$, whereas the critical value under zero temperature, as we discovered in section 4.2, is $\bar{\epsilon}+A=0.172$. Even when there is a non-zero temperature, the polaron still spends the majority of its lifetime oscillating around its localisation site by small amounts. However, occasionally the random forces on the lattice sites in the vicinity of the electron causes a large distortion, such that the effective potential barrier for the electron is significantly lowered, and the electron can escape the well. Once it does that, it is propelled towards one end of the lattice by $\bar{\epsilon}$. But before the electron has time to move far, the random forces may have further distorted the local lattice sites in such a way that a high potential barrier is restored. This then traps the electron again, giving the polaron time to recover its integrity, before the next random time at which the electron jumps out of its potential 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 smooth and regular polaron motion exhibited in fig. 6.

6 Discussions and conclusions

In this study we have presented a new mathematical model describing polaron dynamics in linear peptide chains. The model is dependent on a symmetry parameter, $\beta$, which measures the extent to which the interaction between the polaron’s electron and phonon components is spatially symmetric. We have shown that when $\beta$ takes its extreme values, 0 and 1, the model reduces to existing ones for which it was assumed that the electron-phonon interaction was, respectively, symmetric and antisymmetric. We have justified the physical neccessity of including $\beta$ in the model, in that one should not simply assume the electron to be coupled equally strongly to lattice points on either side, or to be coupled only to the lattice point on one side. Instead, the spatial symmetry should be determined by the adjustable parameter $\beta$.

Apart from $\beta$, we have also identified two composite parameters which are most vital to the intrinsic properties of the polaron. Firstly there is the adiabaticity parameter, $\rho$, measuring the characteristic time scale separation between the electron and phonon, which we justifiably fixed throughout the study. Then there is the effective coupling parameter, $\lambda$, measuring the strength of the electron-phonon interaction. The combination of $\beta$ and $\lambda$ determines the two aspects of the stationary polaron: its maximum localisation probability, and its binding energy. We have computed both of these quantities as functions of $\beta$ and $\lambda$. Moreover, in the infinite lattice limit, we have obtained stationary polaron solutions by analytically integrating the system, 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 which causes minimal energy loss and which, crucially, is directed. Such polaron dynamics could be achieved only if the electron is dislodged from its self-trapping potential well, and the local lattice distortions propagate coherently with the electron, and some mechanism exists which ensures the electron always moves towards one end of the lattice. If the second condition is not met, then over time the electron probability density function would become broader, leading to delocalisation of the polaron. We have found that a constant external force, $\bar{\epsilon}$, on the electron is insufficient for displacing the polaron, unless $\bar{\epsilon}$ is larger than some threshold value, but then the forcing causes rapid energy loss and delocalisation. We have also found that a sinusoidal force, $A\sin(2\pi\tau/T)$, on the electron is never sufficient for displacing the polaron, throughout the range of $A$ that we tested. We then combined the constant and sinusoidal forces, resulting in the mean-shifted periodic forcing (MSPF), $\epsilon=\bar{\epsilon}+A\sin(2\pi\tau/T)$. We have discovered that, for each $\bar{\epsilon}$ which is insufficient on its own to displace the polaron, there is some critical value $A_{\textnormal{c}}$, such that the polaron is displaced if and only if $A\geq A_{\textnormal{c}}$. There is also an optimal value $A_{\textnormal{m}}$, such that the polaron attains maximum velocity at $A=A_{\textnormal{m}}$. The value of $A_{\textnormal{c}}$ is irrespective of the period $T$, whilst $A_{\textnormal{m}}$ is negatively correlated with $T$. As $\bar{\epsilon}$ is decreased, $A_{\textnormal{c}}$ increases, in such a way that the combined magnitude $\bar{\epsilon}+A_{\textnormal{c}}$ remains constant. This suggests that there is a certain amount of extra energy that the electron needs in order to overcome the polaron binding, and how much of it comes from the constant or sinusoidal part is inconsequential, as long as the two parts combine to a large enough overall amplitude. Nevertheless, the split between $\bar{\epsilon}$ and $A$ does determine the manner in which the polaron propagates, specifically its velocity and stability. The velocity is predominantly determined by $\bar{\epsilon}$, and positively correlated with it; but the stability of the polaron is negatively correlated with $\bar{\epsilon}$. By comparing three sets of $\{\bar{\epsilon},A\}$ with the same combined $\bar{\epsilon}+A$, namely $\{0.005,0.167\},\{0.03,0.142\},\{0.1,0.072\}$ (while keeping all other parameters fixed), we found that $\{0.03,0.142\}$ produces optimal balance between polaron velocity and stability.

We have examined how the aforementioned phenomena depends upon $\beta$. To do so, we needed a way of isolating the effect of varying $\beta$. This posed a difficulty, because if we were to fix $\lambda$ and vary $\beta$ then both the maximum localisation probability and binding energy of the stationary polaron would change. We would then be comparing dynamical behaviours of dissimilar polarons. We therefore decided to vary $\lambda$ with $\beta$, in a way that allowed us to generate a set of stationary polarons, one for each combination of $\{\beta,\lambda(\beta)\}$, such that they all had the same maximum localisation probability. Then we launched these polarons using the same external forcing and compared the results. We have found that $\beta=1$, representing a spatially antisymmetric electron-phonon interaction, produces a polaron which is easiest to move, in the sense that the least amount of forcing is required. We have also found that the symmetric model, $\beta=0$, does not make the polaron most difficult to move ($\beta\approx 0.6$ does that). This hints at the existence of some intrinsic mechanism in the $\beta=0$ model which pushes the electron towards one end of the lattice, despite it being coupled to the other end equally strongly.

We have also studied the MSPF under non-zero temperatures, $\theta>0$. The manifestation of thermal effects is random forces on the lattice points. We have found that a non-zero $\theta$ facilitates polaron propagation, in the sense that it lowers the critical amplitude $A_{\textnormal{c}}$, for any given $\bar{\epsilon}$. Moreover, a non-zero $\theta$ results in a gradual onset of polaron motion, meaning the rate of change of polaron velocity with respect to $A$ near $A=A_{\textnormal{c}}$ is small, compared to the onset under $\theta=0$. Our results have also shown that, whenever there is polaron propagation, whether $\theta=0$ or $\theta>0$, the relative displacements between neighbouring lattice points remain under $\mathcal{O}(10^{-2})$. This is a necessary condition which allows us to model the lattice points as point dipoles.

Some of the choices of parameters in the MSPF may be justified physically as follows. It is well known that across the plasma membrane of a living cell, a resting membrane potential is maintained by intercellular chemical processes [28]. It is also well known that within the plasma membrane there exist highly stable transmembrane regions of proteins, for instance the human prolactin receptor 2N7I [29], and the rat monoamine oxidase A 1O5W [30], both of which are $\alpha$-helical structures spanning the entire membrane width. Given a constant potential difference of $\Delta V$ across a linear, homogeneous, isotropic dielectric medium with constant width $d$, the effective electric field inside the medium is given by $E_{0}=\Delta V/(\kappa d)$, where $\kappa$ is the dielectric constant (a.k.a. relative permittivity) of the medium [31]. Assuming the plasma membrane is such a medium, we can then attribute the physical origin of $\bar{\epsilon}$ to the resting membrane potential, and calculate $\bar{\epsilon}$ using eq. 21, namely $\bar{\epsilon}=qE_{0}R/(\hbar\Omega)$. For many cells the values of the $E_{0}$ 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_{0}=-8.4$mV, as well as [33, 34] for $d=78\textnormal{\r{A}}$. However, the value of $\kappa$ for a membrane is highly contentious, due to the fact that it depends sensitively upon a large variety of biophysical attributes of the membrane, such as hydration [35], pH value [36], and structural stability [37]. In a recent review, it was reported that the value of $\kappa$ in literature ranges from 1 to 40 [38]. Feeding these values of $E_{0},d$ and $\kappa$ into the equation for $\bar{\epsilon}$, we find that $\bar{\epsilon}$ ranges from 0.0033 to 0.13. We have taken care to ensure that in this study the values of $\bar{\epsilon}$ falls strictly within this range. For a physical origin of the periodic term, $A\sin(2\pi\tau/T)$, one could look to common electromagnetic radiations which fill the environment around us in the modern age, such as the radiation from telecommunication 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 Ku band frequencies for satellite communications and broadcasting, and the UHF band frequencies for cellular communications [39, 40]. However, the amplitudes of the aforementioned radiations 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 peak power $P$, we can estimate the amplitude $\tilde{A}$ of its output waves at operational distance $d$, by using the well-known formula $P/(4\pi d^{2})=\epsilon_{0}c\tilde{A}^{2}/\kappa$, where $\epsilon_{0},c,\kappa$ are the vacuum permittivity, speed of light, and relative permittivity of the medium, respectively. Feeding $P=1\textnormal{W}$ [41] and $\kappa\leq 40$ [38] into the formula, we obtain dimensionless $A\leq 6\times 10^{-6}/(d/\textnormal{metres})$. This means that in order to obtain $A=0.1$, one needs the operational distance $d$ to be $\mathcal{O}(10^{-5})$ metres, which is unrealistic. It is therefore clear that, in a real cell environment, the effect on a polaron due to a combination of resting membrane potential and random thermal forces are dominant over any external electromagnetic radiation that may commonly be present. This highlights the importance of our observation that, given a constant electric field $\bar{\epsilon}$, there exists some critical temperature $\theta_{\textnormal{c}}$, such that the polaron undergoes directed drift if and only if $\theta\geq\theta_{\textnormal{c}}$. In other words, a combination of $\bar{\epsilon}$ and $\theta$ can be sufficient for displacing the polaron, in a manner which causes minimal energy loss and which is directed, just as a combination of $\bar{\epsilon}$ and $A\sin(2\pi\tau/T)$ can. It has been reported that stationary polarons formed on 3-dimensional lattices, such as an $\alpha$-helix, are more strongly bound and therefore can be stabler during propagation [42]. It is our hope that our model can be adapted to study such 3-D systems, and that the stabilising effect of the helical geometry will enable us to raise $\theta$, from our current values of $\mathcal{O}(10^{1})$K to physiological temperatures.

Acknowledgement

We would both like to thank Larissa Brizhik for kindly answering some questions.