1. Introduction

When considering the properties of crystalline materials, the impact of defects is essential. Point defects such as vacancies and dopants are the defects most commonly considered in both computational and experimental studies of material properties. Furthermore, the modelling of materials at an atomic level is often confined to bulk systems which contain these point defects (1–4). Considerably less is known about extended defects which appear in polycrystalline systems such as surfaces, dislocations and GBs. As nanostructuring of materials is becoming more prevalent, the behaviour of these extended defects is becoming significantly more important (5–8). GBs give rise to structural discontinuities within materials which result in specific structures and potential non-stoichiometry and can lead to the segregation of point defects to varying degrees, depending on the specific structure (9–13). This can significantly affect the macroscopic properties, for example: ionic conductivity, electronic conductivity, thermal conductivity, thermal expansion, elasticity and strength – all of which are crucial for many applications. Therefore, the understanding of interfaces in these materials is key to optimising their performance. Despite this, relatively little is known about the structure and even less of the effects of these interfaces on material properties due to the inherent complexity of the issues.

An example of the importance of interfaces and polycrystallinity is in the fluorite structured fast oxide ion conductors (14), such as yttria stabilised zirconia (YSZ or ZrO₂-Y₂O₃) or trivalently doped ceria (CeO₂), used in solid oxide cells, oxygen membranes and oxygen sensors (6, 15, 16). It is reported that the ionic conductivity within the GBs of these materials is several orders of magnitude lower than the bulk (17–20) with the effect attributed to a wide range of causes including impurities, dopant segregation, defect cluster formations and space charge layers (7, 8, 20–25). In contradiction it has been observed that other materials, such as Bi₂O₃ (26) and nanostructured YSZ, that the ionic conductivity is enhanced (27, 28). Much of the experimental data is based on average effects observed in impedence spectroscopy, where all GBs are treated equally as an average effect (29, 30). In fact, GBs can take on specific structures, an example of this is shown in Figure 1. That is, the atomistic description of what is happening at the GB is incomplete when obtained from macroscopic observations. Another issue which may arise when only considering average effects is that it is likely that different specifically defined interfaces will behave in different ways. As it is difficult to isolate and study the effects of GBs experimentally, computational studies are invaluable to further our knowledge of these defects and their impact on material properties.

The failure to understand the basis of material properties in polycrystalline samples is a significant impediment to the development of new materials and the application of inexpensive processing methods to existing materials. An enhanced understanding of the impact of GBs and polycrystallinity on the properties of materials would allow us to explore alternative routes to optimise their properties and ultimately enhance devices. In order to model the properties of these interfaces we first require a method for accurate prediction of interfacial structures. In this paper we present a computational method for accurately predicting the structure of low angle mirror tilt GBs which can be applied to other interfaces and even heterointerfaces. This method utilises both atomistic simulation and classical molecular dynamic simulation with sophisticated, polarisable force fields derived from ab initio data. Previous theoretical studies of GB structures generally utilise static lattice simulations with empirical force fields and structures based on experimental structures (31, 32). These results often have to be validated using first principles due to the quality and limitations of the force field. In this work the structures are predicted and validated using high-quality force fields derived from ab initio data, this is discussed further in the methodology.

Two fluorite-structured materials are investigated in this study: calcium fluoride (CaF₂) and CeO₂. CaF₂ is the prototypical fluorite material which is a super-ionic conductor at high temperatures (>1100 K) (33). CeO₂ (usually doped) is a highly technologically significant material which is both an ionic and electronic conductor with a wide range of applications including catalysis, solid oxide fuel and electrolysis cells and oxygen sensing (6, 15, 16). We compare their predicted GB structures to experimental structures from the literature obtained via transmission electron microscopy (TEM).

2. Grain Boundary Structures: Generation and Definition

All GBs simulated here were generated using the minimum energy techniques applied to dislocation, interface and surface energies code (METADISE) (34). The most stable GB structures were found by carrying out optimisation scans of the GB. Surfaces with specific Miller indices were first cut and then reflected to form an interface. A potential energy surface (PES) was then calculated using a forcefield by scanning one surface relative to the other. From this scan, a two-dimensional (2D) PES for the boundary was calculated which allowed the minimum energy structure to be identified. The minimum energy GBs were then optimised and the most stable boundary was selected to investigate using molecular dynamics. The 2D potential energy scan along with the GB structure (before and after optimisation) for the Σ9(221) GB in CeO₂ is shown in Figure 2.

GBs are defined by a number of parameters: the crystallographic directions of the axes of the two grains which come together to form the interface (h_i, k_i, l_i ), the rotation axis o = (h_o, k_o, l_o ), the misorientation angle θ around the axis o and the normal axis to the GB plane n. When n is parallel to o the boundary is defined as a twist GB and when n is perpendicular to o the boundary is defined as a tilt boundary. The GBs which are studied in this work are high-angle mirror tilt GBs (n ⊥ o) and the rotation axis is (001).

The geometric definition of the GBs used in this work is the coincidence site lattice model (35). A coincidence lattice site can be defined when there exists a finite fraction of coinciding lattice sites between the two lattices (grains). This model is based on the assumption that when the energy of the GB is low, the coincidence of the atomic sites between the two grains is high, i.e. there are few bonds which are broken across the boundary. The reciprocal density of coincidence lattice sites is known as Σ and is used to characterise the geometry of the GB, as given in Equation (i):

(i)

For cubic lattices, the Σ value can be given by the sum of the squares of the Miller indices of the symmetrical tilt boundary, given by Equation (ii):

(ii)

where δ = 1 if is odd and δ = 0.5 if is even, thus in cubic systems Σ is always an odd number (35, 36). For example, the Σ9(221) GB shown above is defined by the (221) Miller index of the surfaces which are scanned to give this boundary, i.e. (2²+2²+1²) = 9, which is odd so δ = 1 and thus this is written as Σ9(221). The other GBs studied here are defined in the same way.

3. Methodology

Initial GB structures were generated as outlined above using METADISE with shell model interaction potentials for both CaF₂ (37, 38) and CeO₂ (39). These structures were then expanded to at least 30 Ångström (Å) in the x-direction, 22 Å in the y-direction (parallel to the GB) and 76 Å in the z-direction (perpendicular to the GB). Each simulation cell contained two identical GBs as illustrated in Figure 3, with each grain having a depth of at least ~35 Å.

Molecular dynamics simulations were then carried out to determine the average GB structures. The interaction potential used for the molecular dynamics simulations is known as the dipole polarisable ion model (DIPPIM) (40), implemented in the polarisable ion model aspherical ion model (PIMAIM) code (41). The DIPPIM consists of four elements: charge-charge interactions, short-range repulsion, dispersion interactions and polarisation. This is a highly accurate, polarisable, potential, in which the dipoles are solved self consistently at each molecular dynamics step. This leads to a highly accurate description of the dipoles on ions in the simulation which is of particular importance when simulating highly polarisable ions such as F^– and O^2–. The data used to fit the DIPPIM potentials used in this work were calculated using ab initio methods (2, 42, 43). The use of ab initio data allows for non-equilibrium details on the PES to be accounted for which leads to a highly accurate, transferable interatomic potential.

Often interatomic potentials for fluorite materials are derived from equilibrium experimental data or are formed using interatomic potentials from a range of different sources resulting in inconsistent, non-transferable potentials which may have difficulties taking effects of different coordination environments into account, i.e. surfaces and interfaces. Such interatomic potentials are usually better suited to static lattice simulations as opposed to molecular dynamics simulations. In previous work on the surfaces of CeO₂ (44) we have shown that the DIPPIM provides an accurate description of extended defects and the effect of such defects on ionic transport.

The simulation cells were heated to 1473 K for 500 ps in order to simulate annealing of the GB structures, they were then cooled to 573 K for 500 ps and finally simulated at 300 K for 500 ps. Temperature scaling was carried out (at all three temperatures) every 0.025 ps before data collection for analysis began. Final GB structures were generated by averaging over the frames of the trajectory at 300 K. The DIPPIM potential parameters used for CaF₂ were previously derived by Pyper and Wilson et al. (42, 43) and those for CeO₂ were obtained by Burbano et al. (2). All steps of the simulation were carried out using the isothermal-isobaric ensemble (NPT). CaF₂ simulations utilised a timestep of 5 fs and a short-range cut-off of 14 Å and in the case of CeO₂ simulations had a timestep of 4 fs and a short-range cut-off of 11 Å. The GBs which were simulated for CaF₂ were Σ3(111), Σ5(210), Σ5(310), Σ9(221), Σ11(332), Σ13(320) and Σ13(510); and for CeO₂ are Σ3(111), Σ5(210) and Σ9(221).

4. Results and Discussion

Here we present the average predicted structures obtained for GBs in fluorite structured materials and compare these with TEM images obtained from experimental studies. As the F^– and O^2– ions present in CaF₂ and CeO₂ are difficult to image due to their low atomic masses we only compare the cation structures obtained with the aforementioned TEM images. First, we discuss the CaF₂ structures followed by those found for CeO₂. To the authors’ knowledge there are no experimental studies of GB structures in CaF₂ so those found in this study are compared to those of other fluorite materials (CeO₂, ZrO₂ and YSZ).

5. Conclusions

We have presented a computational method for the prediction of the structure of mirror tilt GBs in fluorite structured materials. This method utilises interatomic potentials which are derived from first-principles data meaning the process is entirely predictive. The excellent level of agreement with existing experimental data on the structures of fluorite GBs highlights the power of the method. The ability to accurately predict these structures is an important first step into the computational investigation of the properties of these materials, which is key to future materials and device optimisation. The method presented here can be extended to the prediction of interfaces in different materials, interfaces of different types (i.e. twist GBs) and even heterointerfaces.