Unraveling the Ordered Phase of the Quintessential Hybrid Perovskite MAPbI3─Thermophysics to the Rescue.

Hybrid perovskites continue to attract an enormous amount of attention, yet a robust microscopic picture of their different phases as well as the extent and nature of the disorder present remains elusive. Using specific-heat data along with high-resolution inelastic neutron scattering and ab initio modeling, we address this ongoing challenge for the case of the ordered phase of the quintessential hybrid-perovskite MAPbI3. At low temperatures, the specific heat of MAPbI3 reveals strong deviations from the Debye limit, a common feature of pure hybrid perovskites and their mixtures. Our thermophysical analysis demonstrates that the (otherwise ordered) structure around the organic moiety is characterized by a substantial lowering of the local symmetry relative to what can be inferred from crystallographic studies. The physical origin of the observed thermophysical anomalies is unequivocally linked to excitations of sub-terahertz optical phonons responsible for translational-librational distortions of the octahedral units.

In the light of the growing need for improved photovoltaic and optoelectronic materials, hybrid organic–inorganic perovskites (hereafter HOIPs) stand out as exceedingly promising candidates. They offer exceptional conversion efficiencies and coherent photocarrier transport typical of classic crystalline inorganic semiconductors,[1,2] yet at the same time they also exhibit rather unusual nuclear dynamics reminiscent of that seen in liquids.[3,4] While understanding their atomic structure and dynamics is critical to both rationalize their exceptional photovoltaic performance and improve their environmental and operational stability,[5] this task remains a formidable endeavor that has to date met with limited success.

Three-dimensional HOIPs commonly exhibit a classic ABX3 perovskite architecture, whereby a cuboctahedral cavity is occupied by an organic cation at the A-site. These materials reveal a rich and complex phase behavior largely dictated by cation mobility,[6] typically leading to the emergence of an ordered (orthorhombic) phase at low temperatures, followed by an intermediate tetragonal phase and a high-temperature cubic structure.[7] However, even in the case of prototypical methylammonium lead halides like MAPbI3 or MAPbBr3, a consistent picture of their short- and long-range structures remains a challenge to both experiment and computation. This is even the case well within the ordered phase at low temperatures, where the effects of dynamical disorder associated with either the organic cation or the surrounding inorganic framework are minimized.[8−13] Moreover, the A-site organic molecules and corner-linked PbX6 octahedra exhibit strong dynamical couplings giving rise to preferred short-range ordering of the organic cations, accompanied by distortions of the perovskite framework in the vicinity of the guest molecules. It is common for these local distortions to be observed by probes of local structure rather than by crystallographic techniques sensitive to the long-range order.[6,14] In addition to the above, a robust structural characterization of HOIPs is further hampered by their high propensity to exhibit crystal twinning[15,16] and nanoscale domains of different origins,[14,17−19] once more calling for the use of probes sensitive to local structure. On the experimental front, NMR and radiation-scattering techniques (photon, neutron) can be deployed in a complementary manner to obtain much-needed input into the above, focusing on either the organic or inorganic constituents.[20] At the same time, modern computational methodologies, particularly those relying on first-principles, have become increasingly successful at predicting the local atomic environment and its associated spectroscopic response, thereby enabling the development of models increasingly closer to experimental observation. In our previous works, such a combined experimental and computational strategy was applied for the first time to the case of the ordered phase of MAPbI3.[21] High-resolution inelastic neutron scattering (INS) along with first-principles modeling served to provide substantial evidence for a reduction of the local symmetry around the organic cation relative to the average Pnma crystallographic structure. In this case, the disruption of a fully fledged hydrogen-bonding network is accompanied by short-range octahedral-tilting distortions around the cations. In subsequent studies,[22,23] this improved model (hereafter denoted P1) was used with success in the study of both MA+ and formamidinium (FA+) cations in more complex MA1–FAPbI3 solid solutions. In particular, it was possible to quantify the extent of disorder as well as to shed new light onto the mechanism of physical stabilization of the perovskite framework in these materials. Building upon the above, this work introduces the use of specific heat data to scrutinize further the validity of these two models of MAPbI3. Using these thermophysical data alongside phonon calculations obtained with plane-wave density functional theory (PW-DFT), we introduce a joint experimental–computational protocol allowing for a quantitative validation of candidate structural models of HOIPs. Through its connection to the Vibrational Density of States (VDoS), our quantitative specific heat predictions are further used as a stringent test of the performance of modern density functional approximations (DFAs) in terms of their cost and accuracy to describe MAPbI3.

To set the scene, Figure provides a summary of the current state of affairs concerning the ordered structure of MAPbI3, based on new computational predictions using harmonic lattice dynamics (HLD) calculations and extensive benchmarking with several DFAs: PBE-TS,[24] PBE-D3(BJ),[25,26] PBEsol,[27] and rSCAN.[28,29] These calculations capitalize from our previous works[21,22] and those of Bokdam et al.[30,31] As shown by the INS data in Figure a, sharp spectral features reflect a well-defined local structure, with no discernible inhomogeneous broadenings that would otherwise signal the emergence of static disorder or the presence of distinct domains. Below 20 meV, these features are associated with librational modes of MA+, and they reflect a tendency for it to arrange itself in preferential directions. These observations are not compatible with the formation of an orientational glass.[22,23] Similar considerations apply to the sharp feature observed at ca. 37 meV, associated with a disrotatory internal mode of the cation, This particular mode is a very sensitive probe of the geometry of the (rather weak) N–H···I hydrogen bonds.[22] As extensively discussed in ref (22), the Pnma model cannot reproduce quantitatively a number of INS features, particularly those associated with librational motions of the organic cation. Absolute deviations from the INS data shown in Figure a amount to as much as 10 meV. Furthermore, these results depend far more weakly on the DFA employed, as highlighted in Figure b. In particular, the absence of a distinct trend in predicted transition energies between van der Waals (vdW) corrected and pure semilocal DFAs indicates that this type of interaction does not have a substantial impact on the vibrational properties of MAPbI3. This finding is in line with the results obtained by Pérez-Osorio et al.,[32] which further indicates that vibrational transition energies in the ordered phase of MAPbI3 are not strongly dependent on other factors such as anharmonicity or spin–orbit coupling. These considerations are at odds with those Zhang et al., suggesting that vdW interactions between the organic cation and its neighboring inorganic sublattice contribute to the total electronic energy between these two components by as much as 13.4%.[33] However, we note that this latter work was limited to the use of a highly specific, nonlocal vdW-DF2 functional, hardly benchmarked for hybrid perovskites. As such, further verification is in order. In the Pnma structure, MA+ electric dipoles are arranged in a head-to-tail fashion so as to maximize the number of hydrogen-bonding interactions with the surrounding iodine atoms. The observed overestimation in predicted transition energies relative to experiment thus indicates overbinding of MAPbI3 to the negatively charged inorganic framework.[23] A substantial reduction in the local symmetry around the organic cation leading to the P1 model becomes necessary to remove these discrepancies.[22] This structural model provides a superb description of vibrational transition energies, with absolute deviations below 2 meV (cf. the right panels in Figure b). The phonon band structures shown in Figure S1 provide additional insights into the performance of different DFAs, particularly at low energies. These data confirm that both structural models are mechanically stable at 0 K because no imaginary modes are observed. Interestingly, we find a high degree of consistency in the predicted phonon band structures for all DFAs considered. In particular, the phonon dispersion relations obtained with PBEsol are virtually identical with those obtained with the more elaborate PBE-D3(BJ) functional.[22] Furthermore, we find a similar performance of the celebrated meta-GGA SCAN functional in its numerically regularized incarnation (rSCAN), yet at a considerably higher cost by ca. 1 order of magnitude. In line with the more recent work of Lahnsteiner et al.,[31] we conclude that PBEsol emerges as the most versatile DFA: it is sufficiently cost-effective to describe phonon properties and finite-temperature cation dynamics with reasonable accuracy, and it represents a good alternative relative to the more demanding SCAN approach. On the basis of these results, the remainder of this work will focus on PBEsol outputs.

Figure 1

(a) INS spectrum of MAPbI3[22] along with characteristic modes depicted as colored insets. The lowest-lying molecular internal (disrotatory) mode is marked as τ(C–N). The highest- and lowest-energy molecular librations are marked as (i) and (vi), respectively. The most intense librational feature marked as (v) arises from the reorientation of the whole MA+ about the C–N axis. (b) Fully optimized DFT structures of MAPbI3 (left panels) according to Pnma (middle) and P1 (bottom) models. The two plots on the right show the corresponding errors in predicted transition energies for modes τ(C–N), (i), (v), and (vi), using several DFAs (see main text for further details).

Figures a and 2b show the phonon band structure for the structural models considered in this work. These predictions are also compared to the INS data of MAPbI3 reported by Ferreira and coauthors.[4] Given the large incoherent scattering from hydrogen, the INS data are primarily sensitive to the position of vibrational features associated with hydrogen motions in MAPbI3,[34,35] thereby limiting to some extent a direct comparison with the phonon band structure obtained in the calculations. Nonetheless, at a semiquantitative level, the lowering of symmetry in going from Pnma to P1 leads to an increase in spectral congestion between 2.5 and 4 meV which is more in accord with the experimental data. In this context, we note that the possibility of using perdeuterated specimens to circumvent this limitation is not a viable option, given the preponderance of isotope effects in this material.[36] Our calculated phonon band structures grant us, nonetheless, with the opportunity to provide quantitative predictions for the temperature dependence of the heat capacity within the harmonic approximation and to compare these with the experimental data (see also the Supporting Information). Specific predictions for the Pnma and P1 models using PBEsol are shown and compared with experimental data in Figure c,d over the range 2–75 K. Because the computed heat capacities have been obtained from the cumulative contributions of all normal modes, it is important to note that the present comparison is therefore performed in absolute terms, without any ad-hoc normalization of computed thermophysical data. In line with previous works,[37] we make use of a “Debye-reduced” representation of the form C(T)/T3 to highlight departures from Debye’s law. Additional calculations using a wider set of DFAs are further presented in Figure S2, yet we anticipate that our conclusions as explained below are not dependent on the specific DFA used to obtain the phonon band structure. The experimental specific heat data shown in Figure are quantitatively in line with the results of Fabini et al.,[38] exhibiting a broad peak centered at 5.9 K in C(T)/T3. Further inspection of these data also shows that this low-temperature feature can be clearly reproduced by our HLD calculations, without any assumption on glass formation or disorder. It is, instead, an intrinsic feature of a fully ordered structure, arising from the presence of a dense manifold of normal modes in the VDoS in the millielectronvolts range. The P1 model also provides a better description of the experimental data both above and below the maximum in C(T)/T3 observed in the laboratory. Qualitatively, this closer agreement arises from a broader distribution of low-energy phonon features in the VDoS, which is also visible in the neutron data in Figure a,b as a result of the lower symmetry of P1 relative to Pnma. These considerations are shared by all the DFAs tested in this work, as reported in Figure S2.

Figure 2

(left panels) Calculated phonon dispersion relations for the Pnma (a) and P1 (b) models. The resulting phonon branches are presented in the range 0–4.5 meV and overlaid on the color maps corresponding to INS data from a hydrogenous single-crystal specimen.[4] (right panels) Experimental (black dots) and theoretical (red dots) C(T)/T3. Partial contributions to this quantity from different phonon energy intervals are represented by the shaded curves, as defined in panel c. For the P1 model, these contributions are decomposed further by using 0.5 meV intervals, as shown by the broken lines in (d).

To scrutinize the different phonon contributions further, we have decomposed the predicted C(T)/T3 into partial contributions over selected energy intervals (see Figure c,d). With the benefit of hindsight, we note that the thermophysical data exhibit a good and (largely unforeseen) degree of “resolving power” in the energy domain, via its mapping onto the temperature axis. This dissection exercise also tells us that the dominant contributions come from excitations of phonons at energies of 1.5–4.5 meV, i.e., predominantly optical lattice modes. The INS data in this regime (see Figures and S1) also show an increase in the underlying VDoS at around 2, 3, and 4 meV, and these modes are characterized by a considerable projection onto the hydrogen atoms, reflecting a strong coupling between the cation and the inorganic framework. The Γ-point optical modes involving the highest contributions to the INS intensity are shown in Figure . Interestingly, these modes are linked to distortions of the P1 structure ultimately leading to Pnma. The lowest-energy contributions at around 2 meV can be ascribed to octahedral distortions driven by Pb displacements, which affect the mirror-plane symmetry. The contributions at around 3 and 4 meV are driven by I atom dynamics, reflecting librations and shearing of the PbI6 octahedra, respectively. This quantitative analysis provides further and rather unequivocal evidence on the specific physical origin of the strong departure from Debye’s law, which is related to the excitation of the lowest-energy optical phonons. This picture does not require the presumption of any static disorder, and therefore we anticipate that it will be of importance and relevance in the context of identifying its possible emergence in HOIPs with more complex cation compositions.

Figure 3

Representative Γ-point optical modes of MAPbI3 according to HLD PBEsol calculations using the P1 model. These modes correspond to the most prominent contributions from the hydrogen atoms. Modes (a) and (b) are the first and the third lowest-energy optical modes in the range 2–3 meV, associated with octahedral shearing. Mode (c) is a representative octahedral distortion contributing to the dispersion branch observed at ca. 4 meV.

As reported in Figure , our results can also be used to provide a quantitative assessment of the sensitivity of the heat capacity as an experimental observable to discriminate among plausible structural models as well as different DFAs. As defined in the Supporting Information, we have performed this exercise by looking at the differences between experiment and computational predictions both in terms of their temperature dependence and their associated (cumulative) mean absolute errors (MAEs) (see Figures a and 4b, respectively). Figure a displays the temperature dependence of these deviations. For Pnma, these are as high as 0.6–0.7 K–3 and are least pronounced for PBE-D3(BJ). PBEsol represents an interesting case with positive (negative) deviations of a similar overall magnitude below (above) the peak. This behavior indicates an overall and distinct shift of the computational predictions relative to observation. This shift is primarily caused by an overestimate of the number of normal modes in the range 3–4.5 meV. For the case of P1, these deviations are reduced significantly over the entire temperature range and are predominantly positive. For PBE-D3(BJ), they are below 0.2–0.3 K–3, about a 10% underestimate of the experimental values around the maximum. These considerations are further reflected in the MAEs shown in Figure b. The Pnma model performs consistently worse than P1 for all DFAs considered, by a factor of at least 1.5 in terms of the MAE. For the particular case of MAPbI3, the consistency of this result across all DFAs also tells us that these have reached a sufficiently high level of accuracy, thereby enabling a meaningful comparison between competing structural models.

Figure 4

Quantitative comparison between predictions for the Pnma and P1 models. (a) Deviations from experiment using the DFAs indicated in the figure legends. For further details, see the main text. (b) Mean absolute errors (MAEs) in the predicted C(T)/T3 using different DFAs. For details on the mathematical definitions of these deviations and MAEs, see the Supporting Information.

In summary, we have shown that a detailed scrutiny of heat capacity data constitutes a powerful approach to unravel pending questions concerning the properties of this important material. Our results also serve to demonstrate that a pseudo-orthorhombic P1 structure previously proposed on the basis of INS data[21] continues to provide an improved description of the structure of the ordered phase of MAPbI3 at ambient pressure, relative to what has been inferred from crystallographic studies to date. We underline that this gratifying conclusion has been reached in an entirely independent manner relative to previous studies by focusing on how low-energy modes primarily associated with the inorganic framework (and not the organic cation) manifest themselves in the low-temperature behavior of a hitherto unexploited experimental observable—the heat capacity. The above program of work has been certainly facilitated by the level of quality of existing DFAs, making it possible to survey their relative performance and to ascertain that their potential limitations are now sufficiently minor to enable robust and quantitative model selection. To the best of our knowledge, this is the first time that such model selection protocol has been implemented with success in MAPbI3 or related materials. We also anticipate that its realm of applicability could be further extended to the study of other more complex HOIPs and their mixtures,[39−41] where other experimental probes might also fall quite short at providing the requisite level of physical insight.

Experimental and Computational Details. Specific heat measurements were performed across the low-temperature phase of MAPbI3 down to 2 K by using a Quantum Design PPMS system. The measurements were performed on a pressed powder pellet fixed with Apiezon-N grease by using the same hydrogenous MAPbI3 sample as characterized elsewhere.[21] The calculations of phonon dispersion relations were performed with CASTEP[42] by using the finite-displacement method combined with a nondiagonal supercell approach to solve for the dynamical matrix.[43] Further computational details are provided in the Supporting Information.