Computer simulations of localized small polarons in amorphous polyethylene

Computer simulations of localized small polarons in amorphous polyethylene David Cubero and Nicholas Quirke Citation: The Journal of Chemical Physics 120, 7772 (2004); doi: 10.1063/1.1667471 View online: http://dx.doi.org/10.1063/1.1667471 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/120/16?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Vibron-polaron critical localization in a finite size molecular nanowire J. Chem. Phys. 122, 014701 (2005); 10.1063/1.1828031 Novel method to estimate solubility of small molecules in cis-polyisoprene by molecular dynamics simulations J. Chem. Phys. 115, 6258 (2001); 10.1063/1.1398590 Molecular dynamic simulation of the interaction, at high energy, between the N 2 molecule and polyethylene J. Chem. Phys. 113, 8187 (2000); 10.1063/1.1317228 Clustering of water in polyethylene: A molecular-dynamics simulation J. Chem. Phys. 109, 6476 (1998); 10.1063/1.477293 Subglass chain dynamics and relaxation in polyethylene: A molecular dynamics simulation study J. Chem. Phys. 108, 9912 (1998); 10.1063/1.476430 This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 150.214.182.194 On: Mon, 01 Jun 2015 15:17:42 Computer simulations of localized small polarons in amorphous polyethylene David Cubero and Nicholas Quirke Department of Chemistry, Imperial College, London, SW7 2AY, United Kingdom 共Received 20 October 2003; accepted 14 January 2004兲 We use a simple mean field scheme to compute the polarization energy of an excess electron in amorphous polyethylene that allows us to study dynamical properties. Nonadiabatic simulations of an excess electron in amorphous polyethylene at room temperature show the spontaneous formation of localized small polaron states in which the electron is confined in a spherically shaped region with a typical dimension of 5 Å. We compute the self-trapping energy to be ⫺0.06⫾0.03 eV, with a lifetime on the time scale of a few tens of picoseconds. © 2004 American Institute of Physics. 关DOI: 10.1063/1.1667471兴 I. INTRODUCTION Polyethylene is the simplest organic insulator, playing a very important role in a number of technological applications such as high tension insulation. Despite a vast literature1 concerned with the experimental characterization of its electrical properties, very little is known about the details of electronic transport at the molecular level. An understanding of the mechanisms of charge transport in these materials is important in determining the electronic and optical properties and in the development of new materials with more reliable insulating and other properties. Recent femtosecond spectroscopy experiments have shown the existence of shallow self-trapped polarons in ultrathin alkane layers on a silver共111兲surface.2In addition, recent Car-Parrinello simulations3have shown the spontaneous formation of self-trapped polarons in bulk crystalline polyethylene 共modeled using four chains of seven methylene units in periodic boundary conditions, each one initially alltrans兲. This shallow polaron has been linked to the formation of two opposite trans-gauche defects in a single chain, with the electron trapped in the rotated portion of the chain. In this paper, we will demonstrate the formation of localized polarons in amorphous polyethylene. As in bulk crystalline polyethylene we find very small polarons, with self-trapping energies comparable with the thermal energy kBT, but the geometry of the self-trapped state is very different, being essentially isotropic, reflecting the underling symmetry of the dielectric phase. However, with such small self-trapping energies one could question the physical meaning of such localized polaron states. We will show that the lifetime of these states is large enough to allow them to play a very important role in electron transport. In our previous work we developed a new pseudopotential for electron–polyethylene interactions, which has been used to study electronic states in crystalline and amorphous configurations of polyethylene. These states are found to be in good agreement with the experimental data.4,5 We also have been able to locate the mobility edge 共separating localized and extended states兲in amorphous polyethylene at about the vacuum level. The localized states extend to about ⫺0.3 eV. All the simulations presented in this article start from levels below the mobility edge. The article is organized as follows: in Sec. II we describe the methods and simulation details. In Sec. III we describe the mean-field approach we have used to carry out dynamical simulations. In Sec. IV we describe the self-trapped polaron. Finally, Sec. V provides a short summary and conclusions. II. METHODS We have used a mixed quantum-classical approach to study the dynamics of an excess electron in amorphous polyethylene at room temperature. The electron is treated quantum-mechanically while each methylene group is treated as a classical particle. A fast Fourier transform block Lanczos diagonalization algorithm6was used to compute the adiabatic electronic states of each classical configuration every time step. The time evolution of the system was generated on a Born–Oppenheimer potential energy surface using the molecular dynamics code DLPOLY.7The Hamiltonian of the classical subsystem contains all the standard ingredients used to simulate polyethylene 共bond stretch, valence angles, and dihedrals terms兲plus the interaction energy with the excess electron, which is computed using the Hellman– Feynman theorem.8Transitions between different electronic energy surfaces are allowed and computed by means of the Tully’s fewest switches surface hopping algorithm.9 We have simulated polyethylene systems employing a single chain in a cubic box with periodic boundary conditions. The chain was made of 360 CH2units, though we have also performed simulations with a chain of 1215 CH2units for comparison purposes. The initial equilibrium configurations of amorphous polyethylene at room temperature were generated using the procedure described in Ref. 5. The electronic wavefunctions are represented on a grid of 323or 463 points, depending on the system size. We estimate the uncertainty in the results due to grid size as 0.02 eV. All interactions were truncated at rc⫽9 Å, with all electronic energies including a long-range correction based on the polarizability interaction.5 JOURNAL OF CHEMICAL PHYSICS VOLUME 120, NUMBER 16 22 APRIL 2004 77720021-9606/2004/120(16)/7772/7/$22.00 © 2004 American Institute of Physics This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 150.214.182.194 On: Mon, 01 Jun 2015 15:17:42 III. MEAN FIELD APPROACH The interaction energy between the excess electron at r and the CH2units at Rican be decomposed into repulsive and attractive parts: V共r,Ri兲⫽兺 jVj r共 兩 r⫺Rj 兩 兲)⫹Vp共r,Ri兲,共1兲 where Vj r(r) is a short-range repulsive pair potential5and Vp is the electrostatic energy that accounts for the point chargeinduced-dipole polarization interaction: Vp共r,Ri兲⫽⫺ 1 2兺 jpj•Ej 共0兲Sj共 兩 r⫺Rj 兩 兲,共2兲 with Ej 共0兲⫽⫺eRj⫺r 兩 Rj⫺r 兩 3共3兲 the direct electric field due to the excess electron, and pj⫽ ␣ jEj共4兲 the dipole moment of the united atom jwith a polarizability tensor ␣ j.Sj(r) is the switching function that vanishes when r→0, accounting for the finite size of the methylene units. The local electric field at each atom Ejis the solution of the set of equations Ej⫽Ej 共0兲⫹兺 k⫽jT共Rj,Rk兲• ␣ kEk,共5兲 with T(Rj,Rk)⫽(3R ˆjkR ˆjk⫺1)/Rjk 3,Rjk⫽Rj⫺Rk, and R ˆjk ⫽Rjk /Rjk . Knowledge of the local electric field Ejallows one to compute the force exerted on each atom due to the polarization interaction as 具 ␾ 兩 Fj p 兩 ␾ 典 , where ␾ is the electronic wavefunction and 共Appendix A兲 Fj p⫽eT共r,Rj兲pj⫹兺 k⫽j 3 Rjk 4共共pj•pk兲Rjk⫹共pj•Rjk兲pk ⫹共pk•Rjk兲pj⫺5共pj•Rjk兲共pk•Rjk兲Rjk兲.共6兲 The first term in 共6兲is the force due to the direct field created by the electron 共which because of Newton’s third law is minus the force exerted by the dipole jon the electron兲and the second is the result of the interaction with all other induced dipoles. In Refs. 4 and 5 we solved self-consistently the set of equations 共5兲using an iterative approach for several configurations of bulk polyethylene. However, the calculation is very expensive and utterly prohibitive in a simulation of the dynamics of an excess electron in polyethylene since the interaction energy has to be computed every time step. In 1967 Lekner proposed10 a mean field approach to compute the polarization interaction of a point charge in a system of isotropic dipoles 共with a scalar polarizability ␣ 兲.In this approximation, the local field at each atom is replaced by Ej⫽Ej 共0兲f共 兩 Rj⫺r 兩 兲,共7兲 where f(r) is a function that accounts for the screening of the direct electric field exerted by the electron at the dipole j due to the other induced dipoles. Assuming 共7兲and summing the contribution of all induced dipoles in a statistical way, we obtain an equation for f(r):10 f共r兲⫽1⫺ ␲ n ␣ 冕 0 ⬁dsg共s兲s⫺2 冕 兩 r⫺s 兩 兩 r⫹s 兩 dt f共t兲t⫺2 ␪ 共r,s,t兲 共8兲 with ␪ 共r,s,t兲⫽3 2s2共s2⫹t2⫺r2兲共s2⫹r2⫺t2兲⫹共r2⫹t2⫺s2兲. 共9兲 In Eq. 共8兲nis the number density and g(r) is the pair correlation function of the dipole system, usually taken as the equilibrium pair correlation function in absence of the excess electron. Strictly speaking g(r) should be modified due to the presence of the electron; however, neglecting this perturbation has been shown to be very successful when applied to simple polarizable fluids.11 With this approach f(r) has to be solved numerically using an iterative self-consistent method, but only once and as a result we obtain a pairwise additive interaction energy that can be computed efficiently every time step. The first difficulty in applying this mean-field approach to polyethylene is that the polarizability for the methylene units is not a scalar but a tensor. In fact we will show later that the mean-field approximation breaks down for strong anisotropy in the polarizability tensor. Another difficulty arises when the system is made of dipoles of different polarizability, which again is our case, since our pseudopotential contains two different polarizability centers: one at the center of the CH2units 共in Å3兲, ␣ CH2⫽ 冉 0.5 1.12 1.64 冊 共10兲 共in a coordinate system along the symmetry axes of the methylene unit5兲, and another at the center of the C–C bonds ␣ C–C⫽ 冉 2.103 0 0 冊 共11兲 共along the C–C bond兲. These complications are difficult to overcome theoretically. We show below, however, that this is not necessary, since a naive application of Lekner’s mean-field method to polyethylene already provides good results. In Fig. 1 we have plotted the screening function f(r) 共solid line兲that results when we use Eq. 共8兲with a total average polarizability for the methylene units, ␣ ¯ CH2⫽1 3关Tr共 ␣ CH2兲⫹Tr共 ␣ C–C兲兴⫽1.788 Å, 共12兲 and the pair correlation function g(r) of the CH2units in amorphous PE at room temperature computed from a molecular dynamics simulations in absence of the excess electron. This radial distribution function only contains the pairs connected by a van der Waals interactions, which is consis7773J. Chem. Phys., Vol. 120, No. 16, 22 April 2004 Localized polarons in amorphous polyethylene This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 150.214.182.194 On: Mon, 01 Jun 2015 15:17:42 tent with the pseudopotential itself, since the polarizabilities were computed by ab initio methods for alkanes ignoring the intramolecular dipole–dipole interaction.5 In order to probe this approximation we have computed the effective screening functions from simulations of amorphous PE at room temperature using the full many-body selfconsistent approach 关Eq. 共5兲兴, defined for each dipole species as the average over all dipoles of the quantity pj•Ej (0)/Ej (0) • ␣ j•Ej (0) . The results corresponding to a single configuration of a PE chain containing 1215 CH2units are shown in Fig. 1. The screening function of the ␣ CH2centers 共dashed line in Fig. 1兲is very smooth, suggesting that a mean field approach would succeed for this type of dipole, where the anisotropy of the polarizability tensor is relatively weak. However, the results for C-C dipoles 共crosses兲are not smoothly varying and, even though the data seem to be scattered around the CH2screening function, the dispersion is large enough to reject a mean field treatment. All screening functions share common features. In the limit r→0, f(r) tends to unity, since the direct Coulomb interaction with the electron becomes dominant, and in the opposite limit r→⬁the screening function should approach 共in isotropic systems like amorphous PE兲the Lorentz local factor fL⫽1/(1⫹(8/3) ␲ ␣ ¯ )共see, for example, Ref. 12兲, since at long distances the direct electric field is viewed locally as a constant field. It can be proven that all solutions of the mean field Eq. 共8兲tend to this limit 共see Appendix B兲,as indeed does our numerical solution of the mean field equation in Fig. 1 共solid line兲. However, the screening function measured from the full many-body simulation appears to tend to a higher value 共dashed line兲. The reason for this deviation is that the local electric field was computed solving Eq. 共5兲for all dipoles inside a cutoff sphere of radius rc ⫽9 Å centered at the excess electron. Since the main contribution to the screening function at each dipole comes from the surrounding dipoles, the screening of the direct electric field due to the dielectric itself is not properly accounted for at dipoles close to the surface of the cutoff sphere. A simulation of the same amorphous configuration but with a larger cutoff radius rc⫽9⫹3 Å shows that the proper Lorentz limit is recovered 共dotted line in Fig. 1兲. This deviation is of little practical importance, since the wavefunctions of both systems are practically indistinguishable and the error in the energy levels is lower than 0.01 eV 共provided the same long range correction based on the Lorentz factor fLis added,5 which in this case is as large as ⫺0.41 eV兲. Apart from the large fluctuations in the C–C-bond screening function, the full-many body screening functions are in qualitative agreement with the mean field function, with the former clearly higher near the first minimum. This deviation is responsible for an increase in the absolute value of the polarization energy and thus a decrease of the energy levels. In Table I we present a comparison between the energy levels obtained with the full-many body simulation and the mean-field approach using the Lekner’s screening function plotted in Fig. 1 for both dipole species and the same configuration. It can be seen that all energy levels are shifted down by approximately the same amount of 0.15⫾0.02 eV 共with the error taken from our grid resolution error兲, and that the deviation in the first six wavefunctions is very small 共the projection of one normalized eigenstate on the other is close to unity兲. We have found the same pattern for different configurations and different system sizes. Since our interest is in the lowest energy levels, we will use this mean-field method in our simulations of the dynamics, correcting all energy levels by ⫺0.15 eV, which can be included in the long-range correction constant. Then, the force exerted by the excess electron in each atom due to the polarization interaction is best computed directly from the interaction energy, Vp⫽⫺ 1 2兺 jEj 共0兲• ␣ j•Ej 共0兲f共 兩 r⫺Rj 兩 兲共13兲 关where the switching factor S(r) has been included in the screening factor f(r)]. With Fj p⫽⫺ⵜjVp, where ⵜjis the gradient with respect to Rj, we obtain Fj p⫽eT共r,Rj兲pj⫹1 2共Ej 共0兲• ␣ j•Ej 共0兲兲ⵜjf共 兩 r⫺Rj 兩 兲.共14兲 IV. THE LOCALIZED SMALL POLARON An excess electron introduced into amorphous polyethylene is expected to interact with the dielectric and be trapped in a distorted region. In Fig. 2 we show the time evolution of the three components of the center of mass and FIG. 1. Screening function f(r). The solid line is the mean-field result. The dash and dotted lines correspond to the screening function for the electron– CH2pair computed from a full many-body simulation with a cutoff radius of rc⫽9 and 12 Å, respectively. The crosses correspond to the electron–共C–C bond兲pairs. TABLE I. Comparison of the energy levels between a full many-body selfconsistent and a mean field calculation. nEn MB En MF⫺En MB 具 ␾ n MB 兩 ␾ n MF 典 0⫺0.224 0.148 0.996 1⫺0.191 0.153 0.969 2⫺0.171 0.154 0.985 3⫺0.148 0.139 0.972 4⫺0.104 0.140 0.990 5⫺0.083 0.147 0.920 6⫺0.036 0.138 0.198 7⫺0.025 0.139 0.119 7774 J. Chem. Phys., Vol. 120, No. 16, 22 April 2004 D. Cubero and N. Quirke This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 150.214.182.194 On: Mon, 01 Jun 2015 15:17:42 the localization length of the excess electron in a microcanonical simulation starting from the ground state of an amorphous PE system of length L⫽32.12 Å at room temperature unperturbed by the presence of the excess electron. These quantities were computed under periodic boundary conditions following the procedure sketched in Appendix C. Note that at the beginning the electron is in a state below the mobility edge5and thus localized. However, this state is not stable and the electron explores part of the system, ending up after about 1 ps in a localized state with a smaller length. We present in Fig. 3 the time evolution of the energy spectrum for the same simulation. In this case the electron remains in the ground state, showing an adiabatic relaxation. We have observed similar behavior in all our simulations, though in some cases or if the initial state is not the ground state we witnessed larger relaxation times, of the order of 5–10 ps. In Fig. 4 we show the time evolution of the energy levels in a simulation where the excess electron was initially in the first excited state. The relaxation involves a series of energy hops, indicating a strong nonadiabatic interaction, eventually forming the self-trapped state after about 7 ps. We present in Fig. 5 the results with a smaller system box of length L⫽21.41 Å starting from the ground state. In this case the initial relaxation is not optimally mimicked, since the initial state, though localized, extends through most of the system, but the final self-trapped state is well characterized, showing the same localization length of about 5 Å. These smaller systems are much easier to simulate, reducing drastically the computer time used. In the following we will present results computed using this smaller system size. It is clear from Fig. 3 or 4 that, despite the fact that the kinetic energy is increased due to a reduction of the localization length, the interaction energy becomes more negative. This can only be explained in terms of a smaller repulsion interaction and/or a larger polarization interaction Vp. This is the reason these states are usually called self-trapping polarons. Part of the energy gained goes to the dielectric to create the distortion that makes possible this polaron state. The difference between the energy gained and the distortion energy 共both in absolute value兲should be still positive, since FIG. 2. Components of the excess electron center of mass 共lines, in Angstroms兲as a function of time from an initially unperturbed configuration of amorphous PE at room temperature. Crosses are the localization length of the electron at each time. System boundary edges are included as a reference 共horizontal dotted lines兲, though the actual position is arbitrary because of the periodic boundary conditions. FIG. 3. Electronic energy levels versus time for the same simulation shown in Fig. 2. FIG. 4. Electronic energy levels versus time for a simulation starting from an excited state showing a nonadiabatic relaxation. The black squares denote the current energy level. FIG. 5. The same as Fig. 2 but for a simulation with a smaller system size. 7775J. Chem. Phys., Vol. 120, No. 16, 22 April 2004 Localized polarons in amorphous polyethylene This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 150.214.182.194 On: Mon, 01 Jun 2015 15:17:42 we observe this self-trapped state in every simulation. However, this energy 共or more precisely minus this energy兲, normally called the self-trapping energy,13 is difficult to measure from a microcanonical simulation, since the total energy is conserved and thus it is transferred to the dielectric in the form of thermal energy. In order to compute the self-trapping energy it is more convenient to perform a simulation in the canonical ensemble. We have used the Andersen thermostat, which generates trajectories in the canonical ensemble.14 In this case, the quantity that is minimized in the self-trapping process is not the total energy but the free energy ⌬Gst . We have computed this free energy using the acceptance ratio method for classical systems, which is described in detail in Ref. 15. To sum up, we need to compute the constant Cthat satisfies the equation 兺 mh共U0⫺U1⫹C兲⫽兺 m⬘ h共U1⫺U0⫺C兲,共15兲 where h(x)⫽(1⫹exp ␤ x)⫺1is the Fermi function ( ␤ ⫽1/kBT) and U0and U1are the total energies of the classical systems U0and U1at the configurations m共in a simulation of system 1兲or m⬘共equilibrium configurations of system 0兲. The free energy difference ⌬Fis then given by ␤ ⌬F⫽⫺ln共n1/n0兲⫹ ␤ C,共16兲 where n1and n0are the total number of configurations sampled in each system. Since we are interested in the energy difference between the unperturbed ground state and the self-trapped state, we chose U0as the Hamiltonian of the polyethylene system without the excess electron, and U1as the full Hamiltonian with the electron hold at the ground state. Therefore, the self-trapping energy is given by ␤ ⌬Gst⫽F1⫺共F0⫹ 具 E0 典 0兲,共17兲 where 具 E0 典 0is the average of the electronic ground state energy over the equilibrium configurations of the unperturbed polyethylene system. Averaging over 5⫻2 simulations of 10 ps of adiabatic dynamics we obtained ⌬F ⫽⫺0.27⫾0.01 eV and 具 E0 典 0⫽⫺0.21⫾0.02 eV, which implies ⌬Gst⫽⫺0.06⫾0.03 eV. Moreover, the averaged localization length of the self-trapped polaron was 5.0⫾0.1 Å and of spherical shape, in agreement with the microcanonical simulations. It is worth noting that a direct calculation of the selftrapping energy as the total energy difference of the two states 共apart from neglecting the entropy change兲does not produce an acceptable value due to the large thermal fluctuations of each variable 共⬃0.1 eV兲. It is much more efficient to use an approach such as the acceptance ratio method, in which the energy difference of both systems is computed for each configuration of both simulations. We have observed a correlation between the selftrapping polaron and the local density of the dielectric. In all cases we found that the center of mass of the excess electron was within2Åoftheglobal minimum of the local atomic density 共defined as the number of atomic centers in a sphere of radius equal to the electron localization length divided by the volume of the sphere兲. Furthermore, this minimum local density was always smaller in the presence of the electron than the minimum local density in the unperturbed system. These changes, which in some systems amounted to a 50% reduction in local density, clearly indicate the formation of a cavitylike region by means of which the electron lowers its energy. In order to compare our results in amorphous PE with the self-trapped state found in crystalline PE3we have monitored the dihedral angle distribution and found no significant difference between the gauche population before and after the polaron is formed, reflecting the fact that in the disordered phase the electron has many more ways to distort the dielectric than in the crystalline phase, in which the electron can only create a cavitylike region by abruptly creating two gauche defects in the otherwise all-trans chains. In addition, we have observed that the self-trapped polaron is sensitive to the temperature of the system. Simulations at 350 K show a larger polarization energy, with a smaller localization length of about 4.2 Å, while a system at 200 K shows a more extended excess electron with a typical length of 6.4 Å. The calculations presented above predict that the selftrapping energy is of the same order as the thermal energy at room temperature (kBT⫽0.026 eV), and as a result we would expect rapid trapping and detrapping giving rise to hopping conduction assisted by phonons. This is in fact what we observe in the simulations. In Fig. 6 we present the components of the center of mass of the electron as a function of time for two typical long simulations showing a hop between self-trapping states at different positions. The lifetime of each state is about one order of magnitude larger than the relaxation time to the self-trapping states 共⬃ps兲, demonstrating the physical relevance of these self-trapped states, regardless of their small self-trapping energies. The time evolution of the electronic energy in these simulations shows that this is an adiabatic process, the electron staying in the ground state surface during the hop. The observed lifetime is in agreement with the widely accepted model for adiabatic hopping conduction in polaron theory,16 in which the lifetime is estimated as ␶ ⫽ ␯ ⫺1exp共Ea/kBT兲,共18兲 with ␯ being a typical phonon frequency responsible for the detrapping, and Eaan activation energy, here to be identified with the self-trapping free energy ⌬Gst⫽⫺0.06⫾0.03 eV. Since the optical frequencies in polyethylene are in the range 700–1600 cm⫺1,17 the corresponding phonon energies h ␯ are about one order of magnitude larger than Eaand thus only the acoustic or torsional 共transversal兲phonons, with frequencies in the range 0⭐ ␯ ⭐250 cm⫺1,18 can be responsible for the destruction of the self-trapping state. As a result, ␯ ⫺1 is of the order of a picosecond 共or larger兲, which is about the same time scale for the relaxation to these states, and the exponential factor in 共18兲provides the factor of 10 observed in the simulations. Furthermore, this analysis and the simulation results show that the dynamic behavior at the temperature considered is dominated by the low-laying energy states and thus 7776 J. Chem. Phys., Vol. 120, No. 16, 22 April 2004 D. Cubero and N. Quirke This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 150.214.182.194 On: Mon, 01 Jun 2015 15:17:42 providing a posteriori justification of the use of the Lekner’s mean field method, which describes correctly only the lowest energy states. V. CONCLUSIONS We have shown that even when a complex nonpolar dielectric like amorphous polyethylene with two species of anisotropic polarizable dipoles is considered, a simple mean field approach neglecting these details produces good results when one is only interested in the lowest energy levels of the excess electron. This method makes dynamical simulations feasible. Nonadiabatic simulations of an excess electron in amorphous PE at constant energy show the formation of a small self-trapped polaron with a reduced localization length of 5 Å on picosecond timescales. The wavefunction is isotropic and centered around a small cavitylike region created by the electron. In order to calculate the self-trapping energy we have performed adiabatic simulations in the canonical ensemble using the Andersen thermostat. The self-trapping energy is estimated to be ⌬Gst⫽⫺0.06⫾0.03 eV from the free energy change between an unperturbed dielectric system and a system perturbed by the presence of the electron. This is the energy required to remove the self-trapped state and is small when compared to the activation energy for excitation to extended states from unperturbed configurations, i.e., ⬃0.3 eV, thereby providing a justification of the method employed in Ref. 5, in which the effect of the excess electron on the dielectric was neglected. From Ref. 5, we expect electrons thermally excited to energy levels above the mobility edge to provide a contribution to the zero-field mobility of order of 10⫺3cm2/Vs, which is in good agreement with the highest values found in experiments. However, the smallness of the self-trapping energy suggests a hopping mechanism assisted by phonons, which is in fact observed in the simulations. The lifetime of each selftrapping state is observed to be on the time scale of a few tens of picoseconds. We show that this timescale is consistent with an adiabatic model for detrapping and the computed value of the self-trapping energy. Since the electron is localized and nondegenerate, the corresponding contribution to the mobility can be computed using the Einstein formula13 by measuring the diffusion coefficient from very long time simulations. Preliminary calculations show that the contribution to the mobility due to hopping between these self-trapped states may well be about the same order of magnitude as the mobility due to excited electrons above the mobility edge. These simulations are reported elsewhere.19 APPENDIX A: POLARIZATION FORCES In principle, the force on each dipole jcould be computed from the gradient of the polarization energy Vp. However, it is easier to calculate it from the local electric field. Assume a simple dipole pat r, then the electrostatic force on the dipole will be F⫽lim d→0,qd→p 关⫺qE共r兲⫹qE共r⫹d兲兴⫽共p•ⵜ兲E共r兲.共A1兲 If we now consider the full system with dipoles at Rjand an electron at r0, the electric field everywhere is given by E共r兲⫽E共0兲共r兲⫹兺 kT共r,Rk兲•pk,共A2兲 with E(0)(r)⫽⫺e(r⫺r0)/ 兩 r⫺r0 兩 3. Using these equations, together with ⵜE共0兲共r兲⫽eT共r,r0兲共A3兲 and ⳵ T共r,r⬘兲m,n ⳵ xl ⫽3 兩 r⫺r⬘ 兩 4 冉 共xl⫺xl ⬘兲 兩 r⫺r⬘ 兩 ␦ mn ⫹ 共xm⫺xm ⬘兲 兩 r⫺r⬘ 兩 ␦ ln⫹ 共xn⫺xn ⬘兲 兩 r⫺r⬘ 兩 ␦ lm ⫺5共xl⫺xl ⬘兲共xm⫺xm ⬘兲共xn⫺xn ⬘兲 兩 r⫺r⬘ 兩 3 冊 ,共A4兲 we obtain 共6兲. FIG. 6. The same as Fig. 2 but for two different long simulations showing a hopping mechanism between the self-trapped states. 7777J. Chem. Phys., Vol. 120, No. 16, 22 April 2004 Localized polarons in amorphous polyethylene This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 150.214.182.194 On: Mon, 01 Jun 2015 15:17:42 APPENDIX B: LONG DISTANCE LIMIT OF THE SCREENING FUNCTION The first step to find the long distance limit of the mean field f(r) is to prove ␦ 共t⫺r兲⫽3 8t⫺2 冕 兩 t⫺r 兩 t⫹rdss⫺2 ␪ 共r,s,t兲,共B1兲 where ␪ is given by Eq. 共9兲. Since for t⫽r 冕 兩 t⫺r 兩 t⫹rdss⫺2共r,s,t兲 ⫽共s⫺r⫺t兲共s⫹r⫺t兲共s⫺r⫹t兲共s⫹r⫹t兲 2s3 冏 兩 t⫺r 兩 t⫹r ⫽0, 共B2兲 it is only left to prove that it is properly normalized. Reversing the order of integration 冕 0 t⫹␧dt 冉 3 8t⫺2 冕 兩 t⫺r 兩 t⫹rdss⫺2 ␪ 冊 ⫽ 冕 0 ␧ds 冕 兩 r⫺s 兩 r⫹sdt䊐 ⫹ 冕 ␧ rds 冕 r⫺s r⫹␧dt䊐⫹ 冕 r 2r⫹␧ds 冕 s⫺r r⫹␧dt䊐 ⫽3 8 冉 共␧⫺r兲2共␧⫹5r兲 6r3⫺共␧⫹r兲共␧2⫹2␧r⫺11r2兲 6r3 冊 ⫽1, where 䊐⫽3t⫺2s⫺2 ␪ /8. If we reverse the order of integration in 共8兲, take the limit r→⬁, and use 共B1兲and g(⬁)⫽1, we obtain f共⬁兲⫽1⫺8 3 ␲ n ␣ f共⬁兲,共B3兲 and therefore f(⬁)⫽fL. APPENDIX C: CALCULATION OF THE ELECTRONIC CENTER OF MASS AND LOCALIZATION LENGTH WITH PERIODIC BOUNDARY CONDITIONS The center of mass of the excess electron and the localization length are defined from the wavefunction ␾ 共r兲as 具 ␾ 兩r兩 ␾ 典and 冑 具 ␾ 兩 r"r 兩 ␾ 典 ⫺ 具 ␾ 兩 r 兩 ␾ 典 • 具 ␾ 兩 r 兩 ␾ 典 , respectively. However, under periodic boundary conditions these quantities are not well defined mathematically even for localized states, and different values are obtained if the origin of the system is chosen so that the boundaries cross a region where the electron has a significant density probability. Following the localization criterion presented in Ref. 5, we have overcome this problem by choosing each origin component at the grid points so that the absolute value of the flux ⌽across the corresponding boundary is minimized, where ⌽⫽ 冕 dS ␾ n•ⵜ ␾ 共C1兲 and nis a unit vector perpendicular to the boundary surface. This procedure guarantees that wave function has decayed at the boundaries and the localization length is well defined. 1See, e.g., L. A. Dissado and J. C. Fothergill, Electrical Degradation and Breakdown in Polymers 共Peregrinus, London, UK, 1992兲and references therein. 2N. H. Ge, C. M. Wong, R. L. Lingle, J. D. McNeill, K. J. Gaffney, and C. B. Harris, Science 287, 288 共1998兲. 3S. Serra, S. Iarlori, E. Tosatti, S. Scandolo, M. C. Righi, and G. E. Santoro, Chem. Phys. Lett. 360, 487 共2002兲. 4D. Cubero, N. Quirke, and D. F. Coker, Chem. Phys. Lett. 370,21共2003兲. 5D. Cubero, N. Quirke, and D. F. Coker, J. Chem. Phys. 119, 2669 共2003兲. 6M. H. Gutknecht, Acta Numerica 共Cambridge U.P., Cambridge, 1997兲. 7W. Smith and T. Forester, J. Mol. Graphics 14, 136 共1996兲. 8K. Drukker, J. Comput. Phys. 153,225共1999兲. 9J. Tully, in Classical and Quantum Dynamics in Condensed Phase Simulations, edited by B. Berne, G. Ciccotti, and D. Coker 共World Scientific, Singapore, 1998兲. 10J. Lekner, Phys. Rev. 158, 130 共1967兲. 11 B. Space, D. F. Coker, Z. H. Lui, B. J. Berne, and G. Martyna, J. Chem. Phys. 97, 2002 共1992兲. 12N. W. Ashcroft and N. D. Mermin, Solid State Physics 共Saunders, Fort Worth, 1976兲. 13See A. L. Shluger and A. M. Stoneham, J. Phys.: Condens. Matter 5,3049 共1993兲and/or J. Appel. Solid Stat. Phys. 21, 193 共1968兲. 14H. C. Andersen, J. Chem. Phys. 72,2384共1980兲. 15D. Frenkel and B. Smit, Understanding Molecular Simulation 共Academic, San Diego, 2002兲. 16N. F. Mott and E. A. Davis, Electronic Processes in Non-crystalline Materials 共Clarendon, Oxford, 1979兲. 17S. F. Parker, J. Chem. Soc., Faraday Trans. 92,1941共1996兲. 18A. V. Savin and L. I. Manevitch, Phys. Rev. B 67, 144302 共2003兲. 19D. Cubero, N. Quirke, and D. F. Coker, in preparation. 7778 J. Chem. Phys., Vol. 120, No. 16, 22 April 2004 D. Cubero and N. Quirke This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 150.214.182.194 On: Mon, 01 Jun 2015 15:17:42