\section*{Manuscript}

\section{Introduction}

Metal halide perovskites (MHPs) with the \ch{ABX3} formula (\ch{A} = organic or inorganic cation; \ch{B} = metal cation; \ch{X} = halide anion) have attracted a great deal of attention as low-cost and high-performance semiconductors for applications in photovoltaics~\cite{kePerovskiteSolarCell2014, kimHighEfficiencyPerovskiteSolar2020} and light-emitting diodes~\cite{vanleRecentAdvancesHighEfficiency2018, shenEfficientPureBlue2022}. MHPs are typically synthesized using a facile solution processing deposition techniques. Halide salts are mixed in solution with metal halide precursors, where solvent evaporation causes the formation of colloidal particles in solution, with further evaporation resulting in the growth of these particles in a polycrystalline perovskite film~\cite{yanHybridHalidePerovskite2015, nayakMechanismRapidGrowth2016, mcmeekinCrystallizationKineticsMorphology2017, flatkenSmallangleScatteringReveal2021}. This simple fabrication procedure offers wide tunability in compositions and dimensions. However, it generally introduces a wide variety of defects in the material, ranging from point defects to crystal surfaces and grain boundaries~\cite{ballDefectsPerovskitehalidesTheir2016, onoReducingDetrimentalDefects2020}. For the majority of these defects it is found that they have limited electronic effect on the MHP, with point defects~\cite{kangHighDefectTolerance2017, chuSoftLatticeDefect2020, chuLowfrequencyLatticePhonons2020, xueIntrinsicDefectsPrimary2022}, surfaces~\cite{haruyamaTerminationDependenceTetragonal2014, haruyamaSurfacePropertiesCH3NH3PbI32016, seiduAtomicElectronicStructure2021} and grain boundaries~\cite{guoStructuralStabilitiesElectronic2017, shanSegregationNativeDefects2017, mckennaElectronicProperties1112018, caiAtomicallyResolvedElectrically2022} mainly resulting in electronically benign defect levels. While it is common knowledge that defects induce instability problems in MHPs, for instance through a defect-driven accumulation of defects at grain boundaries~\cite{parkAccumulationDeepTraps2019, phungRoleGrainBoundaries2020, digirolamoIonMigrationInducedAmorphization2020} or material degradation at perovskite interfaces~\cite{apergiDecompositionOrganicPerovskite2022}, the understanding of the dynamics of the degradation reactions at an atomistic scale and their impact on the long-term stability of the materials and devices is limited.

Until now, the common understanding in the literature has been that during the degradation of MHPs the material disintegrates into its precursors. Specifically, various experiments have demonstrated the formation of \ch{PbX2} and amorphous \ch{PbX_{2-x}} (with \ch{X} = \ch{I} or \ch{Br}) using X-ray diffraction (XRD) measurements in the degradation of MHPs~\cite{coningsIntrinsicThermalInstability2015, fanLayerbyLayerDegradationMethylammonium2017, yangComprehensiveUnderstandingHeatinduced2018, albertiPbClusteringPbI22019, luDecisiveInfluenceAmorphous2021}. Moreover, metallic lead has also been found as a degradation product in MHPs~\cite{albertiPbClusteringPbI22019, luDecisiveInfluenceAmorphous2021}. Using electron microscopy it has been established that the degradation of halide perovskites tends to occur at surfaces and grain boundaries~\cite{albertiPbClusteringPbI22019, manekkathodiObservationStructuralPhase2020, luCompositionalInterfaceEngineering2020}. For that reason, a variety of efforts have been devoted to minimizing the occurrence of such grain boundaries in perovskite films to enhance the stability of the materials and devices~\cite{yenLargeGrainedPerovskite2016, seoIonicLiquidControl2016, yadavalliMechanismsExceptionalGrain2020, aliEnlargingCrystalGrains2020}. However, despite this knowledge on the perovskite films and the decomposition products, a detailed understanding of the degradation pathways of halide perovskites is still lacking.

In this work, our objective is to provide insights into the role of surfaces and grain boundaries in the long-term stability of MHPs. To do so, we perform reactive molecular dynamics simulations, using a reactive force field (ReaxFF) we developed for CsPbI3~\cite{polsAtomisticInsightsDegradation2021}. This ReaxFF force field makes use of a dynamical bond order to describe the breaking and creation of bonds~\cite{vanduinReaxFFReactiveForce2001, chenowethDevelopmentApplicationReaxFF2008, senftleReaxFFReactiveForcefield2016}. The simulations allow us to characterize the evolution of the atomic species under thermal stress and to establish a stability trend of a variety of surfaces. Furthermore, based on the simulation trajectories, we establish what structural features make a surface stable or unstable, and in case of decomposition, through which atomistic processes it proceeds. Finally, we find that additional point defects are generally detrimental to perovskite stability, whereas grain boundaries can have either positive or negative effects on the stability of the perovskite lattice.

\section{Results and Discussion}

In the following sections, we present our findings on the stability of the surfaces of inorganic \ch{CsPbI3}. In section~\ref{sec:structural_models} we show the structural models of the surfaces and grain boundaries that we investigate in this work. We specifically highlight the equivalence between the orthorhombic and cubic surface models. In section~\ref{sec:phase_stability_surfaces}, before we show the defect-induced chemical instability, we first illustrate the effects of surfaces on the phase stability of \ch{CsPbI3} by comparing the structural details of the surfaces with those in the bulk. Section~\ref{sec:surface_stability} evaluates the stability of different perovskite surfaces under thermal stress and elaborates on the mechanism of defect-induced degradation. We then explore the effects of additional defects located at surfaces, by analyzing the stability of the surfaces of perovskites with defects in section~\ref{sec:surface_defect_effects}. Finally, in section~\ref{sec:grain_boundaries}, we investigate the effects of grain boundaries on the stability of the perovskite lattice.

\subsection{Structural models} \label{sec:structural_models}

To model the surfaces, we use slab models, shown in Figure~\ref{surface_terminations}a-c. These include \ch{CsPbI3} slabs with the (110), (020) and (202) planes of orthorhombic \ch{CsPbI3} exposed. Of these surfaces, the (110) and (020) planes of orthorhombic \ch{CsPbI3} are those most prevalent in XRD experiments~\cite{suttonCubicOrthorhombicRevealing2018}. Despite the limited occurrence of the orthorhombic (202) plane in experiments, we include a slab with this orientation for completeness. This allows for the investigation of slabs with corners, edges, and faces of the \ch{PbI_{x}} octahedra exposed, which is the case for the (110), (020) and (202) orthorhombic surfaces, respectively. Additionally, each surface is created with varying terminations. We differentiate between the following slab terminations: stoichiometric, \ch{Pb}-poor, \ch{Pb}-rich, an explanation for these names can be found in Supporting Note 1. We emphasize that although we here refer to the slabs and their exposed planes in the orthorhombic form, all slabs attain cubic structures at elevated temperatures due to thermal fluctuations. For comparison, the equivalent cubic structure found when heating each orthorhombic surface is also shown in Figure~\ref{surface_terminations}. Additional details of the surface models can be found in Supporting Note 1.

\begin{figure*}[htbp!]
    \includegraphics{images/Figure_1_Surface_Terminations_V3-2.pdf}
    \caption{\label{surface_terminations} Structural models of the different surface orientations for orthorhombic \ch{CsPbI3} perovskite slabs. (a) (110) / (100), (b) (020) / (110) and (c) (202) / (111) surfaces are shown for the orthorhombic (left) / cubic (right) phases. For each of the slabs the exposed surface is shown at the top.}
\end{figure*}

A variety of \ch{CsPbI3} grain boundaries is studied to find their effects on the stability of halide perovskites. We focus on the $\mathrm{3 \Sigma (112) (0.4, 0)}${}, $\mathrm{5 \Sigma (210) (0.4, 0)}${} and $\mathrm{3 \Sigma (111) (0, 0)}${} grain boundaries, which are created from the cubic phase of \ch{CsPbI3}~\cite{guoStructuralStabilitiesElectronic2017}. We specifically investigate these grain boundary models because they do not contain unsaturated atoms with dangling bonds, the effects of which are probed with the above surface models. Additional details on the creation and naming of the grain boundaries can be found in Supporting Note 2.

\subsection{Phase stability near surface} \label{sec:phase_stability_surfaces}

To assess the effects of surfaces on the bulk structure, we simulate \ch{CsPbI3} slabs of the most commonly encountered (110) orthorhombic surface~\cite{suttonCubicOrthorhombicRevealing2018}, which is equivalent to the (100) cubic surface~\cite{eperonInorganicCaesiumLead2015, huBismuthIncorporationStabilized2017, haqueEffectsHydroiodicAcid2018, geSilverIodideInduced2019, liuModifyingSurfaceTermination2021}. Slabs with a thickness of 4, 6, 8 or 10 octahedral cages and a \ch{Pb}-poor or \ch{Pb}-rich termination are simulated at a constant temperature of \SI{300}{\K}. We use the \ch{Pb}-\ch{I}-\ch{Pb} valence angle ($\theta$) to probe the effects of the surfaces on the bulk structure of the perovskite during the simulation. The time-averaged values of the \ch{Pb}-\ch{I}-\ch{Pb} angles oriented in the direction perpendicular to the surfaces are shown in Figure~\ref{surface_effects} together with the structural model of the slab with 10 layers. The structural models of slabs with a smaller thickness can be found in Supporting Note 3.

\begin{figure*}[htbp!]
    \includegraphics{images/Figure_2_Structural_Effects_V4.pdf}
    \caption{\label{surface_effects} Time-averaged values of the \ch{Pb}-\ch{I}-\ch{Pb} valence angle perpendicular to the surfaces grouped together along the depth of the slabs, with $z = 0$ indicating the center of the slab for the (a) \ch{Pb}-poor and (b) \ch{Pb}-rich terminations. The red, orange, green and blue lines indicate the slabs of 4, 6, 8 and 10 octahedral cages thick, respectively. The time-averaged structure of the 10 layer slabs of \ch{CsPbI3} are shown next to the graphs. The values of the valence angle as obtained from DFT calculations for the orthorhombic (\SI{148}{\degree}) and cubic (\SI{180}{\degree}) phases are indicated with the dotted and dashed lines, respectively.}
\end{figure*}

We find that regardless of the termination, the \ch{Pb}-\ch{I}-\ch{Pb} valence angle near the center of the perovskite slab has an on average lower value ($\theta_{\mathrm{MD}} = \SI{155}{\degree}$) than close to the surface ($\theta_{\mathrm{MD}} = 160-175\SI{}{\degree}$). A comparison of these valence angles to those found for bulk structures using density functional theory (DFT) shows us the center of the slab is orthorhombic ($\theta_{\mathrm{DFT}} = \SI{148}{\degree}$), while the regions near the surface tend to exhibit a more cubic structure ($\theta_{\mathrm{DFT}} = \SI{180}{\degree}$). We hypothesize that this cubic-like surface layer can explain the enhanced stability of nanostructured \ch{CsPbI3} in experiments~\cite{shpatzdayanEnhancingStabilityPhotostability2018, zhaoSizeDependentLatticeStructure2020} against the conversion into the nonperovskite yellow phase~\cite{stoumposSemiconductingTinLead2013, stoumposRenaissanceHalidePerovskites2015, marronnierAnharmonicityDisorderBlack2018}. It has been shown that the large distortion in the orthorhombic phase of the perovskite results in poor \ch{Cs}-\ch{I} contacts in the structure which destabilizes the perovskite lattice~\cite{strausUnderstandingInstabilityHalide2020}. As such, we reason that as a result of nanostructuring, these cubic-like surface layers provide structural stability. Furthermore, we find that the size of this cubic-like surface layer is finite, approximately \SI{2}{\nm} thick. Notably, this observation agrees with experimental findings for \ch{CsPbI3}, where progressively larger nanocrystals undergo a transition from a cubic-like phase to an orthorhombic phase~\cite{zhaoSizeDependentLatticeStructure2020}.

\subsection{Surface stability} \label{sec:surface_stability}

To investigate the thermal stability of perovskite surfaces, we subject a collection of perovskite slabs (Figure~\ref{surface_terminations}) to a range of different temperatures (\SI{300}{\K} to \SI{700}{\K}). From these simulations we determine the onset temperature of material degradation, which we use as an indication for the stability of the perovskite surfaces. In the definition of material degradation we include all processes that result in a crystal lattice that deviates from the pristine form, which include but are not limited to the formation of defect complexes, the clustering of atomic species and the breakaway of atoms from the surface. The onset temperatures for the degradation of all slabs can be found in Table~\ref{onset_temperatures}.

\begin{table}[htbp!]
    \caption{\label{onset_temperatures} Onset temperatures for lattice degradation of perovskite slabs with a different orientations and terminations. For each surface the octahedral feature that protrudes the surfaces is indicated. The (*) denotes a stable surface up to and including \SI{700}{\K} and (-) the absence of such a surface in the investigations.}
    \begin{tabular}{ccccc}
        \toprule
        Orientations        & Surface feature   & Stoichiometric    & \ch{Pb}-poor      & \ch{Pb}-rich      \\ \midrule
        (110)               & Corner            & -                 & *                 & \SI{550}{\K}      \\
        (020)               & Edge              & \SI{550}{\K}      & \SI{300}{\K}      & \SI{300}{\K}      \\
        (202)               & Face              & -                 & \SI{400}{\K}      & \SI{300}{\K}      \\ \bottomrule
    \end{tabular}
\end{table}

We find that (110) oriented orthorhombic slabs possess the highest resistance to thermal stress. The \ch{Pb}-poor surface is stable for temperatures up to \SI{700}{\K}, with the \ch{Pb}-rich surface only showing decomposition of the lattice from \SI{550}{\K} and higher. This stability matches the trend in surface formation energies observed for the equivalent cubic surfaces in the ReaxFF force field validation (see Supporting Note 4). A lower thermal stability is observed for the (020) orthorhombic surface. The stoichiometric (020) surface of \ch{CsPbI3} is stable up to \SI{450}{\K} with the rearrangement of surface iodine atoms occurring at \SI{500}{\K} (details in Supporting Note 5) and degradation of the surface occurring at \SI{550}{\K} and higher. Both the \ch{Pb}-rich and \ch{Pb}-poor termination of the (020) orthorhombic surface are unstable, exhibiting the clustering of atomic species near the surface at \SI{300}{\K}. Finally, the (202) orthorhombic surface is the least stable, with the \ch{Pb}-poor termination degrading from \SI{400}{\K} and up and the \ch{Pb}-rich termination already decomposing at \SI{300}{\K}. On the basis of these results, we rank the surface orientations from most stable to least stable as: (110) $>$ (020) $>$ (202). We emphasize that this observed stability trend correlates well with the appearance of these surfaces in XRD experiments, in which the most stable surfaces appear predominantly~\cite{suttonCubicOrthorhombicRevealing2018, eperonInorganicCaesiumLead2015, huBismuthIncorporationStabilized2017, haqueEffectsHydroiodicAcid2018, geSilverIodideInduced2019, liuModifyingSurfaceTermination2021}.

For more insights into the degradation dynamics, we take a closer look at the degradation of a (110) orthorhombic slab simulated at \SI{600}{\K}. In Figure~\ref{degradation_dynamics}a-c we show snapshots of the lattice degradation during the simulation. Consistent with our earlier observations, the snapshots show that under these conditions the (110) \ch{Pb}-poor surface remains intact, while the \ch{Pb} rich surface shows decomposition, leading to the formation of a \ch{Pb_{x}I_{y}} complex. To quantify the degradation processes, we analyze the radial distribution functions (RDFs) between the \ch{Pb}-\ch{Pb} atom pairs in the simulated system in Figure~\ref{degradation_dynamics}d-e and for the remaining atom pairs in Supporting Note 6. On the \ch{Pb}-rich surface, the \ch{Pb}-\ch{Pb} RDF shows a simultaneous increase in peak at \SI{4.2} {\angstrom} and a decrease in peak at \SI{6.3} {\angstrom}, the timing of which coincides with the onset of degradation seen in Figure~\ref{degradation_dynamics}a-c. On the contrary, the \ch{Pb}-poor surface does not show any change in peak intensity and therefore shows no signs of decomposition. Based on these observations, we conclude that the degradation of the perovskite lattice can be characterized by the clustering of \ch{Pb} atoms, which has been observed in experiments at perovskite grain boundaries before~\cite{albertiPbClusteringPbI22019}. By comparing the position of the emerging \ch{Pb}-\ch{Pb} peak (\SI{4.2}{\angstrom}) to the interatomic distances of \ch{Pb} species from DFT in layered \ch{PbI2} (\SI{4.55}{\angstrom}) or yellow phase \ch{CsPbI3} (\SI{4.65}{\angstrom}), we can conclude that the decomposition product is not a pure form of either of these two material phases, but a more amorphous \ch{Pb_{x}I_{y}} domain. The formation of these domains matches observations from several experimental reports in which the formation of \ch{PbI2} domains is observed in electron microscopy experiments at grain boundaries~\cite{albertiPbClusteringPbI22019, manekkathodiObservationStructuralPhase2020, luDecisiveInfluenceAmorphous2021}.

\begin{figure*}[htbp!]
    \includegraphics{images/Figure_3_Surface_Degradation_V2.pdf}
    \caption{\label{degradation_dynamics} Degradation of a (110) oriented orthorhombic perovskite surface at \SI{600}{\K}, with the \ch{Pb}-rich surface on the top and the \ch{Pb}-poor surface at the bottom of the slab. (a-c) Structural snapshots of the slab are shown from \SI{0.0}{\ns} to \SI{2.0}{\ns}. (d-e) Time evolution of the \ch{Pb}-\ch{Pb} radial distribution function (RDF) for the \ch{Pb}-rich and \ch{Pb}-poor surfaces of the perovskite slab, demonstrating that the degradation of the \ch{Pb}-rich surface starts from a clustering of \ch{Pb} species.}
\end{figure*}

From all the decomposing \ch{CsPbI3} surfaces reported above, we observe that all surfaces degrade through a similar mechanism. To illustrate the key steps of this mechanism, we show the degradation of the \ch{Pb}-rich orthorhombic (110) surface at \SI{600} {\K} as an example in Figure~\ref{degradation_mechanism}. In the first step of the degradation process, an iodine Frenkel defect is formed in the perovskite lattice (Figure~\ref{degradation_mechanism}a-c). As a result of this, two \ch{PbI_{x}} octahedra form an edge-sharing complex, as opposed to the corner-sharing geometry in a regular perovskite lattice. Here, two \ch{Pb} atoms are bound together by two \ch{I} atoms. This edge-sharing complex acts as a metastable state in our simulations, exhibiting lifetimes of up to \SI{50}{\ps}. In the next step, an additional \ch{I} atom is added to the edge-sharing complex, transforming it into a face-sharing complex, in which the interatomic distance of the \ch{Pb} atoms involved in the complex significantly decreases (Figure~\ref{degradation_mechanism}d). We regard the formation of this face-sharing complex as the starting point of the degradation of the surface. Shortly after its formation, this face-sharing complex breaks away from the surface, initiating the decomposition of the perovskite lattice near the surface (Figure~\ref{degradation_mechanism}e-f).

We highlight that the critical steps of the degradation mechanism of perovskite surfaces (Figure~\ref{degradation_mechanism}g-i) resemble the degradation near iodine vacancies in the bulk of \ch{CsPbI3}~\cite{polsAtomisticInsightsDegradation2021}. However, we note that two distinct features make perovskite surfaces more prone to degradation: 1) an abundance of dangling bonds and 2) a lack of steric hindrance. Using DFT calculations it has been shown that the density and type of dangling bonds affects the formation energies of perovskite surfaces~\cite{wangDensityFunctionalStudies2015}. Here, we posit that such dangling bonds impact not only the thermodynamic stability of these surfaces but also the dynamical stability. In particular, we find that undercoordinated \ch{Pb} species readily form new bonds that then progressively degrade the halide perovskite. Additionally, the stability trend of the perovskite surfaces found earlier can be classified according to the protruding surface features as: corner $>$ edge $>$ face. We stress that in this ranking the \ch{PbI_{x}} octahedra, and specifically the \ch{I} species, experience an increasingly smaller steric hindrance from the \ch{Cs} species at the surface, allowing for the movement of these octahedra and thus a larger tendency for perovskite degradation. Finally, these two factors also explain the high thermal stability of the \ch{Pb}-poor (110) orthorhombic surface. In addition to the absence of \ch{Pb} dangling bonds, the \ch{Cs} atoms also sterically hinder the movement of the undercoordinated \ch{I} species at the surface. Together, this inhibits the decomposition, making the surface very stable.

\begin{figure*}[htbp!]
    \includegraphics{images/Figure_4_Surface_Degradation_TwoColumn_V2-1.pdf}
    \caption{\label{degradation_mechanism} Snapshots of the degradation of a \ch{Pb}-rich (110) orthorhombic perovskite surface. (a) Non-degraded perovskite surface. (b-c) A relatively long-lived iodine Frenkel defect at the perovskite surface, resulting in an edge-sharing complex. (d-e) A face-sharing complex at the surface that starts to break away from the surface. (f) Decomposed perovskite lattice at the surface. (g-i) Schematic representation of the critical steps in the degradation mechanism where neighboring \ch{PbI_{x}} octahedra change from corner- to edge- to face-sharing. The species highlighted in pink are the \ch{I} atom involved in the formation of the surface defect. The figure makes use of a shifted time axis, with $t' = \SI{0}{\ps}$ corresponding to $t = \SI{0.75}{\ns}$ in simulation time.}
\end{figure*}

\subsection{Effect of additional point defects on surfaces} \label{sec:surface_defect_effects}

To assess the effects of additional point defects on the stability of \ch{CsPbI3} perovskite surfaces, we use the \ch{Pb}-poor (110) orthorhombic surface as a model system because this perovskite surface has shown high thermal stability (Section~\ref{sec:surface_stability}). Since the surface is dominated by \ch{Cs} and \ch{I} species, with theoretical support for the occurrence of vacancies of both species~\cite{longEffectSurfaceIntrinsic2019}, we look at the effect of such vacancies at \SI{600}{\K}, individually and when they are clustered. In the remainder of this work we refer to the \ch{Cs} and \ch{I} vacancies as $\mathrm{V_{\ch{Cs}}}${} and $\mathrm{V_{\ch{I}}}${}. The structural models used for the investigation of these defects can be found in Supporting Note 7.

In their isolated form, the defects do not significantly affect the stability of the perovskite lattice. The $\mathrm{V_{\ch{Cs}}}${} defect remains on the surface of the slab, where it migrates across the surface as shown in Figure~\ref{defect_effects}a, without resulting in the degradation of the lattice. In contrast, the motion of $\mathrm{V_{\ch{I}}}${} is not bound to the surface of the perovskite. As shown in Figure~\ref{defect_effects}b, a $\mathrm{V_{\ch{I}}}${} defect can migrate into and out of the bulk of the perovskite. Moreover, in some cases the $\mathrm{V_{\ch{I}}}${} defect causes the perovskite surface to degrade; however, this process is not restricted to perovskite surfaces and also occurs in the bulk of inorganic perovskites~\cite{polsAtomisticInsightsDegradation2021}. Whenever a $\mathrm{V_{\ch{I}}}${} and a $\mathrm{V_{\ch{Cs}}}${} defect are created in a closely spaced pair, as shown in Figure~\ref{defect_effects}c-d, the stability of the perovskite surface is significantly affected. During simulation, we observe that this defect cluster tends to stay together at a fixed position in the lattice, only occasionally splitting into isolated $\mathrm{V_{\ch{Cs}}}${} and $\mathrm{V_{\ch{I}}}${} defects. When the defect pair remains clustered, the pair is particularly detrimental for the stability of the perovskite lattice. Specifically, we observe \ch{Pb} species readily moving away from their original position the lattice in close proximity to the defect pair ($t = \SI{0.5}{\ns}$ in Figure~\ref{defect_effects}), initiating the degradation of the lattice by allowing \ch{Pb} and \ch{I} species to cluster and form an amorphous \ch{Pb_{x}I_{y}} domain. We connect the tendency for a defect pair to act as a degradation center to the earlier-established factors for perovskite degradation: the presence of \ch{Pb} dangling bonds and limited steric hindrance for \ch{I} species locally around the defect pair.

\begin{figure*}[htbp!]
    \includegraphics{images/Figure_5_Defect_Migration_V2-3.pdf}
    \caption{\label{defect_effects} Snapshots of different point defects in \ch{CsPbI3} slabs. (a) Top view of a surface with a $\mathrm{V_{\ch{Cs}}}${} defect, showing the motion of the defect across the surface. (b) Side view of a surface with a $\mathrm{V_{\ch{I}}}${} defect, showing the migration of the defect away from the surface into the bulk. A top (c) and side view (d) of a perovskite surface with a $\mathrm{V_{\ch{Cs}}}${} and $\mathrm{V_{\ch{I}}}${} defect pair, showing the degradation of the perovskite surface caused by the presence of this defect at the surface. The $\mathrm{V_{\ch{Cs}}}${} and $\mathrm{V_{\ch{I}}}${} defects are indicated with a blue and red sphere, respectively.}
\end{figure*}

\subsection{Grain boundaries} \label{sec:grain_boundaries}

We assess the stability of grain boundaries from simulations at \SI{600}{\K}. The structure of the model systems after \SI{200}{\ps} of simulation is shown in Figure~\ref{gb_evolution}. The $\mathrm{3 \Sigma (112) (0.4, 0)}${} (Figure~\ref{gb_evolution}a) and $\mathrm{5 \Sigma (210) (0.4, 0)}${} (Figure~\ref{gb_evolution}b) grain boundaries exhibit clustering of atomic species in the grain boundary region. This clustering results in the formation of amorphous \ch{Pb_{x}I_{y}} domains, an observation that is consistent with the surfaces presented above (Section~\ref{sec:surface_stability}) and the experimentally observed degradation of perovskites at grain boundaries~\cite{manekkathodiObservationStructuralPhase2020, luDecisiveInfluenceAmorphous2021, albertiPbClusteringPbI22019}. As a result of the absence of dangling bonds at the grain boundary, we attribute the material instability to the lack of steric hinderance, that is the lack of \ch{Cs} species that could block the clustering of closely-spaced \ch{Pb} and \ch{I} species. Contrary to the other two grain boundaries, the $\mathrm{3 \Sigma (111) (0, 0)}${} grain boundary (Figure~\ref{gb_evolution}a), also known as a twinning plane~\cite{mckennaElectronicProperties1112018}, does not show any degradation throughout the duration of the simulation (\SI{2}{\ns}). Interestingly, the surface analogues of this grain boundary, the (202) orthorhombic surface or the equivalent high temperature (111) cubic surface, are unstable from \SI{400}{\K}. Unlike the surface, which has dangling bonds, the $\mathrm{3 \Sigma (111) (0, 0)}${} grain boundary features intact face-sharing \ch{PbI_{x}} octahedra at the grain boundary. In addition, the cavity filling \ch{Cs} species sterically hinder the movement of these face-sharing \ch{PbI_{x}} octahedra, further stabilizing the $\mathrm{3 \Sigma (111) (0, 0)}${} grain boundary.

\begin{figure*}[htbp!]
    \includegraphics{images/Figure_6_Grain_Boundary_Degradation_V3-1.pdf}
    \caption{\label{gb_evolution} Dynamical evolution of the \ch{CsPbI3} grain boundaries at \SI{600}{\K} after \SI{200}{\ps} from their initial structure. (a) $\mathrm{3 \Sigma (112) (0.4, 0)}${} grain boundary. (b) $\mathrm{5 \Sigma (210) (0.4, 0)}${} grain boundary. (c) $\mathrm{3 \Sigma (111) (0, 0)}${} grain boundary. The blue areas highlight the grain boundaries in the structures.}
\end{figure*}

\subsection{Conclusion}

In summary, using a ReaxFF force field, we study structural and thermal stability effects of surfaces and grain boundaries in the inorganic halide perovskite \ch{CsPbI3}. We show that surfaces affect the crystal phase close to the surface, which attains a cubic-like structure that is approximately \SI{2}{\nm} in size. We believe that this surface region is responsible for the enhanced structural stability of nanostructured \ch{CsPbI3}. Under a thermal stress that ranges from \SI{300}{\K} to \SI{700}{\K}, we find a stability trend for orthorhombic \ch{CsPbI3} of: (110) $>$ (020) $>$ (202). This trend matches the occurrence of these surfaces in experiments, in which the most stable surfaces are most predominantly observed. Comparing all investigated structures, we propose two important factors that are responsible for the degradation of perovskites: 1) the presence of dangling bonds, particularly for \ch{Pb} species and 2) a lack of steric hindrance, especially for \ch{I} species. These two factors explain the high thermal stability observed for the \ch{Pb}-poor (110) orthorhombic surface. The surface only has \ch{I} dangling bonds, not \ch{Pb} dangling bonds, and the motion of the \ch{I} species is blocked by the \ch{Cs} species on the surface, resulting in a stable surface. All other surfaces, which do not satisfy these two conditions, eventually decompose through a common mechanism. In this mechanism first an iodine Frenkel defect is formed in close proximity to the surface; after some time this defect grows to form a complex of two face-sharing \ch{PbI_{x}} octahedra; finally the complex breaks away from the surface, growing larger by accumulating more \ch{Pb} and \ch{I} species, leading to the formation of a large amorphous \ch{Pb_{x}I_{y}} domain on the surface.

The stability of aforementioned surfaces deteriorates when additional point defects ($\mathrm{V_{\ch{I}}}${} and $\mathrm{V_{\ch{Cs}}}${}) are created on the surfaces, which results in the formation of \ch{Pb} dangling bonds and a decrease of steric hindrance locally. The latter factor also explains the stability of the investigated grain boundaries. For example, although the $\mathrm{3 \Sigma (112) (0.4, 0)}${} and $\mathrm{5 \Sigma (210) (0.4, 0)}${} grain boundaries do not contain any dangling bonds, the lack of steric hindrance for \ch{Pb} and \ch{I} species facilitates the clustering of these species, leading to the degradation of the material. On the contrary, the $\mathrm{3 \Sigma (111) (0, 0)}${} grain boundary contains cavity filling \ch{Cs} species that sterically hinder the typical degradation-inducing movement of the face-sharing octahedra, making it stable up to \SI{600}{\K}.

Based on the above, we propose that strategies to stabilize halide perovskites can include the following aspects: i) passivating defects, primarily those leading to unsaturated bonds, i.e. halide vacancies, through passivating agents such as halogens like \ch{F}~\cite{liCationAnionImmobilization2019a} and \ch{Cl}~\cite{gaoHalideExchangePassivation2022} or carbonyl~\cite{xieDecouplingEffectsDefects2021} and azo~\cite{maStablePerovskiteSolar2021} containing ligands; ii) the grafting of surfaces with sterically hindering groups, such as phenylalkylammonium~\cite{wangAllinorganicCesiumLead2018, guoPhenylalkylammoniumPassivationEnables2021, songSurfaceStabilizationFormamidinium2021} and even bulkier organic groups~\cite{wuEfficientStableCsPbI32019}; (iii) optimizing the synthesis conditions, through the type of precursors, solvents, additives etc., to stimulate the growth of non-detrimental surface orientations~\cite{foleyControllingNucleationGrowth2016}, reduce the formation of grain boundaries~\cite{peteramalathasControlledGrowthLarge2020} and suppress the formation of defects altogether~\cite{yangDefectSuppressionOriented2021}.

\section{Computational details}

\subsection{Molecular dynamics}

ReaxFF molecular dynamics simulations were performed in AMS2021~\cite{rugerAMS2021}. All ReaxFF simulations were done using the earlier developed \ch{CsPbI3} ReaxFF force field~\cite{polsAtomisticInsightsDegradation2021}, a validation of this force field for \ch{CsPbI3} surfaces can be found in Supporting Note 5. Before the molecular dynamics runs, all structural models were optimized with the ReaxFF force field. The dynamical simulations used a simulation timestep of \SI{0.25}{\fs}. The thermostat and barostat use a damping constant of $\tau_{T} = \SI{100}{\fs}$ and $\tau_{p} = \SI{2500}{\fs}$, respectively. Whenever slabs were simulated, the vacuum layer used was at least \SI{50}{\angstrom} in the $z$-direction. In the case of slabs, the barostat was only allowed to scale the non-vacuum directions of the model system ($x$- and $y$-directions). The initial velocities of the particles were assigned according to a Maxwell-Boltzmann distribution of the initial temperature. Simulation snapshots and structural models were all visualized using OVITO~\cite{stukowskiVisualizationAnalysisAtomistic2009}.

In the simulations investigating the structural effects of surfaces on the bulk structure of perovskites, we equilibrated the systems to the target temperature of \SI{300}{\K} in an NPT-ensemble for \SI{200}{\ps}. During the equilibration stage we employed a Berendsen thermostat and Berendsen barostat~\cite{berendsenMolecularDynamicsCoupling1984} to control the temperature and pressure, respectively. Production runs were started from the final frame of the equilibration run and took \SI{200}{\ps}. In the production runs, the temperature and pressure were controlled with a NHC-thermostat~\cite{martynaNoseHooverChains1992} with chain length of 10 and an MTK barostat~\cite{martynaConstantPressureMolecular1994}. The time-averaged values of the \ch{Pb}-\ch{I}-\ch{Pb} valence angles were extracted from the time-averaged structure we obtained by averaging the atomic positions over the full duration of the production simulation.

The stability of the perovskite surfaces and grain boundaries was investigated using a three-step approach. In the first two steps, the equilibration stage, we used a Berendsen thermostat and Berendsen barostat to control the temperature and pressure. The first step was used to slowly heat the system from \SI{300}{\K} to the desired target temperature during \SI{100}{\ps}, with the second step maintaining the system at its constant target temperature for \SI{100}{\ps}. The full equilibration of the system was run in the non-reactive mode of ReaxFF in AMS2021, in which the bonds can only be updated but not newly formed, to prevent unwanted reactions during the equilibration of the system. The final stage, the production simulation, was run in an NPT-ensemble for which the starting point was the final frame of the equilibration. The temperature and pressure were controlled, respectively, with an NHC-thermostat with a chain length of 10 and an MTK-barostat. Each of the production runs was \SI{2}{\ns} long.

\subsection{Density Functional Theory}

Density functional theory (DFT) calculations were performed with the projector augmented wave (PAW) method as implemented in the Vienna Ab-Initio Simulation Package (VASP)~\cite{kresseInitioMolecularDynamics1993, kresseInitioMoleculardynamicsSimulation1994, kresseEfficiencyAbinitioTotal1996, kresseEfficientIterativeSchemes1996, kresseUltrasoftPseudopotentialsProjector1999}.  The electron exchange-correlation interaction was described using the Perdew, Burke and Ernzerhof (PBE) functional~\cite{perdewGeneralizedGradientApproximation1996, perdewGeneralizedGradientApproximation1997} with long-range dispersive interactions accounted for by the DFT-D3(BJ) dispersion correction~\cite{grimmeEffectDampingFunction2011}. We treated the outermost electrons of \ch{Cs} (5s\textsuperscript{2}5p\textsuperscript{6}6s\textsuperscript{1}); \ch{Pb} (5d\textsuperscript{10}6s\textsuperscript{2}6p\textsuperscript{2}) and \ch{I} (5s\textsuperscript{2}5p\textsuperscript{5}) as valence electrons with the plane wave basis set expanded to an energy cutoff of \SI{500}{\eV}.

The geometries were optimized by allowing all: the ionic positions, cell shape and cell volume to change until convergence of \SI{1E-3}{\meV} and \SI{10}{\meV\per\angstrom} was reached in energy and forces, respectively. The Brillouin zones were sampled using a Monkhorst-Pack mesh~\cite{monkhorstSpecialPointsBrillouinzone1976}, with the following $k$-space grids resulting in energy convergence to within \SI{1}{\meV}/atom: \ch{CsI}: $12 \times 12 \times 12$; \ch{PbI2}: $11 \times 11 \times 7$; cubic \ch{CsPbI3}: $10 \times 10 \times 10$; orthorhombic \ch{CsPbI3}: $7 \times 7 \times 5$; yellow phase \ch{CsPbI3}: $13 \times 6 \times 4$.

\begin{suppinfo}

Structural models of the surfaces; structural models of grain boundaries; time-averaged structures of 4, 6 and 8 octahedral layer thick perovskite slabs; ReaxFF force field validation using surface energies; illustration of the rearrangement of surface iodine for the stoichiometric (020) orthorhombic surface; explanation of RDFs and overview of the RDFs for remaining atom pairs; structural models of the defective surfaces

\end{suppinfo}

\section*{Supporting Notes}

\subsection{Structural models for surfaces}

\subsubsection{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~\ref{surface_ratios}. The naming scheme is based on the relative ratio of the \ch{I} and \ch{Pb} species at each surface, resulting in the following names: stoichiometric (\ch{I}:\ch{Pb} $=$ 3), \ch{Pb}-poor (\ch{I}:\ch{Pb} $>$ 3) and \ch{Pb}-rich (\ch{I}:\ch{Pb} $<$ 3)

\begin{table}[htbp!]
    \caption{\label{surface_ratios} Relative ratio of the occurrence of atomic species as determined for perovskite slabs with the same termination at both sides of the perovskite slab.}
    \begin{tabular}{ccccc}
        \toprule
        Orientation             & Termination       & \ch{Pb}:\ch{Cs}   & \ch{I}:\ch{Cs}    & \ch{I}:\ch{Pb}    \\ \midrule
        \multirow{2}{*}{(110)}  & \ch{Pb}-poor      & 0.89              & 2.78              & \textbf{3.13}     \\
                                & \ch{Pb}-rich      & 1.13              & 3.25              & \textbf{2.89}     \\ \midrule
        \multirow{3}{*}{(020)}  & Stoichiometric    & 1.00              & 3.00              & \textbf{3.00}     \\
                                & \ch{Pb}-poor      & 0.93              & 3.15              & \textbf{3.15}     \\
                                & \ch{Pb}-rich      & 1.07              & 2.86              & \textbf{2.86}     \\ \midrule
        \multirow{2}{*}{(202)}  & \ch{Pb}-poor      & 0.94              & 3.00              & \textbf{3.20}     \\
                                & \ch{Pb}-rich      & 1.06              & 3.00              & \textbf{2.81}     \\ \bottomrule
    \end{tabular}
\end{table}

\clearpage

\subsubsection{Surface models}

In this work the surface models are created in AMS2021~\cite{rugerAMS2021}. An orthorhombic unit cell is used as the starting point for the surface models. Each model slab is at least \SI{50}{\angstrom} thick, with the lateral dimension spanning at least \SI{45}{\angstrom}. An overview of the full structural models of the doubly terminated (\ch{Pb}-poor and \ch{Pb}-rich) slabs is shown in Figure~\ref{surface_models}.

\begin{figure*}[htbp!]
    \includegraphics{si-images/SI_1_Surface_Models_V1.pdf}
    \caption{\label{surface_models} Structural models of the \ch{Pb}-poor surfaces as created from orthorhombic \ch{CsPbI3}. (a) (110) orthorhombic surface. (b) (020) orthorhombic surface. (c) (202) orthorhombic surface. The \ch{Pb}-rich and \ch{Pb}-poor terminated surfaces are indicated with the red and blue color, respectively.}
\end{figure*}

\clearpage

\subsection{Structural models for grain boundaries}

\subsubsection{Grain boundary names}

The names of the grain boundaries follow from the work of Guo et al.~\cite{guoStructuralStabilitiesElectronic2017}. The grain boundaries in this work are created with the coincidence site lattice (CSL) model. In this approach, two grains are tilted with a certain angle, until their surface planes coincide. An example of such a grain boundary structure is the $\mathrm{3 \Sigma (112) (0.4, 0)}${} grain boundary. In this naming scheme the $3$ in front of $\Sigma$ denotes that every third point on the surface planes coincides. The number that follows $\Sigma$, $(112)$, denotes the lattice plane along which the intersection occurs. The final two numbers indicate the overall fractional shift of the grains with respect to each other. For $(0.4,0)$ this means that one of the grains is shifted along the $[100]$ direction by $40\%$ of the length of the unit cell vector and without any shift along the $[010]$ direction. Note these shift directions hold because the grain boundaries are formed perpendicular to the $z$-axis.

\clearpage

\subsubsection{Grain boundary models}

The grain boundary models in this work are created with the \texttt{aimsgb} package~\cite{chengAimsgbAlgorithmOpensource2018}. A cubic \ch{CsPbI3} 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~\ref{gb_models}. In each model the grain boundaries are separated by at least \SI{25}{\angstrom}. The lateral dimensions of the grain boundary models are at least \SI{30}{\angstrom} in size.

\begin{figure*}[htbp!]
    \includegraphics{si-images/SI_2_Grain_Boundary_Models_V4.pdf}
    \caption{\label{gb_models} Structural models of the grain boundaries as created from cubic \ch{CsPbI3}. (a) $\mathrm{3 \Sigma (112) (0.4, 0)}${} grain boundary. (b) $\mathrm{5 \Sigma (210) (0.4, 0)}${} grain boundary. (c) $\mathrm{3 \Sigma (111) (0, 0)}${} grain boundary. The blue areas highlight the grain boundaries in the structures.}
\end{figure*}

\clearpage

\subsection{Structural effects of surfaces}

\begin{figure*}[htbp!]
    \includegraphics{si-images/SI_3_Figure_Averaged_Slabs_V1.pdf}
    \caption{\label{averaged_slabs} Time-averaged perovskite slab structures for both \ch{Pb}-poor and \ch{Pb}-rich slabs of 4, 6 and 8 layers of octahedral cages.}
\end{figure*}

\clearpage

\subsection{Force field validation}

Slab geometries were created from optimized bulk structures of the cubic phase of \ch{CsPbI3} 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 \times 6 \times 1$ $k$-points mesh was used during slab optimizations. At least \SI{15}{\angstrom} 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~\ref{validation_structures}.

\begin{figure*}[htbp!]
    \includegraphics{si-images/SI_7_Validation_Structures_V1.pdf}
    \caption{\label{validation_structures} Structural models of the slabs used in the validation of the force field. (a) \ch{Pb}-poor terminated slab. (b) \ch{Pb}-rich terminated slab. The simulation cell is indicated with the black lines.}
\end{figure*}

Surface energies ($\gamma$) were calculated as
\begin{equation}
    \gamma = \frac{E_{\mathrm{slab}} - \left( n \cdot E_{\ch{CsPbI3}} \pm E_{\mathrm{precursor}} \right) }{2 A}
\end{equation}
where $E_{\mathrm{slab}}$ is the energy of the perovskite slab, $n$ the number of formal units of \ch{CsPbI3}, $E_{\mathrm{\ch{CsPbI3}}}$ is the energy of cubic \ch{CsPbI3}, $E_{\ch{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 \ch{CsI} as the precursor, resulting in $n = 6$ and an addition of the precursor energy ($+ E_{\mathrm{precursor}}$) for the \ch{Pb}-poor terminated slab and $n = 7$ and a subtraction of the precursor energy ($- E_{\mathrm{precursor}}$) for the \ch{Pb}-rich terminated slab.

An overview of the surface areas and calculated surface energies of the perovskite slabs with the \ch{Pb}-rich and \ch{Pb}-poor termination is shown in table~\ref{validation_values}. 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~\cite{seiduAtomicElectronicStructure2021, longEffectSurfaceIntrinsic2019}. In both computational methods the \ch{Pb}-poor surface termination is found to be significantly more stable than the \ch{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.

\begin{table}[htbp!]
    \caption{\label{validation_values} The areas ($A$) and surface energies ($\gamma$) for the different \ch{CsPbI3} slab terminations as calculated with DFT and ReaxFF.}
    \begin{tabular}{cccc}
        \toprule
        Surface termination     & $A$ ($ \SI{}{\angstrom\squared}$)  & $\gamma_{\mathrm{DFT}}$ ($\SI{}{\joule\per\meter\squared}$) & $\gamma_{\mathrm{ReaxFF}}$  ($\SI{}{\joule\per\meter\squared}$) \\ \midrule
        \ch{Pb}-poor            & 39.58                             & 0.095                                                     & -0.087                        \\
        \ch{Pb}-rich            & 39.58                             & 0.222                                                     & 0.144                         \\ \bottomrule
    \end{tabular}
\end{table}

\clearpage

\subsection{Surface rearrangements of iodine}

For the stoichiometric orthorhombic (020) surface it was observed that \ch{I} species would migrate across the surface from \SI{500}{\K} and up. To illustrate this behavior, some simulation snapshots from \SI{500}{\K} simulations are shown in figure~\ref{surface_rearrangements}. We underline that the jump of an \ch{I} atom to a neighboring \ch{Pb} atom (Figure~\ref{surface_rearrangements}a) forms local \ch{Pb}-rich and \ch{Pb}-poor octahedra, which become the degradation centers at elevated temperatures.

\begin{figure*}[htbp!]
    \includegraphics{si-images/SI_4_Surface_Rearrangement_V1.pdf}
    \caption{\label{surface_rearrangements} Snapshots of the rearrangement of surface iodine atoms for the orthorhombic (020) stoichiometric surface at \SI{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 \ch{Pb}-poor and \ch{Pb}-rich surface, respectively. A shifted time axis $t'$ is used.}
\end{figure*}

\clearpage

\subsection{Radial distribution functions}

The RDF between the atoms of the species $a$ and $b$ is calculated using
\begin{equation}
    g_{ab} \left( r \right) = \frac{1}{N_{a} N_{b}} \sum^{N_{a}}_{i = 1} \sum^{N_{b}}_{j = 1} \langle \delta \left( \lvert \vec{r}_{i} - \vec{r}_{j}  \rvert - r \right) \rangle,
\end{equation}
which sums over Dirac delta functions $\delta$ that are offset by the interatomic distance between the $i$th atom of species $a$ located at $\vec{r}_{i}$ and the $j$th atom of species $b$ located at $\vec{r}_{j}$. The quantity is averaged over the full length of the simulation indicated by the ensemble averaging ($\langle \cdots \rangle$).

In this work, we employ the radial distribution function (RDF) as implemented in the \texttt{MDAnalysis} package~\cite{michaud-agrawalMDAnalysisToolkitAnalysis2011, gowersMDAnalysisPythonPackage2016}. Surface specific RDFs are obtained by only including those atoms that are within 4 layers of \ch{Pb} atoms from the surface. The time information is obtained by dividing the trajectory of the degradation simulation into five blocks of equal length, each corresponding to \SI{400}{\ps} of data. Subsequently, the RDF is determined for each of these spatial and temporal blocks, resulting in RDFs with space- and time-resolved information on the degradation in our system. A demonstration of the partitioning into surface-specific regions and the remaining RDFs (atom pairs: \ch{Pb}-\ch{Pb}, \ch{Pb}-\ch{Cs}, \ch{Pb}-\ch{I}, \ch{Cs}-\ch{Cs}, \ch{Cs}-\ch{I} and \ch{I}-\ch{I}) of the \ch{Pb}- and \ch{Pb}-poor (110) orthorhombic surface of \ch{CsPbI3} are shown in Figure~\ref{rdfs_combined}.

\begin{figure*}[htbp!]
    \includegraphics{si-images/SI_5_RDF_Atom_Pairs_Combined_V1.pdf}
    \caption{\label{rdfs_combined} Time-resolved radial distribution functions (RDFs) for the (110) surface of orthorhombic \ch{CsPbI3}. The figure shows the RDFs of the \ch{Pb}-\ch{Pb}, \ch{Pb}-\ch{Cs}, \ch{Pb}-\ch{I}, \ch{Cs}-\ch{Cs}, \ch{Cs}-\ch{I} and \ch{I}-\ch{I} atom pairs. The RDFs in the top half correspond to the \ch{Pb}-rich surface (red) and the figures at the bottom half to the \ch{Pb}-poor surface (blue).}
\end{figure*}

\clearpage

\subsection{Defective surfaces}

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~\ref{defect_slabs}a, the top surface has a $\mathrm{V_{\ch{Cs}}}${} defect and the bottom surface a $\mathrm{V_{\ch{I}}}${} defect. For closely spaced pairs of $\mathrm{V_{\ch{Cs}}}${} and $\mathrm{V_{\ch{I}}}${} defects, Figure~\ref{defect_slabs}b, we create the pair of defects at the top and bottom of the slab. The slabs in the simulations are at least \SI{50}{\angstrom} thick. The lateral dimensions are also at least \SI{50}{\angstrom} wide.

\begin{figure*}[htbp!]
    \includegraphics{si-images/SI_6_Defect_Slabs_V1.pdf}
    \caption{\label{defect_slabs} Structural models used to study the effect of point defects on the stability of the \ch{Pb}-poor (110) orthorhombic perovskite surface. The blue and red circles are used to indicate the position of the \ch{Cs} and \ch{I} vacancy, respectively. (a) A slab with a single $\mathrm{V_{\ch{Cs}}}${} at the top surface and a single $\mathrm{V_{\ch{I}}}${} at the bottom surface. (b) A slab with closely-spaced defect pairs of $\mathrm{V_{\ch{Cs}}}${} and $\mathrm{V_{\ch{I}}}${} at both the top and bottom surface.}
\end{figure*}

\clearpage


\section{Introduction}

The study of temporal evolution of social structure is one of the significant application domains of complex collective systems research. In particular, social fragmentation transition, i.e., transition of social states between many disconnected communities with distinct opinions and a well-connected single network with homogeneous opinions, is a timely research topic with high relevance to various current societal issues \cite{holme2006,zanette2006,kozma2008,bohme2011,sayama2020b,blex2020,levin2021,sasahara2021}. We had previously studied this problem using numerical simulations of adaptive social network models \cite{sayama2020} and found that two individual behavioral traits, homophily (i.e., tendency to strengthen connections to similar individuals and weaken those to dissimilar ones) \cite{mcpherson2001,kossinets2009,bakshy2015} and attention to novelty (i.e., tendency to strengthen connections to individuals whose opinions stand out compared to others), had the most statistically significant impact on the outcomes of social network evolution \cite{sayama2020}. Specifically, when homophily was strong, the social network evolved into fragmented states of many disconnected clusters with diverse opinions, but when attention to novelty was strong, the social network evolved to well-connected yet informationally homogeneous states. However, the previous study was rather limited in terms of the range of parameter values examined, and possible interactions between multiple behavioral traits were largely ignored, especially about the other behavioral trait, social conformity (i.e., how strongly individuals assimilate themselves to social neighbors).

In this study, we examined a broader spectrum of social network dynamics through a larger-scale numerical experiment of the same model with expanded parameter sweep ranges by an order of magnitude in each parameter dimension, resulting in a total of 116,640 simulation runs. To capture nontrivial interactions among behavioral parameters, we modeled and visualized the dependence of outcome measures on the model parameters using artificial neural networks. Results show that, while the competition between homophily and attention to novelty is still the primary determinant of social fragmentation when social conformity behavior of individuals is weak, another transition plane emerges at an intermediate homophily level when individuals have strong social conformity behavior, which was not previously known.

In what follows, we describe the model of adaptive social networks and the settings and results of the large-scale numerical experiments. We further discuss implications of the results for social evolution and potential future research directions.

\section{Model}

Our original model \cite{sayama2020} describes distributed opinion dynamics on an adaptive
social network made of $n$ nodes. Adaptive networks \cite{gross2009,sayama2013camwa} are a class of dynamical network models in which node states and edge connectivities co-evolve in adaptation to each other. In our model, node $i$ has its own opinion state $x_i \in \mathbb{R}$. Nodes are connected through
weighted directed edges that represent the information flow from source to target nodes. The edge weight is denoted as $w_{ij} \in
\mathbb{R}_{\geq 0}$, where $i$ is the target node and $j$ is the source node.

The adaptive network dynamics, i.e., the co-evolution of node states and edge weights, are governed by the following differential equations:
\begin{align}
\frac{d x_i}{d t} &= c \left( \langle x \rangle_i - x_i \right) + \epsilon \label{dxdt}\\
\frac{d w_{ij}}{d t} &= h F_h(x_i, x_j) + a F_a(\langle x \rangle_i, x_j) \label{dwdt}\\
\langle x \rangle_i &= \frac{\sum_j w_{ij} x_j}{\sum_j w_{ij}} \label{localav}
\end{align}
Here, $\langle x \rangle_i$ (Eq.~(\ref{localav})) represents the weighted local average of neighbors' opinions (i.e., social norm) perceived by node $i$.
Parameter $c$ and noise term $\epsilon$ in Eq.~(\ref{dxdt}) represent the strength of social conformity and stochastic fluctuation of node states, respectively.
Parameters $h$ and $a$ in Eq.~(\ref{dwdt}) represent the strength of homophily and attention to novelty, respectively.
$F_h$ and $F_a$ in Eq.~(\ref{dwdt}) are functions that describe the increase/decrease of edge weights because of homophily and attention to novelty, respectively. $F_h$ and $F_a$ can be any functions that monotonically decrease (for $F_h$) or increase (for $F_a$) as the distance between the two arguments increase. In this study, we used the following simple functions for $F_h$ and $F_a$:
\begin{align}
F_h(x_i, x_j) &= \theta_h - |x_i - x_j| \\
F_a(\langle x \rangle_i, x_j) &= |\langle x \rangle_i - x_j| - \theta_a
\end{align}
Here $\theta_h$ and $\theta_a$ are the default values of $F_h$ and $F_a$, respectively, when the two given arguments are equal. These functions describe that the edge from node $j$ to node $i$
tends to become strengthened when $j$'s state is similar to $i$'s (i.e., homophily) and distant from the local average (i.e.,
attention to novelty), or weakened otherwise. We restricted $w_{ij}$
to be always nonnegative, and any negative values resulted from numerical simulation of Eq.~(\ref{dwdt}) would be rounded up to zero.

Simulating this adaptive social network model from a random initial condition produces a sequence of social network configurations in which node states (opinions) spread through social ties and edge weights (connection strengths) also change due to node states (an example is shown in Fig.~\ref{screenshot}). We implemented the numerical simulator of the model in Python 3.7
with NetworkX \cite{hagberg2008} and PyCX \cite{sayama2013}\footnote{The simulator code is available from the
  author upon request.}.

\begin{SCfigure}[][t]
\centering
\includegraphics[width=3in]{screenshot.png}
\caption{A snapshot of the proposed adaptive social network model visualized in the middle of a simulation. Colors of nodes represent their states (opinions), while the shade of edges represent their weights (connection strengths).}
\label{screenshot}
\end{SCfigure}

This model is known to exhibit social fragmentation transition, i.e., phase transition between fragmented and homogenized social network states, as the individuals' behavioral parameters are varied (Fig.~\ref{transition}) \cite{holme2006,zanette2006,kozma2008,bohme2011}. Our previous study \cite{sayama2020} showed that, when homophily ($h$) is stronger or attention to novelty ($a$) is weaker, the social network is more inclined to become fragmented into many disconnected small clusters with various opinion states (Fig.~\ref{transition}, left), and in the opposite settings social homogenization is more likely to occur (Fig.~\ref{transition}, right). Meanwhile, the potential effect of social conformity ($c$) was unclear in the previous analysis, which is the main focus of the present study.

\begin{figure}[tp]
\centering
\includegraphics[width=\columnwidth]{transition.pdf}
\caption{Examples of final states of the adaptive social network simulation with $n=100$ that demonstrate social fragmentation transition. Left: Fragmented state with large $h$ and small $a$. Right: Homogenized state with small $h$ and large $a$. Specific parameter settings are shown beneath each panel.}
\label{transition}
\end{figure}

\section{Experiments}

\subsection{Settings}

We conducted numerical simulations of the above adaptive social network model to systematically investigate the effects of individual behavioral parameters ($c$, $h$, $a$, $\theta_h$, and $\theta_a$) on the course of social network evolution. The parameter values used are as follows:
\begin{itemize}
\item Network size: $n \in \{30, 100, 300\}$
\item Behavioral parameters: $c, \; h, \; a, \; \theta_h, \; \theta_a \in \{0.003, 0.01, 0.03, 0.1, 0.3, 1.0\}$
\end{itemize}
The above range of values for behavioral parameters was an order-of-magnitude larger in each dimension than what was examined before \cite{sayama2020}. Each parameter value combination was simulated 5 times with independently generated random initial conditions. This resulted in a total of $3 \times 6^5 \times 5 = 116,640$ simulation runs, taking a substantial amount of computational time and resource. Simulations were thus conducted in parallel on four designated PCs for over a few months.

In each simulation, the initial configuration of the network was randomly generated so that every pair of nodes were connected by two directed edges (in both directions) with a randomly generated weight ($w_{ij} \in [0, 1]$) in each direction and each node had a random node state
sampled from the normal distribution ${\cal N}(0, 1)$. Equations (\ref{dxdt}) and (\ref{dwdt}) were numerically simulated using a simple Euler forward method with
time step size $\Delta t = 0.1$ for $t \in [0, 100]$. The stochastic behavior of node states represented by $\epsilon$ in Eq.~(\ref{dxdt}) was simulated by adding a random number
sampled from ${\cal N}(0, 0.1^2)$ to $x_i$ at each discrete time step $\Delta t$.

\subsection{Outcome measures}

At the end of each simulation run ($t=100$), we converted the final network configuration into an undirected network by replacing the two directed edges between each pair of nodes with a single undirected edge whose weight was the average of the original two edges' weights. Then the Louvain modularity maximization method \cite{blondel2008fast} was applied to the undirected network to detect community structure in the final network configuration. Within each detected community, we calculated the average node state (called
``average community state'' hereafter). Using the results of these steps, we calculated the following five
network metrics as final outcome measures:
\begin{enumerate}
\item Average edge weight (= arithmetic average of all the edge
  weights in the network)
\item Number of communities
\item Modularity of the community structure
\item Range of average community states (= difference between largest
  and smallest average community states)
\item Standard deviation of average community states
\end{enumerate}
These outcome measures were averaged over five independent simulation runs for each combinations of parameter values. The first three outcome measures capture the structural properties of the social network, while the last two capture the opinion diversity in the social network. When the social network is fragmented, the average edge weight takes a small value, while all the other measurements takes large values. The opposite pattern is realized when the social network is homogenized. This allows us to easily detect which state the adaptive social network evolved into in quantitative ways.

\section{Results}

In order to capture and visualize the effects of the five behavioral parameters on the five outcome measures (including possible nonlinear interactions among those behavioral parameters), we modeled the parameter-outcome mapping using artificial neural networks with Wolfram Research Mathematica 12's artificial neural network predictor \cite{mathematica}. Natural logarithms of the five behavioral parameter values were used as five-dimensional input vectors, and the five outcome measures obtained from simulation results under those parameter settings were used as five-dimensional output vectors. The combinations of these input and output vectors were gathered for the whole simulation runs for each network size ($n \in \{30, 100, 300\}$) and used as the data set to train an artificial neural network model for specific $n$. The performance goal of training was set to maximizing the accuracy of outcome prediction \cite{mathematica}.

Illustrative results with $n=300$ are shown as heat maps of each outcome measure in Fig.~\ref{phase-spaces-structure} (for final network structure) and Fig.~\ref{phase-spaces-state} (for final node states).
The competition between homophily ($h$) and attention to novelty ($a$) is still observed as the primary determinant of social fragmentation in a low-conformity ($c$) regime (top rows of all panels in Figs.~\ref{phase-spaces-structure} and \ref{phase-spaces-state}), seen as the diagonal phase transition plane in the plots. However, another vertical transition plane emerges at an intermediate homophily level in a high-conformity regime (bottom rows of all panels in Figs.~\ref{phase-spaces-structure} and \ref{phase-spaces-state}), which was not previously known. Similar patterns were observed for other outcome measures and network sizes. This new result shows that, when individuals' social conformity ($c$) is sufficiently strong, homogenization of the social network can occur even without attention to novelty. This implies that social conformity and attention to novelty, while very different in their intentions and actions at microscopic individual levels, have similar effects of promoting connections among individuals in an adaptive social network.

\begin{figure}[tp]
\includegraphics[width=\columnwidth]{phase-spaces-structure.pdf}
\caption{Phase diagrams of adaptive social network evolution in terms of network structure outcome measures. Each plot shows outcome dependence on homophily ($h$, horizontal axis), attention to novelty ($a$, vertical axis) and conformity ($c$, varied from top-left to bottom-right) modeled using artificial neural networks. (a) How the number of communities depends on $h$, $a$ and $c$. (b) How the modularity of the community structure depends on $h$, $a$ and $c$. Red and blue regions correspond to fragmented and homogenized network states, respectively. $n = 300$, $\theta_h = 0.1$, and $\theta_a = 0.1$. Similar patterns were observed for other outcome measures and network sizes.}
\label{phase-spaces-structure}
\end{figure}

\begin{figure}[tp]
\includegraphics[width=\columnwidth]{phase-spaces-state.pdf}
\caption{Phase diagrams of adaptive social network evolution in terms of node state outcome measures. Each plot shows outcome dependence on homophily ($h$, horizontal axis), attention to novelty ($a$, vertical axis) and conformity ($c$, varied from top-left to bottom-right) modeled using artificial neural networks. (a) How the range of average community states depends on $h$, $a$ and $c$. (b) How the standard deviation of average community states depends on $h$, $a$ and $c$. Red and blue regions correspond to fragmented and homogenized network states, respectively. $n = 300$, $\theta_h = 0.1$, and $\theta_a = 0.1$. Similar patterns were observed for other outcome measures and network sizes.}
\label{phase-spaces-state}
\end{figure}

The result shown above also reveals a previously unrecognized competition between social conformity ($c$) and homophily ($h$) when attention to novelty is weak (i.e., low-$a$ regions; near the bottom edge of each heat map in Figs.~\ref{phase-spaces-structure} and \ref{phase-spaces-state}). Namely, when $c$ is low social fragmentation dominates, but when $c$ is high social homogenization becomes possible for smaller values of $h$. This is quite intriguing because these two behaviors (social conformity and homophily) have very similar effects at an individual level (i.e., they both make ego and alter similar to each other). In fact, their differences are often very vague and undetectable in empirical social network studies \cite{shalizi2011}. Meanwhile, these two behaviors are mechanistically distinct, because social conformity is about node dynamics while homophily is about edge dynamics. This finding, that their competitive balance may lead to very different societal outcomes down the road, offers a lot of implications for how we should consider our social interactions and behaviors in this highly interconnected world.

\section{Conclusions}

In this study, we conducted large-scale numerical simulations of our adaptive social network model to investigate the transition points between fragmentation and homogenization of social networks in a multidimensional behavioral parameter space. Artificial neural network-based modeling and visualization of the parameter-outcome mapping revealed a new phase transition plane for strong social conformity ($c$) and weak attention to novelty ($a$) regimes, which was previously unrecognized. The overall multidimensional phase space structure shows a nonlinear interaction among the three key behavioral mechanisms (social conformity, homophily, and attention to novelty). Within the range of parameter values tested so far, it appears that social homogenization (blue regions in Figs.~\ref{phase-spaces-structure} and \ref{phase-spaces-state}) occupied a greater volume in the log-scale parameter space than social fragmentation did.

This study presents a concrete example of complex collective systems research to study {\em Artificial Society}, i.e., study of hypothetical models of {\em society-as-it-could-be}. Such theoretical/mathematical/computational exploration of social systems can play valuable roles complementary to more empirical social science research, in the same spirit of Artificial Life research \cite{langton1992} that complements traditional biology. Computational examination of hypothetical scenarios, such as changing individual behaviors in our model, allows for exploration of various possible forms of our society and may lead to a discovery of novel possible social states which would not be realized just by analyzing empirical data obtained from real society \cite{sayama2020b}. Such exploratory endeavor is becoming increasingly important and relevant in today's highly automated, interconnected society, as our daily interactions are moving away from traditional, ``natural'' forms and becoming more and more mediated by artificially designed, ``engineered'' communication platforms. This has become even more manifested because of the recent COVID-19 pandemic (think about Zoom, YouTube, Slack, and other social media/collaboration platforms). We hope that studies like ours presented here may help re-evaluate and re-design the algorithms and interfaces of online human communications and interactions for the betterment of our social network evolution.

This study is still limited in several aspects. First, we did not explore variations of the amplitude of stochastic fluctuations ($\epsilon$) or functional shapes of homophily and attention to novelty ($F_h$ and $F_a$). Second, phase transition planes were identified only by numerical simulations while analytical estimate of transition conditions is not accomplished yet. Third, we assumed that the behavioral parameter values would apply uniformly to all individuals in society with zero behavioral diversity. Future research directions are naturally to address each and all of these limitations in the current model. In particular, introducing individual behavioral diversity within a collective complex system is known to produce unexpected, nontrivial macroscopic outcomes \cite{sayama2020b,sayama2009}. Such behavioral heterogeneity should be represented in future models to gain more nuanced, more realistic collective outcomes. Finally, quantitative comparison and validation of model behaviors with actual social network evolution data will ultimately be needed. However, obtaining such empirical data of social network evolution has been extremely difficult, and this will remain one of the major challenges in adaptive social network modeling research.