Surface reconstruction of tetragonal methylammonium lead triiodide

We present a detailed first-principles analysis of the (001) surface of methylammonium lead triiodide (MAPbI3). With density-functional theory we investigate the atomic and electronic structure of the tetragonal (I4cm) phase of MAPbI3. We analysed surfaces models with MAI- (MAI-T) and PbI2-terminations(PbI2-T). For both terminations, we studied the clean-surface and a series of surface reconstructions. We find that the clean MAI-T model is more stable than its PbI2-T counterpart. For the MAI termination,reconstructions with added or removed units of nonpolar MAI and PbI2 are most stable. The corresponding band structures reveal surface states originating from the conduction band. Despite the presence of such additional surface states, our stable reconstructed surface models do not introduce new states within the band gap.


I. INTRODUCTION
Perovskite solar cells (PSCs) have attracted immense attention within the photovoltaic community due to their rapidly rising power conversion efficiency (PCE): it reached 25.5% 1 only nine years after the invention of the state-of-the-art PSC architecture in 2012 (PCE ∼10%) 2,3 . The hybrid (organic-inorganic) halide perovskite (HP) methylammonium (MA) lead triiodide (CH 3 NH 3 PbI 3 or MAPbI 3 ) has been the most common PSC photoabsorber for a long time, and it is still a major focus of both experimental and theoretical studies, along with the rising isostructural material based on formamidinium (FA). HPs have also received significant recognition in luminescence and light detection [4][5][6][7][8][9] .
The application of these proposed solutions requires a) Electronic mail: azimatu.seidu@aalto.fi an understanding of the surface properties and possible surface reconstructions of HPs. This includes several aspects, such as morphology control during the growth of the HP thin films, HP-interlayer interface engineering, and the passivation of intrinsic defects at the interfaces and grain boundaries. A comprehensive understanding of the atomic and electronic structure of MAPbI 3 surfaces would advance the development of this class of novel materials and their applications. The surfaces of MA- [44][45][46][47] and FA-based 48-50 perovskites have been investigated theoretically and experimentally 4,51-58 , but understanding of the non-pristine surfaces is still lacking. Most of the MA-perovskite surface studies are focused on the stability of the two main terminations of HP (001) surfaces: MAI-and PbI 2 -terminated (shortened as MAI-T and PbI 2 -T hereafter, respectively) with little to no consideration of possible surface reconstructions. In our previous work, we investigated the atomic and electronic structure of (001) surfaces of cesium lead triiodide (CsPbI 3 ) using first-principles density functional theory (DFT) calculations and surface-phase-diagram (SPD) analysis 59 . For both cubic (α) and orthorhombic (γ) phases, we found that the CsI-termination is more stable than PbI 2 -termination, and the former class features a series of stable surface reconstructions with added or removed valence-neutral CsI and PbI 2 units. Our previous study established a systematic method for understanding stable surface reconstructions with a representative HP and motivates the present work in regard to both materials engineering and theoretical methodology.
In this work, we present a comprehensive DFT study of the (001) surface of the room-temperature tetragonal phase of the more popular HP, MAPbI 3 . Haruyama et al. have carried out preliminary studies for this surface and identified some stable surface terminations dependent on growth conditions, with some of them beyond the regular "clean surface" models 44,45 . Nevertheless, an extensive exploration of surface terminations and reconstructions with the addition or removal of constituent elements CH 3 NH 2 (MeNH 2 ), Pb, I, and their complexes is lacking.
We aim to establish such a systematic theoretical description by means of DFT, ab initio thermodynamics 60-62 , and SPD analysis.
It is worth noting that the unique charge state of the organic MA cation introduces additional complexity into this DFT study compared to our work on CsPbI 3 . Simply separating MAPbI 3 into its constituents MA, Pb, and I in a way similar to the decomposition of CsPbI 3 into Cs, Pb, and I 2 is not thermodynamically sensible. The charge-neutral CH 3 NH 3 radical is not stable on its own and far less suitable as a thermodynamic reference system than Cs is for CsPbI 3 . We therefore use the neutral CH 3 NH 2 (MeNH 2 ) and H 2 molecules in this work. Similar to Ref. 59, we will classify the thermodynamic stability of considered MAPbI 3 surfaces for different growth conditions and analyze their electronic structure.
The remainder of this paper is organized as follows. In Sec. II, we briefly outline the computational details of our DFT calculations and summarize the thermodynamic constraints for the growth of bulk MAPbI 3 , as well as the MAI-T and PbI 2 -T surfaces. In Sec. III, we first analyze the stability of the clean-surface models (MAI-T and Pb 2 -T) and the reconstructed models with missingand add-atoms and complexes. We then discuss the impact of surface reconstruction on both the atomic and electronic structure. Finally, we conclude with a summary in Sec. IV.

II. COMPUTATIONAL DETAILS
All DFT calculations were performed using the Perdew-Burke-Ernzerhof exchange-correlation functional for solids (PBEsol) 63 implemented in the all-electron numeric-atom-centered orbital code fhi-aims [64][65][66] . We chose PBEsol because it describes the lattice constants of MAPbI 3 well at moderate computational cost 67,68 . In our previous study on CsPbI 3 surfaces 59 , we also tested the PBE functional, but found only negligible changes in the surface phase diagram. We expect the same to be true for MAPbI 3 . Scalar relativistic effects were included by means of the zeroth-order regular approximation 69 . As with the PBE test, the inclusion of full spin-orbit coupling did not affect the conclusions of our CsPbI 3 study 59 and we expect the same for MAPbI 3 . Standard fhi-aims tier-2 basis sets were used in combination with Γ-centered 4×4×4 (bulk) and 4×4×1 (surfaces) k-point meshes. The bulk structures were optimized with the analytical stress tensor 70 until forces were below 5.0 × 10 −3 eV · Å −1 . For the surface slab models, we fixed the lattice constants and all atomic positions except for atoms in the top and bottom MAPbI 3 units (the surface atoms). A surface-dipole correction 71 was applied in all surface calculations.
In the interest of open science 72 , we made all relevant calculations included in this work available on the Novel Materials Discovery (NOMAD) repository 73 .

Bulk and surface structures
As experimentally reported, the tetragonal phase of MAPbI 3 is stable from ∼ 160 to ∼ 330 K, including room temperature 74,75 . The structure belongs to the polar space group I4cm (No. 108) as a result of its intrinsic polarization along the principal axis 74 . Considering several possible disordered MA alignments 76,77 , we constructed a series of 2 × 2 × 2 supercells with different MA orientations and optimized their structures with DFT. We then take the structure with the lowest energy. The lattice parameters of this structure (Fig. 1)  In this work, we focus on the (001) surfaces, which are the major facet of HPs 44,45,78 and the most relevant surfaces of MAPbI 3 . Due to the polar bulk structure, it is not possible to build a surface supercell by repeating several bulk layers along the [001] direction as this would result in a polar surface model (see Figure 2a left). Such a model will induce artefacts into the calculated properties of the system such as the unphysical removal of band degeneracies and reduction of the band gap (Fig. 2a right) and would ultimately lead to a polar catastrophe, in which the valence band at one end of the slab lie higher in energy than the conduction bands at the other end.
To circumvent these artefacts, we constructed a symmetric slab model by introducing a "domain wall" in the slab. As sketched in Figure 2b left, such a domain wall is a PbI 2 -containing (001) plane located at the center of the slab. The atomic structures on opposite sides of the domain wall is mirrored with respect to this plane, so that the [001] components of the MA dipole moments on opposite sides cancel each other, giving rise to a nearly vanishing overall dipole moment. As a result, the polar artefacts vanish and the surface band structure (Fig. 2b right) exhibits a proper band gap and the right degeneracies. Similar approaches have been successfully employed in previous studies 44,79 . With the approach illustrated in Fig. 2b, we constructed symmetric clean surface models in a way similar to our previous work for CsPbI 3 . Specifically, the MAI-T surface model consists of 6 MAI and 5 PbI 2 layers alternately stacked along the [001] direction. Similarly, the PbI 2 -T surface model has 7 PbI 2 and 6 MAI alternating layers. By inserting a 40 Å-thick vacuum layer to separate neighboring slabs along [001] and including surface-dipole correction 71 in the DFT calculations, we minimized the interaction between neighboring slabs. Figure 3 depicts the optimized structures of both clean MAI-T and PbI 2 -T surfaces. The top views of both phases show a similar in-plane tilting pattern of PbI 6 octahedra and in-plane alignment of MA dipoles as in the bulk. The side views demonstrate that the mirror symmetry of both slab models with respect to the domain wall is maintained after geometry optimization. We note that in MAI-T, MA dipoles at both top and bottom surfaces point inwards (sketched in Fig. 2) as a result of hydrogen bonding with the surface I ions.
We studied various add-and missing-atom surface models based on both MAI-T and PbI 2 -T clean surfaces. All add-atom models (i X ) were constructed by adding the atoms or atom-complexes X to the surface, while for missing-atom models (v X ), atoms or complexes X were removed from the topmost X-containing layers. For MAI-T surfaces as an example, v MeNH 2 , v H , v MA , v I , and v MAI were constructed by removing atoms from the topmost MAI layer, while v Pb and v PbI 2 indicate the removal of atoms from the PbI 2 layer below the topmost MAI layer. For models with double missing-or addatoms (i.e., v 2X or i 2X ), we considered both line and diag- Only the more stable model will be presented and discussed in Sec. III. For instance, we find the line modes to be more stable in both v 2MAI and i 2PbI 2 .

B. Thermodynamic constraints for stable MAPbI 3 bulk and surfaces
We applied the grand potential analysis to investigate the stability of a variety of different surface reconstructions. Neglecting finite temperature contributions, the grand potential (Ω) is Here, ∆H indicates the standard formation energy of the model system, E the total energy, µ i the chemical potential of species i in its most stable form, x i the number of atoms of this species in the structure, and ∆µ i the change in the chemical potential away from its value in the element's most stable phase, µ i . ∆µ i represents the control of experimental growth conditions and is both a meaningful and convenient parameter to vary in phase diagrams. The relative stability between two structures is determined by comparing their grand potentials, with the structure lower in grand potential considered more stable. Details of the grand potential analysis are described in our previous work on surface reconstruction of CsPbI 3 59 . We first consider conditions for stable MAPbI 3 in the bulk. In order to avoid the formation of elemental Pb and I, molecular MA (as a whole instead of elemental C, N, and H for simplicity), as well as bulk MAI and PbI 2 , the region of the phase diagram for stable MAPbI 3 is determined by the inequalities, The inequalities in Eq. (3) can be rearranged as Inequalities in Eq. (2) define the domains of variables µ MeNH 2 , µ H , µ Pb , and µ I , and the inequalities in Eq. (4) define the region for growth of "stable-bulk MAPbI 3 " in the phase diagram. µ MeNH 2 , µ H , µ Pb , and µ I can be calculated for the stable reference structures of MeNH 2 (molecule), H (H 2 molecule), Pb (P 6 3 /mmc), and I (I 2 molecule) with DFT, respectively. Formation energies ∆H in Eq. (4) can be calculated with DFT, too.
Equations (2) and (4) only serve to determine the bulk stability. For the stability of (clean and reconstructed) surface models, we need to solve Eq. (1) to obtain the SPDs. Note that the bulk and surface are not in isolation from each other. The final surface stability is determined by the intersection of the SPD and the stable-bulk region.

A. Thermodynamic stability analysis of bulk and surface terminations
The PBEsol-calculated formation energies of bulk The energy required for tetragonal MAPbI 3 to decompose into MAI and PbI 2 , i.e., the difference between the left and the right values of either inequality, is as small as 0.06 eV. Such a narrow stability region reflects the general instability of tetragonal MAPbI 3 . SPD analysis helps identify the stability of the two considered surface terminations. Figure 4 shows that, at ∆µ MeNH 2 = ∆µ H = 0, the MAI-T and PbI 2 -T clean surfaces are stable in the Pb-poor and Pb-rich limits, respectively. We consider MAI-T the more stable surface since the region for stable MAI-T covers a wider range of  that claimed the stability of MAI-T over PbI 2 -T and is similar to the CsPbI 3 surface properties that we reported earlier 59 . Our discussions will therefore focus on MAI-T surfaces from here on. Data for PbI 2 -T surfaces, including the relaxed surface-reconstructions and the SPDs, are given in the Supplementary Material (SM) .

B. Identification of stable reconstructions of MAI-T surfaces
SPDs for the considered surface reconstructions of the MAI-T surfaces are shown in Figure 5 (SPDs for the PbI 2 -T counterparts are available in Fig. S1 of SM). It is not surprising that the ∆µ I /∆µ MeNH2 and ∆µ I /∆µ H SPDs display similar features, as MeNH 2 and H are closely related to each other through the organic MA component of the material. In these two SPDs, which are given in the Pb-rich limit (∆µ Pb = 0), we observe the following stable surface structures: i 4MAI (in the MeNH 2and H-rich limit), v 4MAI (in the MeNH 2 -and H-poor limit), clean surface, v PbI 2 , v 2PbI 2 , i PbI 2 , i 2PbI 2 , i 2MAI , and v MAI . The major difference in the appearance of these two phase diagrams lies with the v 4H surface, which is observed in the H-poor and I-rich limit. With our choice of 2D slices, this surface reconstruction appears in one quadrant of only one of these two 2D phase diagrams.
The MeNH 2 -and H-rich (thus MA-rich) limit (∆µ MeNH 2 = ∆µ H = 0) creates a third 2D slice of the total phase diagram, shown on the right side of Figure 5. In this SPD, we find i 4MAI and v 2PbI 2 to be stable. Except for v 4H , all the observed stable reconstructions are valence-neutral, i.e., with addition or removal of MAI or PbI 2 units, net charges are not induced in the system, which is similar to what we previously found for the CsPbI 3 surfaces 59 . We notice that in the MeNH 2 -, H-, and Pb-rich limit i 4MAI dominates over PbI 2 -derived reconstructions. That is, on the MAI termination layer at the MAPbI 3 surface, the tendency for growing an extra MAI layer is greater than for growth of PbI 2 units, which would eventually transform the system into PbI 2 -T. This finding again verifies that MAI-T is more stable.
We are particularly interested in the most relevant reconstructions, which we define as those regions in the SPDs that intersect the stable bulk region. It is these overlapping regions of bulk and surface stability that are viable standalone surfaces in the laboratory. These relevant models are the clean surface, v PbI 2 , v 2PbI 2 , i PbI 2 , and i 2PbI 2 in the Pb-rich limit, and i 4MAI in the MA-rich limit. Different from CsPbI 3 , for which we observe a relatively broad range in chemical potential for the clean surface, the range for its stability on MAI-T (001) at ∆µ Pb = 0 is very narrow in terms of ∆µ MeNH 2 and ∆µ H . In addition, it is only stable in I-deficient growth conditions.
The optimized geometries of the most relevant surface models are given in Figure 6 (relevant reconstruction models of PbI 2 -T are presented in Fig. S2 of SM). Because the surface atomic structure varies mainly to accommodate the absence or addition of atoms, our discussion of geometric rearrangement will be with reference to the clean surface in the following. We observe that PbI 2 removal causes noticeable atomic structure changes in the reconstructed surfaces. The topmost PbI 2 layer of v PbI 2 displays (PbI 6 ) 2 (PbI 5 ) polyhedra, while in v 2PbI 2 there are two isolated PbI 5 polyhedra (Fig. 7). Interestingly, no migration of surface I anions occurs in v 2PbI 2 , which is the main characteristic of the equivalent removal on α-CsPbI 3 59 . This is very likely due to the different A-site cations: the hydrogen bonding between the ammonium group and the surface I anions would stabilize the latter, so that the surface Pb-I units are relatively regularly distributed.
In i PbI 2 and i 2PbI 2 , each added PbI 2 unit is linked to a surface I atom via Pb-I bonding, giving rise to a PbI 3 tetrahedron that contains one Pb and three I atoms as its vertices (Figure 7). Notably, the i 2PbI 2 reconstruction shows characteristic PbI 5 PbI 3 polyhedra (Fig. 7), as previously reported by Haruyama et al. 44 . This asymmetric distribution of surface I atoms results from the removal of two linearly-aligned PbI 2 units, which is very different to the same surface reconstruction of α-CsPbI 3 where the diagonal mode is more stable.
Finally, we find a relatively regular alignment of the added MAI units, forming a uniform sheet on the MAI-T surface in i 4MAI . The C-N bonds of all added MA + cations point towards the surface to form hydrogen bonds with the topmost I anions. The average shortest H(N) · · · I distance is 2.68 Å, a typical value for hydrogen bonding 84,85 . The fact that extra MAI units can readily grow above the already existing MAI surface-termination layer, as also indicated by the SPDs in surfaces in MAI-rich situations. This would be detrimental to device performance since MAI is very poor in transporting charge carriers.

C. Electronic properties of bulk and most relevant reconstructed MAPbI3 surface models
In this section, we focus on the electronic properties of bulk MAPbI 3 and the most relevant reconstructed MAI-T surface models. The band structures of the most relevant PbI 2 -T surface reconstructions are presented in Fig. S3 of SM. Figure 8 depicts the band structures of the bulk and the pristine MAI-T surface of MAPbI 3 . For the bulk, we adopt the high-symmetry k-point path of a simple-cubic lattice in the 2 × 2 × 2 supercell model for simplicity. Our plots show the band structure along M-X-Γ-M with M = 1 2 , 1 2 , 0 , X = 0, 1 2 , 0 , and Γ = (0, 0, 0), i.e., within the a * b * plane of the Brillouin zone (identical to the ab = (001) plane in real space in our cases).

Electronic properties of the bulk and the clean MAI-T surface
The bulk band structure exhibits a direct band gap at the Γ point. The element projected density of states (PDOS) in Fig. 9 reveals that the valence band (VB) is dominated by I-5p orbitals. The VB maximum (VBM) exhibits a noticeable contribution from Pb-6s which gives rise to the well-known antibonding character and thus introduces the noticeable band dispersion at Γ (Fig. 8). The conduction-band minimum (CBM) consists mainly of Pb-6p orbitals. The MA cation shows no significant contributions at the band edges.
Next, we investigate, if our MAI-T surface model introduces surface or mid-gap states. The bands of the clean MAI-T surface are shown in the middle panel of Fig. 8. Both the VB and CB edges of the surface nearly align with the bulk bands at the M point. The band gap of the surface at M is slightly larger than the bulk, which could be a quantum confinement artefact of the slab model, but the bands themselves agree. At Γ and X, however, the CB of the clean surface extends below the CBM of the bulk. At these points, the shapes of the bulk and surface bands at Γ are different.
To understand the nature of these states, we first analyse the PDOS in Fig. 9. The PDOS verifies that the apparent band gap of the clean surface is ∼ 0.2 eV less than the bulk, as already indicated by the band structure. We plot the charge densities of the lowest four CB states of our surface model in Fig. 10. We find that these states are surface states that come in nearly degenerate pairs, with the partners of each pair on opposite sides of the slab. The slight degeneracy lift in each pair is caused by a small relaxation induced structural asymmetry in the two bulk halves that make up our surface slab model (see Section II A 1). The characters of these band edge wave functions are the same as the ones we observed for the α phase of CsPbI 3 59 . The DOS of the clean surface in Fig. 9 is consistent with the interpretation of conduction band derived surface states. The right panel exhibits small bumps as shown by the arrow in Fig. 9 (right panel) in the CB at approximately 1.25 and 2.10 eV compared to the bulk, indicating the rearrangement of bands in the slab model.
The behaviour of the VBM is a little different. It is pinned to the domain wall at the mid-plane of the slab model. This quasi-two-dimensional state is still dispersive and its band very closely matches that of the bulk, which provides good evidence for the quality of our slab model. The pinning of the state to the mid-plane is reasonable. The dipoles on either side of the domain wall point away from the mid-plane, raising the energy of an electron residing on the mid-plane and positioning it as the maximum of the occupied states (VBM). Again, the lack of artefacts in the slab electronic structure and good agreement with the bulk support the quality of the model, even with the pinning of the VBM to the domain wall.

Electronic properties of the most relevant MAI-T surface models
The band structure and PDOS of the most relevant surface models observed in Fig. 5, (v PbI 2 , v 2PbI 2 , i PbI 2 , i 2PbI 2 , and i 4MAI ), are shown in Figures 11 and 12, respectively. Similar to Fig. 8, the bulk band structure is included as background for comparison in Fig. 11. We find that the band structures of two of the most relevant reconstructed surface models, v PbI 2 and i 4MAI , resemble the band structure of the clean MAI-T surface shown in Fig. 8.
For the others, there are flat bands near and below the bulk CB edges (CBEs), which are more pronounced at the M-point. They are most visible in i PbI 2 . In the i PbI 2 and i 2PbI 2 surfaces, we observe increased intensities within the CB of the PDOS at ∼ 1.6 and ∼ 1.9 eV, which corresponds to the flat bands at M-point in Fig. 11. These small peaks come from Pb states, suggesting that the flat bands are indeed due to the added PbI 2 units. To confirm this, we plot the charge distribution of CBE at M for these two reconstructions in Figure 13.
Even though the CBE at M in these two reconstructions and in v 2PbI 2 belongs to a band which is flat across the entire Brillouin zone, the wave functions of these states (Fig. 13) still resemble the surface states of the pristine slab, especially for v 2PbI 2 and i 2PbI 2 . This is somewhat surprising since the CBE states of the cleansurface model belong to a dispersive band. The CBE states of the reconstructed surfaces in Fig. 13 still come in nearly degenerate pairs with each partner appearing on opposite sides of the slab, as for the clean surface. However, we can ignore the state at the bottom of the slab, since it corresponds to the clean and not the reconstructed surface. The states in i PbI 2 and i 2PbI 2 at the top of the slab, on the reconstructed surface side, have considerable weight on the added PbI 2 units. This wave-function localization explains the flat character of the band.

IV. CONCLUSION
In summary, we have investigated the stability and electronic structure of MAPbI 3 surfaces in the tetragonal phase from first principles. To circumvent the polar catastrophe in our supercell calculations, we build a slab geometry from two MAPbI 3 bulk segments with opposite polarity, effectively introducing a domain wall in the middle of the slab. Our surface science study reveals that the methylammonium-iodine (MAI) termination is more stable than the PbI 2 termination. We further observe that the removal or addition of polar units that induce zero net charge in the system lead to more stable sur-face reconstructions. MAI-terminated surfaces introduce conduction-band derived surface states near the conduction band edge, which result in a surface band gap that is slightly smaller than the bulk band gap. The stable reconstructions do not introduce further surface states in the band gap, which bodes well for the transport properties across interfaces with these reconstructions. Our study opens up future work on surface adsorbates, defects and interfaces.

V. SUPPLEMENTARY MATERIAL
See Supplementary Material for surface phase diagrams of PbI 2 -terminated models, crystal and electronic band structures of the most relevant reconstructed surfaces of PbI 2 -T models .

VI. AUTHORS CONTRIBUTIONS
All authors contributed equally in this work.

ACKNOWLEDGMENTS
We acknowledge the computing resources from the CSC-IT Center for Science, the Aalto Science-IT project, and Xi'an Jiaotong University's HPC Platform. We further acknowledge funding from the Väisälä Foundation and the Academy of Finland through its Key Project Funding scheme (305632) and postdoctoral grant no. 316347.

VII. CONFLICT OF INTEREST
The authors have no conflicts to disclose.

VIII. AIP PUBLISHING DATA SHARING POLICY
The data that supports the findings of this study will be openly available in Novel Materials Discovery (NO-MAD) repository at 73