The excess electron at polyethylene interfaces †

This work investigates the energy and spatial properties of excess electrons in polyethylene in bulk phases and, for the first time, at amorphous vacuum interfaces using a pseudopotential single-electron method (Lanczos diagonalisation) and density functional theory (DFT). DFT calculations are made employing two approaches: with pseudopotentials/plane waves and the local-density approximation; and with all-electron Gaussian basis functions at the B3LYP level of theory, supplemented with a lattice of ghost atoms. All three methods predict similar spatial localisation of the excess electron, but a reliable comparison of its energy can only be made between the Lanczos and DFT using Gaussian bases. While Lanczos predicts that an excess electron would preferentially localise in nanovoids with diameters smaller than 1 nm, DFT suggests that it would localise on surfaces in nanovoids larger than 1 nm. Overall we conclude that in DFT studies of polyethylene/vacuum interfaces at the current level of theory, orbital-based methods provide a useful representation of excess electron properties. which represents the simulation cell with a vacuum gap of 1.0 nm in the x -direction. The chains are replicated above and below to illustrate the bonds across the periodic boundaries in the y direction. To build the lamellar chains we first place four chains of 20 carbons (–CH 2 ) inside the central cell. Then, we bond two neighbouring chains with four –CH 2 – groups in the y -direction at each of the two chain endings. At the y axis periodic boundaries we also add two 4 –CH 2 – groups to create a semi-infinite lamella. An equivalent bulk crystalline polyethylene system comprises four chains of 24 carbons as shown in panels (c and d) bonded across the periodic boundaries in the z direction.


I. Introduction
Electron trapping processes are central to energy transfer in nature 1 and to a wide range of technologies including photocatalysis, 2 photovoltaics, 3,4 organic thin film transistors, lightemitting diodes, 5 radiation damage and repair, 6,7 and in DNA based molecular electronics, 8,9 and electrical insulation. 10 In particular, electron trapping processes have been linked to dielectric breakdown in polymeric insulator shielding, 11 which is significant for high-voltage cables used for power transmission. Around 9% of the generated electricity is lost through transmission and distribution, 12 amounting to an estimated financial loss of 3 trillion U.S. dollars a year worldwide. In the light of these financial and technological challenges, the present work investigates the nature of interface trap states in polymeric materials 13 and their representation in single and multi-electron models. Our focus is on polyethylene, chemically the simplest organic insulator which is used in a wide range of applications from home electric devices to high-voltage cables. 14 Although dry low-density polyethylene has a very low electrical conductivity (10 À16 S cm À1 ), 15 this material still traps electrons (and holes) that are injected through contact with the conductor, forming a space charge of excess electrons. The traps are gap states caused by physical 16 and chemical disorder, 17 including interfaces and impurities. On the one hand, the electrostatic forces arising from trapped electrons, are thought to lower the energy barriers to local conformational change producing microvoids, 18 which eventually initiate the failure of polyethylene insulating materials through aging and dielectric breakdown. 11,19 On the other hand, the creation of new interfaces through the mixing of nanoparticles with polymers to produce nanocomposites has been reported to improve dielectric properties. 20 Therefore, we need a much better understanding of the behaviour of excess electrons at polyethylene interfaces to address key technological phenomena such as the aging and breakdown of insulators and nanodielectrics.
As a first step, this work studies the energy and localisation of excess electrons in polyethylene in bulk, and at interfaces with vacuum. These excess-electron properties have been calculated using Lanczos diagonalisation and single-electron semi-empirical pseudopotentials, which we refer to hereinafter as the 'Lanczos method'. The Lanczos method has been used to study pure systems such as alkanes, 21,22 rare gases, 21,23 and water. 24 In previous work we have used this pseudopotential approach to predict the density of states for excess electrons in bulk polyethylene, including the mobility edge between delocalised states and those states representing electrons localised (trapped) at nanovoids in bulk polyethylene with radii of up to 0.4 nm. 25 The Lanczos polyethylene pseudopotential was fitted to experimental data for the energy of excess electrons in alkanes (from a difference in work functions of electrons in a metal in vacuum and the same metal in the alkane fluid, 26 so with respect to a zero of potential in vacuum away from the metal surface) and ab initio data. 22 For materials with complex chemistries the Lanczos method is more difficult to apply since we would need a semi-empirical pseudopotential for each component. This limitation could be overcome by using, for example, all-electron density functional theory (DFT) methods. However, it is yet unclear how to determine excess-electron properties in DFT; an important question given importance of electron localisation physics, [27][28][29] chemistry, 30,31 and biology, 32 and the focus of the present work.
In this work, we use DFT and Lanczos calculations to investigate the behaviour of excess electrons at vacuum interfaces of polyethylene. This behaviour is evaluated as function of polyethylene morphology (lamellar and amorphous) and the vacuum gap d. We then compare the results obtained from two DFT methods and to Lanczos'. The first uses the CRYSTAL14 code 33 at the B3LYP 34 level of theory and LCAO Gaussian basis sets (in what follows the LCAO basis is assumed to comprise Gaussian functions). The second uses Quantum Espresso 35 with Becke88 and Perdew86 for the exchange-correlation (LDA + BP) and a plane wave basis. This comparison of Lanczos and DFT approaches provides insight into the question of how best to represent the excess electron in DFT.
In previous work using DFT, two Kohn-Sham orbitals have been routinely chosen to represent the excess-electron states. The more common choice is the lowest-unoccupied molecular orbital (LUMO) of the neutral N-electron system, 36,37 which assumes that this orbital's energy is that of an electron added to a neutral N electron system. However, in principle at least, the LUMO in DFT is unsuitable to describe any polarisation effects exerted by an added electron. The second choice is the highest-occupied molecular orbital (HOMO) of the charged N + 1 system. 38 This HOMO contributes to the ground state total energy and is correlated with the other N electrons (see Discussion in ref. 31). The HOMO of the N + 1 electron system can be rigorously identified in exact DFT with the electron affinity (A) of the N-electron system (cf. Fig. 1 in ref. 31 and Discussion therein). We however employ a hybrid functional so that it remains to be seen to what extent the HOMO(N + 1) provides a good approximation to the true excess electron energy.
Our manuscript is organised as follows. In Section II, we briefly describe the methods employed to build the bulk and interfacial systems and the parameterisation of the DFT and Lanczos calculations (the full methodology is given in the ESI †). In Section III, we present the degree of localisation of excess electrons in lamellar systems and compare with previous studies. We next show the degree of localisation and energies for the bulk and interfacial systems of amorphous polyethylene. Finally, Section V concludes with our main findings and implications for future work.

II. Methodology
We study two models of bulk polyethylene: a bulk crystal (C) composed of 4 chains with 24 carbons bonded across the periodic boundaries, and an amorphous polyethylene (A) phase built with four rings of 20 (or 40) carbons each with no terminal methyl groups. These amorphous configurations are prepared with Materials Studio using the COMPASS2 39 force field with no periodic boundaries in the x direction to avoid splitting the chains when the vacuum is imposed. The cubic simulation cell has a length of 1.29 nm and a density 40 of 0.86 g cm À3 . Our lamellar (L) polyethylene systems are composed of four chains with initially 20 carbon atoms in each one. Chains are folded in the x direction parallel to the surface normal vector of the yz plane by adding four methylene groups including across the periodic boundaries (see Fig. 1). The lateral dimensions of the systems are set to 0.986 and 0.74 nm in the y and z directions, 41 respectively, with a length of 2.68 nm in the x direction. Surfaces are created by separating both folded faces with vacuum of thickness 0.1, 0.5, 1.0, 1.5, 2.0, and 3.0 nm. The geometry of all systems is optimised with the ab initio package CP2k 42 in its QUICKSTEP 43 implementation. 44 The optimised configurations are then used as input for calculations with the codes Quantum Espresso and CRYSTAL14 which both use periodic boundary conditions. The Quantum Espresso plane-wave code uses local-density approximation (LDA) supplemented by Becke88 and Perdew86 for the exchangecorrelation (LDA + BP). Plane-wave calculations are restricted to neutral systems due to the difficulty with convergence for charged interfaces. 45 CRYSTAL14 was used to perform all-electron calculations at the B3LYP theory level. This exchange-correlation functional has given excellent estimates of band gaps 46 and optical gaps 47 of extended systems. We use a standard 6-31G** 48,49 basis set to represent the local atomic orbitals as primitive Gaussian functions. The calculations are made for systems with N and N + 1 electrons (with background positive charge) using the same simulation cell parameters and nuclear configuration. The optimised configurations from CP2k are also used as an input to run the Lanczos method. 22 Full details of our computational methodology are given in the ESI. † Representative configurations of these phases are shown in Fig. 1.

III. Results and discussion
We first calculate the electronic (fundamental) band gap for bulk polyethylene phases before reporting the effect of vacuum gaps in Fig. 2. For bulk crystalline polyethylene using CRYSTAL14 (B3LYP and Gaussian basis) we obtain a band gap of 8.79 eV. This value is in excellent agreement with the experimental result of 8.8 eV. 37,50 Other authors have also found good agreement with Gaussian basis sets. For example, in their study of all trans polyethylene with Gaussian sets and two local-density functionals, Miao et al. 51 obtained band gaps of 7.7 eV with Gaspar-Kohn-Sham (GKS) and 8.0 with Perdew-Zunger (PZ). In contrast, the band gap in our crystalline polyethylene system decreases to 6.22 eV using LDA + BP and plane waves (PW) in the Quantum Espresso package. This value is in agreement with that of 6.46 eV reported by Serra and colleagues 36,52 who employed LDA (BP and BLYP) and plane waves. More recently, Chen et al. 53 obtained a gap of 6.9 eV using plane waves and the PBE functional in the VASP code. In addition to these crystalline values, in the present work we provide gaps for the lamellar and amorphous phases. Forming a lamellar phase lowers the CRYSTAL14 band gap slightly to 8.2 eV while an amorphous phase lowers it to 8.09 eV (and 7.56 eV with ghost atoms, see below). In conclusion, we find that the band gap in non-crystalline polyethylene; lamellar and amorphous, is reduced by 6.7% and 8.0% with respect to the pure crystalline structure.
Creating free space above the polyethylene surface is similar to generating voids and defects in studies of semiconductors. Previous studies have found it necessary to use so-called ghost atoms in an augmented basis set to ensure that the electrons can penetrate this space to lower the total energy. [54][55][56] In our study we have added a lattice of ghost hydrogen atoms 57 throughout the simulation cell. These atoms represent points in space where the electronic orbitals are defined by Gaussian functions but have no electronic or nuclear charges. The construction of ghost-atom lattices is detailed in the ESI. † Fig. 2 shows that ghost atoms have a dramatic effect on the band gaps, lowering them by over an eV, and changing their variation with the slab separation. The band gap now falls with increasing slab separation rather than increasing. Clearly not using a basis augmented with ghost atoms restricts the excess electron in an unphysical way even for bulk amorphous polyethylene (see Fig. 3). The excellent agreement with the experimental band gap for the basis without ghost atoms must then be seen as fortuitous. We note in passing that an augmented basis is also required for fluid methane. In previous work 58 we showed that neither the LUMO(N) nor the HOMO(N + 1) DFT (CRYSTAL14) energy levels could represent the experimental  methane excess electron energy as a function of density, except for the highest liquid density. However, using an augmented basis DFT gives good agreement with the experiment over a wide range of densities. 59 Introducing a lattice of ghost atoms requires a new parameter to be set: the distance d g between ghost atoms. This approach is in contrast to using ghost atoms  in vacancy calculations, where a real atom at a position r is replaced by a ghost atom. One method of choosing d g would be to find the minimum in the total energy as a function of this distance. We find that the total energy decreases marginally (o0.001%) with d g but that for d g o 0.141 nm the code fails to converge, so total energy is not a useful function in this respect. However, the LUMO(N) and HOMO(N + 1) energies for the polyethylene/vacuum systems are relatively insensitive to d g for d g r 0.185 nm (see Fig. 3). We therefore have some freedom in choosing this parameter and we choose d g = 0.185 nm.
In Fig. 4 we plot the excess electron density for a representative bulk amorphous polyethylene configuration using each method: Lanczos ground-state, PW LUMO(N), LCAO HOMO(N + 1), LCAO LUMO(N), LCAOG HOMO(N + 1), and LCAOG LUMO(N), where G signifies the augmented basis with ghosts. These images show that the LCAO method without ghost atoms misplaces the excess electron. All other DFT methods localise the electron in agreement with Lanczos. An analysis of the overlap of the excess electron wave function C with the ground state Lanczos wave function C l , measured as the dot product hC|C l i for five configurations (see Fig. S1, ESI †) indicates that the LCAO LUMO(N) and the PW LUMO(N) show the best agreement with averages of hC|C l i = 0.65 and 0.67, respectively. Fig. 5(a) illustrates the average profiles of the charge density of the excess electron for the vacuum gaps formed by lamellar surfaces from the single electron Lanczos calculation. The extra electron prefers to localise away from the slab (bulk-like) regions for all separations, showing a transition in its crosssectional profile from small to larger separations. For separations larger than 2.0 nm the excess electron prefers to be located near the surface on the vacuum side with the wave function decaying exponentially towards the centre of the gap (normal to the surface). In the plane of the interface the charge density is delocalised. For slits smaller than 2.0 nm the charge density on opposing surfaces overlaps to give a parabolic profile across the gap. Note that the average profiles at d = 3.0 nm are asymmetric with respect to the origin of coordinates since the optimised geometry is not perfectly symmetric. In addition, Fig. 5(b) shows that the first three excited states localise near the interface for d = 3.0 nm, whereas the fourth excited state lies at the middle of the vacuum slab. Fig. 5(c) illustrates the average profiles of the plane wave LUMO(N) charge density for the neutral lamellar surfaces separated by d = 0.1, 0.5, 1.0, 2.0, and 3.0 nm. These profiles are very similar to those from Lanczos shown in Fig. 5(a) and (b), except the electron density peaks further away from the surface (around 0.5 nm at d = 3 nm). A significant difference is detected at d = 3.0 nm, where LUMO+2 and LUMO+3 are observed to form states extended throughout the vacuum, as shown by Fig. 5(d). Our LUMO+3 disagrees with that predicted by Righi et al. 36 since theirs resembles the conduction band of bulk lamellar polyethylene with the electron density residing in the bulk phase. Fig. 6(a) shows that for amorphous polyethylene the excess electron (from Lanczos) forms a surface state in a similar fashion to the lamellar phase at d = 3.0 nm. We have 65.3% of the excess electron in a region between 0.2 nm inside each surface and 0.7 nm into the vacuum, the peak in the electron density is at 0.2 nm above the polyethylene surface. In contrast to the lamellar phase, it is now the second excited state that is extended into the vacuum as shown in Fig. 6(b). Fig. 6(c) presents the cross-sectional averages of the LUMO of neutral systems calculated with plane waves. Despite the change in the surface morphology from lamellar to amorphous and the increase of the cross sectional area of the simulation cell from 0.729 (= 0.986 Â 0.740) nm 2 in the former to 1.674 (= 1.294 Â 1.294) nm 2 in the latter, the charge densities look similar with a small fraction inside the slab as d increases. The only qualitative difference is the complete loss of symmetry of the LUMO at 3.0 nm. However, looking at the LUMO+1, this state is only 0.03 eV more positive in energy (i.e. of the order of k B T, where k B is the Boltzmann's constant) and it is likely they both represent the same state, restoring the symmetry. Fig. 6(e and f) compares the cross-sectional averages, from CRYSTAL14 with ghost atoms, of the electron density of LCAOG LUMO(N) and HOMO(N + 1) for the amorphous surfaces at the B3LYP level of theory. The density of the LUMO(N) peaks around 0.2 nm above the interface for a gap of 3 nm in agreement with Lanczos and plane wave DFT while for the HOMO(N + 1) the peak is further out at 0.5 nm or more above the interface. The decay of the excess electron density into the vacuum is exponential for Lanczos and DFT 60,61 (see Fig. S2, ESI †) yielding localisation lengths l loc for the 4 nm vacuum gap of 0.24 nm for Lanczos, 0.31 nm for the PW LUMO(N), 0.31 nm for LCAOG LUMO(N), and 0.14 nm for the LCAOG HOMO(N + 1). We note that while the excess electron is delocalised in the directions parallel to the interface in our semi-infinite flat periodic surface model, in a real material the electron may well localise in y and z as well due to irregularities or the finite size of the surface.
We now focus on the variation of the kinetic, potential, and total energies of the excess electron in amorphous polyethylene with the vacuum gap d calculated using the Lanczos method, as illustrated in Fig. 7. In the bulk material the excess electron energy is À0.37 eV, agreeing with the much larger systems of bulk polyethylene studied in our previous work, 25 dropping as we add vacuum to a minimum of À0.53 eV at d = 0.5 nm and then increasing to plateau at around À0.2 eV for d 4 2.0 nm.  This implies that an excess electron would prefer to localise in the bulk phase or a small planar void, rather than in a vacuum gap wider than 1.0 nm. For the larger system of four chains with 40 carbons each, Fig. 7(a) shows that the excess electron energy is more negative for surface states than in the bulk with a value of À0.09 eV at d = 0 nm. The excess electron energies of the two systems are in good agreement except for the bulk where there are clearly system size effects which is to be expected as we have shown previously that the electron trapping is in naturally occurring nanovoids whose distribution will depend on system size and nanostructure. The minimum in the excess electron energy for planar vacuum gaps between 0.1 and 1.0 nm is also consistent with our previous study of the excess electron energy of spherical nanovoids in bulk polyethylene, although the energy of the planar gaps falls off more rapidly with d. The behaviour of the kinetic and potential energies, which plateau rapidly as the gap increases as expected as the excess-electron charge density is localised at the polyethylene surfaces on the vacuum side. The excess electron is not present away from the surface and hence, the kinetic energy is not a decaying function of Bd 2 , which would be the case if the wave function was extended throughout the vacuum, whilst the potential energy becomes much more positive with increasing d as the localised electron interacts only with the closest surface. To summarise, the Lanczos method implies that an excess electron in polyethylene containing nanovoids either formed naturally in the amorphous phase or due to cracks or other imperfections that are nanoscale in at least one dimension, would localise preferentially in spaces with gaps less than 1.0 nm either as a surface state (localised normal to the plane) on a planar surface or inside a nanovoid localised in three dimensions. Fig. 8 shows the energy of the LUMO(N), LCAOG excess electron compared to the Lanczos ground state in amorphous polyethylene. In order to make a meaningful comparison with Lanczos and with DFT data for different vacuum gaps, the LUMO(N) energies must be corrected to a common energy scale. We make this correction with respect to the zero of the electrostatic potential at the centre of the vacuum gaps, as it has been done previously in metal-organic frameworks (MOF). 62 The average Hartree potential at the centre of the vacuum is observed to approach a negligible value, both for LUMO(N) and HOMO(N + 1) LCAOG, as shown in Fig. S3 (ESI †). In the case of the HOMO(N + 1) the decay takes place further from the polyethylene slab, which is expected given that the Hartree potential in this case includes the contribution due to the excess electron. This implies that the electrostatic potential drop at the polyethylene surfaces is negligible, in full agreement with our DFT calculations reported in ref. 22 and with simple electrostatic considerations: there should be no dipole correction at all as long as the molecules are not drastically truncated at the surfaces, since the atoms that make up polyethylene have no permanent dipoles. In contrast to the Lanczos data, using the LUMO(N) energy as an approximation to the excess-electron energy E ee , implies that an excess electron would prefer to localise in vacuum gaps d Z 1 nm, rather than in smaller gaps. The LUMO(N) energy is small and positive for d o 0.2 nm where Lanczos is small and negative. Nonetheless, there is excellent agreement between the two approaches for d Z 2 nm with E ee = À0.2 eV within the accuracy of both calculations. The HOMO(N + 1) energies, agree very well with Lanczos for small gaps, but are around 0.5 eV more negative for large gaps. If we take the HOMO(N + 1) as representative of the excess electron then there is a slight preference for localisation in gaps in the range 1 to 2 nm. For the plane-wave data. The Hartree potential shows huge variation across the vacuum gaps and we cannot correct it to a zero of potential. Previous work using this approach has reported energy differences between energy levels, not absolute values. 36 Nevertheless the data are shown as a dotted line in Fig. 8, falling as approximately 1/d.

IV. Conclusion
We have compared various DFT methods of estimating the properties of an excess electron with data from the single electron pseudopotential approach using Lanczos diagonalisation. All methods localise the excess electron similarly (with an extended basis for CRYSTAL14) with more of less equivalent localisation sites and localisation lengths. Only the LUMO(N) and HOMO(N + 1) obtained from CRYSTAL14 can be corrected to a reference zero potential energy at infinity. The plane wave method used here employs core pseudopotentials which make it impossible to calculate the average Hartree potential reliably however higher order methods employed with plane waves such as GW (see for example 63 ) may give more reliable predictions.
Although the Lanczos method predicts a minimum in the excess electron energy for vacuum gaps between 0.1 and 1.0 nm at around À0.6 eV, DFT suggests that the excess electron would prefer to localise on surfaces with vacuum gaps larger than 1 nm. Both LUMO(N) and Lanczos agree that the energy of an excess electron in planar vacuum gaps larger than 2 nm is around À0.2 eV. The HOMO(N + 1) energies differ from those of the LUMO(N) and the ground state of Lanczos by a shift of nearly 0.5 eV for larger gaps. On the other hand, LUMO(N) energy differs from Lanczos and HOMO(N) for bulk polyethylene by a similar amount. While the differences between the two approaches are significant at this level of theory, the current work has not considered polaronic effects that could distort the local nuclear configuration and change the energy of excess electrons in the vacuum gaps. We note though that for the bulk polyethylene system the effect on localisation energies has been found to be small 64 of less than 0.1 eV. In Fig. 8, it can be seen that the Lanczos energies lie within the predictions of both methods.
Overall, the similarity between the LUMO(N) with the HOMO(N + 1) orbitals using a Gaussian basis with ghosts atoms and their agreement with Lanczos are encouraging and we conclude that in DFT studies of polyethylene/vacuum interfaces at the current level of theory, these orbital-based methods provide a useful representation of excess electron properties.

Conflicts of interest
There are no conflicts to declare.