What Happens at Surfaces and Grain Boundaries of Halide Perovskites: Insights from Reactive Molecular Dynamics Simulations of CsPbI3

The commercialization of perovskite solar cells is hindered by the poor long-term stability of the metal halide perovskite (MHP) light-absorbing layer. Solution processing, the common fabrication method for MHPs, produces polycrystalline films with a wide variety of defects, such as point defects, surfaces, and grain boundaries. Although the optoelectronic effects of such defects have been widely studied, the evaluation of their impact on the long-term stability remains challenging. In particular, an understanding of the dynamics of degradation reactions at the atomistic scale is lacking. In this work, using reactive force field (ReaxFF) molecular dynamics simulations, we investigate the effects of defects, in the forms of surfaces, surface defects, and grain boundaries, on the stability of the inorganic halide perovskite CsPbI3. Our simulations establish a stability trend for a variety of surfaces, which correlates well with the occurrence of these surfaces in experiments. We find that a perovskite surface degrades by progressively changing the local geometry of PbIx octahedra from corner- to edge- to face-sharing. Importantly, we find that Pb dangling bonds and the lack of steric hindrance of I species are two crucial factors that induce degradation reactions. Finally, we show that the stability of these surfaces can be modulated by adjusting their atomistic details, by either creating additional point defects or merging them to form grain boundaries. While in general additional defects, particularly when clustered, have a negative impact on the material stability, some grain boundaries have a stabilizing effect, primarily because of the additional steric hindrance.

1 Structural models for surfaces 1

.1 Surface termination names
The different surface terminations are named according to the relative concentration of species at these surfaces. To determine the surface concentration of the species, we create perovskite slabs with the same termination at both ends and consequently calculate the ratios of the atoms that make up these slabs. An overview of the ratios found for the different slabs is shown in Table S1. The naming scheme is based on the relative ratio of the I and Pb species at each surface, resulting in the following names: stoichiometric (I:Pb = 3), Pb-poor (I:Pb > 3) and Pb-rich (I:Pb < 3).

Surface models
In this work the surface models are created in AMS2021. 1 An orthorhombic unit cell is used as the starting point for the surface models. Each model slab is at least 50Å thick, with the lateral dimension spanning at least 45Å. An overview of the full structural models of the doubly terminated (Pb-poor and Pb-rich) slabs is shown in Figure S1.

Grain boundary models
The grain boundary models in this work are created with the aimsgb package. 3 A cubic CsPbI 3 unit cell is used as the starting point for the grain boundary models. An overview of the initial structures of the grain boundary models is shown in Figure S2. In each model the grain boundaries are separated by at least 25Å. The lateral dimensions of the grain boundary models are at least 30Å in size.

Structural effects of surfaces
The time-averaged structures of the 10 layer thick slabs of CsPbI 3 are shown in Figure 2 in the manuscript. For a complete overview of the slabs with varying thickness, the timeaveraged structural models of the slabs of 4, 6 and 8 octahedral cages thick are shown in Figure S3, which demonstrates that thinner slabs attain more cubic-like structures.

S7
Slab geometries were created from optimized bulk structures of the cubic phase of CsPbI 3 in both DFT and ReaxFF. During relaxation of the slabs, only the ionic positions were allowed to change until convergence was obtained in energy and forces. In the DFT calculations a 6 × 6 × 1 k-points mesh was used during slab optimizations. At least 15Å of vacuum was included to minimize the interaction between periodic images of the slab. The structural models of the perovskite slabs are shown in Figure S4. Surface energies (γ) were calculated as where E slab is the energy of the perovskite slab, n the number of formal units of CsPbI 3 , E CsPbI 3 is the energy of cubic CsPbI 3 , E precursor is the energy of the perovskite precursor and A the surface area of the perovskite slab. To make use of a common point of reference for the validation, we pick CsI as the precursor, resulting in n = 6 and an addition of the precursor energy (+E precursor ) for the Pb-poor terminated slab and n = 7 and a subtraction of the precursor energy (−E precursor ) for the Pb-rich terminated slab.

S8
An overview of the surface areas and calculated surface energies of the perovskite slabs with the Pb-rich and Pb-poor termination is shown in Table S2. From this comparison we find that the used ReaxFF force field matches the stability trend of our DFT simulations, as well as those found in literature. 4,5 In both computational methods the Pb-poor surface termination is found to be significantly more stable than the Pb-rich surface termination. Although the relative magnitude for DFT and ReaxFF does not exactly match, the comparison validates that the force field is transferable for the simulation of surfaces.  Figure S5. We underline that the jump of an iodine atom to a neighboring Pb atom ( Figure S5a) forms local Pb-rich and Pb-poor octahedra, which become the degradation centers at elevated temperatures.

a) b)
t' = 0 ps t' = 10 ps Figure S5: Snapshots of the rearrangement of surface iodine atoms for the orthorhombic (020) stoichiometric surface at 500 K. (a) Stuctural arrangement before the migration of surface iodine. (b) Surface geometry after surface iodine rearrangment, with the blue and red circle indicating a Pb-poor and Pb-rich octahedra, respectively. A shifted time axis t ′ is used.

S10
The RDF between the atoms of the species a and b is calculated using which sums over Dirac delta functions δ that are offset by the interatomic distance between the ith atom of species a located at ⃗ r i and the jth atom of species b located at ⃗ r j . The quantity is averaged over the full length of the simulation indicated by the ensemble averaging (⟨· · · ⟩).
In this work, we employ the radial distribution function (RDF) as implemented in the

S12
Defective perovskite surfaces are created with vacancies. The defects in these structural models are created by removing the atoms of certain species, thus resulting in vacancies in the lattice. For the isolated defects, we create slabs with a two vacancies in total, one of each kind. As shown in Figure S7a, the top surface has a V Cs defect and the bottom surface a V I defect. For closely spaced pairs of V Cs and V I defects, Figure S7b