⇦Chinese Physics Letters, 2022, Vol. 39, No. 4, Article code 046101Express Letter Defects in Statically Unstable Solids: The Case for Cubic Perovskite $\alpha$-CsPbI$_3$ Xiaowei Wu (吴晓维)1†, Chen Ming (明辰)1†, Jing Shi (石晶)2†, Han Wang (王涵)3, Damien West4, Shengbai Zhang (张绳百)4, and Yi-Yang Sun (孙宜阳)1* Affiliations 1State Key Laboratory of High Performance Ceramics and Superfine Microstructure, Shanghai Institute of Ceramics, Chinese Academy of Sciences, Shanghai 201899, China 2Department of Physics, Jiangxi Normal University, Nanchang 330022, China 3Chemical Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, California 94720, USA 4Department of Physics, Applied Physics, and Astronomy, Rensselaer Polytechnic Institute, Troy, New York 12180, USA Received 24 December 2021; accepted 7 March 2022; published online 10 March 2022 †These authors contributed equally to this work.
*Corresponding author. Email: yysun@mail.sic.ac.cn Citation Text: Wu X W, Ming C, Shi J et al. 2022 Chin. Phys. Lett. 39 046101    Abstract High-temperature phases of solids are often dynamically stable only. First-principles study of point defects in such solids at 0 K is prohibited by their static instability, which results in random structures of the defect-containing supercell so that the total energy of the supercell is randomly affected by structural distortions far away from the defect. Taking cubic perovskite $\alpha$-CsPbI$_3$ as an example, we first present the problem incurred by the static instability and then propose an approach based on molecular dynamics to carry out ensemble average for tackling the problem. Within affordable simulation time, we obtain converged defect ionization energies, which are unattainable by a standard approach and allow us to evaluate its defect tolerance property. Our work paves the way for studying defects in statically unstable solids.
cpl-39-4-046101-fig1.png
cpl-39-4-046101-fig2.png
cpl-39-4-046101-fig3.png
cpl-39-4-046101-fig4.png
cpl-39-4-046101-fig5.png
 DOI:10.1088/0256-307X/39/4/046101 © 2022 Chinese Physics Society Article Text In recent years, inorganic halide perovskites have demonstrated high performance in photovoltaic[1–4] and light-emitting[5–8] applications. Among these materials, CsPbI$_3$ has been most intensively studied for making solar cells. Currently, the highest power conversion efficiency based on CsPbI$_3$ has reached 20.4%.[9] Compared with its counterparts containing organic cations (e.g., CH$_3$NH$_3$PbI$_3$), CsPbI$_3$ has higher thermal and chemical stability.[10,11] The main issue of CsPbI$_3$ is that, at room temperature, it stabilizes in a non-perovskite phase,[12] instead of the $\alpha$-phase. The latter is the desired phase for making solar cells. Several breakthrough studies have successfully stabilized CsPbI$_3$ in the $\alpha$-phase[1,2,13] or other distorted perovskite phases, such as $\beta$-and $\gamma$-CsPbI$_3$,[9,14–17] which are all shown to be promising for making high-efficiency solar cells. There appears a question whether CsPbI$_3$ can take over CH$_3$NH$_3$PbI$_3$ and its organic family if the stability issue can be favorably solved. Defect tolerance of CH$_3$NH$_3$PbI$_3$ and similar lead halide perovskites is among the most attractive properties of these materials.[18–25] It is perceived that they could be free of intrinsic deep-level defects that serve as carrier recombination centers, the most undesired defects in solar cells. This finding is an important contribution to the defect physics of semiconductors beyond the materials themselves. This property also allows the materials to be prepared without demanding control on the purity and crystallinity, which facilitates the low-cost production of solar cells. However, the static instability of $\alpha$-CsPbI$_3$ presents a serious issue for theoretical study of its defect properties. As we will show, the commonly used approach based on the density functional theory (DFT) fails to produce credible results on the positions of defect levels in the band gap and, therefore, cannot be used to study the defect tolerance of $\alpha$-phase CsPbI$_3$. The dynamic behaviors of defects in halide perovskites were discussed recently from different perspectives. For example, Rakita et al. noted that the fast diffusion of intrinsic defects is in accord with the “self-healing” phenomenon observed in these materials.[26] Cohen et al. found that the eigenvalues of defect states in halide perovskites are dynamic by exhibiting substantially larger fluctuation within the band gap than that in traditional semiconductor materials,[27] and noted that the local and instantaneous arrangement of atoms causes the large fluctuation. The structural instabilities have been discussed for bulk halide perovskites.[28,29] However, how static instability of bulk materials affects defect calculations has not yet been studied even though this may not be a serious issue for studying the stable phases of CsPbI$_3$.[30–33] In this Letter, we propose a first-principles approach to calculating defect transition levels (or ionization energies) at finite temperature, critically important for statically unstable systems such as $\alpha$-CsPbI$_3$. The calculation is carried out by performing molecular dynamics (MD) simulations, in which the structure is highly distorted from the high-symmetry structure of $Pm\bar{3}m$ space group at any given time. The distortions are strong enough to break the harmonic approximations for the lattice dynamics so that no stable phonon spectrum can be obtained from the high-symmetry structure. In the presence of a defect, the distortions as visualized by the [PbI$_6$] octahedron rotations are at random. This fails the static calculations as commonly used for studying defects. As we will show, statistical ensemble average in the time period of typical MD simulations is sufficient to determine the defect transition levels, which in turn allows us to evaluate the defect tolerance of $\alpha$-CsPbI$_3$. Computational Method. Our calculations were carried out with the Vienna ab initio simulation package (VASP)[34] using projector augmented wave (PAW)[35,36] potential to describe the interaction between ion cores and valence electrons and using plane waves as basis set. The generalized gradient approximation (PBEsol)[37] was used for the exchange-correlation functional. A $3\times3\times3$ supercell including 135 atoms was used for modeling the intrinsic defects. The Brillouin zone was sampled by the special $k$-point (0.5, 0.5, 0.5), where the band gap is located. The cutoff energy of the plane waves was set to be 238 eV. This setup of $k$-point and cutoff energy yields a band gap of 1.17 eV for the $3\times3\times3$ supercell without structural distortions. As a convergence test, we checked the unit cell using a $10\times10\times10$ $k$-grid and 500 eV cutoff. The obtained band gap is 1.19 eV. MD simulations were conducted with the $NVT$ ensemble using the Nosé–Hoover thermostat.[38,39] A time step of 2.5 fs was used. Other computational details will be described in the following. Stability of Bulk $\alpha$-CsPbI$_3$. To illustrate the problem, we first show in Fig. 1(a) the variation of the band gap of $\alpha$-CsPbI$_3$ calculated by DFT-based MD simulation at 650 K. The averaged band gap is 1.67 eV, as marked by a dashed line through the fluctuating MD results. For comparison, we also show the band gap from the high-symmetry $Pm\bar{3}m$ structure as a black dashed line, which is 1.17 eV only. Such a small band gap in a high-symmetry structure has also been noted in early theoretical studies on the halide perovskites.[40] The large difference between the band gaps obtained from the two ways mentioned above reflects the strong electron-phonon coupling in this material. As the defect transition levels are all discussed with respect to the band gap, an underestimated band gap will directly affect the calculation of defect properties. Figure 1(b) shows the snapshot at the last step from a 25-ps MD simulation at 650 K. The structure is visibly distorted from the perfect $Pm\bar{3}m$ structure. However, we note that this does not affect the fact that a diffraction measurement will yield a high-symmetry structure. Figure 1(c) shows the averaged structure based on the MD trajectory from the last 20 ps of the whole simulation. The averaged structure reasonably recovers the $Pm\bar{3}m$ symmetry with only slight distortion due to the finite averaging time. In practical diffraction measurements, such as x-ray diffraction, the measuring time is several orders of magnitude longer than 20 ps. Therefore, a $Pm\bar{3}m$ structure is obtained.
cpl-39-4-046101-fig1.png
Fig. 1. (a) Variation of the band gap of $\alpha$-CsPbI$_3$ in DFT-based MD simulations at 650 K. Black dashed line shows the band gap of the structure with the $Pm\bar{3}m$ symmetry. (b) Snapshot of the atomic structure of $\alpha$-CsPbI$_3$ at the last step of the MD simulation. (c) Averaged structure based on the MD trajectory from the last 20-ps MD simulation. (d) Phonon spectra of $\alpha$-CsPbI$_3$ at 0 K from finite-difference calculation (left) and at 650 K from TDEP calculation (right).
To further show the static instability and dynamic stability of $\alpha$-CsPbI$_3$, we compare the phonon spectra obtained at 0 K and 650 K in Fig. 1(d). The spectrum at 0 K obtained by standard finite-difference method[41] indicates severe instability by exhibiting large imaginary vibrational modes across the whole Brillouin zone, which cannot be described by common distortions in perovskite structures, namely, imaginary modes at zone-center indicating proper ferroelectricity, at zone-boundary indicating rigid octahedron titling, and Jahn–Teller distortion.[29] In contrast, the spectrum at 650 K obtained by temperature-dependent effective potential (TDEP) approach[42] confirms the dynamical stability of $\alpha$-CsPbI$_3$ at this temperature. To obtain the spectrum at 650 K, the MD trajectory during the last 20-ps simulation (8000 MD steps) was used. In the following, we will use “static instability” to describe the case that, at 0 K, the material does not possess a stable static structure.
cpl-39-4-046101-fig2.png
Fig. 2. (a) Schematic of the potential energy surface (PES) of the $v_{_{\scriptstyle \rm I}}$ defect in $\alpha$-CsPbI$_3$. The red arrow marks the unstable-balance point on the PES corresponding to the middle structure, which is obtained by directly removing an I atom from the supercell with the $Pm\bar{3}m$ symmetry and then relaxing to meet the zero-force criterion. The left black arrow marks the lowest-energy structure obtained by following the imaginary modes from a frequency calculation on the middle structure. The right black arrow marks the structure after six iterations of calculations following the imaginary modes. (b) Energy evolution of the $v_{_{\scriptstyle \rm I}}$ defect in six iterations of imaginary mode-following calculations (see the text).
Issue due to Static Instability in Defect Study. We next show the problem in a defect calculation incurred by the static instability of $\alpha$-CsPbI$_3$. Take the iodine vacancy ($v_{_{\scriptstyle \rm I}}$) as an example. If starting from a supercell based on the $Pm\bar{3}m$ structure, we end up with a highly symmetric structure for the $v_{_{\scriptstyle \rm I}}$ defect after structural relaxation. However, the structure is located at an unstably balanced point in the configuration space, as marked by an upward arrow in Fig. 2(a). From a frequency calculation, it is found that this structure exhibits 39 imaginary frequencies. By following each frequency, we further relax the structures. As shown in Fig. 2(b) (step 1), all the imaginary modes lead to structures of lower energy, and the energies for several modes are significantly lower than the high-symmetry structure by up to 0.34 eV. By picking the lowest-energy structure, we further conduct a frequency calculation. Again, 33 imaginary frequencies are found [step 2 in Fig. 2(b)]. By iterating this process several steps (steps 3–6), the total energy seemingly reaches a plateau with total energy about 0.57 eV lower than the high-symmetry structure. However, it should be noted that the structure still exhibits 8 imaginary modes. The cause of the problem above is not difficult to identify. By inspecting the structures of the two low-energy structures as marked by two downward arrows in Fig. 2(a), it can be seen that the imaginary modes lead to rotations of the iodine octahedra as they are close to the instantaneously stable structure of $\alpha$-CsPbI$_3$. In the perovskite structure, the octahedron rotation exists universally in bulk and low-dimensional perovskites without defects.[43,44] The structures in Fig. 2(a) on the left and right are taken from the end structures of step 1 and step 6 in Fig. 2(b), respectively. Each individual low-energy structure from the calculations in Fig. 2(b) could be associated with one or more octahedron rotations. These distorted structures are randomly distributed in the configuration space. Falling into any one of them purely relies on the initial perturbation. In the examined case above, the perturbation is realized by following one of the imaginary modes. It should be emphasized that none of these low-energy structures are metastable. They are just reached by zero-temperature structural relaxation. Even worse, the possibility of being trapped in such low-energy structures could increase with larger supercells. Apparently, defect properties from such calculations are purely random and hard to explain physically. MD-Based Ensemble Average for Defect Calculations. We present a straightforward approach to sampling the configuration space by using MD simulations. To study the defect tolerance, one is often concerned with the defect transition level as defined by $$ \epsilon(q/q{^\prime})=\frac{E_{\scriptscriptstyle{\rm D}}^q-E_{\scriptscriptstyle{\rm D}}^{q^{\prime}}}{q^{\prime}-q}-E_{\rm VBM},~~ \tag {1} $$ which indicates whether the defect is deep or shallow. $E_{\scriptscriptstyle{\rm D}}^q$ and $E_{\scriptscriptstyle{\rm D}}^{q^{\prime}}$ are the total energies of the supercells containing a point defect of charge $q$ and $q^{\prime}$, respectively. $E_{\rm VBM}$ is the energy of valence band maximum (VBM) of the host material. In common defect calculations on $\alpha$-CsPbI$_3$,[45] the total energies are obtained by optimizing the supercell structures at zero temperature, which suffer from the randomness problem as illustrated in Fig. 2. Figure 3(a) shows the MD simulation results for the $v_{_{\scriptstyle \rm I}}$ defect demonstrating an effective way of sampling the configuration space. The left panel shows the total energies as a function of simulation time for the neutral and 1$^+$ charge states. The difference between the two average energies over the last 20 ps gives rise to the transition level $\epsilon(0/+)$ with respect to the VBM. Here, $\epsilon(0/+) = 1.68$ eV indicates that $v_{_{\scriptstyle \rm I}}$ is a shallow donor. Figure 3(b) shows the results for an acceptor defect, I-on-Cs antisite (I$_{\rm Cs}$). In CH$_3$NH$_3$PbI$_3$, the similar defect I$_{\rm MA}$ (MA = CH$_3$NH$_3$) was found to have a deep $\epsilon(0/-)$ transition level due to the formation of I trimer,[18] which is a stable species commonly found in electrochemistry.[46] Our results show that I$_{\rm Cs}$ also has a deep $\epsilon(0/-)$ level inside the band gap. By inspecting the atomic structure, I trimer is also observed, which will be further discussed below.
cpl-39-4-046101-fig3.png
Fig. 3. (a) Total energies of the $v_{_{\scriptstyle \rm I}}$ donor defect in $\alpha$-CsPbI$_3$ as a function of the MD simulation time for the neutral and 1$^+$ charge states. (b) Total energies of the I$_{\rm Cs}$ acceptor defect as a function of the MD simulation time for the neutral and 1$^-$ charge states. Dashed lines show the averaged values over the last 20-ps simulations. The right panels show the histogram statistics of the total energy from each step of the MD simulations, where the smooth solid lines show the Gaussian fitting to the histograms. The arrows mark the Gaussian peaks.
The right panels in Fig. 3 show the histograms of the total energies. By fitting the histogram with a Gaussian function, we obtain the peak positions and standard deviations. The peak positions from the Gaussian fitting typically deviate from the direct ensemble averages in the corresponding left panels in Fig. 3 by 0.01 eV or smaller. However, in several simulations we observe larger deviations up to 0.07 eV. Another indicator of the reliability of simulation is the standard deviation ($\sigma$) obtained from the Gaussian fitting. In principle, $\sigma$ should be the same as that from the simulation of the perfect supercell if the supercell size is large enough. Our results show that for the perfect supercell $\sigma = 0.72$ eV, while for the defects, the calculated $\sigma$ deviates from this value by $\pm 0.1$ eV, which provides an estimate on the error bar of the calculated transition levels. Defect Transition Levels. For the material to be defect tolerant, the key indicator is that there exist no efficient recombination centers. In general, the $\epsilon(0/+)$ and $\epsilon(0/-)$ transition levels need to be considered, while transitions between higher charge states are carrier traps only.[47] Deep $\epsilon(0/-)$ and $\epsilon(0/+)$ levels serve as the most efficient non-radiative recombination centers. In a recent work,[48] it was noted that I$_{\rm Pb}$ in the negative charge state could still have a sizable capture coefficient for electrons. We do not consider such possibility in this work. Figure 4(a) shows the defect transition levels for the intrinsic defects, where the levels shown by solid bars correspond to the normal transitions of charge states, while those by dashed bars (not the shaded bars) correspond to the abnormal transitions. Here, by “normal” transitions, we mean that a cation vacancy or an anion interstitial is an acceptor and a cation interstitial or an anion vacancy is a donor, while the “abnormal” transitions are opposite to the normal transitions. To evaluate the normal transition for an antisite defect, it can be considered as a complex of an interstitial and a vacancy. Comparing with previous calculations employing the standard approach with the static high-symmetry structure,[45] we obtain different results. In the previous study, a half of the defects exhibit transition levels about 1.5 eV above the conduction band minimum (CBM) or below the VBM. Also, I$_{\rm Pb}$ was predicted to be a deep acceptor and Pb$_{\rm Cs}$ with a transition level 0.1 eV below the CBM was considered to be a deep donor. In contrast, here we find four deep-level defects (more than 0.1 eV away from the band edges) having the normal transitions (four thick solid bars in the gap) and four deep-level defects having the abnormal transitions (four dashed bars in the gap), while all other transition levels are shallow. All the deep levels can be understood based on the strong covalency in lead halide perovskites,[18] namely, they always involve the formation of either a Pb dimer or an I trimer. Figure 4(b) shows the averaged structures of the deep-level defects over 20-ps simulations, where the Pb dimers and I trimers can be clearly seen. Understanding on the Deep Levels. A general understanding on the condition for forming the Pb dimers and I trimers and the associated deep transition levels are derived from the above results. If the defect in a particular charge state has two or more electrons in the defect band or conduction band, a Pb dimer will form. Similarly, if the defect in a particular charge state has two or more holes in the defect band or valence band, an I trimer will form. The deep-level transitions are always between two charge states with and without the Pb dimer or I trimer formed. The formation of a Pb dimer or an I trimer is associated with electronic energy gain,[18] which is to be competed with the elastic energy loss due to the structural distortion.[33] When there are two extra electrons (or holes), the electronic energy gain is doubled so that it can overcome the elastic energy loss. This is the case for all the eight deep transition levels in Fig. 4(a). For Pb$_{\rm I}$ and I$_{\rm Pb}$, their $\epsilon(0/+)$ and $\epsilon(0/-)$ transitions are between the states having two and three electrons (or holes). Therefore, both charge states involved in the transition have the Pb dimer or I trimer formed, explaining why the $\epsilon(0/+)$ and $\epsilon(0/-)$ transitions of Pb$_{\rm I}$ and I$_{\rm Pb}$ are shallow. For Cs$_i$ and $v_{_{\scriptstyle \rm Cs}}$, the formations of Pb dimer or I trimer are difficult, even in the abnormal charge state, because they need to overcome large energy barriers to escape from a perfect [PbI$_6$] octahedron. Therefore, these two defects always exhibit shallow transition levels. It should be emphasized that having the deep levels does not mean that the material cannot be defect tolerant. It has been proposed that,[18] in order for the deep levels serving as recombination centers, the Fermi level of the material needs to be close to the band edges. It is known that the halide perovskites usually have negligible amount of equilibrium carriers, rendering the Fermi level always located at the mid-gap. In such a case, the charge states associated with deep-level transitions are not thermodynamically stable. To check this condition, we also calculate the $\epsilon$($-/-$2) and $\epsilon$($+/+$2) transitions for $v_{_{\scriptstyle \rm Pb}}$, Pb$_i$, Cs$_{\rm I}$ and I$_{\rm Cs}$ defects, which are all shallow transitions as shown in Fig. 4(a). These results are consistent with the previous study.[18] Discussion on the MD-Based Approach. We calculated the defect formation energy in neutral charge state by $$ E_{\scriptscriptstyle{\rm D}}-E_0+\sum_i n_i(E_i+\mu_i),~~ \tag {2} $$ where the total energies of the defect-containing supercell ($E_{\scriptscriptstyle{\rm D}}$) and defect-free supercell ($E_0$) were calculated by ensemble averages from MD simulations. The total energies of elemental phases ($E_i$) were determined by calculations at 0 K and amended by $\frac{3}{2}k_{\scriptscriptstyle{\rm B}}T$ for each atom (see the following discussion for details). The chemical potentials ($\mu_i$) here were measured with respect to the corresponding $E_i$. The $n_i$ in Eq. (2) is the number of incorporated or removed atoms when forming the defect. The allowed region for $\alpha$-CsPbI$_3$ to be thermodynamically stable with respect to the binary compounds is shown in Fig. 4(c). The calculated defect formation energies of all the intrinsic defects as a function of chemical potentials are shown in Fig. 4(d). Under the I-rich condition (point B), all the nominal acceptor defects ($v_{_{\scriptstyle \rm Cs}}$, $v_{_{\scriptstyle \rm Pb}}$, I$_{i}$, Cs$_{\rm Pb}$, I$_{\rm Cs}$ and I$_{\rm Pb}$) exhibit low formation energies below 1 eV. Therefore, this condition should be avoided. Under the I-poor condition (point A), $v_{_{\scriptstyle \rm I}}$ and Pb$_{i}$ have relatively low formation energies, while the other four nominal donor defects have large formation energies above 1 eV. By varying the chemical potential, it is possible that all the defects have formation energy higher than 1 eV.
cpl-39-4-046101-fig4.png
Fig. 4. (a) Defect transition levels in $\alpha$-CsPbI$_3$ calculated by MD-based ensemble average. Red bars show the acceptor levels and blue bars show the donor levels. The four deep “normal” transition levels are shown by thick solid bars. The four deep “abnormal” transition levels are shown by dashed bars. The $\epsilon$($-/-$2) and $\epsilon$($+/+$2) transitions are shown by shaded bars. (b) Atomic structures of the deep-level defects obtained by averaging the trajectories from the last 20-ps MD simulations. The chosen charge states for these defects exhibit either a Pb dimer or an I trimer. For clarity, Cs atoms are not shown here except for the cases of Cs$_{\rm I}$ and Cs$_{\rm Pb}$ defects. (c) Determination of chemical potentials for calculating defect formation energies in $\alpha$-CsPbI$_3$. The gray area shows the allowed chemical potential region for $\alpha$-CsPbI$_3$ to be thermodynamically stable. (d) Defect formation energies of intrinsic defects in neutral charge state at A and B points in the chemical potential space corresponding to I-poor and I-rich conditions, respectively.
As we employed MD simulations to obtain defect properties at a finite temperature, it is necessary to describe the physical meaning of the quantities that were calculated above and the approximations that were made. The key quantity here is the defect formation energy, from which the defect transition levels are calculated. We consider a system as shown in Fig. 5 illustrating the case of doping a binary host material with an element A, where the element B in the host is to be replaced. The reservoirs A and B contain A and B atoms with chemical potentials of $\mu_{\scriptscriptstyle{\rm A}}$ and $\mu_{\scriptscriptstyle{\rm B}}$, respectively. The defects could be charged and carry a charge $q$, where the electrons are retrieved from or released to the reservoir of electrons, i.e., the Fermi level $E_{\scriptscriptstyle{\rm F}}$.
cpl-39-4-046101-fig5.png
Fig. 5. Schematic showing the thermodynamics of defect formation. The case of doping a binary host material is taken as an example, where an element A (gray balls) is to substitute for the element B in the host (brown balls). Three A-on-B defects are shown in the illustration of the defect system so that three holes (white balls) appear in the reservoir A and three extra B atoms appear in the reservoir B. Six little holes appear in the electron reservoir and represent the case of $q=2e$.
Consider the defect formation process as a chemical reaction with the reactants and products of the reaction represented by the upper and lower panels of Fig. 5, respectively. To reach thermal equilibrium, it is required that $$ G_{\rm host}+N_{\scriptscriptstyle{\rm D}}\mu_{\scriptscriptstyle{\rm A}}-N_{\scriptscriptstyle{\rm D}}qE_{\scriptscriptstyle{\rm F}}=G_{\rm defect}+N_{\scriptscriptstyle{\rm D}} \mu_{\scriptscriptstyle{\rm B}},~~ \tag {3} $$ where $G_{\rm host}$ and $G_{\rm defect}$ are the Gibbs free energies of the host and defect systems, respectively, and $N_{\scriptscriptstyle{\rm D}}$ is the total number of A-on-B defects formed during the reaction. Using $G=H-TS$ and separating the entropy into configurational entropy ($S_{\rm c}$) and vibrational entropy ($S_{\rm v}$), we have $$ -k_{\scriptscriptstyle{\rm B}}T{\ln}\Big(\frac{N_{\scriptscriptstyle{\rm D}}}{N_{\rm site}}\Big)=\Delta H-T\Delta S_{\rm v}-\mu_{\scriptscriptstyle{\rm A}}+\mu_{\scriptscriptstyle{\rm B}}+qE_{\scriptscriptstyle{\rm F}},~~ \tag {4} $$ where $\Delta H$ and $\Delta S_{\rm v}$ are the changes in enthalpy and vibrational entropy per defect, respectively, and $N_{\rm site}$ is the total number of lattice sites that can accommodate the defects. The relation $S_{\rm c} \approx k_{\scriptscriptstyle{\rm B}}N_{\scriptscriptstyle{\rm D}}{\ln}(N_{\rm site}/N_{\scriptscriptstyle{\rm D}})$ has been used under the approximation that $N_{\scriptscriptstyle{\rm D}} \ll N_{\rm site}$, i.e., the dilute limit of defect concentration. Now, defining the right hand side of Eq. (4) as the defect formation energy ($\Delta E_{\rm f}$), we have $$ N_{\scriptscriptstyle{\rm D}}=N_{\rm site}{\exp}\Big(-\frac{\Delta E_{\rm f}}{k_{\scriptscriptstyle{\rm B}}T}\Big).~~ \tag {5} $$ This is the relation generally used to calculate defect concentration in a material. To be general for a more complex defect involving the exchange of more atomic species, the chemical potential terms can be grouped into $\sum_{i}n_i\mu_i$, where $i$ is the index for the atomic species and $n_i$ is negative (positive) if the atomic species is to be doped into (removed from) the host. Then, the defect formation energy is given by $$ \Delta E_{\rm f} \equiv \Delta H-T\Delta S_{\rm v}+\Sigma_{i}n_i\mu_i+qE_{\scriptscriptstyle{\rm F}}.~~ \tag {6} $$ The defect formation energy defined in this way is not simply the change in Gibbs free energy between the defect and host systems because the configurational entropy due to the presence of defects has been considered separately. Also, note that in practical calculations the chemical potentials $\mu_i$ are measured with respect to the corresponding $E_i$ and written as $E_i+\mu_i$ as in Eq. (2), while the Fermi level $E_{\scriptscriptstyle{\rm F}}$ is measured with respect to $E_{\rm VBM}$ and written as $E_{\rm VBM}+E_{\scriptscriptstyle{\rm F}}$. In common calculations on $\Delta E_{\rm f}$, the contributions from the $PV$ term and the vibrational entropy are often ignored. We have followed these approximations in this work. In such a case, $H$ is replaced by the internal energy $U$. Moreover, as the movements of nuclei are treated classically, the contribution from kinetic energies is canceled according to the equipartition theorem so that the internal energy $U$ is further replaced by the potential energy $E$ (or total energy in the common context of DFT calculations). As a result of the treatment mentioned above, $\Delta H$ in Eq. (6) is replaced by $E_{\scriptscriptstyle{\rm D}}-E_0$ in Eq. (2). As $E_{\scriptscriptstyle{\rm D}}$ and $E_0$ are both evaluated at finite temperature, the $E_i$ terms in chemical potentials are also supposed to be obtained at the same temperature by MD simulations. However, because the $E_i$ terms are all from calculations on single atomic species and the averaged potential energy equals the averaged kinetic energy, it is valid to simply amend $E_i$ by $\frac{3}{2}k_{\scriptscriptstyle{\rm B}}T$ because the temperature-dependent part of the internal energy is $3k_{\scriptscriptstyle{\rm B}}T$. The focus of this paper is to discuss the issue and solution on defect study of statically unstable materials. The proposed approach was to demonstrate that the MD-based calculations are able to remove the randomness in defect structures and their associated total energies. Nevertheless, it is still worthwhile to discuss the possible errors that could exist in our results. The calculation of band edge positions of lead halide perovskites suffers from both the well-known DFT errors and the strong spin-orbit coupling (SOC) effect.[49,50] Hybrid functional calculations with SOC employing a sufficiently large supercell and full structure relaxations would be desired for the study of lead halide perovskites.[33] However, based on our previous studies[18,33] and the results in this work, there is an important conclusion which is consistent, i.e., the appearance of deep-level defects is always accompanied by the formation of a Pb dimer or an I trimer suggesting that the results from semilocal functionals are still worth of reference. Nevertheless, whether the Pb dimers and I trimers observed in our PBEsol calculations can exist in high-level calculations will need further examination. Before concluding, we comment that the method proposed above should be generally applicable to the study on defect properties in other statically unstable semiconductor materials. One important application would be on the study of thermoelectric materials, which often exhibit rich phase transitions at elevated temperatures.[51–53] The high-temperature phases of such materials are mostly statically unstable similar to $\alpha$-CsPbI$_3$. Even more complicated is the possibility of liquid-like behavior of cations in some thermoelectric materials, e.g., cuprous chalcogenides.[54] Even though defect properties are of central importance for the electronic properties of thermoelectric materials, commonly used first-principles approaches cannot reliably describe such liquid-like materials. Our approach is expected to be straightforwardly applied to such systems. In summary, by MD simulation and phonon calculation we show the static instability of $\alpha$-CsPbI$_3$ with the $Pm\bar{3}m$ symmetry. We demonstrate the band gap issue and, more seriously, the randomness issue in the calculation of defect transition levels incurred by the static instability. We propose to use MD simulation to sample the configuration space in such defect calculations. This approach is able to obtain sufficiently accurate defect transition levels for evaluating the defect tolerance of this material. The errors in the results can be estimated by the standard deviation of the total energies from the ensemble average through a Gaussian fitting. By the MD-based approach, we observe eight deep defect transition levels, which are always associated with the formation of Pb dimers or I trimers. It is found that two electrons in the conduction band will induce the formation of a Pb dimer and two holes in the valence band will induce the formation of an I trimer. According to the previous proposal,[18] such deep-level transitions can always be avoided if the Fermi level is located at the mid-gap, which is actually the case for most lead halide perovskites. This MD-based approach is expected to be straightforwardly generalized to other statically unstable materials. Acknowledgments. Y.-Y. Sun acknowledges supports from Shanghai International Cooperation Project (Grant No. 20520760900) and the Major Science and Technology Programs of Yunnan (Grant No. 202002AB080001-1). C. Ming acknowledges supports from the Natural Science Foundation of Shanghai (Grant No. 19ZR1421800). S. Zhang acknowledges supports from the US National Science Foundation (Grant No. CBET-1510948). Part of the calculations was conducted at Stampede supercomputer resources at TACC made available by Extreme Science and Engineering Discovery Environment (XSEDE) through Allocation (Grant No. TG-DMR180114), which was supported by US National Science Foundation (Grant No. ACI-1548562). References Inorganic caesium lead iodide perovskite solar cellsQuantum dot–induced phase stabilization of α-CsPbI 3 perovskite for high-efficiency photovoltaicsAll-Inorganic Perovskite Solar CellsBifunctional Stabilization of All-Inorganic α-CsPbI 3 Perovskite for 17% Efficiency PhotovoltaicsFast Anion-Exchange in Highly Luminescent Nanocrystals of Cesium Lead Halide Perovskites (CsPbX 3 , X = Cl, Br, I)Quantum Dot Light-Emitting Diodes Based on Inorganic Perovskite Cesium Lead Halides (CsPbX 3 )Enhancing the Brightness of Cesium Lead Halide Perovskite Nanocrystal Based Green Light-Emitting Devices through the Interface Engineering with Perfluorinated IonomerNanocrystals of Cesium Lead Halide Perovskites (CsPbX 3 , X = Cl, Br, and I): Novel Optoelectronic Materials Showing Bright Emission with Wide Color GamutSurface Engineering of Ambient-Air-Processed Cesium Lead Triiodide Layers for Efficient Solar CellsChemically Stable Black Phase CsPbI 3 Inorganic Perovskites for High‐Efficiency PhotovoltaicsA review on the stability of inorganic metal halide perovskites: challenges and opportunities for stable solar cellsHigh-temperature structural evolution of caesium and rubidium triiodoplumbatesSolvent-controlled growth of inorganic perovskite films in dry environment for efficient and stable solar cellsThermodynamically stabilized β-CsPbI 3 –based perovskite solar cells with efficiencies >18%Thermodynamically Stable Orthorhombic γ-CsPbI 3 Thin Films for High-Performance PhotovoltaicsChlorine doping for black γ-CsPbI3 solar cells with stabilized efficiency beyond 16%Stabilizing γ‐CsPbI 3 Perovskite via Phenylethylammonium for Efficient Solar Cells with Open‐Circuit Voltage over 1.3 VStrong Covalency-Induced Recombination Centers in Perovskite Solar Cell Material CH 3 NH 3 PbI 3Defect Tolerance in Methylammonium Lead Triiodide PerovskiteUnusual defect physics in CH 3 NH 3 PbI 3 perovskite solar cell absorberIdentifying defect-tolerant semiconductors with high minority-carrier lifetimes: beyond hybrid lead halide perovskitesInstilling defect tolerance in new compoundsIodine chemistry determines the defect tolerance of lead-halide perovskitesPoint defect engineering in thin-film solar cellsHigh Defect Tolerance in Lead Halide Perovskite CsPbBr 3When defects become ‘dynamic’: halide perovskites: a new window on materials?Breakdown of the Static Picture of Defect Energetics in Halide Perovskites: The Case of the Br Vacancy in CsPbBr 3Dynamic shortening of disorder potentials in anharmonic halide perovskitesAssessment of dynamic structural instabilities across 24 cubic inorganic halide perovskitesIntrinsic Point Defects in Inorganic Cesium Lead Iodide Perovskite CsPbI 3All-inorganic halide perovskites as candidates for efficient solar cellsBenign Deep-Level Defects in Cesium Lead Iodine PerovskiteDefect tolerance in CsPbI 3 : reconstruction of the potential energy landscape and band degeneracy in spin–orbit couplingEfficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis setProjector augmented-wave methodFrom ultrasoft pseudopotentials to the projector augmented-wave methodRestoring the Density-Gradient Expansion for Exchange in Solids and SurfacesA molecular dynamics method for simulations in the canonical ensembleCanonical dynamics: Equilibrium phase-space distributionsStructural and electronic properties of hybrid perovskites for high-efficiency thin-film photovoltaics from first-principlesFirst principles phonon calculations in materials scienceTemperature dependent effective potential method for accurate free energy calculations of solidsThe classification of tilted octahedra in perovskitesOctahedron rotation evolution in 2D perovskites and its impact on optoelectronic properties: the case of Ba–Zr–S chalcogenidesIntrinsic point defects in inorganic perovskite CsPbI 3 from first-principles predictionInverting the Triiodide Formation Reaction by the Synergy between Strong Electrolyte Solvation and Cathode Adsorption for Lithium–Oxygen BatteriesDefect tolerance in chalcogenide perovskite photovoltaic material BaZrS3Correctly Assessing Defect Tolerance in Halide PerovskitesDensity Functional Calculations of Native Defects in CH 3 NH 3 PbI 3 : Effects of Spin–Orbit Coupling and Self-Interaction ErrorDiscovering lead-free perovskite solar materials with a split-anion approachThermal Conductivity during Phase TransitionsNanoscale Behavior and Manipulation of the Phase Transition in Single-Crystal Cu 2 SeGeTe ThermoelectricsRecent Advances in Liquid‐Like Thermoelectric Materials
[1] Eperon G E, Paternò G M, Sutton R J, Zampetti A, Haghighirad A A, Cacialli F, and Snaith H J 2015 J. Mater. Chem. A 3 19688
[2] Swarnkar A, Marshall A R, Sanehira E M, Chernomordik B D, Moore D T, Christians J A, Chakrabarti T, and Luther J M 2016 Science 354 92
[3] Liang J, Wang C, Wang Y, Xu Z, Lu Z, Ma Y, Zhu H, Hu Y, Xiao C, Yi X et al. 2016 J. Am. Chem. Soc. 138 15829
[4] Wang Y, Zhang T, Kan M, and Zhao Y 2018 J. Am. Chem. Soc. 140 12345
[5] Nedelcu G, Protesescu L, Yakunin S, Bodnarchuk M I, Grotevent M J, and Kovalenko M V 2015 Nano Lett. 15 5635
[6] Song J, Li J, Li X, Xu L, Dong Y, and Zeng H 2015 Adv. Mater. 27 7162
[7] Zhang X, Lin H, Huang H, Reckmeier C, Zhang Y, Choy W C H, and Rogach A L 2016 Nano Lett. 16 1415
[8] Protesescu L, Yakunin S, Bodnarchuk M I, Krieg F, Caputo R, Hendon C H, Yang R X, Walsh A, and Kovalenko M V 2015 Nano Lett. 15 3692
[9] Yoon S M, Min H, Kim J B, Kim G, Lee K S, and Seok S I 2021 Joule 5 183
[10] Wang Y, Chen Y, Zhang T, Wang X, and Zhao Y 2020 Adv. Mater. 32 2001025
[11] Xiang W, Liu S F, and Tress W 2021 Energy & Environ. Sci. 14 2090
[12] Trots D M and Myagkota S V 2008 J. Phys. Chem. Solids 69 2520
[13] Wang P, Zhang X, Zhou Y, Jiang Q, Ye Q, Chu Z, Li X, Yang X, Yin Z, and You J 2018 Nat. Commun. 9 2225
[14] Wang Y, Dar M I, Ono L K, Zhang T, Kan M, Li Y, Zhang L, Wang X, Yang Y, Gao X et al. 2019 Science 365 591
[15] Zhao B, Jin S F, Huang S, Liu N, Ma J Y, Xue D J, Han Q, Ding J, Ge Q Q, Feng Y et al. 2018 J. Am. Chem. Soc. 140 11716
[16] Wang K, Jin Z, Liang L, Bian H, Wang H, Feng J, Wang Q, and Liu S F 2019 Nano Energy 58 175
[17] Ye Q, Ma F, Zhao Y, Yu S, Chu Z, Gao P, Zhang X, and You J 2020 Small 16 2005246
[18] Agiorgousis M L, Sun Y Y, Zeng H, and Zhang S 2014 J. Am. Chem. Soc. 136 14570
[19] Steirer K X, Schulz P, Teeter G, Stevanovic V, Yang M, Zhu K, and Berry J J 2016 ACS Energy Lett. 1 360
[20] Yin W J, Shi T, and Yan Y 2014 Appl. Phys. Lett. 104 063903
[21] Brandt R E, Stevanović V, Ginley D S, and Buonassisi T 2015 MRS Commun. 5 265
[22] Walsh A and Zunger A 2017 Nat. Mater. 16 964
[23] Meggiolaro D, Motti S G, Mosconi E, Barker A J, Ball J, Perini C A R, Deschler F, Petrozza A, and De Angelis F 2018 Energy & Environ. Sci. 11 702
[24] Park J S, Kim S, Xie Z, and Walsh A 2018 Nat. Rev. Mater. 3 194
[25] Kang J and Wang L W 2017 J. Phys. Chem. Lett. 8 489
[26] Rakita Y, Lubomirsky I, and Cahen D 2019 Mater. Horiz. 6 1297
[27] Cohen A V, Egger D A, Rappe A M, and Kronik L 2019 J. Phys. Chem. Lett. 10 4490
[28] Gehrmann C and Egger D A 2019 Nat. Commun. 10 3141
[29] Yang R X, Skelton J M, da Silva E L, Frost J M, and Walsh A 2020 J. Chem. Phys. 152 024703
[30] Huang Y, Yin W J, and He Y 2018 J. Phys. Chem. C 122 1345
[31] Zhang X, Turiansky M E, and Van de Walle C G 2021 Cell Rep. Phys. Sci. 2 100604
[32] Zhang J, Zhong Y, and Li G 2021 J. Phys. Chem. C 125 27016
[33] Ming C, Wang H, West D, Zhang S, and Sun Y Y 2022 J. Mater. Chem. A 10 3018
[34] Kresse G and Furthmüller J 1996 Comput. Mater. Sci. 6 15
[35] Blöchl P E 1994 Phys. Rev. B 50 17953
[36] Kresse G and Joubert D 1999 Phys. Rev. B 59 1758
[37] Perdew J P, Ruzsinszky A, Csonka G I, Vydrov O A, Scuseria G E, Constantin L A, Zhou X, and Burke K 2008 Phys. Rev. Lett. 100 136406
[38] Nosé S 1984 Mol. Phys. 52 255
[39] Hoover W G 1985 Phys. Rev. A 31 1695
[40] Brivio F, Walker A B, and Walsh A 2013 APL Mater. 1 042111
[41] Togo A and Tanaka I 2015 Scr. Mater. 108 1
[42] Hellman O, Steneteg P, Abrikosov I A, and Simak S I 2013 Phys. Rev. B 87 104111
[43] Glazer A M 1972 Acta Crystallogr. Sect. B 28 3384
[44] Ming C, Yang K, Zeng H, Zhang S, and Sun Y Y 2020 Mater. Horiz. 7 2985
[45] Li Y, Zhang C, Zhang X, Huang D, Shen Q, Cheng Y, and Huang W 2017 Appl. Phys. Lett. 111 162106
[46] Zhang X P, Li Y N, Sun Y Y, and Zhang T 2019 Angew. Chem. Int. Ed. 58 18394
[47] Wu X, Gao W, Chai J, Ming C, Chen M, Zeng H, Zhang P, Zhang S, and Sun Y Y 2021 Sci. Chin. Mater. 64 2976
[48] Zhang X, Turiansky M E, and Van de Walle C G 2020 J. Phys. Chem. C 124 6022
[49] Du M H 2015 J. Phys. Chem. Lett. 6 1461
[50] Sun Y Y, Shi J, Lian J, Gao W, Agiorgousis M L, Zhang P, and Zhang S 2016 Nanoscale 8 6284
[51] Chen H Y, Yue Z, Ren D, Zeng H, Wei T, Zhao K, Yang R, Qiu P, Chen L, and Shi X 2019 Adv. Mater. 31 1806518
[52] Chen L, Liu J, Jiang C, Zhao K, Chen H, Shi X, Chen L, Sun C, Zhang S, Wang Y et al. 2019 Adv. Mater. 31 1804919
[53] Zhang X, Bu Z, Lin S, Chen Z, Li W, and Pei Y 2020 Joule 4 986
[54] Zhao K, Qiu P, Shi X, and Chen L 2020 Adv. Funct. Mater. 30 1903867