Effect of disorder on transport properties in a tight-binding model for lead halide perovskites

The hybrid organic-inorganic lead halide perovskite materials have emerged as remarkable materials for photovoltaic applications. Their strengths include good electric transport properties in spite of the disorder inherent in them. Motivated by this observation, we analyze the effects of disorder on the energy eigenstates of a tight-binding model of these materials. In particular, we analyze the spatial extension of the energy eigenstates, which is quantified by the inverse participation ratio. This parameter exhibits a tendency, and possibly a phase transition, to localization as the on-site energy disorder strength is increased. However, we argue that the disorder in the lead halide perovskites corresponds to a point in the regime of highly delocalized states. Our results also suggest that the electronic states of mixed-halide materials tend to be more localized than those of pure materials, which suggests a weaker tendency to form extended bonding states in the mixed-halide materials and is therefore not favourable for halide mixing.

amount of disorder can lead to a phase transition from a metallic to an insulating state. This raises the question of why the mixed-halide perovskites have good transport properties in spite of being disordered. The answer to this question is unknown, and there are ongoing attempts to elucidate the origin of the unusual transport properties of these materials, as can be seen in ref. 17. In particular, it has been suggested that the nature of the atomic orbitals that form the valence and conduction bands is unfavourable for the formation of deep trap states where non-radiative recombination takes place 18 . In addition, a high mobility is important to extract carriers faster than they can recombine or be trapped, which is possible in a material where states remain largely delocalized despite the large degree of disorder caused by the random orientations of the organic molecules, compositional variations (iodide and bromide mixing, formamidinium [FA] and MA mixing, etc.), and structural imperfections (vacancies, etc.). Interestingly, ref. 12 suggested that disorder could in fact be helping in electron-hole separation, which is an important step in the operation of a solar cell.
The effects of disorder on the electronic properties of these materials can be investigated using simulations of disordered structures. First-principles calculations with more than a few unit cells are in principle possible but computationally demanding 12 . One can instead use a tight-binding model, which reduces the number of degrees of freedom per unit cell drastically and allows the simulation of larger systems 19,20 . We take the model developed in ref. 21, including the parameters obtained there, as the starting point for our calculations. This simplification allows us to perform systematic simulations on systems containing 512 unit cells (effectively containing ~6000 atoms) at a rather low computational cost.

Results
Tight-binding model. We consider a model where electrons can occupy and hop between the Pb and halide atomic sites, and at each site there are four orbitals available to it, one s and three p orbitals. More specifically, for MAPbI 3 the relevant orbitals are Pb 6s, Pb 6p, I 5s and I 5p. We do not include any orbitals associated with the MA molecules, because these orbitals are far away from the valence band maximum (VBM) and conduction band minimum (CBM). We focus on the cubic phase, where the structure exhibits symmetry between the three axes x, y and z. In this case, the halide atoms are arranged to form perfect octahedra that are all aligned with each other, which keeps the number of system parameters (i.e. orbital energies and orbital overlaps) at the reasonable value 9. Spin-orbit coupling (SOC) is also included in the model, which adds two more parameters to the model, namely the SOC strengths for Pb and I. Treating the orbital energies and overlaps as fitting parameters, Boyer-Richard et al. obtained system parameters that reproduce the main features in the electronic band structure of MAPbI 3 21 .
It is worth mentioning here that a recent study used a somewhat similar tight-binding model to investigate SOC effects in a different family of halide perovskite materials 22 .
The model is described by the Hamiltonian: Figure 1. Structure of the perovskite material MAPbI 3 . The I atoms are located at the edges of the octahedra, with the Pb atoms at the centers of the octahedra and the MA molecules at the centers of the spaces between the octahedra. The disorder in the systems has two sources: (1) the random, and in fact temporally fluctuating, orientations of the MA molecules producing a random electrostatic potential landscape in the lattice and (2) the substitutional halide atoms, i.e. Br atoms replacing some of the I atoms.
SCIEnTIFIC RepoRts | 7: 8902 | DOI: where  i,α is the energy of orbital α (i.e. s or p) at site i, t iαjβ is the hopping strength between orbital α at site i and orbital β at site j, α a i, is the particle annihilation operator for orbital α at site i, and ˆ † α a i, is its hermitian conjugate. Spin-orbit coupling, which couples different p orbitals and different spin states, for atom A (i.e. Pb or I) is described by the Hamiltonian: In the absence of disorder, we set  i,α , α β t i j and λ A to the values given in ref. 21. Two types of disorder are present in the problem that we wish to study. The disorder in the orientations of the MA molecule can be simulated by introducing random fluctuations in the on-site energies α i,  , which we take to depend only on the location of the orbital and therefore to be independent of α. A substitution of I by Br in general results in changes in the orbital energies as well as the hopping strengths between orbitals at neighboring sites 23 . We consider both of these mechanisms as possible causes of the changes that accompany halide substitution. Consequently, we explore two different ways to investigate the effect of Br substitution. In one set of calculations, we simulate mixed-halide materials by setting the Br orbital energies to lower values than the I orbitals, while keeping the hopping strengths fixed. We find that in order to obtain the correct band gap for MAPbBr 3 we must set the Br orbital energies 0.86 eV lower than those of I, and this is the value that we use in our simulations. In another set of of calculations, we fix the orbital energies and set the hopping strength between Pb and Br sites to be 0.73 of the value for Pb-I hopping, where as above the value 0.73 is chosen because it gives the correct band gap for MAPbBr 3 .
An interesting question here relates to the definition of the ordered state. Since each MA molecule contributes a potential energy to the different orbitals depending on its orientation, the ordered state would have well-defined orientations of the molecules, leading to a nonzero potential felt by electrons at neighbouring atoms. Recent molecular dynamics simulations suggest some correlation between neighbouring MA molecules 8,9,[24][25][26][27] . The correlations are such that in one direction the MA molecules tend to be anti-parallel to each other and the other two directions neighbouring molecules tend to be perpendicular to each other. Such correlations imply that the MA molecules will not produce a macroscopic polarization that will lead to a large contribution to the electric field inside the crystal. In our model, we ignore any net contribution to the potential from the MA molecules in the ordered state and effectively treat the molecules as electron donors in that case.
Another question relates to the logic behind taking experimental parameters that are measured for a disordered system and using them to extract the parameters of the disorder-free model. The idea of this approach would be that there are two possible outcomes: if we find that the physical properties (and in particular the band structure) of the ordered and disordered states are only slightly different from each other, then the results justify our use of ignoring the difference between the ordered and disordered states when calculating the parameters of the model based on the properties of realistic, and hence disordered, systems. If on the other hand we find that the physical properties of the ordered and disordered states are drastically different, we would have to go back and reevaluate our basic assumptions. As we shall see below, one of our main conclusions is that realistic levels of disorder do not significantly modify the electronic states of the materials under study.
Inverse participation ratio. Disorder typically has the effect of turning the extended energy eigenstates of perfectly ordered potentials into localized states. This localization in the energy eigenstates then translates into the transport properties, e.g. resulting in a reduced mobility. A good quantity that can be used to analyze such localization properties is the inverse participation ratio (IPR), which for a given quantum state ψ i ( ) is given by: the material. We shall also analyze the IPR for the states at the VBM and CBM, which can serve as indicators for the transport properties (e.g. the mobility) of holes and electrons. For the VBM and CBM calculations, we take the mean values for the states that lie within 0.1 eV from the band edges.

Computational details. For disorder of strength 
∆ in the on-site energies α i,  , we include a Gaussian-distributed random contribution to the orbital energies with the standard deviation of such an ensemble set to ∆. Our use of the Gaussian distribution is motivated by the central limit theorem. Each I atom is surrounded by four equally distant MA molecules in its immediate vicinity followed by larger numbers of additional MA molecules when we consider neighbouring unit cells. Each Pb atom has eight equally distant MA molecules serving as its nearest-neighbour molecules. These numbers are far from approaching infinity, but they are sufficiently large to suggest that regardless of the exact distribution of the potential energy fluctuation contributed by each molecule the sum of these contributions will have an overall shape resembling a Gaussian distribution. As mentioned above, we shift all four orbitals at a given site by the same amount. We apply the same rules for calculating the fluctuations to both Pb and I atoms. When simulating mixed-halide materials, we first generate configurations with the two halide species randomly distributed in the lattice (with a specified average concentrations x and 1 − x) and then we either lower the orbital energies at the Br sites by 0.86 eV or reduce the coupling strengths between the Br and Pb sites by the factor 0.27, as explained above.
We investigate the properties of the system by considering a supercell composed of 8 × 8 × 8 cubic unit cells. Since each cubic unit cell contains one Pb atom and three halide atoms and we are including four orbitals and two spin states for each atomic site, our supercell has 16,384 states that an electron can occupy. The size of the Hilbert space is therefore 16,384. Each basis state in the 16,384-dimensional Hilbert space corresponds to the electron being in a given orbital at a given atomic site in the supercell with spin state up or down. We set the quasi-momentum to zero in our calculations.
Calculation results. We now present the results of the numerical calculations explained above. For reference and in order to have some insight about the system in the tight-binding model, we start with some general considerations. In Fig. 2(b) we show the band structure in the absence of disorder, calculated using a single unit cell and a variable quasi-momentum. In Fig. 2(a) we show the band structures that one would obtain if SOC is ignored. As can be seen from the figure the most obvious effect of including SOC is to create a splitting in the conduction band and reduce the band gap, compared to the case where SOC is ignored. In Fig. 2 we show the density of states (DOS) as a function of energy, along with the contributions from the different orbitals, for the case including SOC. The DOS shows that the s orbitals contribute mostly to states at the bottom of the band structure in our model. To investigate whether these orbitals need to be kept in the model, we calculate the band structure for a model where the s orbitals are ignored and we plot the results in Fig. 2(d). One can see clearly that ignoring the s orbitals strongly modifies the valence band, including the states that determine the band gap of the material. The s orbitals therefore need to be kept in the model.
It is interesting to note that the band structure in Fig. 2(b) does not exhibit a Rashba splitting, although evidence of such a splitting has been observed in recent experiments and it is believed to be related to SOC 32 . The reason is that the Rashba splitting requires both SOC and inversion symmetry breaking 22,33 . Our model does not contain any inversion symmetry breaking. In the supplementary material we break the symmetry and show that the Rashba splitting appears in that case. The exact microscopic mechanism for inversion symmetry breaking in the materials under study is presently not well understood. Once this mechanism is established, the tight-binding model can be modified to represent it accurately. The Rashba splitting should not affect our results on the localization of electronic states presented below.
The energy levels in Fig. 2(a-c) can be divided into three bands: a deep valence band (to which we shall refer as VB −1) that lies a few eV below the valence band, the valence band (VB) the conduction band (CB).
Disorder in on-site energies. We would like to investigate the localization properties of the energy eigenstates. First, to visualize the localization process, we plot in Fig. 3 the probability distribution for a few different disorder settings. As the disorder strength increases, the states become increasingly localized and for any given state the occupation probability becomes concentrated in an increasingly small number of unit cells. In this figure and below, we shall go up to a disorder strength of 10 eV, which is higher than what we can expect in a realistic perovskite material, in order to show the behaviour in the limit where the electronic states becomes localized to a single or very few unit cells.
In Fig. 4 we plot the IPR as a function of ∆. The IPR for both the VB and CB remains close to zero all the way up to ∆ . 0 2 eV, after which it starts increasing and approaches unity for large values of  ∆ . If we take the mean values of the IPR over entire energy bands, the increase of the IPR is gradual. The VB −1 is the most sensitive band to disorder, and its IPR grows the fastest among the three bands. This feature can be understood by considering that this band is formed from s orbitals, which have a hopping strength that is smaller than the hopping strength between p orbitals by a factor of 3-4. The result that the electronic stats become localized when the disorder level reaches the range 0.2-1 eV can be understood intuitively by considering that each band in the band structure in Fig. 2(b) has a width on the order of 1 eV. As a result, one can intuitively expect that localization will occur when the disorder strength becomes comparable to the energetic width of the bands.
The IPR for the states at the VBM and CBM remains very close to 1/512 up to ∆ ~ 0.3-0.5 eV [ Fig. 4(b)]. The insensitivity of these states to small amounts of disorder can be understood based on the fact that the DOS is very small at the VBM and CBM [see Figs 2(c) and SI1(a)]. The disorder-induced localization can be understood as a hybridization between different basis states caused by the disorder when seen as a perturbation to the Hamiltonian in the zero-disorder limit. Energy levels that are far from other energy levels tend to hybridize less with other energy levels as a result of perturbations in the Hamiltonian. This insensitivity of the VBM and CBM states to SCIEnTIFIC RepoRts | 7: 8902 | DOI:10.1038/s41598-017-09442-4 disorder will be more clearly visible when we discuss halide mixing below. Once  ∆ reaches 0.5 eV or higher, however, the increase in the IPR of the VBM and CBM state is quite abrupt, which might reflect the onset of a localization phase transition. The VBM states become localized for a somewhat smaller value of ∆ compared to the CBM states, which implies that hole transport should be more sensitive to on-site energy fluctuations than electron transport.
Halide mixing. We now consider the case of disorder resulting from halide mixing, specifically the substitution of Br atoms for some of the I atoms in MAPbI 3 to give MAPbI (3−x) Br x . The IPR as a function of Br concentration x is plotted in Fig. 5. We distinguish between five different IPR values: the mean values over all states in the VB −1, VB and CB and the IPR values of the states at the VBM and CBM. The mean values over the occupied states can reflect the bond strengths in the material: delocalized states typically correspond to stronger bonding. In the case where the orbital energies change when substituting Br for I, the mean value of the IPR over the VB −1 states exhibits the largest variations as x is varied, with the dependence being weak for small values of x but with a large peak at x = 0.8. In the model where the hopping strengths change when substituting Br for I, all three bands exhibit some increase in the IPR for mixed-halide materials with the most pronounced increases occurring at x = 0.8. These results suggest that halide mixing with relatively small Br concentrations is easier than mixing with small I concentrations. The IPR values near the band edges reflect the transport properties of electrons and holes. We find that (in both disorder models) the IPR for states near the VBM increases by about 10% in 50-50 mixed-halide materials as compared to pure I or pure Br materials, which suggests that the hole mobility will be almost unaffected by halide mixing. The IPR for states near the CBM is essentially independent of halide concentrations, suggesting that the disorder resulting from halide mixing has a lower effect on the electron mobility than it does on the hole mobility. This result makes sense, because the states near the CBM are composed mainly of Pb orbitals.
We note here that recent experiments suggest that doping with Cl might improve transport properties 34,35 . However, these results are widely believed to be related to defect formation during material synthesis. Such defects are not included in our model, and therefore we would not expect to capture these effects in our calculations.

Discussion
In order to discuss the relevance of our results on on-site disorder to the materials in the family of MAPbI 3 , we need to have an estimate for the parameter ∆ in these materials. Fluctuations related to the fluctuating orientation of the MA molecule were investigated in ref. 9 using molecular dynamics simulations, where the energetic barrier for the MA molecule rotation was found to be around 10-15 meV. Reference 36 calculated the partial charges on the different atoms in the MA molecule in MAPbI 3 . Using these charges and the relative positions of the atoms in the material, Coulomb's law gives that the rotation of the MA molecule would produce energy fluctuations on the order of 0.2-0.3 eV at the locations of the neighbouring atoms. If we include the fact that the dielectric constant is on the order of 10, the charges on the MA molecule will be screened by a similar factor, and we again obtain the similar estimate 0 02 ∆ .  eV. Figure 4 shows that the states start to become localized around  ∆ = .
0 2 eV, which is about one order of magnitude larger than the value of  ∆ that results from the MA orientation fluctuations. We note here that ref. 12 estimated potential fluctuations to be in the range 0.2-0.5 eV. These numbers are consistent with the electronic state localization found there.
Our discussion above has concentrated on the IPR, which quantifies the localization of the energy eigenstates of the system. This localization should be reflected in the transport properties. In particular, the mobility is defined as where v d is the drift velocity obtained upon applying an electric field of strength E app . For a given (incoherent) scattering rate, the drift velocity can be expected to scale as the spatial size of the relevant wave functions. The argument here is that after a given scattering event the charge carrier moves ballistically until the next scattering event, and when the electronic states have a certain characteristic spatial size the electrons will remain within the same length scale between scattering events. As a result, one can expect that the mobility is proportional to  , while blue corresponds to lower probabilities and red corresponds to higher probabilities (such that strong colours, especially with a mostly blue background and a few red squares, imply strongly localized states). Panels (a-c) correspond to on-site disorder with 0 1  ∆ = . eV (a), 1 eV (b) and 10 eV (c). The tendency for the wavefunction to localize with increasing ∆ is clear from these panels. Panels (d) and (e) show the probability distribution for disorder originating from halide mixing. Here we assume disorder in the hopping strengths, and we take a 50-50 mixture of I and Br. In Panel (d) we show the results for the model where the Br orbital energies are different from those of I, while in Panel (e) we show the results for the model where the hopping strength between Pb and Br is different from that between Pb and I. IPR −1/3 . These relations agree with the intuitive expectation that localization of the energy eigenstates, which leads to increased values of the IPR, leads to a reduction in the mobility.
To conclude, motivated by the emerging lead halide perovskite materials and their remarkable transport properties, we have analyzed the localization properties of the energy eigenstates of a tight-binding model that contains disorder stemming from random MA molecule orientations or from halide mixing. The IPR increases as the on-site disorder strength is increased, indicating a degradation of transport properties. However, our results suggest that the disorder strength in the lead halide perovskites puts them in the regime of delocalized states, where the mobility is not significantly reduced because of the disorder. We also find that the disorder originating from halide mixing has a minimal effect on the IPR near the band edges and therefore on the mobility. However, we do observe more dramatic changes deeper in the VB, which could indicate weaker bonds for the mixed-halide perovskites, which in turn is consistent with the phase separation that has been observed experimentally. Combined with other favourable material properties, such as the absence of deep trap states, our results help explain the good transport properties of these materials.  Figure 4. The IPR as a function of the on-site disorder strength ∆. In Panel (a) we plot the mean value of the IPR for the three energy bands: VB −1 (red squares), VB (green circles) and CB (blue triangles). The IPR for all three bands increases gradually as  ∆ approaches and exceeds 1 eV. We note here that although we are plotting the mean values over five instances, instance-to-instance fluctuations were small and almost indiscernible in the majority of cases. In Panel (b) we plot, for five instances, the IPR for the states at the VBM (red symbols) and CBM (green symbols). When we calculate the IPR for the VBM and CBM, we take the averages over the states that lie within 0.1 eV from the band edges. The solid and dashed lines are the respective mean values for the five instances. We can see that both lines increase and decrease again as  ∆ increases. The counter-intuitive decrease occurs in the regime where there is no boundary between the VB and CB, and this delocalization tendency should therefore not be taken as reflecting the transport properties at the VBM and CBM. In (b) we see that the tendency towards localization is more sudden than when taking the mean values for entire energy bands. The localization of the states at the VBM and CBM occurs around ∆ . − 0 2 1  eV, with the VBM states becoming localized for somewhat smaller values of ∆.  Figure 5. The IPR as a function of Br concentration x. Panels (a) and (b) correspond, respectively, to the models where the orbital energies and hopping strengths are modified when an I atom is replaced by a Br atom. The different data sets correspond to the mean values of the IPR over the VB −1 (red squares), VB (green circles) and the CB (blue triangles), as well as the IPR values for for VBM (magenta inverted triangles) and CBM (cyan diamonds).