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.


S1. Computational Details
The Plane-Wave Pseudo-Potential (PW-PP) formulation of Density Functional Theory (DFT) under Periodic Boundary Conditions (PBCs) was used along with the developer´s version of the CASTEP code. 1,2 We used three different semi-local Density Functional Theory Approximations (DFAs), namely, two exchange-correlation functionals within the Generalized-Gradient-Approximation (GGA), i.e. the original incarnation of the Perdew-Burke-Ernzerhof (PBE) 3 functional along with its solid-state version, PBEsol. 4 Furthermore, we employed the regularized version (rSCAN) 5 of the meta-GGA-type SCAN functional. 6 The core electrons were described by a set of hard norm-conserving PPs, while the electronic wave functions were defined using a PW basis set with a kinetic energy cutoff of 900 eV. The pseudopotentials were generated on-the-fly using the parent functional in all cases but in the rSCAN calculations, where PBE pseudopotentials were used instead. For the PBE functional we adopted two different atom-pair-wise approximations to account for van-der-Waals (vdW) T he numerical settings were the same as those used in our previous studies on MAPbI 3

. 11,12
A Monkhorst-Pack grid was used to maintain a constant k-spacing of 0.05Å -1 . All structures were accurately relaxed at atmospheric pressure to minimize residual atomic forces.
The convergence criteria in the variation of the total energy, Hellmann-Feynman forces, external stress, maximum displacement, and the self-consistent field (SCF) were defined as 1 × 10 -12 eV/atom, 1 × 10 -5 eV/A, 0.0001 GPa, 1 × 10 -6˚A , and 1 × 10 -12 eV/atom, respectively. Following geometry optimization, phonon frequencies and eigenvectors were calculated from dynamical matrices constructed by numerical differentiation of the analytical gradients with respect to atomic displacements. For the finite-difference calculations, we used a displacement amplitude of 0.01Å. The non-diagonal supercell method by Lloyd-Williams and Monserrat was employed to reduce the size of the supercell required to obtain the force constants. 13 The calculated phonon band structure gives us the harmonic energy E kj of a given mode j for a given wave-vector k. The Vibrational Density of States (VDoS) is then an explicit function of energy, defined as the number of modes per unit energy.
Likewise, the heat capacity at constant volume is given by the following double-summation where T is the temperature and k B is Boltzmann´s constant. As written, C V is an adimensional quantity that is equal to zero at T =0 and is bounded from above by the Dulong-Petit limit at high temperature. This summation can be restricted arbitrarily, in order to obtain partial contributions to the heat capacity over a given energy range, as illustrated in Figs. 2c and 2d in the main text.
The experimentally accessible specific-heat capacity is obtained at constant pressure, C p (T ). C p (T ) is related to C V (T ) by the following equation: 14 S3 where B 0 (T ) is the bulk modulus, V (T ) is the volume at temperature T, and α V is the volumetric thermal expansion coefficient. The second term proportional to temperature on the right-hand side of Equation (3) is the so-called lattice-dilatation term, which is negligible at the low-temperature conditions of relevance to the present work. Hence, C p (T ) ≈ C V (T ).
To compare with computational predictions, the experimentally obtained molar C p (T ) has been normalized by the universal gas constant R.

S2. Deviations and Statistical Methods
First, a spline interpolation was performed between the (denser) computational data sets obtained for each DFA. This operation allowed us to compare the experimental measurements and the computational predictions at the same temperatures. Consequently, the differences ∆C p /T 3 were obtained as: where y i and x i stand for the predicted and the experimental values, respectively.
Finally, the (cumulative) Mean Absolute Errors (MAE) in the prediction of the 'Debyereduced' heat-capacity values were computed using: where n is the size of the data set (in this case, the number of experimental points).

S3. Additional Computational Results
Figure S1: Phonon-dispersion relations calculated for the Pnma and P1 models of the ordered phase of MAPbI 3 , using a range of DFAs.
x x x x x x x x x xx x x x xx x x x xx x x x xx x x x x S5 Figure S2: Experimental (black dots) and computational (colored empty circles) thermophysical data: C p (left panels) and C p (T )/T 3 (right panels). The calculations have been performed using the DFAs reported in Fig. S1.