Halogen Vacancy Migration at Surfaces of CsPbBr$_3$ Perovskites: Insights from Density Functional Theory

Migration of halogen vacancies is one of the primary sources of phase segregation and material degradation in lead-halide perovskites. Here we use first principles density functional theory to compare migration energy barriers and paths of bromine vacancies in the bulk and at a (001) surface of cubic CsPbBr$_3$. Our calculations indicate that surfaces might facilitate bromine vacancy migration in these perovskites, due to their soft structure that allows for bond lengths variations larger than in the bulk. We calculate the migration energy for axial-to-axial bromine vacancy migration at the surface to be only half of the value in the bulk. Furthermore, we study the effect of modifying the surface with four different alkali halide monolayers, finding an increase of the migration barrier to almost the bulk value for the NaCl-passivated system. Migration barriers are found to be correlated to the lattice mismatch between the CsPbBr$_3$ surface and the alkali halide monolayer. Our calculations suggest that surfaces might play a significant role in mediating vacancy migration in halide perovskites, a result with relevance for perovskite nanocrystals with large surface-to-volume ratios. Moreover, we propose viable ways for suppressing this undesirable process through passivation with alkali halide salts.


I. INTRODUCTION
Halide perovskites are exciting materials with exceptional optoelectronic properties, wide tunability, and a broad range of applications spanning solar cells 1-3 , light-emitting diodes (LEDs) 4,5 , photo-detectors [6][7][8] and X-ray scintillators 9 . Solar cells based on lead-halide perovskites ABX 3 with A=(CH 3 NH 3 ) + (methylammonium, MA + ), (NH 2 CHNH 2 ) + (formamidinium, FA + ), Cs + , B=Pb 2+ , and X=Cl − , Br − , I − , can be processed at low-temperature, and have exceeded power conversion efficiencies of 25% 10 . Yet, commercialisation of perovskite-based solar cells and other devices is hampered by the lack of stability of the perovskite absorbers towards moisture, oxygen, light, heat, and electric fields 11 . Various strategies have been applied to improve the stability of these materials, including encapsulation, (partial) replacement of the A site cation 12 , and passivation 13 . All-inorganic leadhalide perovskites CsPbX 3 have seen their own surge of interest, in particular because colloidal CsPbX 3 nanocrystals can exhibit very high photoluminescence quantum yields, with band gap energies and emission spectra tunable over the entire visible spectral region 14 . However, even all-inorganic halide perovskites can exhibit poor stability under electric fields 15 . Material degradation and phase separation in both organicinorganic and all-inorganic halide perovskites have been attributed to the migration of mobile ionic species 16 .
Ion migration in halide perovskites has been studied since the 1980s 17 . The dominant migrating species in these materials are halogen ions [18][19][20][21] , mediated by the presence of halogen defects. The mechanism of halogen migration has been studied both experimentally 17,[22][23][24][25][26] and using first principles simulation techniques such as density functional theory (DFT) 19,23,[27][28][29][30] . And while reported activation energies for these migration processes span a wide range from ∼0.1 to ∼1.0 eV [18][19][20][21]23,31,32 , there is a consensus that halogen migration is the primary channel for the ionic conductivity observed in halide perovskites. The large spread of the experimental values has been linked to synthesis conditions, experimental techniques, and the role of grain sizes for defect formation in polycrystalline thin films 33,34 .
Prior first principles calculations of ion migration in both organic-inorganic and all-inorganic halide perovskites have, with a few exceptions, focused on ion migration in the bulk. However, in perovskite nanocrystals, ion migration at surfaces is expected to play an increasingly large role with decreasing particle size. Furthermore, in halide perovskite thin films, Kelvin probe force microscopy was used to demonstrate that ion migration is dominant at grain boundaries 35,36 . More recently, the effect of surfaces and grain boundaries has also been explored via first principles calculations 30,37 . However, these computational studies have provided a mixed picture: Meggiolaro et al. reported that migration energy barriers of interstitial iodine vacancies are little affected by surfaces 37 . On the other hand, Oranskaia et al. showed a clear effect of the surface on Br vacancy and interstitial migration in MAPbBr 3 and FAPbBr 3 30 . For both materials, the activation energies of Br vacancy migration were computed to be 0.3 eV lower than in the bulk. The variation in calculated results can have a number of sources such as differences in the applied level of theory, for example the choice of the DFT exchange-correlation functional and the method used for calculating migration barriers. A complication that is particular to the organic-inorganic perovskites, is that in the majority of DFT calculations, the rotational dynamics of the organic cation at room and higher temperatures are not taken into account. Instead structural models with fixed orientations of the molecular moieties are used, leading to significant differences in the potential energy landscape depending on the choice of molecular orientation 38 . Indeed, Oranskaia et al. also showed that activation energies for Br migration significantly depend on the orientation of the organic moiety. The uncertainties associated with the choice of a suitable structural model for organic-inorganic perovskites at elevated temperatures, and the important role of surfaces in all-inorganic halide perovskite nanocrystals, are motivating our first principles study of bromine vacancy migration in the bulk and at a surface of cubic CsPbBr 3 .
Our DFT calculations of vacancy-mediated bromine migration paths and energy barriers in CsPbBr 3 , show a significant dependence on the presence of a surface. We find that the migration barrier in the bulk is about twice as large as the one at the surface. We show that variations of the Br migration barriers are correlated with variations in the Pb-Br bond length of bonds in the vicinity of the vacancy: halide migration at the surface is facilitated by the larger structural flexibility of the surface as compared to the bulk. Furthermore, migration paths considerably differ between the surface and the bulk, which can also be traced back to more flexible bonds at the surface. Finally, we study the effect of surface modification with alkali halide monolayers, demonstrating that a NaCl passivation layer leads to an increase of the migration energy of Br vacancy migration to almost the value in the bulk of the material.

II. METHODS
CsPbBr 3 is orthorhombic with Pbnm symmetry at room temperature and undergoes two successive phase transitions to tetragonal (P4/mbm) at 88°C and to cubic (Pm3m) at 130°C 39,40 . Figure 1 depicts the bulk and surface slab structures of cubic CsPbBr 3 used in this work. For constructing these structural models, we first performed a geometry optimization starting from the experimental high temperature crystal structure of CsPbBr 3 with Pm3m symmetry using DFT within the PBEsol approximation 41 as implemented in the Vienna Ab−initio Software Package (VASP) 42,43 . The resulting optimized lattice parameter of 5.86 Å is in very good agreement with the experimental value from X-ray diffraction 44 . To model the surface of CsPbBr 3 , we designed two (001) surface slab supercells with distinct surface terminations by repeating the primitive Pm3m unit cell with PBEsol-optimized lattice parameters twice along the [100] direction, once along the [010] direction and six times along the out-of-plane [001] direction, with the bottom three layers fixed to bulk positions and the top three layers fully mobile. Unless otherwise specified, all our calculations were performed with this unit cell setup. Surface A is PbBr 2 -terminated and surface B is CsBrterminated. The surface energy is converged to within 25 meV with respect to the slab thickness. To avoid spurious interactions between periodic images, we inserted 30 Å of vacuum along the (001) direction. The bulk system features the same number of layers without vacuum. We used the projector augmented wave (PAW) method 45 , a cutoff energy for the planewave expansion of 300 eV and a k-point grid with 4 × 4 × 4 points for the bulk and 4 × 4 × 1 points for the slab system. For geometry optimizations, we used a convergence criterion of 0.05 eV/Å. In all structural optimizations, the volume and shape of the unit cells was kept fixed.
We computed migration paths and migration energies using the climbing-image Nudged Elastic Band (cNEB) approach 46 , an optimization method for identifying the minimum energy path between a given initial and final state. We used three cNEB images, i.e., intermediate structure snapshots between the initial and the final state of the system, to simulate the migration of the Br vacancy. In the cNEB method, the images along the reaction path are optimized such that the highest energy image is driven up to the saddle point. The migration barrier represents the amount of energy necessary for an ion or a defect to move from the initial to the final state and is calculated as the difference between the energy of the saddle point and the energy of the initial state of the transition.

A. Structural changes upon vacancy formation
We start by performing geometry optimizations of the bulk and A-and B-terminated surfaces shown in Figure 1. As expected, the bulk system remains unaffected by further geometry optimization, while the slab structure is compressed at the surface, with axial Pb-Br bonds by more than 4 % shorter than in the bulk and almost unchanged equatorial bonds. The average relative variation of Pb-Br bonds per layer as compared to the bulk Pb-Br bond length of 2.93 Å is shown in Figure 2(a), where we have averaged over all axial and equatorial Pb-Br bonds (per layer), respectively. A DFT study on intrinsic point defects in CsPbBr 3 showed that under Br-poor growth conditions, bromine vacancies have the lowest formation energy among all possible point defects 47 . In our surface slab structures, a Br vacancy can occupy three symmetry-inequivalent positions in each layer of surface A and B: axial, equatorial along the [100] direction, and equatorial along the [010] direction. Note, that these three vacancy positions also have different energies in our structural model for the bulk, which is an artifact of our asymmetric unit cell. We define the formation energy of a vacancy in the slab (E slab f ) and in the bulk (E bulk f ) as the difference between the energies of the unit cell with and without vacancy. In Table I, we report the binding energy E B = E bulk f − E slab f to quantify by how much a vacancy prefers to bind to the surface as compared to the bulk. E B is largest in L 1 and converges to zero in subsequent layers, in agreement with results for MAPbI 3 by Meggiolaro et al. 37 , suggesting that the surface is more prone to defects. We further find that binding to surface A is preferred over binding to surface B, in line with observations of iodine vacancy clus-  tering at MAI-terminated surfaces in MAPbI 3 48 . For completeness, we also report E B for the two equatorial vacancy positions and note that this value is significantly larger for the vacancy along [010] because of our 2 × 1 × 6 unit cell setup. In the following, we will only discuss ion migration between axial vacancies. In Figure 2(b) and (c) we show the average Pb-Br bond length variation with respect to the undistorted bulk bond length upon introduction of the axial Br vacancy and geometry optimization. The creation of a Br vacancy at surface A leads to severe distortions of the system, featuring axial Pb-Br bonds reduced by up to 20 %. In contrast, introducing a vacancy at surface B leads to smaller variations and a less distorted structure. We have also calculated the variation of the Pb-Br bond length for a Br vacancy generated in the deeper lying surface layers, and find that even though the absolute value of the bond length variation differs for the two surface terminations, in both cases the variation in axial bond length is ∼5 times larger than the variation in equatorial bonds. Furthermore, the deeper the vacancy is created, the less compressed the structure is at the surface. In Figure 2(d) we show that formation of a Br vacancy within the bulk has similar con-sequences, i.e., a bond length compression in the vicinity of the vacancy. However, the distortions in the bulk are highly suppressed due to a more rigid structure, with fewer degrees of freedom in comparison with the surface slabs, leading to a zero average variation of the bond lengths when averaging over all mobile layers. Note that the large compression of more than 20% in L 1 of surface A is an artifact of the asymmetric unit cell. In a 2 × 2 × 6 cell, the compression of axial bonds is smaller than in the 2 × 1 × 6 cell, with a relative bond length compression of 4.6% in L 1 of surface A without a vacancy, 3.4% in L 1 of surface A with a vacancy and hardly any variation with respect to the undistorted bulk for the case of a vacancy in the bulk. However, the trends for subsequent layers are similar to the 2 × 1 × 6 unit cell.

B. Br vacancy migration in CsPbBr 3
Next, we use the cNEB method to determine the energy barrier of Br vacancy migration between two adjacent axial vacancy positions at both surfaces and for vacancy migration in layers L 2 and L 3 of the A-terminated surface. The migration barrier is calculated as the difference between the total energies of the initial state and the saddle point. Table II summarizes our results. For the bulk structure, we compute a migration barrier of 0.65 eV. Our specific unit cell setup is by design not suitable for direct comparison with experimental results or bulk calculations of halide migration in symmetric structural models, such as the one used in Ref. 49 . However, our setup allows us to realize the same defect concentration and in-plane boundary conditions in the bulk and slab unit cells, and hence compare trends. It is worth mentioning that our calculations are in good agreement with the experimental values ranging between 0.72 and 0.66 eV reported in the literature 17,22 . However, our result is 140 meV larger than the calculated migration barrier reported by Zhang et al. 49 . This discrepancy may be explained based on different structures, defect concentra-tions, and approximations for the exchange correlation energy (Ref. 49 uses the orthorhombic phase of CsPbBr 3 and the PBE approximation). Furthermore, we expect that our calculations represent an upper bound on the migration barriers, since we are neglecting the large, anharmonic vibrations reported for CsPbBr 3 and other halide perovskites at room and higher temperatures 50 . Table II and Figure 3(a) and (b) show the activation energy for Br migration at the A-and B-terminated surfaces, and within subsurface layers of surface A. Our first main finding is that the migration energy at both surfaces is substantially lower than that in the bulk, more than a factor of two at the B surface. We have confirmed that our finding of a significantly lower migration energy at the surface also holds in a 2 × 2 × 6 unit cell setup, where we calculate migration energies of 0.48 eV in the bulk and 0.28 eV at the surface with values of 0.20 eV and 0.26 eV in surface layers L 2 and L 3 ). Furthermore we find that migration barriers at surface A and B are very similar; the migration barrier at surface A is only 90 meV larger than that at surface B. Interestingly, the migration barrier in subsurface layer L 2 is ∼110 meV lower than directly at the surface, and only slowly increases in subsequent subsurface layers. We show in Figure 3(a) that the variation of the migration barrier with the layer number is correlated with relative bond lengths variations of axial and equatorial bond lengths with respect to the bulk. Smaller migration energies are associated with a significant axial compression of the surface. The slightly higher migration barrier at the surface of A as compared to subsequent surface layers is correlated with a subtle interplay between longer equatorial and shorter axial bond lengths as compared to the bulk. However, note that trends in migration barriers as a function of surface depth should be viewed with caution: In our calculations, the bottom layers of the surface slab are fixed to the bulk atomic positions, a constraint that might affect the magnitude of the calculated migration barriers in the layers adjacent to the fixed layers. We therefore calculated migration barriers for a surface slab with four mobile layers as well. The migration barriers for this system, also shown in Figure 3(a) and Table II, are slightly lower but follow the same trends. Observation of bulk-like migration barriers deeper into the surface would likely require structural models with more surface layers. Our results are also in line with observations by Oranskaia et al. showing that larger lattice distortions lead to smaller migration energies for the through-cell migration of a Br vacancy in organic-inorganic perovskites 30 .
A detailed analysis of the migration paths of the Br ion, associated with the energy profiles shown in Figure 3(b), reveals significant qualitative differences between the migration paths in the bulk and at the two surfaces, see Figure 3(c). In the bulk, the halide ion moves along an almost straight line from one axial vacancy position to the other -the shortest possible path. Contrary to that, the migration path is curved at both surfaces, with the saddle point deviating from the straight line. Analogous curved paths for vacancy migration between an equatorial and an axial position have previously been reported for oxide perovskites based on neutron diffraction 51  paths have been reported for Pb-based halide perovskites as well 23,49 . We quantify the curvature of the migration path, δ , by computing the perpendicular distance of the Br ion in the cNEB saddle point configuration to the straight line between the initial and final positions of the Br ion, schematically represented in Figure 3(c). As reported in Table II, we find that in the bulk, the Br ion follows an almost straight line, with a deviation more than 7 times lower than at the A and B surfaces. This finding highlights the more flexible nature of the surface, which can deform and accommodate a defect more easily, explaining the lower energy of Br vacancy migration as compared to the bulk. Finally, we observe that at both surfaces, the saddle point is bowed away from the surface, with δ at the A surface almost double of what it is at the B surface. We find that δ is correlated with the compression of axial Pb-Br bond lengths and can be traced back to the structural symmetry of the two surfaces: Formation of a Br vacancy is associated with the breaking of one bond at the B surface, leading to less restructuring as compared to surface A, where two bonds are broken and both the Pb-Br layer above and below the vacancy adjust to vacancy formation and migration. Motivated by the correlation between migration barriers and surface restructuring, we investigate the effect of surface modification on Br vacancy migration energies. Surface modification is a common strategy for passivating surface and interfacial defect states in halide perovskites 55 . Chemical surface treatment with organic ligands has been shown to increase photoluminescence lifetimes and quantum yields 56,57 . However, organic ligands may lead to problems with stability. Therefore, alkali halides have recently been suggested as interface modifiers between the halide perovskite absorber and the electron-or hole-transport layers in solar cells, with some studies showing that they lead to enhanced stability and device performance 58,59 . Moreover, a first principles study by Apergi et al. demonstrated that alkali halide surface modifiers allow for improved electronic level alignment between the halide perovskite absorber and NiO hole transport layer with wide tunability to match those of various perovskite compositions 60 .
Here, we investigate four alkali halide monolayers, NaBr, NaCl, KBr and KCl, and their effect on Br vacancy migration at the surface of CsPbBr 3 . We construct our passivated systems by placing the monolayer on top of surface A and relaxing the structures. In Figure 4(a) we show the particular case of a slab structure passivated with a NaCl monolayer. Upon geometry optimization we find that the axial Pb-Br bonds of the surface slab are significantly less compressed than those of the unpassivated structure for NaCl-and NaBr-passivated surfaces. In Figure 4(b), we show that the variation of Pb-Br axial bonds at the A-surface is less than 2 % and that of equatorial bonds is negligible. In comparison with the undistorted CsPbBr 3 bulk structure, K-based monolayers feature longer bonds, leading to larger distortions induced by the lattice mismatch between the perovskite and the passivation layer. In fact, passivating the surface with a KBr monolayer does not reduce distortions and yields very similar bonds as compared to the unpassivated system. In contrast, NaCl and NaBr monolayers have bond lengths that differ by only 0.08 Å from the bulk. We therefore hypothesize that passivation with NaCl and NaBr should lead to an increase of the Br vacancy migration barrier at the surface.   Following the same approach as before, we introduce a Br vacancy in the first A-surface layer of the NaCl-passivated and NaBr-passivated surface slab structure, respectively, and calculate the energy of Br vacancy migration in the surface layer. For the NaCl-passivated system, we find a migration barrier of 0.57 eV, only 80 meV lower than that computed for migration in the bulk. Interestingly, however, the vacancy follows a curved migration path, with a larger deviation from the straight line (δ = 1.68 Å) than in the unpassivated system. Furthermore, we find that passivating the surface with NaBr leads to a migration barrier of 0.48 eV, slightly larger than that computed for the migration at the surface of the unpassivated slab system. This result, reinforces our hypothesis that larger variations of Pb-Br bond lengths lead to smaller migration energies and indicates the potential of simple alkali halide salts for suppressing halogen vacancy migration at surfaces of halide perovskites.

IV. CONCLUSIONS
In conclusion, we performed a first principles DFT study of Br vacancy migration in CsPbBr 3 and showed that the migration barrier within the close-packed bulk structure of cubic CsPbBr 3 is roughly twice as large as that at either of the distinctly terminated (001) surfaces of the system. Our calculations suggest that the significant reduction of the migration barrier at the surface is due to the "softer" structure of the surface which allows for significant bond lengths variations as compared to the bulk. Motivated by this observation, we studied the effect of surface modification with alkali halide monolayers and demonstrated that passivation with NaCl significantly decreases the structural distortions seen in the unpassivated surface, in particular the compression of the axial Pb-Br bonds. Consequently, NaCl passivation leads to an increase of the Br vacancy migration barrier at the surface back to almost the value it has in the bulk. Our results highlight the important role of surfaces in determining perovskite stability by facilitating ion migration. The dependence of vacancy activation barriers on Pb-Br bond lengths, in particular the importance of axial bond length compression, suggests that strain engineering, for example via epitaxial growth, could be another viable route for suppressing ion migration in halide perovskites 61 . We believe that future computational studies should be directed towards elucidating the role of grain boundaries, in particular in polycrystalline MAPbI 3 . With the advent of machine-learning force fields with DFT accuracy, reliable structural models of large supercells of organic-inorganic halide perovskites and the inclusion of temperature effects in large-scale molecular dynamics simulations have become computationally feasible 62 .