Molecular statics simulation of CdTe grain boundary structures and energetics using a bond-order potential

The structure and energetics of coincidence site lattice grain boundaries (GB) in CdTe are investigated by mean of molecular statics simulations, using the Cd–Zn–Te bond-order potential (second iteration) developed by Ward et al (2012 Phys. Rev. B 86 245203; 2013 J. Mol. Modelling 19 5469–77). The effects of misorientation (Σ value) and interface plane are treated separately, complying with the critical need for full five-parameter characterization of GB. In addition, stoichiometric shifts, occurring between the inner interfaces and their adjacent atomic layers, are also predicted, revealing the energetic preference of Te-rich boundaries, opening opportunities for crystallography-based intrinsic interface doping. Our results also suggest that the intuitive assumption that Σ3 boundaries with low-indexed planes are more energetically favorable is often unfounded, except for coherent twins developing on {111} boundary planes. Therefore, Σ5, 7 or 9 boundaries, with lower interface energy than that of twin boundaries lying on different facets, are frequently encountered.


Introduction
During the past decade, cadmium telluride (CdTe) solar cells made the transition from an originally promising material into a real and mature commercial technology. Record laboratory cell efficiencies are still increasing, reaching 22.1% in 2016 [1], whereas full size modules overpassed 18%. Nevertheless, a lot of performance improvement paths can still be explored, one of them being in the direction of the open circuit voltage. Recently, the 1 V barrier has been overcome in doped single crystals [2]. This value remains, however, relatively low when compared to the CdTe bandgap (1.5 eV). This is mostly due to the low carrier concentrations achievable in the absorber layer, resulting from the difficult doping of CdTe, and due to parasitic recombinations of the photo-generated carriers. The tight link existing between carrier density and crystallographic defects emphasizes the crucial role of the latter on cell performance. In particular, grain boundaries (GB) have been proven to significantly alter the charge transport in CdTe layers [3][4][5], although their actual effects remain equivocal. It is likely that a better understanding of their structure and formation mechanisms will help uncovering their role in this context. Molecular dynamics and molecular statics (MS) simulations are pertinent tools to resolve the structural and energetic properties of defects, as well as to capture the growth features of thin films. The accuracy of such simulations relies on a careful parameterization of the considered atomistic potentials and an adequate choice of the probed time scales. In the case of CdTe, Ward and Zhou conducted detailed studies in this direction [6][7][8] and developed a bond-order potential (BOP) that has been found to yield particularly accurate results in term of vapor deposition [8] and defect energy trends, such as dislocation line energies [9], or vacancy and interstitial energies [7]. BOPs encompass a wider range of analytic interatomic potentials which adopt Pauling's [10] and Coulson's [11] bond-order concept. This means that they express the binding energy of a system as a sum over all atomic pairs considering a so called bond-order term which is a function of the local environment of an atom [12][13][14].
Using such a bond-order term introduces the advantage over conventional interatomic force field formulations in that they are capable of describing different bonding states of an atom which allows them to access chemical reactions and to include effects associated with directional, topological and energetic bond anisotropy.
The development of BOPs for the CdTe system marks a major breakthrough, enabling to conduct meaningful atomistic simulations without the need for computationally very costly ab-initio, i.e. density functional theory (DFT), calculations. BOPs thus allow performing static and dynamic simulations involving large numbers of atoms. In the present work, these CdTe BOP potentials are employed to model the structure of a specific set of coincidence site lattice (CSL) GB, which were experimentally observed in CdTe layers [15], at the atomic scale.

MS simulation setup and parametrization
MS simulations were performed using the LAMMPS package [16] and the Cd-Zn-Te potential (second iteration) developed by Ward et al [6,17]. The Ovito software [18] was used for visualizing the simulation results and for first order analysis of the data. GB structures were created by successive relaxations of two adjoining crystallites, with controlled orientations. In every case, the total number of atoms included in the simulation was in the order of ∼5×10 4 . The rotational parameters of the interface, that is, boundary misorientation and interface planes, were thus initially defined, whereas the energy minimization step was in charge of relaxing the translational parameters. In order to probe the complete energy space, and not only one specific local energy minimum, the translational parameters were actually pre-constrained prior to each relaxation step. This was performed by additionally translating one side of the boundary into a plane parallel to the interface. Therefore, one relaxed structure was created as output per translation vector. The energy minimization step was done in two steps, comprising a standard energy relaxation procedure, by adjusting atom coordinates, combined with a relaxation of the simulation box to release the stress accumulating in the periodic directions. Periodic boundary conditions were only enforced in the directions coplanar with the interface, thus creating a free surface standing opposite to the GB. This prevents interatomic stresses from building up in this direction and circumvents the issue of GB polarity reversal, which would have had otherwise occurred along this direction. Both crystals were created large enough so that this free surface does not affect the GB energetics. This was confirmed by verifying that the atomic potential energy reaches a plateau after the first ≈10-20 Å away from the surface, thus mimicking a bulk value, as visible in figure 1 and as will be later discussed in greater detail.
The polarity issue originates from the fact that CdTe crystallizes in the non-centrosymmetric F m 43 space group. This leads, for example, to the existence of two chemically distinct 111 ( ) and 111 ( ) planes, the former being exclusively composed of Cd atoms, whereas the latter has a monoatomic Te composition. As a result, the {111} plane family is said to be polar. Therefore, for a given set of rotational parameters, the initial number of atomic configurations to be considered before relaxation strongly depends on the polarity of the boundary planes, as emphasized in table 1. Coupling the possible starting configurations, summarized in table 1 and the atomic translation grid used for relaxation, resulted in up to around 800 different starting configurations which were relaxed for every set of rotational parameters. In the following, the reported quantitative values were arbitrarily averaged over the ten configurations with the lowest energies. Since the function minimizing the energy of different starting configurations is surjective, that is, different initial structures can output similar relaxed boundaries, the reported standard deviations can be zero in the case where all of these ten configurations are identical.
In order to extract quantitative information from the GB region, it is necessary to provide a clear definition of the interface at the atomic level. In the present work, the differentiation between bulk and GB atoms was established based on their coordination geometry. This was performed by adopting the diamond structure identifier available in the OVITO code [18], derived from the conventional common neighbor analysis method, and described in [19]. As such, all particles having at least a first or a second neighbor not positioned on a zincblende lattice site were assigned to the GB. Based on this distinction, three different parameters, namely the GB energy γ, boundary stoichiometry σ, and boundary atomic density ρ were attributed to each interface. The first one, γ, quantifies the extra energy of the interface, relative to the bulk. It is usually obtained by comparing the energy of a system containing a GB, with the energy of a bulk single crystal. However, in the case of zincblende structures, the issue of polarity inversion prevents the use of periodic boundary conditions in the direction perpendicular to the interface, as aforementioned. Therefore, only the energy of the GB plus the energy of the free surfaces thus formed can be directly obtained. As such, it is important to evaluate the contribution of the free surfaces to the relaxed system energy. The energy profiles of a system containing a GB and two free surfaces, on the one hand, and of a single crystal, on the other hand, are displayed in figure 1, in the direction perpendicular to the Plane interfaces. It clearly appears that the long range elastic contributions of the free surfaces are negligible, and only affect the first ∼5 Å next to the crystal edges. The effect of the GB itself is, however, much more pronounced, and covers around 20 Å. In between these interfaces, the average energy per atom has a constant 'bulk' value, which is found to be the same as what is obtained in a single crystal simulation. The GB energy γ can thus be obtained in a conventional way by comparing the energy of the system encompassing the GB, but excluding the first atomic layers next to the surfaces, with an identical volume in a single crystal. Normalizing this value by the interface area provides a GB energy in J m −2 .
Since distorted bonds and deviations from the zincblende coordination strongly affect the electronic band structure of semiconductors, the atomic number density of the boundary ρ is introduced. It quantifies the absolute number of boundary atoms, normalized by the interface area. This leads to values of ρ in units of mol m −2 (for the sake of convenience, most of the results are actually given in unit of density, u.d., with 1 u.d. being equal to 10 −5 mol m −2 ). It may be useful to also introduce a relative indicator, i.e. a parameter quantifying only the excess number density of distorted atoms, as compared to the initial atomic density of the crystallographic plane composing the boundary. To this end, the dimensionless variable P is introduced and defined as follows in equation (1): with ρ the atomic number density of the relaxed boundary, as-defined above, and , r ¢ the density of the atomically sharp A, respectively B, crystallographic planes forming the ΣX A||B interface. The use of the ρ′ denomination is employed to stress that these values correspond to non-relaxed structures, whereas ρ requires an energy minimization step. Therefore, P can be envisaged as a measure of the sharpness of the relaxed interface, where higher values of P indicate boundaries extending over several parallel atomic planes. Moreover, the critical case, namely, when P equals 1, corresponds to GB formed by pure juxtaposition of bulk crystals with minimum relaxation needed. The atomic number densities of low-indexed boundary planes of the cubic sphalerite space group are tabulated in table 2.
Finally, the boundary stoichiometry σ measures the Cd/Te atomic ratio in the distorted area. It is, therefore, a dimensionless variable. It should be noted that, although a stoichiometric shift can be observed in the boundary, it is in any case compensated by an opposite enrichment/depletion in the boundary surroundings, since the overall structures have a perfect one to one Cd/Te ratio. As such, σ only measures the chemical variation between the inner interface and its vicinity.
Even though GB are not planar defects, but rather have an associated defective volume, values are normalized by the interface area to compare simulation data with experimental measurements.

Σ3 interfaces
Due to the remarkably low stacking fault energy of CdTe of around 9 mJ m −2 [20], a significant fraction of GB observed in these layers are Σ3 twin boundaries, namely, around 60% of the total interface network length [15]. The Σ value refers to the reciprocal fraction of coinciding atomic sites in the superlattice formed by the two adjoining crystallites with respect to all atoms in the lattice. Therefore, GB exhibiting this misorientation relationship, thereby defining three out of five rotational parameters, are expected to be composed of small Quantitative data obtained on a broader variety of Σ3 boundaries are presented in figure 3. As can be observed, the energy of the coherent Σ3{111}||{111} interface calculated with the BOP is almost zero. This is due the tendency of the potential to underestimate the stacking fault energy of CdTe as already reported by Zhou et al [9]. The energy of coherent twin boundaries can be generally well approximated by taking half of the stacking fault energy, leading to an expected value of around 5×10 −3 J m −2 . This is, in any case, well below the values calculated for any other boundary planes. The reported data are hardly comparable with other studies, as there is a lack of publications dealing with GB energy calculations in the CdTe literature. One of the only incoherent Σ3 boundaries, being extensively studied, is the {112}||{112} boundary. In this work, the calculated energy of this interface is around 0.53 J m −2 which is slightly lower, but still in good agreement, with the value of 0.6 J m −2 reported by Park et al [21] from first-principle DFT calculations. Furthermore, both, the present study and the work of Park et al [21] emphasize the lowest energy of the Te-rich boundary, supporting the numerical values reported in this work.
A general trend observed in figure 3 corresponds to the increase of the interface energy γ with the boundary atomic number density ρ. One exception nevertheless arises in the case of the Σ3{112}||{112} boundary, where the calculated atomic number density is dramatically lower than for other interfaces with similar energy. It is likely that the crystallographic symmetry of the boundary planes plays a significant role here. The influence of symmetry appears more clearly when correcting for the initial density of the planes involved in the boundary as displayed in table 3. As expected, the P value of the coherent Σ3{111}||{111} boundary is equal to 1, since the interface is formed by simple juxtaposition of sharp {111} atomic planes. In addition, symmetric Σ3 boundaries tend to exhibit the lowest relative densities, implying a higher interface sharpness.
As visible in figure 3, the lowest energy configuration is actually obtained for the Σ3{112}||{255} interface. This result supports the previous observation, namely, that   figure 4 and is found to be composed of three sets of Shockley partial dislocations, with 1/6〈112〉 Burgers vector. These dislocations share a common 〈110〉 line direction in both grains and are of the edge (1 set) and 30°F Interestingly, this dislocation arrangement is the only one that can be resolved by the dislocation analyzer available in the Ovito software [18,22]. In this case, the GB is recognized as a well-ordered structure across which Burgers circuits can be easily mapped. In contrast, for all other boundaries the software returns a 'non-crystalline' area, which is clearly more disordered. The higher order of the former boundary results in a lower energy of this one compared to all others. In figure 3, bars are color coded accounting for the stoichiometry of the boundary. It is apparent that for most of the incoherent twin boundaries studied, lowest energy configurations are either Te-rich (σ lower than 1) or stoichiometric (σ equals to 1). The only exception is the Σ3{112}{127} interface, which is slightly Cd-rich. Since all different possible starting configurations had been tested, it is possible to compare the different features among similar boundaries with reversed chemical composition. The main difference between Cd and Te lies in their electronic configurations. While Cd possesses two valence electrons, Te has six. Therefore, while two doubly-bonded Cd atoms will repel each other due to the Coulomb interaction of their nuclei, the same is mitigated in the case of Te, due to the possible electronic interactions between their remaining valence electrons. This trend seems to be particularly well-captured in the BOP and it explains quite well the energy difference between the Te-or Cd-rich Σ3{112}{122} interface as reported in this work (Δγ=0.2 J m −2 ) and also observed by Park et al [21] (Δγ=0.1 J m −2 ). The structure of the Te-rich Σ3{112} {112} GB is displayed in figure 2(a). Another example of the critical role of stoichiometry on the boundary energy is visible in the case of Σ3{111}{115} GB as displayed in figure 5. Since both boundary planes are polar, four different chemical variants of the same structure exist. It appears once again that Te-rich cores exhibit the lower energies, whereas higher Cd concentration in the core leads to higher σ values. As can be observed in figure 5, the σ value does not fully capture this chemical trend, as it averages the Cd/Te ratio over the full thickness of the boundary. Therefore, it distinguishes between bulk and boundary, rather than between GB core and GB outermost. The larger thermodynamic stability of Te-rich boundaries may be a beneficial factor for CdTe as Zhang et al [23] pointed out that these interfaces are more likely to be efficiently passivated upon chlorine annealing.
When comparing the properties of the atomically sharp GB (plane polarity, ρ′) with their relaxed counterparts (σ, ρ), it appears that no prediction can be easily made a priori, solely based on the features observed on unrelaxed interfaces. This is clearly emphasized by the overall large P values reported in table 3.

Higher-order Σ interfaces
The study of CSL boundaries becomes increasingly complex with increasing values of Σ. Indeed, as mentioned in the previous section, the volume of the CSL unit cell linearly scales with the Σ value. As such, if the boundary planes are kept identical, a Σ9 boundary would require a simulated volume cell three times larger than its Σ3 counterpart, to keep the simulated CSL period number constant. As a result, computation times severely increase and only few high-order Σ boundaries with low-indexed boundary planes can technically be studied.
Second order Σ9 twin boundaries. In materials with a low stacking-fault energy, Σ9 interfaces usually represent a non-negligible fraction of the total GB network length (commonly around one-fifth of the Σ3 fraction [15,24]). This is due to the fact that their density is enforced by successive twinning between adjacent grains according to the CSL multiplication rule [24].
A comparison between the interface properties of first order, Σ3, and second order, Σ9, twin boundaries, both with identical interface planes, is displayed in table 4. As could be expected, the energy of the Σ9 boundaries is systematically higher than that of Σ3 interfaces with the same boundary planes. Nevertheless, the magnitude of this increase strongly varies, depending on the crystallography of the facets considered. Indeed, as displayed in table 4, Δγ can be as low as +2% to as high as +86% in the extreme case of the {112}||{255} GB. In addition, in this case, the sharp increase in Δγ is concomitant to a similar increment of the boundary atomic number density difference Δρ. Nevertheless, for moderate energy changes, as observed for the {111}||{115} and {110}||{110} interfaces, Δρ is actually negative. This means that the extent of the defected area is smaller for the Σ9 interface than for its Σ3 counterpart, but it also means that the average energy per boundary atom still increases, resulting in a positive value of Δγ. Finally, the interface stoichiometry σ does not seem to show any correlation with the boundary misorientation. Therefore, considering the four different boundary sets displayed in table 4, it appears that it is hardly possible to infer, a priori, the structure and energetics of a Σ9 boundary based on the knowledge of a Σ3 interface with the same boundary planes. In addition, when combining all facets, the energy range of Σ9 interfaces, varying between 0.45 and 0.73 J m −2 , overlaps with the one of Σ3 interfaces, ranging from 0.35 to 0.68 J m −2 as visible in figure 3.  As such, when considering the misorientation parameters only, Σ9 boundaries are not clearly differentiable from incoherent twin boundaries from an energetic point of view. This, once again, stresses the need for full-rotational parameters characterization of GB, when assessing their influence on recrystallization behavior or functional properties. A special case is observed for the {111}||{115} interface which exhibits extremely close structural parameters for its Σ3 and Σ9 variants. Both structures are presented in figure 6 and they closely resemble each other. This is a particular feature of boundaries composed of a {111} interface plane. Indeed, Σ3 n type boundaries relate to successive twinning, {111} being the twinning plane. Hence, as visible in figure 6, the second order twin can be obtained by taking the Σ3 variant, fixing the left-hand side of the boundary and performing a pure twist of 180°in the plane of the interface. Therefore, both repeating units present recognizable features and first and second order twins are closely related.
Other Σ boundaries. As emphasized in the previous sections, boundary structures are highly dependent on the interface planes considered. Therefore, in order to extract the influence of the Σ value only, it is necessary to compare GB with identical crystallographic planes. This is Figure 6. Atomic structure of a Σ3{111}{115} boundary (left) and of the corresponding Σ9{111}{115} interface (right). Cadmium atoms are drawn in red whereas tellurium atoms are displayed in blue. Note that bonds are arbitrarily drawn based on the distance between nuclei and not based on actual electronic interactions between atoms. the case described in table 5 for {112}{112} interfaces with different Σ values ranging from 3 to 11. Due to the lack of boundaries studied, it is not possible to draw any general conclusions regarding interface energy trends with the Σ parameter. Some interesting features can, nevertheless, be noted from table 5. First of all, as mentioned in the previous section, the energy of the Σ3 interface is always lower than their higher Σ counterparts. When comparing the boundary energy of the Σ5 and Σ7 boundaries, it becomes apparent that γ does not simply increase with Σ as could be expected. Indeed, the number density of CSL lattice sites intuitively relates to a notion of ordering, which can naturally be assumed to also exist at the boundary. Nevertheless, although a sharp energy difference is observed between Σ3 boundaries and other misorientations, the energy differences between Σ5, Σ7, and Σ11 are distributed in a much narrower range. As in the previous case, the boundary energy is once again found to increase with the Σ value, whereas the stoichiometry of the boundary does not seem to follow any particular trend. Nevertheless, as mentioned in the case of Σ9 boundaries, when the interface plane is disregarded, it appears that no sharp energetic gap is found to separate Σ3 interfaces from higher-order boundaries. As such, the predominance of the former may primarily relate to their ease of nucleation. In addition, GB mobility may also be invoked as a possible explanation for the low occurrence of these low-energy configurations in the annealed thin films. Indeed, highly coherent sigma GB have practically zero mobility. Therefore, they cannot be easily removed during any competitive and capillary coarsening, leading to their prevalence after heat treatments.
A surprisingly low GB energy is measured in the case of Σ7 {111}||{111} interface (not shown here), with a GB energy of only 0.28 J m −2 . This value is even lower than the energy of most stable incoherent Σ3 twin boundaries. This can nevertheless be mitigated by stressing that the energy of symmetric {111} interfaces tends to be underestimated by the BOP, as already previously mentioned in the case of the coherent first order twin.

Conclusions
The structure and energy of CSL boundaries was investigated by mean of MS. As expected, the coherent Σ3 {111}||{111} interface is found to be significantly more stable than any other GB. Nevertheless, when considering other interface planes, a broad spectrum of energies and boundary stoichiometries is observed. Overall, Te-rich boundaries are found to be more energetically stable than their Cd-rich counterparts. This is an important point since Te-rich boundaries are predicted to be less detrimental for the overall cell efficiency upon processing. In addition, symmetric GB tend to exhibit a lower relative atomic number density, indicating a lower degree of relaxation. This, however, does not necessary lead to lower energy configurations, and high-indexed planes can lead to surprisingly stable GB, as it is the case for the Σ3 {112}||{255} interface. In our simulations this is achieved through the formation of a particular dislocation arrangement in the core of the boundary. For a given set of interface planes, the lowest energy configuration does not simply scale with the Σ value, even though the lowest configuration seems to be always obtained for the Σ3 case. When considering the misorientation solely, i.e. Σ only, a significant overlap in terms of interface energy is found to exist between first order twin boundaries and higher Σ ones. Therefore, we found Σ5, 7 or 9 boundaries with lower interfacial energy than incoherent Σ3 ones formed on other interface planes.
These results stress the shortcoming of GB structure/property correlations, when confining the analysis only to the misorientation space. In addition, the overall predominance of Te-rich CSL, and more detailed knowledge about the magnitude of the stoichiometric shift as a function of the interface crystallography may be useful when interpreting the electronic response of such defects or when designing passivation strategies.