Effect of vacancy concentration on the lattice thermal conductivity of CH3NH3PbI3: a molecular dynamics study

Hybrid halide perovskites are drawing great interest for photovoltaic and thermoelectric applications, but the relationship of thermal conductivities with vacancy defects remains unresolved. Here, we present a systematic investigation of the thermal conductivity of perfect and defective CH3NH3PbI3, performed using classical molecular dynamics with an ab initio-derived force field. We calculate the lattice thermal conductivity of perfect CH3NH3PbI3 as the temperature increases from 300 K to 420 K, confirming a good agreement with the previous theoretical and experimental data. Our calculations reveal that the thermal conductivities of defective systems at 330 K, containing vacancy defects such as VMA, VPb and VI, decrease overall with some slight rises, as the vacancy concentration increases from 0 to 1%. We show that such vacancies act as phonon scattering centers, thereby reducing the thermal conductivity. Moreover, we determine the elastic moduli and sound velocities of the defective systems, revealing that their slower sound speed is responsible for the lower thermal conductivity. These results could be useful for developing hybrid halide perovskite-based solar cells and thermoelectric devices with high performance.


Introduction
Halide perovskites (HPs) have emerged as exciting new materials for photovoltaic applications owing to their unique properties of strong light absorption, high carrier mobility and relatively simple manufacturing process. 1 In fact, the power conversion efficiency of perovskite solar cells (PSCs), which use HPs as a light absorber, has rapidly increased from an initial value of 3.8% in 2009 (ref. 1) to a certied value of 25.5% in 2020. 2, 3 In addition to photovoltaics, HPs draw great interest in applications such as thermoelectrics, 4-6 light emitting diodes, 7 photodetectors 8 and lasers. 9 For these applications, the lattice thermal conductivity of the HPs is the most important quantity to consider for obtaining the best performance under the conditions these devices are normally operated. 10 Moreover, as the temperature of devices using HPs could increase rapidly during operation, HPs can be decomposed into their constituent components, leading to degradation of device performance. Meanwhile, a temperature difference between the sandwiched layers of HP-related devices can cause mechanical stress, resulting in device destruction. Therefore, it is necessary to gain insight into the thermal conductivity of HPs for guaranteeing the normal operation of device and designing thermoelectric materials.
The chemical formula of HPs is written as ABX 3 , where A is the monovalent organic or inorganic cation, such as CH 3 NH 3 + (or MA + ), CH(NH 2 ) 2 + (or FA + ) or Cs + ; B is the divalent metallic cation, such as Pb 2+ , Sn 2+ or Ge 2+ , and X is the halogen anion such as I À , Br À or Cl À . Among them, a prototype HP material used in most applications is methylammonium lead iodide CH 3 NH 3 PbI 3 (or MAPbI 3 ). In experiments, MAPbI 3 was found to have an ultralow thermal conductivity k of 0.3-0.5 W m À1 K À1 at room temperature, 5,6,[11][12][13][14] which is useful for thermoelectric and thermal barrier applications. In line with experiments, theoretical calculations with density functional theory (DFT) [15][16][17] and/or classical molecular dynamics (MD) methods [17][18][19][20] have also been performed, revealing that the k value for MAPbI 3 is in the range 0.2 to 1.8 W m À1 K À1 , depending on the crystalline phase and calculation method. It should be noted that the calculations mainly focused on uncovering the mechanism of its ultralow k, demonstrating that the MA + cation plays a critical role in suppressing k, as conrmed in infrared and Raman spectroscopy experiments. [21][22][23] In general, HPs are synthesized from precursor solutions under ambient conditions by simple and low-cost chemical processing, so that various defects are inevitably created in the resultant HP solid. 24,25 Such defects, including vacancies, antisites and interstitials, were found to affect the optoelectronic and thermal properties of HPs. Therefore, numerous investigations have been performed on the unusual defect physics and chemistry of HPs, demonstrating that they are defect-tolerant even though the defect concentration is as high as 10 17 -10 20 cm À3 due to their low formation energies. [26][27][28] It was known that defects usually create deep levels near the mid gap, acting as charge carrier traps and accordingly deteriorating the performance in the Si-or III-V semiconductor-based solar cells. In PSCs, however, defects fortunately act as shallow donors or acceptors, 29,30 thus not severely degrading the optoelectronic properties and charge carrier mobilities of HPs. 31,32 In spite of many investigations on the electronic transport properties of MAPbI 3 , the effect of vacancy defects on its lattice thermal and mechanical properties remains unexplored.
In this work, we investigate the effects of various vacancy defects with different concentrations on the thermal conductivity of MAPbI 3 by performing systematic MD simulations. From the simulation data, we determine the contributions of individual MA + , Pb 2+ and I À ions to the change of heat transport properties. We estimate the elastic constants and sound speeds in perfect and defective systems, revealing their relation with the thermal conductivity. Section 2 describes the MD method used to calculate the thermal conductivity, simulation models with diverse vacancies, and the employed force eld. Section 3 presents the simulation results and discussion, comparing them with the available experimental and theoretical data. Finally, conclusions are provided in Section 4.

Molecular dynamics methods
Thermal conductivity of a crystalline solid is expressed as a second-order tensor, which can be calculated using different kinds of MD methods, such as equilibrium MD (EMD) 33 and approach-to-equilibrium MD (AEMD). 34 The EMD approach is based on the uctuation-dissipation theory at an equilibrium state of system, where the thermal conductivity can be calculated by integrating the heat ux autocorrelation function (HFACF) obtained with a sufficiently long-term run and using the following Green-Kubo relation, 33 where k B is the Boltzmann constant, V is the volume of system, and a and b are the Cartesian coordinates. The second-order tensor f ab (t) ¼ hq a (t) 5 q b (0)i/hq a (0) 5 q b (0)i is the HFACF, and q a (t) is the heat ux written as, where r i is the position vector of the ith atom, and H i is the Hamilton function of the system. The AEMD method yields k by solving the heat ux equation along the ux direction z, 34 where k ¼ k/rC V is the thermal diffusivity with mass density r and volumetric heat capacity C V . For the rst step, the two subregions are equilibrated with a periodic boundary condition at different temperatures as depicted in Fig. 1. Next, an NVE simulation is performed for the entire system, gradually decreasing the difference of average temperatures between the two sub-regions. Then, the temperature difference is tted to the exponential function, giving the thermal diffusivity k (see the ESI for details †).

Modeling and computational details
The MAPbI 3 crystal was found to be stabilized in the orthorhombic phase below 160 K, the tetragonal phase between 160 and 330 K, and the cubic phase over 330 K. 35,36 Since the MAPbI 3 crystal is heated for calculating the k value of MAPbI 3 above room temperature, we adopted the tetragonal phase with I4/ mcm space group at temperatures below 330 K and a pseudocubic phase with Pm 3m space group at higher temperatures. For the EMD simulations, we created 8 Â 8 Â 8 supercells with 6144 atoms for the pseudocubic phase (Fig. 2) and 5 Â 5 Â 5 supercells with 6000 atoms for the tetragonal phase of perfect MAPbI 3 (Fig. S1 †). To create defective MAPbI 3 structures, we remove a certain number of molecular moieties or atoms in the 8 Â 8 Â 8 supercells. With respect to the defect concentration, Walsh et al. 26 pointed out that the concentration of vacancy defects, such as the MA vacancy V MA , Pb vacancy V Pb and I vacancy V I , should exceed 0.4% in cubic MAPbI 3 at room temperature, based on their DFT calculations. Accordingly, the vacancy concentrations in the defective MAPbI 3 were allowed to Fig. 1 Temperature profile during the AEMD simulation. After enough equilibration at temperature (T 1 + T 2 )/2, the simulation box is split into two regions, i.e., hot and cold regions, which are equilibrated at T 1 and T 2 . Here, T 1 > T 2 , and hT 1 (t 1 )i and hT 2 (t 2 )i are the average temperatures of the regions with t 1 < t 2 . L z is the length of the simulation box in the zdirection.
vary from 0.0 to 1.0% with an interval of about 0.2% in this work. The simulation cells for defective systems with these vacancy concentrations were created by randomly removing MA + and Pb 2+ cations with numbers from 1 to 5 with an interval of 1 and I À anions with numbers from 6 to 15 with an interval of 3. Aer removing the ions, we introduced a uniform background charge to neutralize the defective charged MAPbI 3 system as applied in the previous MD simulations. 37 For the AEMD simulations, we adopted the perfect tetragonal phase with simulation supercells with cell lengths of L x ¼ L y x 2.6 nm and L z x 38-152 nm (see Fig. S2 † for L z x 51 nm). In fact, it was conrmed from the previous work 20 that k calculated with L x ¼ L y x 3 nm was changed by only 5% from that with L x ¼ L y x 11 nm.
All the MD simulations were performed using the LAMMPS package. 38 The pppm solver with an accuracy of 10 À5 was used to compute the coulombic interactions, and the velocity Verlet algorithm was used for integrating motion with a time step of 0.5 fs. Periodic boundary conditions were applied in the x, y and z directions. The control of temperature and pressure during the simulation was managed by the Nosé-Hoover thermostat and barostat. In the AEMD simulations, the entire system was initially relaxed by NPT simulations at a temperature of 300 K and pressure of 0 atm for 500 ps. Then, two split regions 1 and 2 were equilibrated by NVT simulations at different temperatures of 325 K and 275 K for 500 ps. Subsequently, NVE equilibration was carried out for sufficient simulation time, which was different depending on the cell length. In the EMD simulations, meanwhile, NVE simulation was carried out for 5 ns aer NPT and subsequently NVT equilibrations for 500 ps. To reduce statistical error, we repeated 5 and 10 independent runs with different initial atomic velocities for both perfect and defective models. For calculations of elastic constants, we performed molecular mechanic (MM) simulations aer NPT equilibration at 330 K for 10 6 steps.

MYP force eld
Choosing a suitable force eld is a very important factor for determining the accuracy of classical MD simulations. To date, several force elds have been developed for MAPbI 3 , 18,19,39-42 and among them the MYP force eld derived from ab initio DFT calculations 39,40 has been widely used to study the lattice thermal properties of MAPbI 3 . 17,20 In particular, Barboni et al. calculated the jumping rate of iodine vacancies in defective MAPbI 3 using the MYP force eld. 37 In this work, we also employed the MYP force eld to investigate the effect of vacancy defects on the lattice thermal properties of MAPbI 3 .
The MYP force eld is composed of the three terms, where U OO , U II and U OI are the organic-organic, inorganicinorganic and organic-inorganic potentials, respectively. Thus, the U OO part contains intermolecular and intramolecular interactions between the MA moieties as described in the standard AMBER force eld, 43 and the U II part includes the interaction between Pb and I atoms with the Buckingham-type potential. Meanwhile, the U OI part is the sum of Buckingham, coulombic and Lennard-Jones (12-6) terms as follows, 39 where r ij is the distance between the i-th and j-th atoms, Z i is the partial charge of the i-th atom, 3 0 is the dielectric constant of a vacuum, and r ij , s ij , 3 ij , s ij are the empirical parameters. These potential parameters are provided in Tables S1-S6 (see the ESI †) and can be found in the previous work. 39 3 Results and discussion

Thermal conductivity of a perfect system
To verify the accuracy of the selected methods, we calculated the thermal conductivity of perfect MAPbI 3 at 300 K using both the AEMD and EMD methods, comparing with the previous theoretical and experimental data. Firstly, the AEMD method was applied to calculate the k value of perfect MAPbI 3 with a tetragonal structure at a temperature of 300 K. Through the analysis of NPT simulation, we found that the equilibrated lattice parameters and elastic properties of the tetragonal phase agreed well with the previous data (Table S7 †). It was known that the k value of an insulating solid, calculated using MD simulations, depends on the size of the simulation box when the size is smaller than the mean free path of the phonon. 20 The dependence of k on the simulation box size can be expressed as, where k(L z ) and k(N) are the thermal conductivities of samples with nite and actual size, respectively, and a is the tting parameter. This indicates that the actual k(N) can be obtained by applying linear extrapolation of 1/k vs. 1/L z . To estimate k(L z ) with increasing L z , we need to determine the volumetric heat capacity C V and mass density r. The heat capacity was estimated using the linear relation between the total energy E and temperature T as Linear tting of the total energy E of the system at different temperatures (see Fig. S3 †) produced C V as 1.327 Â 10 À5 eV K À1ÅÀ3 in good agreement with the previous value of about 1.3 Â 10 À5 eV K À1 A À3 . 20 In the AEMD simulation, it was found that as the NVE simulation proceeds, the average temperature difference DT between the two sub-regions decreases, nally becoming the nal value of zero, as shown in Fig. 3(a). Initially, DT was 50 K for all cases of cell size L z , which were set from 38 nm to 152 nm with 5 intermediate points. While increasing the cell size L z , the decreasing rate of DT gets slower while the tting accuracy increases. In each case of L z , we performed tting of the obtained DT(t) (orange line shown in the inset of Fig. 3(a)) to an exponential function (see the ESI for details †), yielding the thermal diffusivity k and thus the thermal conductivity k(L z ). We then carried out linear tting of 1/k(L z ) vs. 1/L z , as shown in Fig. 3(b), and thus obtained 1/k(N) as 1.172 m K W À1 . This gives the actual k of tetragonal MAPbI 3 as 0.853 W m À1 K À1 , which is in good agreement with the previous simulation data of 0.8 W m À1 K À1 obtained using the same MYP force eld and the same AEMD method. 20 Next, we applied the EMD method, i.e., the Green-Kubo method, to calculate the k value of perfect MAPbI 3 . Fig. 4 shows the principal k a in each Cartesian direction (a ¼ x, y, z) and the volumetric k v for perfect MAPbI 3 as the EMD simulation proceeds, these were dumped every 50 ps to estimate the HFACF and k vs. correlation time (see insets). The HFACF and k were plotted every 1 fs, and they converged to zero and a certain value aer 10 ps, respectively. The principal thermal conductivities were determined to be k x ¼ 0.328, k y ¼ 0.480 and k z ¼ 0.406 W m À1 K À1 , and their average was determined to be k v ¼ 0.405 W m À1 K À1 . These agree well with the previous simulation data of 0.25-0.40 W m À1 K À1 , obtained using the same MYP atomic potential and the same EMD method, 17 and the experimental values of 0.3-0.5 W m À1 K À1 for perfect MAPbI 3 . 5,11,13 It is worth noting that the anisotropy of the thermal conductivity is associated with the structural anisotropy of the MAPbI 3 crystal with the orientation of the MA cation. 10 Although the same interatomic potential was used, the AEMD and EMD methods produced different thermal conductivities. In principle, these two methods should give similar results. It seems that the AEMD method gives the upper limit of thermal conductivity for the bulk, since it refers to an innitely long sample with perfect crystallinity. In Fig. 3(b), the shaded region indicates the experimental range, for which L z is between 130 nm and 330 nm, being considered as the "effective" mean free path in accordance with the experimental coherence length below 150 nm. 44 The AEMD method has a much larger computational burden for a solid with low k like MAPbI 3 because DT diminishes more slowly to zero. Therefore, the EMD method is used for defective MAPbI 3 .
We then investigate the temperature dependence of lattice thermal conductivity of perfect MAPbI 3 as the temperature increases from 300 to 420 K with an interval of 30 K. Fig. 5 shows k as a function of temperature calculated with the EMD method in comparison with the previous theoretical and available experimental data. In our work, there is a minor difference of about 0.015 W m À1 K À1 between the tetragonal (0.449 W m À1 K À1 ) and pseudocubic (0.464 W m À1 K À1 ) phases at 330 K. Some previous work reported a big jump in the thermal conductivity at 330 K upon phase transition from the tetragonal to pseudocubic phase (blue dotted lines), 13,15,18 whereas more recent works revealed that k does not abruptly increase at 330 K, even upon phase transition (red dotted lines). 5,6,14,17,45 Our work agrees with the latter cases, conrming that the lattice thermal conductivity rarely changes upon increasing temperature. Qian et al. 18 adopted different interatomic potentials for different phases, giving a big jump in thermal conductivity upon phase transition from tetragonal to cubic phases, but this makes it difficult to associate the thermal conductivity with the cation dynamics. Meanwhile, Wang et al. 17 used the same MYP potential for different phases, resulting in no jump in thermal conductivity as in our work. This is connected with the viewpoint that the MA cation has a minor effect on the thermal transport at higher temperatures, since the difference in crystal structures between tetragonal and cubic phases is in the distribution and orientation of the MA cations.

Thermal conductivity of defective system
Next, we calculated the lattice thermal conductivities of the defective MAPbI 3 systems in the pseudocubic phase at 330 K, which contains three kinds of vacancies, V MA , V Pb and V I , with different concentrations ranging from 0.0 to 1.0% with an interval of 0.2%. Fig. 6 shows the calculated k values of defective MAPbI 3 as a function of vacancy concentration. It was found that the thermal conductivity decreases gradually when creating the vacancy defects, and shows a decreasing tendency overall, as the vacancy concentration increases. At the beginning of vacancy creation, between 0 and 0.2% concentration, the reductions in k for V I and V Pb were large compared to those for the case of V MA . On the way to 1%, although there are some slight rises at 0.4% for V MA and V I , and at 0.6 and 0.8% for V Pb , we could say that the thermal conductivities of defective systems decrease overall as the concentration increases from 0% to 1%. These slight rises at the middle concentrations may come from the intrinsic uctuation feature of classical MD simulations, although we repeated the simulations at least 10 times. For comparison among the kinds of vacancy, V I causes the greatest reduction of k from 0.464 AE 0.011 W m À1 K À1 at 0% to 0.371 AE 0.006 W m À1 K À1 at 1%, V Pb shows the next greatest reduction to 0.375 AE 0.008 W m À1 K À1 at 1%, and V MA gives the least reduction to 0.398 AE 0.006 W m À1 K À1 at 1% concentration. In general, it can be thought that the creation of vacancy defects in an insulating crystalline solid causes a disturbance of the (acoustic) phonon propagation, reducing the lattice thermal conductivity as these are the main thermal carriers. From our calculation results, it can be regarded that the atomic vacancies of V I and V Pb have a greater scattering effect on phonon transfer than the molecular vacancy of V MA . This is related to the fact that the acoustic phonon modes carrying the heat are mainly from the Pb-I cage, while the MA molecule plays a role in scattering phonons.
For each case of vacancy defects, we discuss the mechanism behind the reduction of lattice thermal conductivity of MAPbI 3 in comparison with the previous theoretical and experimental data. In the case of V MA , it should be noted that the organic MA + cations are the main contributors to the ultralow thermal conductivity of MAPbI 3 . 11,15,17,19,20 It was found that coupling between the rotational and translational motions of the MA cation plays a key role in suppressing the thermal conductivity of MAPbI 3 . 19 In fact, Hata et al. demonstrated through MD simulations that the real MAPbI 3 clearly exhibits a lower k value (0.185 W m À1 K À1 ) compared to the hypothetical PbI 3 bare inorganic frame model (0.326 W m À1 K À1 ), nding that the rotational motions of MA cations mainly induce the suppression of thermal conductivity. 19 Moreover, coupled rotational and translational motions of MA were found to interact with the Pb-I cages and derive couplings between the lattice vibrations, suppressing the phonon transport in MAPbI 3 . The removal of some molecules leads to a lower number of molecular degrees of freedom and thus could favor the thermal transport. 20 According to our calculations, however, k decreases as the concentration of MA vacancies increases, in contrast to the previous data and general insight. The MA + cation seems to have a signicant suppressing effect on k only at low temperature, while it has a relatively small contribution at room temperature. 12 Instead, the internal phonon-phonon scattering dominantly inuences the thermal transport in halide perovskites at room temperature. In fact, similar ultralow thermal conductivities were observed in the all-inorganic halide perovskites with no organic molecule such as 0.36 W m À1 K À1 for single crystalline nanowire CsPbBr 3 (ref. 46) and 0.45 AE 0.05 W m À1 K À1 for CsPbI 3 (ref. 47) at room temperature, being associated with the metal-halide cluster rattling modes. From such considerations, it can be noticed that, at higher temperature, the MA + cation contributes negligibly to the lattice thermal conductivity of MAPbI 3 . The MA molecular vacancies associated with the local perturbations of the PbI charge and structure can  act as scattering centers for heat transfer and induce Rayleighlike scattering of phonons, resulting in a reduction of thermal conductivity of MAPbI 3 . 20 For the case of Pb vacancies, one can see that k decreases more rapidly than in the V MA case as the vacancy concentration increases, indicating that the effect of V Pb on k is more pronounced than that of V MA . Considering that the acoustic phonon modes are dominated by Pb and I atoms, it is natural to see such reduction of k by the creation of V Pb , which can play the role of phonon scattering center. Neither theoretical nor experimental data for this V Pb effect was reported before. Here, it is worth noting that when substituting a lighter metal cation for Pb, the lattice thermal conductivity of HPs was found to decrease for the pseudocubic phase but increase for the tetragonal phase. For instance, upon replacing Pb with Sn in MAPbI 3 , Mettan et al. found a decrease of k from 0.5 W m À1 K À1 to 0.09 W m À1 K À1 for the pseudocubic phase, but an increase to 0.69 W m À1 K À1 for the tetragonal phase. 51 In the case of allinorganic halide perovskites, Yang et al. revealed a decrease from the 0.45 AE 0.05 W m À1 K À1 of CsPbI 3 to 0.38 AE 0.04 W m À1 K À1 for CsSnI 3 . 47 The k reduction upon this replacement is because the third order force constant in the Sn-I bond is larger than that in the Pb-I bond, resulting in stronger anharmonicity. 52 The role of the iodine vacancy is similar to the case of V Pb as shown in Fig. 6, presenting that k rapidly decreases as the V I concentration increases. We discuss the effect of substituting a lighter halogen anion for the I À anion on the thermal conductivity. From experiments, the thermal conductivity of single crystalline MAPbX 3 solid was found to increase from 0.34 AE 0.12 W m À1 K À1 for X ¼ I to 0.44 AE 0.08 W m À1 K À1 for X ¼ Br and to 0.50 AE 0.05 W m À1 K À1 for X ¼ Cl. 13 A similar tendency was observed for the cases of all-inorganic halide perovskites. 53 Such an increase in thermal conductivity with decreasing the atomic number of the halogen atom is associated with the increasing elastic stiffness constants.
Our simulation result of the thermal conductivity of defective MAPbI 3 , containing vacancies such as V MA , V Pb and V I , indicates that the low thermal conductivity of MAPbI 3 originates from the inorganic Pb-I framework lattice. 14,19,20 From this perspective, it is worth noting that the k value of the PbI 2 solid is also very low, 0.22 W m À1 K À1 at 293 K. 14 In many previous studies of HPs, it has been found that a heavier mass of halide atom could lead to lower thermal conductivity. If the constituent atoms are heavier, i.e., with larger atomic number and larger ionic radius, the phonon vibrations in principle get weaker while strengthening its scattering, thereby reducing thermal transport. In particular, the electronegativity of the X atom decreases with increasing the atomic mass, and thus the B-X bond length increases, resulting in lowering the thermal conductivity. For the case of the B-site atom, there are only a few studies for thermal conductivity, and therefore, more work is needed to elucidate the effect of changing the B-site metal. 10 Such big effects of the B-site metal and X-site halogen can provide the possibility of tuning the thermal conductivity of MArelated HPs, MABX 3 .

Elastic properties of perfect and defective systems
Many previous studies analyzed the thermal transport properties of MAPbI 3 using the phonon group velocity and phonon lifetime based on a picture of particle-like phonon propagation in crystals. However, recent studies have proposed a coexistence of crystal-and glass-like transport mechanisms for heat conduction in HPs, demonstrating that the crystal-like transport mechanism is dominant at low temperatures while the glass-like transport plays a dominant role at higher temperatures. 54 For the case of CsPbBr 3 , for instance, crystal-and glasslike contributions to k were found to be 78% and 22% at 50 K, while 30% and 70% at 300 K. 54 This indicates that description of k only by particle-like phonon propagation is not adequate for HPs at room temperature.
For the description of crystal-and glass-like thermal transport, the elastic moduli and sound velocity are interesting because they are evaluated for polycrystalline solids. In fact, Elbaz et al. scaled the thermal conductivity with the sound velocity to explain the measured low thermal conductivities of 0.34-0.73 W m À1 K À1 at room temperature for MAPbX 3 (X ¼ I, Br, Cl). 50 Moreover, the order of thermal conductivities of HPs was determined from the mechanical properties and sound velocities. 17 Therefore, we evaluated the mechanical properties, including elastic constants and moduli, and the sound velocity Table 1 Elastic stiffness constant (C ij ), bulk modulus (B), shear modulus (G), Young's modulus (E), Poisson's ratio (n), mass density (r), average sound velocity (v) and thermal conductivity (k) for perfect and defective MAPbI 3

System
Elastic stiffness constants (GPa) Elastic moduli (GPa) from them for perfect and defective MAPbI 3 , this can be meaningful in explaining the variation tendency of k upon the creation of vacancies. Table 1 presents the calculated elastic stiffness constants (C ij ), elastic moduli including bulk (B), shear (G) and Young's (E) moduli, Poisson's ratio (n), density (r), and sound velocity (v), with the thermal conductivity (k) for perfect and defective MAPbI 3 in the pseudocubic phase at a vacancy concentration of 1.0%. Our results from the classical EMD method were found to be in good agreement with the previous calculations from the DFT method. It should be noted that, although the pseudocubic phase was initially assigned, the nal phase equilibrated by NPT simulation was found to be orthorhombic with very small differences of lattice constants due to the numerical noise being inherent in MD simulation, leading to 9 independent elastic constants. We should point out that the average sound velocity of 1922 m s À1 by Faghihnasiri et al. 49 was evaluated with a mistaken equation and thus we corrected it to be 1628.6 m s À1 by applying the correct equation (eqn (S14) †), which is in good agreement with our results and other calculations. 48 It was found that the elastic constants, moduli and sound velocities of defective MAPbI 3 are smaller than those of the perfect system, leading to the lower thermal conductivity. As shown in Fig. 7 and Table 1, there is a certain relation between the thermal conductivity and elastic modulus or sound velocity; larger elastic constants and sound velocities lead to higher thermal conductivity.

Conclusions
In this work, we have investigated the lattice thermal conductivity of perfect and defective MAPbI 3 using classical MD simulations with the MYP force eld. With the AEMD and EMD methods, the thermal conductivity of perfect MAPbI 3 was determined to be 0.853 and 0.405 W m À1 K À1 at 300 K using the supercells, and its temperature dependence was studied by increasing the temperature from 300 K to 420 K with the EMD method. We made supercell models of defective systems containing vacancy defects such as V MA , V Pb and V I with different concentrations to investigate the effect of vacancy concentration on the thermal conductivity. Our calculations revealed that, as the vacancy concentration increases from 0 to 1% with an interval of 0.2%, the thermal conductivity decreases overall with slight rises at 0.4% for V MA and V I and at 0.6 and 0.8% for V Pb . We have shown that such vacancies can act as phonon scattering centers, thereby causing the decrease in thermal conductivity. Furthermore, the elastic moduli and sound velocities of defective MAPbI 3 were determined to be smaller than those of the perfect system, revealing that the low sound velocity in relation with the bulk and shear moduli is responsible for the ultralow thermal conductivity of MAPbI 3 . Our results are helpful for designing functional hybrid halide perovskites for photovoltaic and thermoelectric devices with high performance.

Author contributions
Song-Nam Hong and Chol-Jun Yu developed and planned the original project. Song-Nam Hong performed the calculations and draed the rst manuscript. Chol-Jun Yu supervised the work. Song-Hyok Choe assisted with the MD simulation. Un-Gi Jong and Yun-Hyok Kye contributed to useful discussions. All authors reviewed the manuscript.

Conflicts of interest
There are no conicts to declare.