Holstein polaron at fixed temperature: influence of the chain length

Based on the semiclassical Holstein model, the dynamics of a quantum particle in one-dimensional molecular chain with a trapping site is modeled. Numerical simulation is used to investigate the dynamics of a polaron in a chain with small random Langevin-like perturbations which imitate the thermostat. Parameter values are chosen such that the polaron energy at the trapping site is much greater than the energy of temperature fluctuations. Results of modeling demonstrate that temperature decay of a polaron depends on the chain length even at very low temperatures.


Introduction
At present, the attention of researchers is attracted to the possible mechanisms of transfer of charged particle (electron or hole) in quasi-one-dimensional biomacromolecules, such as DNA or proteins. The interest in the charge transfer mechanisms is also associated with the possibility of using DNA in molecular electronics [1,2]. Due to its stability the polaron mechanism has been intensively studied (see for example [2][3][4][5] and references therein). Usually the polaron dynamics along unperturbed chains of molecular sites is investigated, and it is assumed that weak fluctuations of the chain (with energy much less than the characteristic energy equal to the depth of the polaron level) change the polaron properties only slightly.
However, this assumption is not entirely correct. The results of computer simulations for uniform chain demonstrate [6,7] that the existence/disrupting of polaronic state depends not only on the thermostat temperature T , but also on the length of the chain N , i.e. on the thermal energy of the chain.
In this work, we consider a uniform chain with a trap site in the middle, that captures the charged particle. With the use of a numerical experiment we have shown that temperature decay of a polaron depends on the chain length even at very low temperatures.

The model
The model is based on the Holstein Hamiltonian for discrete chain of sites [8]: Here variables b n are the amplitudes of the probability of the charge (electron or hole) occurrence on the n-th site (n = 1, . . ., N , N is the chain length), andũ n are site displacements with respect to its center of mass. Parameters ν mn (m ̸ = n) are matrix elements of the charge transition between m-th and n-th sites (depending on overlapping integrals), ν nn is the electron energy on the n-th site; M n is the n-th site's effective mass, K n is the elastic constant. The last term suppose that the charge energy at the sites depends linearly on the sites displacementsũ n , α ′ is the coupling constant.
The dimensionless equations of motion for Hamiltonian (1) have the form: The relation between dimension and dimensionless parameters is as follows. τ is arbitrary characteristic time, the matrix elements η = ν nn±1 τ / , η kk = ν kk τ / , the frequencies of sites oscillation ω = √ τ 2 K/M , the coupling constant is χ = α ′ √ τ 3 / M . To model Langevin thermostat, subsystem (3) involves the terms with friction (γ = τγ/M is the friction coefficient) and the random force Z n (t) with the properties where the dimensionless temperature is T =T [K]/T * , T * is characteristic value. Choosing T * = 1 K and τ = 10 −14 sec, the characteristic energy value E * = k B T * τ / is E * ≈ 0.001309.

About the choice of parameter values
We performed calculations for system (2,3) with the following parameter values: η = 1, χ = ω = 0.5, η nn = 0 except for the trapping site with η kk = −5. For estimation of small values of temperature T Gershgorin circle theorem [9] was applied. Quantum subsystem (2) in vector form will be iḃ = Ab, where matrix A is symmetric tridiagonal matrix with elements a nn = η nn + χu n on the main diagonal and elements η on first diagonal below/above. The eigenvalues of matrix A corresponding to different energy levels lie in the union of the Gershgorin discs, centered at a nn with the radius 2η. For linear chain of N sites with selected parameter values there are N − 1 discs centered at 0 with the radius R = 2, and one disc for trapping site centered at −5 with the radius R = 2, which do not intersect with others, i.e. the energy value of charge localization on the trapping site is separated from other values. Taking into account the interaction with classical subsystem leads to the fact that the centers of these discs change in time; at the moment t the centers will be located at the points η nn + χu n (t), i.e. formally the discs can intersect and the eigenvalues can be very close to each other.
The dynamics of the chain sites u n is described by harmonic oscillators. Under the temperature fluctuations, averaged ⟨u n ⟩ = 0, ⟨u 2 n ⟩ = E * T /ω 2 . Based on standard deviation σ T = √ E * T /ω 2 , it is possible to estimate the changes of the centers of the Gershgorin discs for u n = 3σ T and the temperature T cr , below which the discs intersect with a very low probability. For the chosen values of the parameters, from the condition |η kk |−4η > 6(χ/ω) √ E * T , we obtain T cr ≈ 20.
Previous results of modeling in an uniform chain [6,7] demonstrate that the existence/disrupting of polaron depends not only on the thermostat temperature T , but also on 3 the length of the chain N , i.e. on the thermal energy of the chain N T . At the same temperature, in short chains a polaron does not decay while in long chains a charge is delocalized. The boundary {N T } * between polaron and delocalized states was estimated numerically. Similar results were obtained for homogeneous chain with a defect [10] for range T ≥ 50.
As in the papers [6,7], we estimated boundary {N T } * for these values of parameters. Calculations for N = 50 at different T lead to the following estimation: {N T } * ≈ 5000.
If in 'small temperature region' the polaron disruption depended only on temperature, then at T < T cr in chains of any length, the charge would be localized at the trapping site. In this work we performed numerical simulations at T = 5 < T cr ≈ 20 for the chain lengths N = 500 (N T = 2500 < {N T } * ) and N = 1100 (N T = 5500 > {N T } * ).

Results
For each N we have computed a set of realizations (trajectories of the system (2), (3) from different initial data and with different pseudo-random time-series) and calculate time dependencies averaged over realizations ("by ensemble"). Calculations of individual realizations were performed using the 2o2s1g-method [11].
In the numerical experiment we have calculated (and then averaged over 50 realizations) the probabilities to find the charge at the n-th site p n (t) = |b n (t)| 2 , and the parameter of delocalization R(t) This value correlate with charge distribution in the chain. If a charge is localized at one site, p n (t) ≈ 1, then R(t) ≈ 1. If a charge is uniformly distributed over a N -site chain, p n = 1/N , then R = N .

Numerical simulation for chains of 1100 sites
Number of the trapping site is k = 550. In chains of 1100 sites polaron initial state is gradually disrupting. Since averaging over realizations obscures the essence of the process, let's consider one trajectory as an example of such destruction. Figs. 1-3 show the probability distribution p n (t) for some three time points and the site displacements of the chain u n (t) at the same moments.
At the beginning of the calculation (Fig. 1) the probabilities in the chain are small, but nonzero, except for the polaron at the 550-th site (the probability maximum is cut off). Due to the interaction with the charge, the displacement of this 550-th site (≈ −1.9) is significantly larger than the displacements of other sites, caused mainly by temperature fluctuations.
Polaron destruction takes a long time. Fig. 2 shows that the charge probability on the trap site is still large compared to the probabilities on other sites p 550 > p n , and the displacement of the trap site is close to the displacements of other sites in the chain |u 550 | ≈ max n̸ =550 |u n |. In our opinion, this is no longer a polaronic state, since the displacement of the trap site is mainly determined by temperature fluctuations. It is possible to estimate the 'polaronic' displacement u p associated with the interaction with the charge as follows: at T = 0 stationary mode of eqs. (2,3) isu n = 0, u n = Const. From classical subsystem (3) u n = −(χ/ω 2 )|b n | 2 . For p 550 ≈ 0.025 (Fig. 2, top) u p ≈ −0.05 and |u p | < σ T ≈ 0.16. Fig. 3 shows that at the end of the calculation the charge is in a completely delocalized state, the displacements are caused by temperature fluctuations at all sites, including the trap site. In this state, the charge is most likely to be localized on any site in the chain, not only on the trapping site; in Fig. 3 (top) max n p n ≈ 0.006 on site n = 985.
In different realizations the dynamics of polaron destruction differs in time, but the disruption pattern is qualitatively the same.

Numerical simulation for chains of 500 sites
Number of the trapping site is k = 250.
When started from polaron initial state, modeling results in deformed polaron, but it is not destroyed. Therefore, we performed a calculation for the initial data 'charge is uniformly distributed along the chain', p n = 1/500 for all n = 1, . . . , 500. In this case, a pattern is observed that is opposite to the dynamics for N = 1100: for some time the charge is in a delocalized state, and then the probability to find a charge begins to increase at the trap site, and a polaron state is formed there.
When computing realizations, we monitor the following time-dependencies: the probability maximum p M (t) = max n p n (t), the number of the site n M (t) at which p M is localized, and the displacement u M (t) of this site n M . Also we monitor maximum modulus of displacements u L in the chain far from n M , which depend mainly on temperature fluctuations: |u L | = max n,|n−n M |>10 |u n |, as in [12]. Fig. 4 demonstrates polaron formation process (averaged by realizations). At first, p M (t) can be found on any site, but over time the charge is localized on the trap site in all realizations: ⟨n M (t)⟩ = 250 for t > 7.5 · 10 6 (Fig. 4, top). At the same time p M (t) increases (Fig. 4, middle) and |u M (t)| > |u L (t)| (Fig. 4, bottom), e.g., polaron state is forming, although it is not as perfect as at T = 0.

Conclusion
Results of modeling demonstrate that temperature disruption of a polaron depends on the chain length even at very low temperatures. Here we suppose that a polaron is a self-consistent state when the charge is localized in a small region and site displacements in this region are significantly larger than the ones caused by temperature fluctuations.
We considered simple Holstein model with Langevin thermostat, without interaction between sites in a classical chain. In the simulations for low temperatures it is necessary to compute very long trajectories. In more realistic models of polymers, e.g. Peyrard-Bishop-Dauxois [13] or Su-Schrieffer-Heeger (SSH) [14], Hamiltonian includes terms with intersite potential, which may accelerate the redistribution of energy in the system 'initial polaron + temperature', so similar computer experiments will take less time. We believe that other semiclassical models will demonstrate some similar features, such as dependence of localized state on the chain length at the same temperature.