Polar rotor scattering as atomic-level origin of low mobility and thermal conductivity of perovskite CH3NH3PbI3

Perovskite CH3NH3PbI3 exhibits outstanding photovoltaic performances, but the understanding of the atomic motions remains inadequate even though they take a fundamental role in transport properties. Here, we present a complete atomic dynamic picture consisting of molecular jumping rotational modes and phonons, which is established by carrying out high-resolution time-of-flight quasi-elastic and inelastic neutron scattering measurements in a wide energy window ranging from 0.0036 to 54 meV on a large single crystal sample, respectively. The ultrafast orientational disorder of molecular dipoles, activated at ∼165 K, acts as an additional scattering source for optical phonons as well as for charge carriers. It is revealed that acoustic phonons dominate the thermal transport, rather than optical phonons due to sub-picosecond lifetimes. These microscopic insights provide a solid standing point, on which perovskite solar cells can be understood more accurately and their performances are perhaps further optimized.

O ver last few years, the inorganic-organic hybrid perovskite solar cells have taken the worldwide research community by storm [1][2][3] .Typically, the energy conversion efficiencies of solar cells based on methylammonium lead iodide (CH 3 NH 3 PbI 3 ) have been improved from the starting 4% in 2009 to higher than 20% in 2015 (refs 4,5). Very recently, the successes in growth of inch-sized high-quality single crystals and the device integration on wafers have paved the route to large-scale practical photovoltaic applications 6,7 . These hybrid compounds exhibit distinct physical properties from conventional semiconductors. Their hot-phonon bottleneck effect of energetic carriers is obviously stronger than that of GaAs and other inorganic perovskites, which has been attributed to the acoustic-optical phonon up-conversion [8][9][10] . The mobility of charge carriers is relatively smaller compared with classical semiconductors [1][2][3]11,12 and the scattering mechanism is still under debate [13][14][15] . Resembling the charge transport, the thermal transport is unusual as well. The thermal conductivity is surprisingly low, B0.5 Wm À 1 K À 1 for single crystals at room temperature 16 , which is tenfold lower than that of GaAs (ref. 12) and is even lower than that of amorphous silicon 17 .
It is evident that atomic dynamics underlies these peculiarities. However, the atomic-level description of CH 3 NH 3 PbI 3 is complicated by the hybrid nature where both molecular jumping rotations and phonon excitations have to be taken into account. Moreover, these two components also interact via hydrogen bonds between H and I (refs [18][19][20]. The best approach for this issue is, without a doubt, inelastic neutron scattering (INS). The density functional theory (DFT) lattice dynamics calculations indicate that low-energy phonons are entirely composed of motions of Pb and I (refs [21][22][23][24][25][26], which are easily excited at relatively lower temperatures and hence take main parts in determining the physical properties. These phonons can be efficiently probed throughout the Brillouin zones, due to the larger coherent scattering cross-sections of Pb and I. Simultaneously, the incoherent scattering cross-section of H assures the individual motions of molecules can be determined by extracting the quasi-elastic broadening underneath the elastic line, which is known as quasi-elastic neutron scattering (QENS) method 27 . In this work, we show the complete phonons and jumping rotational modes across the low-temperature phase transition, which are obtained by measuring a large single crystal at two high-resolution time-of-flight neutron spectrometers that cover a wide energy window ranging from 0.0036 to 54 meV. These results are well supported by the complementary molecular dynamics (MD) simulations. It is revealed that the molecular dipole order plays a dominant role in determining charge transport and thermal transport properties of CH 3 NH 3 PbI 3 .

Results
Jumping rotational modes of [CH 3 NH 3 ] þ . CH 3 NH 3 PbI 3 crystallizes in the common perovskite structure where the organic cation [CH 3 NH 3 ] þ occupies the centre of the PbI 6 octahedral cage, as shown in Fig. 1a. Neutron and X-ray powder diffraction investigations suggest that it undergoes successive phase transitions, from cubic (Pm 3m) to tetragonal (I4/mcm) at 330 K, and then to orthorhombic (Pnma) at 165 K (refs 28,29). Our single crystal exhibits a sharp phase transition between 160 and 165 K, as shown in Supplementary Fig. 2, which is in agreement with these reports. While it is believed that this transition is coupled with orientation of [CH 3 NH 3 ] þ , current reports on rotational dynamics are contradictory 30,31 , probably owing to the lack of results on single crystals. Given that this compound is extremely sensitive to ambient conditions 32,33 , measurements on a single crystal (with minimized surface contribution) are anticipated to deliver more intrinsic information. Here, the accurate jumping rotational dynamics of molecules is obtained by performing ultrahigh energy resolution QENS measurements on the large single crystal at the backscattering neutron spectrometer, DNA.
First, we demonstrate the jumping rotational dynamics in the orthorhombic phase. As shown in Fig. 2a, the spectrum at 50 K is basically described by a delta function convoluted to the instrumental resolution (the full width at half maximum is B0.0036 meV). Around 80 K, the QENS signal described by a Lorentzian function is identified underneath the elastic line, indicative of the activation of rotations. The QENS intensity becomes more noticeable at 140 K. The momentum-transfer averaged spectra along [00l] direction are systematically fitted by including a Lorentzian function, a delta function and a constant background, which are convoluted to the instrumental resolution. An example is given at 140 K in Fig. 2b. The half width at half maximum of the Lorentzian function, G, is examined as a function of momentum transfer and temperature. As shown in Fig. 2e, G is almost independent on l. Since G is inversely related to the relaxation time, such a value gives rise to an average relaxation time of 23(1) ps at 140 K. The temperature dependence is fitted to the Arrhenius relation and the activation energy B47.9(6) meV is obtained (Fig. 2f). The elastic incoherent structure factor (EISF), defined as the ratio between elastic intensity and the sum of QENS and elastic intensity 27 , is shown in Fig. 2g at 140 K. It is described by the threefold jumping rotational (C 3 ) model, in which EISF where j 0 is the zeroth-order spherical Bessel function and d H-H is the average distance between two adjacent protons, B1.72 Å in terms of neutron powder diffraction results 28 (see Supplementary  Fig. 3 for details).
In the tetragonal phase, the jumping rotational dynamics becomes more complicated. The spectrum obtained at 180 K at DNA displays a narrower QENS profile than that at 140 K, as compared in Fig. 2b,c. The relaxation time is determined to be 64(2) ps. At the same time, a much broader Lorentzian function is found by using the cold neutron time-of-flight spectrometer AMATERAS, which provides a much wider energy window (Fig. 2d). Such a large width corresponds to a relaxation time as short as 0.71(3) ps. Likewise, EISF is calculated based on DNA data by following the same procedure as above 27 . However, the derived EISF is overestimated because QENS intensity of the faster mode is missed in the total intensity. The recent group theory analysis on jumping rotational dynamics of CH 3 NH 3 PbI 3 suggests a coupled C 3 and fourfold (C 4 ) mode in the tetragonal phase 31 . Therefore, we compare the derived EISF with this C 4 #C 3 model. Shown in Fig. 2h is its comparison with the ratio between elastic intensity and total intensity excluded the C 3 component. The detailed calculation is given in Supplementary Fig. 4. It can be seen that the experimental data is well reproduced in this way, indicating that a C 4 mode sets in after the phase transition to the tetragonal phase and the C 3 mode becomes B30 times faster. Hereto, we are able to conclude on the jumping rotational dynamics. In the orthorhombic phase, the protons in CH 3 and NH 3 undergo jumping rotations with respect to the C-N axis. The relaxation time is dramatically shortened in the tetragonal phase while the C-N axis starts to rotate with respect to the c axis of the crystal structure. The rotational modes are illustrated in Fig. 1b. Our results about the jumping rotational dynamics are in agreement with the recent QENS measurements on a powder sample by Chen et al. 31  Phonons of the orthorhombic and tetragonal phases. The jumping rotational dynamics of an individual molecule is well understood above and now we move to the collective dynamics, that is, phonons, which are measured by selecting several incident energy of neutrons (E i ) at AMATERAS. The obtained dynamic structure factor S(q, E) is shown in Fig. 3  With smaller E i , more details are revealed. We find the peak at 11.3 meV shown in Fig. 3d is actually a superposition of three peaks at 10.6, 11.5 and 12.6 meV (Fig. 3b,e). There seems a gap located between 6 and 10 meV, but several weaker modes are roughly identified as highlighted in the inset. The data with E i of 7 meV (Fig. 3c,f) suggests four low-energy modes at 2.28, 3.11, 3.79 and 4.35 meV. Thanks to the high-resolution, we are able to identify several modes that are overlooked in the previous Raman scattering 26 and INS measurements 21 . Particularly, the latter and the associated DFT calculations are compared with our experimental results and MD simulations, as listed in Table 1    can be seen that our measurements match well with the simulations, except the modes 21 and 22, which are dominated by the molecular twist motion 24 . This discrepancy is also reported in another MD simulation study using the same potential functions 24 , and might be related to the zero-point-energy fluctuations that lead to more lattice distortions as suggested in the DFT calculations on P1 symmetry 21 .
Detailed acoustic phonons of the Brillouin zone (220) are examined with E i of 4 meV. Shown in Fig. 4a-d are the longitudinal and transverse phonons at 5 and 180 K, respectively. Even though the background level is high, the dispersions can be seen in both directions. At 5 K, the dispersion curve of longitudinal acoustic (LA) phonon emanates from zone centre of (220), undergoes a steep increase and finally approaches the top in the middle of zone centre and zone boundary (h about 0.2), where it is overlapped with longitudinal optical (LO) phonons. The dispersion of transverse acoustic (TA) phonon exhibits a moderate slope. Intense transverse optical (TO) phonons are also present. At 180 K, the intensity of phonons is fairly enhanced and we determine the accurate dispersions of LA and TA phonons by fitting the momentum-transfer (filled circle) or energy-transfer (filled square) averaged spectra. The solid lines represent the slopes of dispersion curves near the zone centres, which determine the group velocities: n LA ¼ 2841 m s À 1 and n TA ¼ 1155 m s À 1 . Shown in Fig. 4e is the fitting of the TA phonon, which yields the lifetime of 3.61 (42) ps in the middle of the zone boundary and zone centre. For the LA case, we determine that the lifetime is B4.39(46) ps at 180 K in the vicinity of the zone centre ( Supplementary Fig. 7). Recently, Whalley et al. 35 have calculated the lifetimes by considering the phonon-phonon interaction. At 300 K, the lifetimes of phonons range around 2 ps, which is very close to our results. The product of lifetime and group velocity defines mean-freepath. At 180 K, short lifetimes and small velocities yield mean-free-paths of B125 and 42 Å for the LA phonon and TA phonon, respectively. Since we have first determined these critical quantities, it is allowed to compare with previous theoretical counterparts. Listed in Table 2 is the comparison of velocities and mean-free-paths of acoustic phonons [35][36][37][38] . It is worth noting that our INS results are well supported by the MD simulations based on the same potential functions as we used here 34,38 .   More strikingly, the LO and TO phonons around 2.28 meV completely vanish at 180 K. The data of the (004) zone ( Supplementary Fig. 8) indicates that they survive at 150 K. Thus, we compare the optical phonons located at 10.7 and 12.6 meV at 5, 150 and 180 K, which are two of strongest optical phonons observed in the spectra. At 5 K, these optical phonons are well-defined, as shown in Fig. 3e. With warming up to 150 K, the lifetimes of TO phonons (Fig. 4f) at 10.7 and 12.6 meV are reduced down to 0.91 (15) and 0.62(7) ps while those of LO phonons (Fig. 4g) are reduced down to 0.82(5) and 0.57(2) ps, respectively. These lifetimes are much shorter than 5.40(88) ps for the TA phonon (inset of Fig. 4e). At 180 K, just above the phase transition, both TO and LO phonons are too broadened to be distinguished. This means the actual lifetimes at 180 K are much shorter than those at 150 K. This dramatic broadening is greatly consistent with the Raman scattering study in which the lifetime is shortened from 0.6 to 0.15 ps (ref. 26). The characteristic timescales of jumping rotational modes and phonons are summarized in Table 3.
MD simulations of phonons through the phase transition. The finite-temperature phonon excitations are considered in classical MD simulations using recently developed MYP potential functions 34 . In Fig. 5a, we show the phonon density of state (DOS) at 100 K. It is clear that motions of the inorganic part are dominant in the low-energy region while the high-energy spectrum is mostly contributed by molecular motions. The energy of each mode is determined in Supplementary Fig. 6 and compared with the experimental counterpart in Table 1. The distinct changes of lifetimes of acoustic and optical phonons through the phase transition are consistent with the simulations where the DOS spectrum of optical phonons is much more broadened than that of acoustic ones. As shown in Fig. 5b, the peaks of optical phonons (between 2 and 16 meV) are remarkably broadened in the tetragonal phase, in agreement with the broadening shown in Fig. 4f,g (for medium-energy optical phonons) and Supplementary  Fig. 8 (for low-energy optical phonons). Instead, acoustic phonons shown in the inset are still distinguishable even if their intensities are much weaker than those of optical phonons. The significantly reduced lifetimes of optical phonons in the tetragonal phase imply an additional decay channel taking effect except the phonon-phonon interaction. Whalley et al. 35 have considered the phonon-phonon interaction and directly calculate the lifetimes. It is shown that there are similar lifetimes for phonons with frequencies below B3.5 THz (B14 meV). Namely, the phonon decay based on anharmonicity in this energy region is not selective. In contrast, the MD simulations including the long-range electrostatic interactions indicate that optical phonons are more broadened than acoustic ones. Note that the experimentally observed dipole order has been reproduced using the same methodology of simulations 34 . Thus, the additional scattering is most likely attributed to the activated orientational disorder of dipoles in the tetragonal phase, which is similar to the phonon scattering by paraelectric centres 39,40 .

Discussion
In the tetragonal phase, the temperature dependence of mobility shows a À 3/2 power-law, which was initially attributed to the scattering by acoustic phonons 13 . Later, both the experimental investigation of photoluminescence spectra and the transport theoretical study manifest the Fröhlich interaction between charge carriers and LO phonons is the major scattering mechanism 14,15 . Indeed, considering the deformation potential and the piezoelectric potential as major mechanisms, the calculated mobility is much larger than experimentally measured one 41,42 , indicating that acoustic phonons are less important. Notwithstanding, it is still unclear why the mobility of hybrid perovskites is quite small compared with other inorganic semiconductors that have comparable Fröhlich coupling strengths 14 . More importantly, such a scenario is not able to reproduce the critical-like temperature dependence of mobility at the orthorhombic-to-tetragonal phase transition 43 .
According to the current structural model, the orthorhombic phase is antiferroelectric 28 . Our QENS data suggests that   the dipoles are disordered in the picosecond temporal scale in the tetragonal phase. Thus, this orthorhombic-to-tetragonal phase transition is indeed an order-to-disorder transition of dipoles, just as the dielectric constant displays a sharp peak at the transition temperature 44 . Namely, CH 3 NH 3 PbI 3 is not only a polar semiconductor like GaAs, but it is also antiferroelectric. Extensive earlier reports on electrical transport properties of ferroelectric perovskite oxides and IV-VI semiconductors indicate the critical-point polarization fluctuations result in a singularity of resistivity at the transition temperatures, that is, a positive resistivity-temperature slope as approaching the transition temperature from the low-temperature side [45][46][47] . Similarly, the dipole scattering is also active in hybrid perovskites with polar cations. The important role of molecular dipoles is manifested at the sharp difference of the magnitude of mobility of hybrid CH 3 NH 3 PbBr 3 and all-inorganic CsPbBr 3 (ref. 43). At B120 K where both compounds crystallize in orthorhombic structures, the mobility of the former is about 5 times lower. In addition, the critical anomaly has been analytically solved for antiferroelectrics, given as dr crit dT / ð T Tc À 1Þ À a at T4T c and as dr crit dT / ð1 À T Tc Þ 2b À 1 at ToT c , where r crit is critical resistivity, T c is transition temperature, a and b are critical exponents 47 . Thus, one obtains m / À ð1 À T Tc Þ À 2b at ToT c . We find the temperature dependence of mobility near the phase transition of CH 3 NH 3 PbI 3 can be well reproduced with an extremely small value of b ( Supplementary Fig. 9), indicative of the first-order nature 48 . The sharp change of mobility at the phase transition well coincides with the activation of dipole disorder (see Table 3).
Even though the molecular dipoles are disordered on atomic level as they undergo ultrafast C 4 jumping rotational mode in the picosecond temporal scale, it is possible that there are nanoscale/mesoscale correlations with slower dynamics 49 . This kind of large-scale structures are compatible with ferroelectric domains predicted in DFT 50 and observed using piezoelectric force microscopy 51 . They are perhaps relevant to the polaron picture proposed for hybrid perovskites 52 , which is the key to understand the extremely low recombination rate of carriers.
It is known that lattice thermal conductivity is equal to 1/3Cv 2 t where C is the lattice specific heat, t is the average group velocity and t is the average lifetime of phonons 53 . Even through the origin of ultralow thermal conductivity is naturally ascribed to either lower bulk modulus (proportional to v 2 ) or shorter t, the profound mechanism is still unclear due to the lack of the complete phonon data, especially acoustic phonons. Our INS study suggests the contribution of optical phonons to thermal transport is marginal because of the significantly shortened lifetimes. On the contrary, the thermal transport is dominated by acoustic phonons. Their nanoscale mean-free-paths together with smaller velocities are responsible for such a small thermal conductivity. Experimentally, the thermal conductivity exhibits a smooth temperature dependence except an abrupt discontinuity at the transition owing to the anomaly of the specific heat 16 (Table 3). This is well consistent with the less varied lifetimes of acoustic phonons across the phase transition. If the optical phonons were major heat carriers, a dramatic drop of thermal conductivity would be expected at the transition. Molecular vibrations have been recognized to be detrimental to thermal transport 35,54,55 . The anharmonic DFT calculations by Whalley et al. 35 and MD simulations with MYP potential functions by Caddeo et al. 54 suggest the scattering of low-energy phonons (mostly dominated by the inorganic part) by molecular vibrations acts as the major decay channel, but in another simulation by Hata et al. 55 molecular vibrations only mediate the interaction between phonons consisting of PbI motions. Our results seem to agree with the first case because the twist mode at 37.6 meV is significantly broadened even in the orthorhombic phase ( Supplementary Fig. 10).
Given that the molecular vibrations (high-energy optical phonons) are very flat modes, their interaction with acoustic phonons is in analogy to the rattling-mode-phonon coupling in the paradigm of electron-crystal phonon-glass 56 , where less dispersive modes  ARTICLE highly scatter acoustic phonons, leading to glassy-like thermal conductivity 57 . It is considerably promising to achieve superior thermoelectric materials in hybrid inorganic-organic perovskites by enhancing charge transport. Simulations indicate the thermoelectric figure-of-merit of CH 3 NH 3 PbI 3 can reach to 3 at 600 K with optimal carrier density of 10 19 cm À 3 (ref. 58). It is also predicted that CH 3 NH 3 SnI 3 may have much higher value of figure-of-merit than Pb case 59 . Moreover, the low thermal conductivity does not only directly enable hybrid perovskites to be appealing thermoelectric materials, but also plays a fundamentally important part in photovoltaic properties. The suppressed thermal conduction promotes the acoustic-optical phonon up-conversion, benefitting the slow cooling rate of energetic carriers 10 . Such an enhanced phonon bottleneck effect facilitates the realization of the hot-carrier solar cell concept with an extremely high power conversion efficiency 60 .
In summary, we have precisely determined jumping rotational dynamics of molecules and phonons in hybrid perovskite CH 3 NH 3 PbI 3 by conducting high-resolution QENS and INS measurements, which are all consistent with the finite-temperature MD simulations. The characteristic timescales of various dynamics are completely revealed. The dipole order of organic cations plays a decisive role in determining the physical properties such as carrier mobility and thermal conduction, distinguishing from inorganic perovskites. Such observations are probably general for all hybrid compounds with polar molecules. These fundamental insights are beneficial to further study of perovskite solar cells.

Methods
Single crystal growth and mounting. The crystal was grown by using the solution method described in ref. 6. The single crystal used in neutron measurements is 2.32 g in weight and well-faceted. The as-grown single crystal was sealed into a plastic package under vacuum for storage and transportation. The crystal was mounted onto an aluminium holder in a glove box under helium atmosphere. They were sealed into an aluminium can using indium wire in the glove box. The whole process was completed without exposure to air to minimize the contamination by humidity and oxygen. The details of the sample are given in Supplementary Supplementary  Table 1. The crystal was aligned at room temperature at DNA. It took 12 h to cool down to 140 K from room temperature at which it crystallizes in the tetragonal structure. At base temperature of 10 K, the width of the elastic peak is B0.0036 meV, which is equal to the instrumental resolution of DNA. The crystal is twined in the orthorhombic phase, as shown in Fig. 4e and Supplementary Fig. 11. The four-dimension S(q, E) data were reduced and visualized by using Utsusemi suite 63  MD simulations. MD simulations were performed using the massively parallelized LAMMPS package 65 and the recently developed classical MYP force field 34 , which reasonably reproduces the physical properties of CH 3 NH 3 PbI 3 (refs 38,66,67). In such a force field, non-bonded pairwise interactions are described by the sum of Buckingham, Lennard-Jones and long-range electrostatic potentials. Intra-and inter-molecular interactions of organic cations CH 3 NH 3 þ are described by the AMBER force field 68 . The cut-off distances for all the interactions are chosen to be 10 Å. A time step of 0.5 fs was used in all simulations. After energy minimization, simulation models were first equilibrated for 250 ps under the NPT ensemble with a constant pressure of 1 bar and at various temperatures (100, 140, 180 and 200 K). Then, they were simulated for another 50 ps under the NPT ensemble to calculate the DOS by taking the discrete Fourier transform of the mass-weighted VACFs of atoms, specifically given as 69 : where f is the phonon vibration frequency; v i (t) and v i (0) are the atomistic velocities of atom i at time t and 0, which are collected every 10 fs; m i is the mass of atom i; N is the total number of atoms belonging to certain species (CH 3 NH 3 þ or PbI 3 À ); t ¼ 50 ps is the auto-correlation time period used for data sampling; the angular brackets denote the time-averaged velocity auto-correlation functions. The Savitzky-Golay filter was employed to eliminate the thermal noise ( Supplementary Fig. 13).
Data availability. The data that support the findings of this study are available from the corresponding authors upon reasonable request.