Perspective: Theory and simulation of hybrid halide perovskites

Organic-inorganic halide perovskites present a number of challenges for first-principles atomistic materials modeling. Such “plastic crystals” feature dynamic processes across multiple length and time scales. These include the following: (i) transport of slow ions and fast electrons; (ii) highly anharmonic lattice dynamics with short phonon lifetimes; (iii) local symmetry breaking of the average crystallographic space group; (iv) strong relativistic (spin-orbit coupling) effects on the electronic band structure; and (v) thermodynamic metastability and rapid chemical breakdown. These issues, which affect the operation of solar cells, are outlined in this perspective. We also discuss general guidelines for performing quantitative and predictive simulations of these materials, which are relevant to metal-organic frameworks and other hybrid semiconducting, dielectric and ferroelectric compounds.


I. INTRODUCTION
The perovskite mineral, CaTiO 3 , is the archetype for the structure of many functional materials. 1 Metal halide perovskites have been studied for their semiconducting properties since the 1950s. 2 Only recently organic-inorganic perovskites, such as CH 3 NH 3 PbI 3 (MAPI), have been applied to solar energy conversion, showing remarkably strong photovoltaic action for a solution processed material. 3 The field has progressed rapidly in the last five years. The increase in the power conversion efficiency is supported by over three thousand research publications. [4][5][6][7][8] Other potential application areas of these materials include thermoelectrics, 9,10 light-emitting diodes, 4,11 and solid-state memory. 12,13 Recently we published a short review on the nature of chemical bonding in these materials 14 and on the multiple time scales of motion. 15 We will not repeat that material here. There is also a recent review from Mattoni and co-workers focusing upon the use of molecular dynamics (MD) simulations. 16 In this perspective, we address recent progress and current challenges in theory and simulation of hybrid halide perovskites. We pay particular attention to predicting properties that assess the photovoltaic potential of a material. Factors to consider include light absorption, charge transport, absolute band energies, defect physics, and chemical stability. The total energy, electronic energy levels, dielectric function, and band effective masses can be calculated with electronic structure methods on a representative (static) crystal structure. The lattice and molecular dynamics can describe a variety of dynamic behaviors at finite temperatures. These perovskites combine a complex crystal structure, modulated by the static a) Electronic mail: a.walsh@imperial.ac.uk and dynamic disorder, with a subtle electronic structure requiring methods beyond density functional theory (DFT) to correctly treat the many-body and relativistic effects. As such, the halide perovskites represent a challenge to predictive materials modeling, in a system of great experimental interest, and where there is considerable motivation to improve on the status quo.

II. CRYSTAL STRUCTURES AND LATTICE DYNAMICS
A. Phase diversity (Hybrid) perovskites of the type ABX 3 form a crystal structure with an (organic) A site cation contained within an inorganic framework BX 3 of the corner-sharing octahedra. Halide substitution on the X site (X = Cl − , Br − , I − ), metal substitution on the B site (B = Pb 2+ , Sn 2+ ), and cation substitution on the A site (A = CH 3 NH + 3 , HC(NH 2 ) + 2 , Cs + , Rb + ) lead to varied chemical and physical properties. 19,20 In addition to isoelectronic substitution (e.g., replacing Pb 2+ by Sn 2+ ), it is possible to perform pairwise substitution to form double perovskites (e.g., replacing two Pb 2+ by Bi 3+ and Ag + ). 21,22 In the first report of CH 3 NH 3 PbI 3 by Weber in 1978, the crystal structure was assigned as a cubic perovskite (space group Pm3m). 23,24 The anionic PbI − 3 network is charge balanced by the CH 3 NH + 3 molecular cation. The symmetry of CH 3 NH + 3 (C 3v ) is incompatible with the space group symmetry (O h ) unless the orientation disorder (static or dynamic) is present. The crystal structure solved from X-ray or neutron diffraction data usually spread the molecules over a number of orientations with partial occupancy of the associated lattice sites. A common feature of perovskites is the existence of phase changes during heating (typically from FIG. 1. The high-resolution powder neutron diffraction pattern of the hybrid halide perovskite CH 3 NH 3 PbI 3 is shown in the left panel (adapted with permission from Ref. 15 based on the data in Ref. 17). This illustrates the low and high temperature phase transitions. While an ordered CH 3 NH + 3 sub-lattice is expected in the orthorhombic phase, the orientational disorder increases with higher temperature. The crystallographic unit cells of the pseudo-cubic and orthorhombic perovskite phases are shown in the right panel (adapted with permission from Ref. 18). The associated structure files can be accessed from https://github.com/WMDgroup/hybrid-perovskites.
lower to higher symmetry) as shown in Fig. 1. In hybrid halide perovskites containing methylammonium, these are orthorhombic (Pnma), tetragonal (I4/mcm), and cubic (Pm3m) phases. 17 For CH 3 NH 3 PbI 3 , the Pnma to I4/mcm phase transition is first-order with an associated discontinuity in physical properties, while the I4/mcm to Pm3m phase transition is second-order with a continuous evolution of the structure and properties. 17,25 The phase transitions are linked to a change in the tilting pattern of the inorganic octahedral cages and the orderdisorder transitions of the molecular sub-lattice. 25-27 X-ray diffraction (XRD) measurements upon cooling (heating) suggest the inclusion of tetragonal in orthorhombic phases (and vice-versa). 28 This is often observed for the first-order solidstate phase transitions. In addition, it has been suggested that the presence of multiple photoluminescence peaks at low T is due to the coexistence of ordered and disordered orthorhombic domains. 29 A similar phase behavior tends to be seen for other compositions, however the transition temperatures vary. In CH 3 NH 3 PbI 3 , the orthorhombic to tetragonal transition temperature is 162 K, becoming cubic by around 328 K. CH 3 NH 3 PbBr 3 is cubic above 237 K. 30 In addition, compounds such as HC(NH 2 ) 2 PbI 3 (FAPI) and CsSnI 3 feature the phase competition between a corner-sharing octahedra perovskite phase (black in appearance) and the edge-sharing octahedra molecular crystals (yellow or white in appearance). 31 Only the corner-sharing perovskite phase is of interest for solar energy applications.

B. Local and average crystal environment
The first electronic structure calculation of hybrid halide perovskites was reported by Chang, Park, and Matsuishi in 2004, 34 in the local density approximation (LDA) of density functional theory (DFT). They modeled a static structure where the CH 3 NH + 3 molecule was aligned along 100 (toward the face of the corner-sharing PbI − 3 framework), but found that the barrier for rotation to 111 was less than 10 meV. This small barrier for cation rotation gave credence to a prior model that the molecular sub-lattice was dynamically disordered. 30 Similar barriers were later found within the generalised gradient approximation (GGA) of DFT. 35 Ab initio molecular dynamics (MD), neutron scattering, 36,37 and time-resolved infra-red 38 data all indicate a 1-10 ps reorientation process of the molecular cation at room temperature. As a result of (by definition) anharmonic molecular rotation and large-scale dynamic distortions along soft vibrational modes, the local structure can deviate considerably from that sampled by diffraction techniques. Bragg scattering does not probe the local disorder if it preserves the longrange order on average. Knowledge of these locally broken symmetries is essential for meaningful electronic structure calculations, where the broken symmetry results in the lifting of degeneracy, and a potentially quite different solution.
In spite of the larger cation, FAPI appears to possess a similar time scale of rotation to MAPI. 31 A lighter halide (and therefore a smaller cage) results in faster rotation, in spite of the greater steric hindrance. 39 Together, these data suggest that molecular rotation is a function of the local inorganic cage tilting. The relatively insignificant mass of the organic cation follows the distortion of the cavity.
Spontaneous distortions can also be observed in the vibrational spectra. The calculated harmonic phonon dispersion for MAPI in the cubic phase is presented in Fig. 2. The acoustic phonon modes soften as they approach the M (q = 1 2 , 1 2 , 0) and R (q = 1 2 , 1 2 , 1 2 ) Brillouin zone boundary points. This zone boundary instability can only be realised in an even supercell expansion, where it corresponds to anti-phase tilting between successive unit cells. This behavior is characteristic of the perovskite structure and can be described by the Glazer tilt notation. 40,41 Within the frozen-phonon approximation, the potential energy surface can be traced along the soft acoustic M and R phonon modes. In both cases, this results in a double well with an energy barrier ∼k B T at the saddle point. 32 At room temperature, the structure is dynamically disordered, with continuous tilting. The structure is locally non-cubic but possesses only cubic Bragg scattering peaks. 42 Indeed, the MD simulations of halide perovskites show continuous tilting of the octahedra at room temperature. 31,43,44 As the temperature decreases, the structural instability condenses via the soft mode at the R point (with an energy barrier of 37 meV) into the lower symmetry J. Chem. Phys. 146, 220901 (2017)

FIG. 2. (Left)
The harmonic phonon dispersion for CH 3 NH 3 PbI 3 from a "pseudo-cubic" structure. The imaginary frequencies of acoustic modes at the M (q = 1 2 , 1 2 , 0) and R (q = 1 2 , 1 2 , 1 2 ) Brillouin zone boundary points correspond to an instability expressible in a supercell as alternate tilting of the octahedra. (Right) The imaginary acoustic mode at the R point in a 2 × 2 × 2 supercell expansion shows a double-well potential in the DFT internal energy. The saddle point corresponds to a 1 × 1 × 1 cubic structure, while the two local minima correspond to a distorted structure of lower symmetry. The energy barrier is small enough to allow both minima that can be accessed at room temperature, so the system is expected to exhibit dynamic disorder rather than static disorder. A similar behavior is found at the M point. The figure is adapted with permission from Refs. 32 and 33. The underlying phonon data are available from https://github.com/WMD-group/Phonons. tetragonal phase. This is followed by condensation of the M point (with an energy barrier of 19 meV) to the orthorhombic phase. 32 While the molecular cation continuously rotates with the inorganic tilts in the cubic phase, and is partly hindered in the tetragonal phase, it can only liberate in the low temperature orthorhombic phase.
In the static picture (as in the case of an electronic band structure calculated for a single ionic snapshot), the organic cation plays no direct role in the optoelectronic properties of the material as the molecular electronic levels lie below that of the inorganic framework. Allowing motion, the electrostatic and steric interaction between the organic molecule and inorganic framework couples tilting and distortion of the octahedra to the organic cation motion. The dynamic structural distortions change the bond lengths, angles, and the atomic orbital overlap, perturbing the band-structure and bandgap. 32,[44][45][46] The electronic structure thus becomes sensitive to the temperature, which will be discussed further in Sec. III.

C. Thermodynamic and kinetic stability
Ab initio thermodynamics has emerged as a powerful tool in materials modeling, with the ability to assess the stability of new materials and place them on equilibrium phase diagrams even before the experimental data are available. [47][48][49] The total energy from the DFT calculations approximates the internal energy of the system. By including lattice vibration (phonon) and thermal expansion contributions, the Gibbs free energy and other thermodynamic derivatives can be evaluated. 50 In context of photovoltaic materials, this has been applied to Cu 2 ZnSnS 4 and used to identify the processing window where a single-phase compound can be grown in equilibrium. 51 For the tin sulfide system, it shows the close competition between SnS, SnS 2 , and Sn 2 S 3 . 52 An issue with hybrid perovskites and other metal-organic frameworks is that the calculated heat of formation is close to zero. The decomposition reaction CH 3 NH 3 PbI 3 → CH 3 NH 3 I + PbI 2 (1) has been predicted to be exothermic. 53 Subsequent calorimetric experiments have supported the prediction that hybrid lead halide perovskites are metastable. 54 It is likely that these materials are only formed due to entropic (configurational, vibrational, and rotational) contributions to the free energy. The concept of metastable materials is attracting significant interest. [55][56][57][58] These are materials that do not appear on an equilibrium phase diagram but can be synthesised with a finite (useful) lifetime. For compounds that do not represent a thermodynamic ground-state, the chemical kinetics become critical, and formation and stability and can be particularly sensitive to local gradients in chemical potential (e.g., compositional, thermal, electronic). Although kinetic factors can be calculated with first-principles techniques, this is a more cumbersome and costly process than the equilibrium bulk thermodynamics, which requires only the total energies of local minimum structures. To our knowledge, there have been no rigorous attempts to model the kinetics of decomposition pathways for hybrid perovskites over complete chemical reactions.

D. Anharmonic lattice vibrations and thermal conductivity
Electronic structure theory is most often carried out in the Born-Oppenheimer approximation where the nuclei are static classical point charges. To consider thermal vibrations, expansion, or heat flow, the theoretical framework of the lattice dynamics can be used. 50 In the harmonic approximation, the small-perturbation lattice dynamics are fully specified by the second-order forceconstants of individual atoms. These are readily constructed into the so-called dynamical matrix. The eigenstates of this matrix are the normal modes of vibration with an associated frequency. The description of collective vibrational excitations in crystals can be simplified with second quantization to the creation and annihilation of phonon quasi-particles, specified by these normal modes. Thermal expansion coefficients, system anharmonicity (e.g., modal Grüneisen parameters), and the temperature-dependence of other properties can be calculated in the quasi-harmonic approximation (QHA). In this formalism, the lattice dynamics are harmonic at a given temperature; however, the cell volume is scaled by thermal expansion to account for the finite-temperature anharmonic effects.
The thermal expansion coefficient of MAPI in the cubic phase has been calculated with the QHA. The value is sensitive to the exchange-correlation functional used. For example, a value of 3.0 × 10 −5 /K is calculated with the Perdew-Burke-Ernzerhof (PBE) functional with the Tkatchenko-Scheffler dispersion corrections. 59 The PBEsol functional produces a value of 12.5 × 10 −5 /K. 18 These compare to finite temperature scattering measures of 1.9 × 10 −5 /K by X-ray 60 and 13.2 × 10 −5 /K by neutron diffraction. 17 Even taking the smallest value above, the expansion coefficient is one order of magnitude greater than silicon. 61 This highlights the strong deviation from the harmonic behavior in halide perovskites.
In the harmonic approximation (and similarly the QHA), the eigenmodes of the dynamical matrix are orthogonal and the resulting phonons are non-interacting. Consequently, the phonon lifetimes are infinite as the phonons do not scatter; the thermal conductivity is ill-defined. To calculate phononphonon scattering, and so its contribution to the finite thermal conductivity, the anharmonic lattice dynamics need to be considered. A computational route is to use the perturbative manybody expansion, e.g., as implemented in PHONO3PY, 62 which includes the third-order force constants. For CH 3 NH 3 PbI 3 , 41 544 force evaluations are required to evaluate the thirdorder force constants, compared to 72 for the second-order (harmonic) force constants. 32 Consequently, these calculations are vastly more computationally expensive. Using this approach, the phonon-phonon scattering rates are calculated to be three times larger in MAPI compared to the standard covalent semiconductors CdTe and GaAs. 32 The phonons barely exist for a full oscillation before they split or combine into another state. The resulting mean free paths are on the nanometer rather than more typical micrometer scale. The lattice thermal conductivity is extremely low, 0.05 Wm 1 K 1 at 300 K. 32 This combination of high electrical conductivity and low thermal conductivity makes these compounds potential thermoelectrics. 9,10 In highly anharmonic systems perturbatively including the third-order force constants may not be sufficient to describe the true dynamics. Yet, going further in the lattice dynamics formalism becomes prohibitive. In addition to the computational cost, it is not obvious whether the fundamental tenant of the lattice dynamics, of expanding in small displacements around a minimum structure, is correct for such soft and highly anharmonic materials. In contrast, MD treats anharmonic contributions to all orders, but as it stochastically explores the phase space, long integration times are required to sample rare events. Finite size effects also mean that only phonon modes commensurate with the supercell are sampled, so size convergence has to be carefully considered.
The classical interatomic potentials derived from the first-principles calculations have been developed for hybrid perovskites. [63][64][65] Such models are able to correctly reproduce crystal structures, along with mechanical and vibrational properties. The calculation of thermal conductivity from molecular dynamics simulated for MAPI predicted values of 0.3 to 0.8 Wm 1 K 1 at 300 K. 66,67 Although still ultra-low, these values are greater than the values calculated by perturbation theory. It should be noted that these values give upper limits for thermal conductivity as they refer to a defect-free isotopically pure bulk sample.

III. ELECTRONIC STRUCTURE
Despite the dynamic disorder just discussed, in many respects, halide perovskites display the characteristics of traditional inorganic semiconductors, with a well-defined electronic band structure and electron/hole dispersion relations. However, subtleties emerge upon closer examination, when the electronic structure is correctly modeled.

A. Many-body and relativistic effects
Perhaps surprisingly, the local and semi-local exchangecorrelation functionals provide a reasonable estimate for the bandgaps of these heavy metal halide materials. This is due to cancellation of errors. For Pb-based perovskites, the conduction band has mainly the Pb 6p character. Due to the large nuclear charge, the electronic kinetic energy requires a relativistic treatment and spin-orbit coupling (SOC) becomes significant. The first-order effect is a reduction in the bandgap by as much as 1 eV, 68 as the degenerate 6p orbitals are split and moved apart in energy. This is shown in Fig. 3 for the bromide compounds. The typical bandgap underestimation of exchange-correlation functionals formed within the local density approximation (LDA) or generalized gradient approximation (GGA) is offset by the absence of relativistic renormalisation.
SOC is not expected to have a large impact on the structural properties of the Pb-based compounds as the (empty) conduction band is mainly affected. By the Hellmann-Feyman theorem, the force on atoms depends only on the electron density, which is provided by the occupied orbitals. Accurate force-constants (as needed in both molecular and lattice dynamics) can be calculated without SOC considerations. 69 There have been a number of electronic structure calculations considering many-body interactions beyond DFT. 68,70,71 Quasi-particle self-consistent GW theory shows that the band dispersion (and so density of states, optical character, and effective mass) is considerably affected by both the GW electron correlation and spin-orbit coupling. 68 Some materials see only a rigid shift of the band structure (retaining DFT dispersion relations), 72,73 but this is not the case for hybrid perovskites. This point has not been fully appreciated, in part, because DFT codes are more widespread and convenient to generate data.
A consequence of SOC when combined with a local electric field is the Rashba-Dresselhaus effect, a splitting of bands in the momentum space. 74 This can be understood as an electromagnetic effect, where the magnetic moment (spin) of the electron interacts with a local electric field, to give rise to a force which displaces it in the momentum space. Up and down spins are displaced in opposite directions, and this displacement is a function (in both size and direction) of the local electric field, which will depend on the local dynamic order. For a static structure, this is demonstrated in Fig. 3 for CH 3 NH 3 PbBr 3 . Neglecting SOC, the cubic phase has band extrema at the R point (a direct bandgap). With SOC, the valence and conduction bands each split into valleys symmetrical around R. The splitting is much more pronounced in the Pb 6p conduction band (compared to the Br 4p valence band), as expected from the Z 2 dependence of spin-orbit coupling. This asymmetry in the band extrema results in direct-gap like absorption and indirect-gap like radiative recombination, which we discuss later. The relativistic spin-splitting can only occur in crystals that lack a centre of inversion symmetry, a prerequisite for generating a local electric field. The cubic representation of CsPbBr 3 has an inversion centre, so, while SOC affects the bandgap through the separation of Pb 6p into p 1 2 and p 3 2 combinations, no splitting of the band extrema away from the high symmetry points is observed (see Fig. 3). This is true only for a static cubic structure. As discussed earlier, hybrid halides will have continuous local symmetry breaking. The calculations based on static high symmetry structures are not representative of the real (dynamic) system and can be misleading.
The calculation of electronic and optical levels associated with intrinsic and extrinsic point defects will be particularly sensitive to the electronic structure method used. Neglect of SOC and self-interaction errors can result in an incorrect position of the valence or conduction band edges, thus introducing spurious errors in defect energy levels and predicted defect concentrations. Du 75 showed how for the case of an iodine vacancy, a deep (0/+) donor level is predicted for GGA-noSOC, while a resonant donor level is predicted for GGA-SOC and hybrid functional with SOC treatments of electron-exchange and correlation.

B. Electron-phonon coupling
Going beyond the static lattice approximation with perturbation theory, we can consider the interaction of the electronic structure with vibrations of the lattice. Electronphonon coupling can perturb the electronic band energies (changing the bandgap) and couple electronic excitations (the hole and electron quasi-particles) into vibrational excitations (phonon quasi-particles). In a semiconductor, charge carrier scattering is often dominated by this electron-phonon interaction. The strength of these processes can set a limiting value on mobility. Electron-phonon coupling is often calculated in a second-order density functional perturbation theory calculated for a static (rigid ion) structure. For normal covalent systems, this term is expected to dominate over the first order contribution from the acoustic deformation potential as vibrations are typically small. These calculations are difficult to converge, as integration is over both electronic reciprocal space and vibrational reciprocal space, and the electron-phonon interaction is often found to be a non-smooth function. 76 In a recent work, 32 we developed a method to calculate the electron-phonon interaction of soft anharmonic phonon modes and applied this to the acoustic zone boundary tilting in hybrid halide perovskites described earlier. We solve a onedimensional Schrödinger equation for the nuclear degree of freedom along the phonon mode, and then combine the resulting thermalised probability density function (which includes zero point fluctuations and quantum tunneling) with a bandgap deformation potential along this mode. The method includes quantum nuclear motion, goes beyond the harmonic regime, but only contains the first-order contribution to the electronphonon coupling of the bandgap deformation. A positive bandgap shift of 36 meV (R point phonon) and 28 meV (M point phonon) is predicted at T = 300 K. Saidi et al. sampled all non-soft harmonic phonons using a Monte Carlo technique, 46 finding significant differences with the (more standard) perturbation theory results. The electron-phonon interactions can be calculated with MD, but as with phonon-phonon scattering, achieving convergence with respect to electronic (k-point sampling and basis set) and vibrational (q-point sampling and supercell size) parameters, while maintaining sufficient integration time to capture rare processes, is costly.
Recently a "one shot" method has been developed to calculate bandgap renormalization and phonon-assisted optical absorption and applied to Si and GaAs. 77 Nuclei positions are carefully chosen as a representative sample from the thermodynamic ensemble, and the electronic structure is needed for this static structure only-a significant increase in the computational efficiency. Such techniques may provide a promising method to calculate the electron-phonon coupling of complex materials, but so far are only valid in the harmonic phonon approximation. They have not yet been tested for the family of hybrid halide perovskites or other more complicated crystal structures.

C. Charge carrier transport
We now consider some aspects of charge carrier transport in hybrid halide perovskites. The minority-carrier diffusion length is defined as the average length a photo-excited (or electronically injected) carrier travels before recombining. In a working photovoltaic device, the diffusion length must be sufficient for photogenerated charges to reach the contacts. The minority-carrier diffusion length is a product of the diffusivity D and lifetime τ of minority charge carriers, L d =

√
Dτ. The minority-carrier diffusion lengths in MAPI are reported to be considerably longer than other solution processed semiconductors. 78 Long lifetimes (large τ) can be partly attributed to the "defect-tolerance" of hybrid perovskites (discussed in Sec. IV C), reducing the rate of ionized-impurity scattering and non-radiative recombination.
The effective mass of both electrons and holes in hybrid halide perovskites is small (though careful calculations including spin-orbit coupling indicate that the band extrema do not show a parabolic dispersion relation, and so the concept of effective mass is ill-defined 68 ). Given effective masses of <0.2m e , the carrier mobility of MAPI (<100 cm 2 V 1 s 1 ) is modest in comparison to conventional semiconductors such as Si or GaAs (>1000 cm 2 V 1 s 1 ). 4 The carrier mobility must be limited by strong scattering.
The low temperature mobility in this material reduces as a function of temperature as T 1.5 , which provides circumstantial evidence for being limited by acoustic phonon scattering. 79,80 However, if we only consider acoustic phonon scattering (which is elastic due to the population of acoustic modes), the calculated mobility is orders of magnitude larger than the experiment. A key realisation is that the soft nature of these semiconductors results in optical phonon modes (see Fig. 2) below thermal energy. 18,69 Optical phonon scattering is inelastic and dominates once the charge carriers have sufficient energy to generate the phonon modes. 81 Through solving the Boltzmann transport equation parameterised by the DFT calculations, at room temperature, the scattering from longitudinal optical phonons is most relevant in limiting mobility. 82,83 The carrier mobility will be further limited by scattering from point and extended lattice defects. 84 The fluctuations in the electrostatic potential resulting from the dynamic disorder provide a macroscopic structure from which carriers will also scatter. 43,85

IV. PHOTOPHYSICS AND SOLAR CELLS
Recent research interest in hybrid halide perovskites is mainly due to their use as the active layer in efficient solar cells. There are areas of the underlying physics which are not yet developed and which may be limiting progress in the field. Ion migration is poorly understood and has been correlated with hysteresis effects 88,89 and device degradation. Defects which act as recombination centres have not been identified and characterised. Additionally, interfaces have not been optimised for the optimal charge carrier extraction. We will now outline these issues where theory and simulation have much to contribute.

A. Ion migration
The charged point defects in the bulk allow for mass transport of ions and can result in spatial fluctuations of the electrostatic potential. For solid-state diffusion to be appreciable in magnitude, there needs to be a high concentration of defects and a low activation energy for diffusion.
The equilibrium concentration of charged vacancy defects is calculated as being in excess of 0.4% at room temperature in MAPI. 14 Low defect formation energies and free-carrier concentrations found across the hybrid halide perovskites indicate that Schottky defects are prevalent across this family of materials. While each point defect is charged, they are formed in neutral combinations so that a high concentration of lattice vacancies does not require a high concentration of electrons or holes to provide charge compensation.
The ion migration rate is given by where ∆H diff is the activation energy for solid-state diffusion and ν is the attempt frequency. In MAPI, the diffusion of methylammonium cations, iodide anions, and protons has been considered in the literature. 88,90,91 Activation energies calculated from first principles show that the predominant mechanism for ion migration is the vacancy assisted hopping of iodide ions. 88 Based on a bulk activation energy of 0.58 eV, 88 a rate of 733 hops/s would be expected at T = 300 K, with an associated diffusion coefficient of 10 12 cm 2 s 1 . Effective activation energies as low as 0.1 eV have been reported experimentally, 92,93 which likely correspond to diffusion along extended defects (dislocations, grain boundaries, surfaces). 94,95 The corresponding diffusion rate of 10 5 cm 2 s 1 is very fast, but comparable to the surface diffusion of iodine observed in other compounds. 96 It is also comparable to the diffusion coefficient of 4 × 10 −6 cm 2 s 1 predicted by the classical molecular dynamics. 97 Modeling ion diffusion at device scales is not yet possible with ab initio methods. Parametrised drift-diffusion modeling of ion and electron density indicate that slow moving ions can explain the slow device hysteresis. 89,98 A vacancy diffusion coefficient of the order of 10 12 cm 2 s 1 is consistent with both predictions and transient measurements. 88 It has been suggested that ion migration within mixedhalide compositions is the result of a non-equilibrium process J. Chem. Phys. 146, 220901 (2017) FIG. 4. Ion transport occurs in halide perovskites: they are mixed ionic-electronic conductors. The vacancy-mediated diffusion of halide anions has been associated with both the current-voltage hysteresis of solar cells and the rapid interchange between iodide, bromide, and chloride materials. The microscopic origin of the reversible ion segregation observed in mixed (Br,I) systems remains unresolved and a subject of debate. Alloyed materials have been found to phase separate upon illumination, but recover their initial state when the light source is removed. The phase separation is associated with a striking red-shift in the photoluminescence spectra. A statistical mechanical analysis of the ground-state DFT calculations suggested a large miscibility gap, 86 while the charge carriers generated upon illumination can provide an additional driving force for the phase separation. 87 The results from a simple thermodynamic model are shown in the right panel, where the free energy of mixing contains contributions from the ground state (∆F groundstate ) with an additional component due to the difference in bandgaps between the mixed (I,Br) and phase separated I-rich phases (∆F photo ). The latter contribution requires local carrier concentrations approaching 10 21 cm 3 to make a substantial contribution to the overall mixing energy. induced by photoexcitation. X-ray diffraction measurements by Hoke et al. 99 show that under illumination the mixed halide perovskite MAPb(I 1x Br x ) 3 segregates into two crystalline phases: one iodide-rich and the other bromide-rich. This segregation leads to reduced photovoltaic performance via charge carrier trapping at the iodide-rich regions. In some reports, after a few minutes in the dark, the initial single phase XRD patterns are recovered. This reversible process is unusual and defies the common assumption made that ion transport and electron transport are decoupled.
A schematic outline of the phase segregation process is shown in Fig. 4. A phase diagram constructed from the firstprinciples thermodynamics found a miscibility gap for a range of stoichiometries at room temperature. 86 This suggests that a mixed-halide material is metastable and will phase segregate after being excited by light, which follow a decreasing free energy gradient toward halide-rich areas formed prior to light excitation (such as grain boundaries). The accumulation of charge carriers increases the lattice strain and drives further halide segregation. Our calculations indicate that the transition between mixing and segregation will occur at a local carrier concentration of 10 21 cm 3 , which would require accumulation into small regions of the material.

B. Electron-hole recombination
The open-circuit voltage (V OC ) of a solar cell is determined by the rate of charge carrier recombination in the material, as no photogenerated charges are being extracted and so all are recombining. When operated to generate power, the rate of recombination competes with the rate of charge extraction, limiting the fill factor of the solar cell. In addition, rates of recombination contribute to the photovoltaic potential of a material.
Recombination is usually separated into three channels: non-radiative, radiative, and Auger. These, respectively, correspond to one-, two-, and three-electron processes. Assuming that the prefactors for the rates of these processes are constant, the carrier density in an intrinsic semiconductor can be modeled as a rate equation where G is the rate of electron-hole generation and n is the density of charge-carriers. While non-radiative recombination is limiting in many inorganic thin-film technologies, hybrid perovskites are not significantly affected. This is surprising for the high density of defects expected for a material deposited from solution at relatively low temperatures, leading to the material being described as "defect tolerant." 100 Radiative (bimolecular) recombination is slower than would be expected for a direct bandgap semiconductor. Recent calculations 101,102 revealed how relativistic Rashba splitting can suppress radiative recombination at an illumination intensity relevant to an operating solar cell. After photoexcitation, electrons thermalise to Rashba pockets in the conduction band minima away from the high symmetry point in the reciprocal space. This leads to an indirect charge recombination pathway as the overlap in k-space between occupied states near upper valence and lower conduction bands diminishes. It has also been suggested that direct recombination is suppressed due to the pockets of minima being spin-protected. 102 Direct gap radiative recombination is reduced by a factor of 350 at solar fluences, as electrons must thermally repopulate back to the direct gap. 101 This is in agreement with the temperature-dependence of the bimolecular rate measured experimentally 28 and calls into question the validity of models where a global radiative recombination rate independent of carrier concentration is used.
Auger recombination is only significant at fluences well above solar radiation, but it is important for understanding laser photophysics.
Ferroelectric effects could contribute to the electron-hole separation due to the electrostatic potential fluctuations in real space. Although the molecular cation plays no direct role in charge generation or separation, it could have a part to play in charge transport through the formation of polar domains. 85,103 Macroscopic ferroelectric order is not necessary to explain the device behavior in a 3D drift-diffusion simulation. 104 A multiscale Monte Carlo code based on a model Hamiltonian parameterized for the inter-molecular dipole interaction in MAPI explored the results of this dynamic polarisation. 43 This predicts the formation of antiferroelectric domains which minimise energy via the dipole-dipole interaction, which work against a cage-strain term preferring ferroelectric alignment. 36 This provides electrostatically preferred pathways for electrons and holes to conduct. Developing more accurate models and measurements of the nature and effects of lattice polarisation in these materials are the subject of on-going research efforts.

C. Defect levels in the bandgap
To understand why the rate of non-radiative recombination is low, we consider the known defect properties of hybrid perovskites. Defects appear to have a minimal impact upon charge carrier mobility and lifetime, 105 which can be attributed to a combination of large dielectric constants and weak heteropolar bonding.
Under the Shockley-Read-Hall model for semiconductor statistics, non-radiative recombination is mediated through deep defect states in the gap. 106 Shallow defect states can act as traps but the carriers are thermally released to the band before recombination can occur. Hybrid perovskites-with high dielectric constant and low effective mass-show a tendency toward benign shallow defects under the hydrogenic model 107 where m * m 0 is the effective mass ratio, 0 is the static dielectric constant, and n is an integer quantum number for the given energy level. Atomic units are used and so energy is given in Hartrees.
In Table I, we give the first hydrogenic defect level for MAPI, Si, and CdTe, where the binding energy for MAPI is only 3 meV. For ionic materials, one would expect a large central cell correction that could result in much deeper levels, for example, as seen for the color centres in alkali halides. 108 It was shown numerically that the on-site electrostatic potentials in the I-II-VII 3 perovskites are relatively weak owing to the small charge of the ions (e.g., Cs + Pb 2+ I − 3 ) compared to other perovskite types (e.g., Sr 2+ Ti 4+ O 2− 3 ), 103,109 which would also support more shallow levels. In addition, arguments based on covalency have also been proposed. 105

D. Beyond the bulk: surfaces, grain boundaries, and interfaces
As perovskite solar cells approach commercial viability, 6 there are considerations to be made beyond the bulk materials. Surfaces, grain boundaries, and interfaces will influence the device performance and long-term stability and become increasingly important as the science is scaled up from lab to production line. Accurate interface modeling requires consideration of halide migration, ion accumulation, charge carrier transport, and charge carrier recombination at the defect states. There has been preliminary work that provides insights, but real systems offer much deeper complexity.
Perovskite films fabricated through solution processing methods are multicrystalline and so the formation of grain boundaries is inevitable. The resulting microstructure provides pathways for ion conduction, electron-hole separation, and recombination. The shallow traps introduced are evidenced through improved device performance with increasing grain size 112 and their thermal activation. Calculations suggest that grain boundaries do not introduce deep defects and consequently have negligible effect upon the rate of nonradiative recombination. 113,114 This is in conflict with spatially resolved photoluminescence 115 and cathodoluminescence 116 measurements which evidence greater non-radiative loss at grain boundaries.
Recent calculations using nonadiabatic MD and timedomain DFT 117 indicate that grain boundaries localize the electron and hole wavefunctions and provide additional phonon modes. This leads to increased electron-phonon coupling which in turn will give a higher rate of non-radiative recombination.
The typical device structure for a perovskite cell is the perovskite absorber layer sandwiched between an electron transport layer (e.g., TiO 2 , SnO) and a hole transport layer (e.g., spiro-OMeTAD, PEDOT-PSS). At the interface, there are two key considerations. One is that the bands should be electronically matched so as to allow efficient charge extraction without large energy loss. The second is that the formation of defects should be minimised as these act as sites for recombination, can lead to mechanical degradation of the device, and have been linked to hysteresis. 118 The commonly used hole transporter spiro-OMeTAD is hygroscopic so that stability in humid air is a concern. 119 This has prompted the development of screening procedures 120,121 to identify alternative contacts. The electroniclattice-site (ELS) figure of merit considers band alignment, lattice match, and chemical viability via the overlap of atomic positions. 120 Using this figure of merit, Cu 2 O is identified as a possible earth abundant hole extractor, while oxide perovskites, such as SrTiO 3 and NaNbO 3 , have been identified as possible electron extractors. As with the majority of screening techniques, the candidate materials meet the necessary but perhaps not sufficient conditions. Further refinements may consider the change in electronic properties as the lattice strain and chemical inhomogeneity at the interface is introduced. II. A collection of common issues that can arise in the simulation of hybrid perovskites. Note that for convergence of supercell size, unusual behavior can be observed due to the fact that octahedral titling modes of perovskites are allowed in even cell expansions (e.g., 2 × 2 × 2) and suppressed in odd cell expansions (e.g., 3 × 3 × 3) of the cubic lattice. The lattice dynamics are particularly sensitive to basis set convergence and plane-wave codes may require an energy cutoff 25% higher than a typical electronic structure calculation. For a cubic halide perovskite, k-point sampling of at least 6 × 6 × 6 is required to give reasonable total energy and electronic structure, so a Γ-point approximation is only valid for very large supercells and should be tested carefully for the property of interest.

V. CONCLUSION
We have outlined many of the physical properties that make hybrid perovskites unique semiconductors, but also challenging for contemporary theory and simulation. A number of practical points relating to the issues we have encountered while running simulations of these materials are summarized in Table II.
The volume of work in this field has not allowed us to address all active areas of research, including that around perovskite-like structures with lower dimensionality (e.g., Ruddleston-Popper phases) 5,122,123 and double perovskites with pairwise substitution on the B site, 21,22,124,125 which are both attracting significant interest. The optoelectronic properties of inorganic perovskites such as CsSnX 3 and CsPbX 3 (X = Cl, Br, I) are also promising and provide fertile ground for future research, especially for applications in solid-state lighting. 126,127 Attempts are also being made to distil our understanding of halide perovskites into computable descriptors for large-scale screening toward the design and discovery of novel earth-abundant non-toxic semiconductors. 105,[128][129][130] There are many opportunities ahead as we pick apart the relationship between organic and inorganic components, electronic and ionic states, as well as order and disorder in this complex family of materials.

ACKNOWLEDGMENTS
The primary research underpinning this discussion was performed by Jarvist M. Frost (molecular dynamic and Monte Carlo investigations), Federico Brivio (crystal and electronic structure), Jonathan M. Skelton (lattice dynamics and vibrational spectroscopy), and Lucy Whalley (bandgap deformations). We are indebted to our large team of external collaborators including the groups of Mark van Schilfgaarde, Saiful Islam, Simon Billinge, Piers Barnes, and Mark Weller. L.D.W. would like to acknowledge support and guidance from the staff and students at the Centre for Doctoral Training in New and Sustainable Photovoltaics. This work was funded by the