{'fraction_non_alphanumeric': 0.027039587646288395, 'fraction_numerical': 0.006421902065993494, 'mean_word_length': 4.851668726823238, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Linear-scaling density functional theory simulations of polar semiconductor nanorods

Binary polar semiconductors in the wurtzite structure can be grown in the form of nanorods aligned along the wurtzite [0001] direction. In such structures, very large dipole moments have been observed experimentally. We have studied the distribution of charge in GaAs nanorods, so as to elucidate the origin of the dipole moments. To make contact with the realistic experimental regime, we need to model systems of several thousand atoms, necessitating the use of a Linear-Scaling formulation of DFT. We use the ONETEP code, which combines the benefits of O(N) computational effort and plane-wave accuracy. We show that both the direction and magnitude of the dipole moment of a nanorod, and its electric field, depend sensitively on how its surfaces are terminated and do not depend strongly on the spontaneous polarization of the underlying lattice. Furthermore, we observe that the Fermi level for an isolated nanorod always coincides with a significant density of electronic surface states at the polar surfaces, which are either mid-gap states or band-edge states. These states pin the Fermi level, and therefore fix the potential difference across the rod. We provide evidence that this effect can have a determining influence on the polarity of nanorods, and has consequences for the way a rod responds to changes in its surface chemistry.


Introduction
Nanostructured materials are becoming ever more important due to the possibilities they enable for control and optimisation of materials properties. Within this broad class of materials, semiconductor nanocrystals have a prominent role, with current and future applications in optoelectronics, thermoelectrics, and photovoltaics, and as pressure sensors and catalysts. The valuable electronic and optical properties of semiconductor nanocrystals are strongly connected to their size and shape [1,2,3,4,5,6,7,8]. Nanocrystals of binary semiconductors, particularly, though not exclusively, those with noncentrosymmetric crystal structures such as wurtzite, have been shown experimentally to exhibit large dipole moments [9,10,11]. These dipole moments strongly affect their internal electronic structure and therefore their optical properties, as well as their interactions with the environment and with other nanocrystals. However, the origin of these dipole moments, and indeed the factors influencing the charge distribution generally, are not well-understood.
Developments in electronic structure theory, particularly the availability of accurate quantummechanical simulation methods with linear-scaling computational cost [12,13,14,15,16,17,18,19,20,21], make it possible to simulate whole nanostructures of realistic size. This raises the possibility of simulating realistic nanoparticles and systematically varying controllable aspects such as sizes, shapes, compositions and surface terminations. The development of linear-scaling DFT methodologies which combine low-order-scaling with systematic, plane-wave accuracy, as in the onetep code [21,22], in particular, make it possible to disentangle the competing factors which determine the distribution of charge in such systems.
We have performed LS-DFT calculations [23], using the onetep code, of the ground-state charge distributions in wurtzite-structure GaAs nanorods of sizes comparable to those found in experiment. We showed that the dipole moment observed is strongly dependent on the surface termination. In particular, there is a correlation between the dipole moment and the energy gap observed between the position of the highest occupied state at one polar surface and the lowest unoccupied state at the other polar surface.
Subsequently [24], we extended these simulations by examining the behaviour of the dipole moment as a function of the length and cross-sectional area of the rod, and the composition of the rod as it varies between different III-V element pairs. We introduced an interpretation of these behaviours in terms of the concept of Fermi level pinning (FLP). FLP is known to occur in semi-infinite semiconductor surfaces when there are partially-occupied states at the Fermi level, and can be expected to occur when the dipole moment of a rod is large enough that the resulting electric field adjusts the Fermi level enough that the potential difference along the rod is comparable to the bandgap. We showed that a finite-surface version of FLP then plays a crucial role in determining the polar characteristics of nanorod systems.
In the current article, we present two new features: firstly, we present a detailed description of the methodology of these simulations (Section 2), which could not be included in the previous works due to limitations of space, and secondly, we revisit the results of our earlier work [23] on variation of dipole moment with surface termination, to examine how they can be interpreted in terms of the concept of Fermi level pinning (Section 3).

Electronic Energy Minimisation
This work uses linear-scaling density-functional theory (LS-DFT) as implemented in the onetep code [21,22]. We represent the single-particle density matrix, which contains a complete description of the Kohn-Sham system, using a set of non-orthogonal, strictly localized orbitals where the density kernel K αβ is a generalisation of occupancies to a nonorthogonal representation. The orbitals {φ α (r)} are referred to as Non-orthogonal Generalised Wannier Functions (NGWFs). Each of these NGWFs is expressed using an underlying basis of periodic cardinal sine, or 'psinc' functions, D i (r), which are effectively representations of a delta function in a bandwidth-limited plane-wave basis. Each function D i (r) = D(r − r i ) is centred on a grid point r i , and the NGWFs are then written φ α (r) = where c iα are coefficients which are non-zero only for grid points r i within a sphere of radius R φ around the atom on which NGWF α is centred. The NGWFs are thus strictly localised within this sphere.
The spacing of the real-space grid is determined by an energy cutoff E cut which determines the kinetic energy of the largest G-vectors, as in plane-wave pseudopotential approaches. In this work, E cut = 400 eV was employed throughout. Four NGWFs per atom were used for Ga and As (initialised to s−and p-like atomic orbitals), and one NGWF per atom for H (initialised to s). Large NGWF radii of R φ = 10 a 0 were employed to ensure well-converged electronic structure and accurate forces. The ion cores were represented using norm-conserving pseudopotentials with three and five valence electrons respectively for Ga and As. The non-linear core correction approach was used to include core-valence interactions [25], which are important for both Ga and As. The local density approximation was used for exchange and correlation. Due to the strict localisation of the NGWFs, the density n(r) = αβ φ α (r)K αβ φ β (r), the overlap matrix S αβ = φ α |φ β and the Hamiltonian H αβ = φ α |Ĥ|φ β , can all be constructed in linear-scaling computational effort, through the use of the 'FFT box' approximation [26]. Minimisation of the total energy E T = αβ K αβ H βα −E dc [n] then proceeds via a double loop: the outer loop optimises the spatial form of the NGWFs φ α (r) via a conjugate gradients algorithm, using the energy gradient ∂E T /∂φ(r) to construct the search direction. The inner loop optimises the matrix elements of the density kernel K αβ . The NGWF optimisation has the effect of producing an in-situ optimised, minimally-sized, local orbital representation with an accuracy demonstrably equivalent to that of a plane-wave basis [27]. Sparse matrix algebra operations are performed using an atom-blocked scheme with sparsity determined by the overlaps of the local orbitals and by a real-space cutoff for the density kernel. Efficiently parallelised sparse matrix algebra routines [28,29] mean that the calculations scale well to hundreds of parellel cores: most of the current results were performed using 256 Intel Core i7 cores (64 processors).
The kernel is optimised using a combination of methods [30] including LNV [31] and penalty functional minimisation [32], designed to ensure that the energy is minimised with respect to K αβ subject to the requirements that the density matrix is correctly normalised, in that Tr[ρ] = N e , the number of electrons, and that it idempotent, in thatρ 2 =ρ. This has the effect of ensuring integer occupancies of zero or one for all single-electron eigenstates of the Hamiltonian, which is appropriate in an bulk insulator with a clearly defined bandgap or in an isolated system with a clear HOMO-LUMO gap. However, in the nanorod systems studied in this work there is often no significant HOMO-LUMO gap in the full density of states, which at first sight might question the appropriateness of a scheme employing integer occupation numbers. Fortunately in these cases there is always a clear gap in the local density of states everywhere in the rod, and the HOMO and LUMO have vanishing overlap due to being localised at opposite ends of the rod.
In these nanorods with widely separated polar surfaces, a density matrix corresponding to integer occupancies of the polar surface states is then more appropriate than one with 'smeared' metallic occupancies. In Ref. [24] we investigated the effect of smearing of occupancies around the Fermi level by performing test calculations on a few smaller rods with explicit diagonalisation of the Hamiltonian, followed by occupation of the states according to a Fermi-Dirac distribution with a self-consistently determined chemical potential. In all cases tried, while there was minor variation of the total energy, due to the slightly different occupations, there was not any major variation of the dipole moments observed (under 5% change in all cases).

Electrostatic Interactions
The strict localisation of the NGWFs means that in a large simulation cell, regions of vacuum at a significant distance from all the atoms (greater than R φ ) are not within the localisation region of any NGWFs. This means that apart from the 'whole-cell' operations such as calculation of the Hartree and XC energies, little extra computational effort is associated with the description of extra vacuum in the simulation cell. In a periodic simulation, it is therefore simple to pad a nanostructure with sufficient vacuum to prevent the wavefunctions of periodic images from actually overlapping.
However, the long-range of the Coulomb interaction, combined with the large dipole moments encountered in these systems, means that despite extra vacuum being free, it is still not feasible to converge electrostatic finite size effects with respect to simulation cell size in a periodic simulation. We therefore truncate the Coulomb interaction with a real-space cutoff [33], of the general form: where Θ(r) is the Heaviside step function. In practice, it makes more sense in these elongated systems to employ a cutoff on the surface of a cylinder rather than a sphere. The method for achieving this has been detailed elsewhere [34]. The truncation of the Coulomb interaction automatically negates all multipole terms in the interaction between different nanorods. Of particular importance is that it prevents dipole-dipole term, which would otherwise act to damp the magnitude of the dipole that would be observed, as seen in Ref. [34].

Geometry Optimisation
The structures of the nanocrystals were determined using geometry optimisation with the BFGS algorithm. Positions were initialised to atomic positions as cleaved from the perfect crystal structure (itself determined by geometry optimisation of the comparable bulk wurtzite structure using the castep plane-wave DFT package [35]). Initial positions of the surface terminating atoms were determined by choosing a typical bond-length perpendicular to the surface. The non-centrosymmetric nature of the Wurtzite structure means that surfaces perpendicular to the c-axis necessarily exposes a plane of Ga atoms on one polar surface and a plane of As atoms on the other. We will refer to these as the 'Ga surface' and 'As surface' respectively, regardless of the ligand species attached (eg H). The rods relaxed in this process are of length ∼ 3.5 nm and width ∼ 1.2 nm.
We employ a shorthand notation to denote the possible types of terminating species. In this notation, 'H' is hydrogen, 'B' is a bare (unterminated) surface, and 'P' is 'pseudo-hydrogen', is a hydrogen-like nucleus with a variable nuclear charge to account for terminating ligands of varying electronegativity: (5/4)e neutralises the Ga surfaces, and (3/4)e neutralises the As surface [36]. The first letter indicates the lateral surface termination, and the letter after the slash indicates the polar surface termination.
For large nanorods this generated initial configurations with relatively low forces which relaxed in around 50-100 BFGS steps depending on the surface termination. Convergence was defined as when all forces were relaxed to under 0.1 eV/Å, positions were converged to within 0.005Å and the energy gain per step was under 0.03 meV/atom.

Local Density of States
One of the most useful means to understand and interpret the information provided by an electronic structure calculation is a description encompassing both the spatial and energetic distribution of the single-particle eigenstates. In the local-orbital framework provided by ONETEP this is achieved by performing a single diagonalisation of the Hamiltonian matrix in the NGWF basis, at the end of the calculation once self-consistency has been achieved. While this comes with a O(N 3 ) computational cost, the prefactor is low because the NGWF basis is so small (four NGWFs per atom for Ga and As). Therefore, this diagonalisation remains feasible up to system sizes of around 10000 atoms, particularly if a parallel eigensolver such as ScaLAPACK is used. The generalised eigenproblem to be solved is of the form The matrix M β n describes the eigenvectors, which take the form |ψ n = β M β n |φ β . Eq. 2 can be solved using the ScaLAPACK routine PDSYGVX, to obtain the eigenvalues { n } and where in practice we replace the delta function with a Gaussian broadening of 0.1 eV. The local density of states in a given region I is calculated by including a projection of each eigenstate onto the local orbitals of region I, as The non-orthogonality of the projectors has been taken into account by raising the index of one of the terms, by implicitly constructing the contravariant dual of the NGWF (as in [37,38]). The functions |φ α need not be explicitly constructed in real space, since the relationship φ α |φ β = δ αβ implies that we can re-write Eq. 4 as follows: For these nanorod systems, a natural choice of regions over which to decompose the DOS is sets of slabs of atoms along the c-axis. Atoms are thus decomposed into groups according to their position along c, and a separate LDOS curve obtained for each group (as in Figure 2).

Calculating Dipole Moments
The total polarisation of an isolated nanostructure can be expressed as the sum of the ionic and electronic polarisation. The ionic polarisation is simply The electronic polarisation d el , meanwhile, is given in an isolated system by To evaluate this in the local orbital framework, we first evaluate the matrix elements of the position operator in the NGWF representation, R αβ = φ α |r|φ β , within the FFT box (adjusting for the origin of the FFT box so as to give a result relative to the origin of the simulation cell). Then we write Eq. 7 as indicating that the polarisation can be calculated directly from the density kernel and the NGWFs, without requiring any knowledge of the individual eigenstates. This allows dipole moments to be calculated easily for extremely large systems. The end result is the total polarisation which for a charge-neutral system is independent of the origin used for the position operator as long as it is the same in both terms.  Figure 1 shows the atomic configurations before and after relaxation for four different choices of surface termination. It can be seen that in the case of H-termination on the polar and lateral surfaces, the unrelaxed structure is not significantly altered, except that the As-H and Ga-H bonds align with where the next atom in the bulk structure would be, to restore the tetrahedral configuration of each Ga and As site, as would be expected. The differing electronegativities of Ga and As also mean that the Ga-As bonds in the surface layer tilt somewhat. In the cases with bare surfaces, there is much more significant relaxation. When the lateral surfaces are unterminated, a major distortion of the surface layer occurs so as to produce flat, threefold coordinated layers of Ga and As atoms, with one ridge occurring per primitive cell in the c−direction. This has the effect of partially 'healing' the local bandgap, restoring overall insulating behaviour. A bare polar surface is seen to be similarly energetically unfavourable: surface atoms of both species are again flattened into planes, so as to result in charge-neutral surface. Excess atoms appear to have a tendency towards being expelled from the structure, although the relaxation scheme employed here does not actually cause them to do so in practice. Clearly, in some cases, this suggests that the structures considered here will not be stable. A semi-infinite polar surface with similar termination would be unstable to reconstructions involving changes in the composition of the surface (addition or removal of atoms). At some size, therefore, one would expect the same thing to happen even for a finite-sized nanorod. This possibility is not accounted for within the current model, but remains an important issue for future work to investigate.   Figure 1, in particular the dipole moments of the table, suggest that to obtain a qualitative understanding of the dependence of the dipole moment on the surface termination it is not necessary to relax the surface. In fact, the variations in electronic structure between rods with different surface species can best be understood by considering very large unrelaxed nanorods with pristine unreconstructed surfaces. Figure 2, presents the local density of states of four simulated nanorods of length ∼ 12 nm and width ∼ 2 nm, varying in the surface species with which dangling bonds are terminated. It shows that, as might be expected, the presence of a strong electrostatic field shifts the peaks of the density of states up in energy at the negatively-charged end and down at the positively charged end. The most notable feature, discussed but not fully explained in Ref. [23], is that the dipole moment itself clearly does not depend strongly on the underlying polarisation of the wurtzite structure, but depends very strongly on the surface termination of both the polar and lateral surfaces. This dependence can be sufficient to reverse the direction of polarisation (hence the  Table 1. Dipole moments (in Debye) of the four types of rod whose LDOS are shown in Fig. 2. oppositely-directed electric fields in Figs 2(b) and 2(d) compared to 2(a) and 2(c)). Table 1 shows the variation of the dipole moment between the four types of rod. This variation cannot be explained simply in terms of the differing ionic charges of the terminating species: arguments based on counting ionic charges and the charge associated with Ga-As, Ga-H and As-H bonds can be used to estimate roughly what change in polarisation might be expected to result from a given change in surface termination, in the absence of electronic or ionic polarisation (i.e. what would be obtained simply by cleaving the bulk structure without allowing electron density to relax). These uniformly predict both that changes in the dipole moments of the rods would be considerably larger than those observed, and do not explain even qualitatively the variations seen in Fig. 2.

Variation of Dipole Moments with Surface Termination
Examining the different figures, the most obvious effect the variation of surface termination has on the electronic structure is through the position and number of surface states in the bandgap. It can be observed that in all cases, the polarisation and the corresponding voltage drop along the rod are such that states adjacent to the valence band edge at one polar surface are aligned in energy with states adjacent to the conduction band at the other polar surface. However, this occurs in a different way depending on the nature of the surface termination: in the rods with H-terminated polar surfaces, the surface-induced states specific to the polar surfaces (i.e. in the top and bottom curves) consist of redistribution of weight in the DOS around the valence band maximum on the Ga surface and around the conduction band maximum on the As surface. Alignment between these states can thus be achieved by an electric field which raises the potential at the Ga surface and and lowers it at the As surface.
Conversely, for the rods with bare (unterminated) polar surfaces, the surface-induced 'dangling-bond' states specific to the polar surface are adjacent to the conduction band on the Ga surface and the valence band at the As surface, and are much more strongly peaked around a specific energy. This means an electric field of opposite sign is necessary to bring the surface states into energetic alignment.
When one examines the rods without lateral-surface termination, Figs. 2(c) and 2(d), one sees a significant density of states associated with the dangling bonds on the lateral surfaces adjacent to the valence band maximum: these have the effect of almost closing the bandgap.
In summary, therefore, one can explain the trends of Table 1 in terms of an argument based on pinning of the Fermi level as follows: exposing layers of Ga on one face and As on the other necessarily creates a charge imbalance between opposite polar surfaces, and thus a dipole forms. However, if this dipole moment were too large, unoccupied electronic states would be pushed below the Fermi level on the positively charged end and raised above it on the negatively charged end. This would have two effects: a higher-energy configuration overall, and a reversal of the dipole moment due to the charge transfer. However, the bulk of the rod is itself polarisable, so instead charge flows to oppose the field by polarising the bonds along the rod axis until the positions of the filled and empty states are aligned at opposite ends of the rod. The dipole moment is thus limited in magnitude by Fermi level pinning, which effectively ensures that the in the final electronic configuration, the voltage drop between the ends of the rod cannot exceed the difference between the energies of the occupied states at one end and the unoccupied states at the other end.
The differences in sign and magnitude of the dipole moments can be explained in terms of the effects of the surface species on the local density of states at the surfaces: bare lateral surfaces significantly close the bandgap, resulting in a lower dipole moment for H/B and B/B compared to H/H and B/H respectively. Meanwhile, going from H-termination to bare surfaces at the polar ends switches the location of the occupied states, moving them from the Ga end to the As end, and vice-versa for the empty states. Therefore the sign of the dipole moment switches from negative to positive, while the smaller gap resulting from the position of the dangling-bond states inside the local gap at each end ensures the magnitude reduces.