Molecular disorder and translation/rotation coupling in the plastic crystal phase of hybrid perovskites

The complexity of hybrid organic perovskites calls for an innovative theoretical view that combines usually disconnected concepts in order to achieve a comprehensive picture: (i) the intended applications of this class of materials are currently in the realm of conventional semiconductors, which reveal the key desired properties for the design of e ﬃ cient devices. (ii) The reorientational dynamics of the organic component resembles that observed in plastic crystals, therefore requiring a stochastic treatment that can be done in terms of pseudospins and rotator functions. (iii) The overall structural similarity with all inorganic perovskites suggests the use of the high temperature pseudo cubic phase as the reference platform on which further re ﬁ nements can be built. In this paper we combine the existing knowledge on these three ﬁ elds to de ﬁ ne a general scenario based on which we can continue the quest towards a fundamental understanding of hybrid organic perovskites. With the introduction of group theory as the main tool to rationalize the di ﬀ erent ideas and with the help of molecular dynamics simulations, several experimentally observed properties are naturally explained with possible suggestions for future work.


Introduction
Having demonstrated the first truly low-cost photovoltaic devices with more than 20% conversion efficiency, hybrid organic perovskites (HOP) are currently a phenomenal scientific and technological revolution for the entire sector. [1][2][3][4][5][6][7] Being brought back to the forefront among the current hot topics, they have also revealed undeniable assets for light emitters 8 and water splitting applications 9 along with potential outcomes in spintronics. 10 Meanwhile and despite a tremendous effort of the DFT community, there is still need for further theoretical work in connection with both fundamental and applied issues. In this context, there is hope for a theoretical framework and generic tools to streamline the analysis of recent but also older experimental and computational data on HOP. To achieve these goals, we believe that HOP must be placed in context and viewed along with other classes of materials (Fig. 1). This had made it possible to exploit the concept and tools developed over decades for conventional semiconductors 11 and to identify key differences such as reverse band edge electronic state ordering. 12 This line of comparison also points out that the performance of devices depends heavily on the quality of the material, for instance the degree of crystallinity in silicon-based solar cells. Given the inherent complexity of HOP as compared to silicon, it also raises the question of a bulk crystal structure, related phase transitions and their impact on the physical properties. 13 In this regard, a comparison with all inorganic perovskites (AIP) has already yielded interesting results, among which is the proposal to consider their cubic high temperature phase as the reference phase for HOP. 14 This prompts the question: to what extent the hybrid character is important for the intended applications? In other words, how it affects materials growth and properties and how it can be accounted for in theoretical approaches or numerical modeling. Regarding the latter, given the symmetry of possible organic cations that fit in the cuboctahedral cage of the 3D inorganic framework, the high temperature reference phase must be thought as pseudo-cubic with static or dynamical orientational disorder of the organic moieties. From early investigations on methylammonium (CH 3 NH 3 + , hereinafter referred to as MA) lead-based perovskites, we know that the organic cation is tumbling between various orientations with a typical timescale of picoseconds. This led to the suggestion that, in HOP, excitons may undergo a dielectric screening due to the rotational degrees of freedom of the MA molecule, which has recently been confirmed experimentally. [14][15][16] From the device point of view, this promotes improved transport properties based on free carriers. 17,18 Such an orientational degree of freedom is also reminiscent of plastic crystals. 14,19,20 Indeed, liquid and plastic crystal phases are possible intermediate stages between ordered crystal and liquid phases. Liquid phases exhibit both translational and orientational/conformational disorder, whereas in a plastic crystal phase atoms and the molecular centers of mass (CM) are located at well-defined positions within the periodic lattice. If the disorder is mainly related to molecular rotations rather than molecular conformations, the phase is frequently referred to as a "rotor phase". In that sense, there is clear experimental evidence of the rotational motion of the organic cations in the high temperature pseudo-cubic phase of HOP, namely from NMR, dielectric measurements, millimeter wave spectroscopy and neutron diffraction. [21][22][23][24] Although a couple of scenarios for the orientation of the molecular cation are considered, based on the early work by Onoda et al., 22 there is yet no theoretical framework that provides an in-depth uniform analyses.
In this paper, after a brief introduction to AIP and HOP, we will first revisit the available experimental data related to phonons in AIP and put them into perspective with HOP for which similar data is still lacking. Next, we will subsequently develop two approaches that allow taking the stochastic molecular degrees of freedom into account. The first simplified approach is based on the concept of pseudospins (PS), recently introduced by one of us, 19 which offers a discrete representation of possible random orientations of the molecular cation accounting for both its elastic and electric multipoles. We then examined how rotator functions, which provide continuous variables, can contribute to rationalize the properties of HOP. To complement the available experimental data, we performed molecular dynamics (MD) simulations at several temperatures and for several super-cells. This work presents for the first time MD simulations [25][26][27][28] combined with group theory analysis 19 applied to HOP. All these different approaches, when implemented with the support of the group theory, are shown to help streamline the couplings between stochastic molecular excitations and lattice vibrations, the structural phase transitions and possible orientational glassy behavior of perovskite structures.

Computational method
First principles molecular dynamics simulations were performed using the CP2K package, which implements a hybrid methodology combining Gaussian basis functions and plane waves. [29][30][31] All simulations were performed under NVT ensemble using periodic boundary conditions with the temperature being controlled by using a Nosé-Hoover thermostat with three chains and a time constant of 50 fs. 32,33 The electronic structure properties were calculated using the PBE functional with the addition of the Grimme D3 scheme to correct the dispersion interactions. 34,35 The Kohn-Sham orbitals were expanded in Gaussian basis sets (DZVP-MOLOPT for Pb, I, C, N and H) 36 and the norm-conservative GTH pseudopotentials were employed. 37,38 The auxiliary plane-wave basis set was defined by an energy cutoff of 300 Ry. The time-step for the integration of the dynamic equation was set to 1 fs.
The preparation of the systems was made using the available crystallographic data, obtained at 400 K for the cubic phase of MAPbI 3 . Different systems were prepared by the replication of the unit cell for a larger quasi-cubic simulation box. System-222 (96 atoms) contains 2 unit cells in each Cartesian direction. Analogously, larger systems having up to 5 unit cells per direction were prepared. The largest case, System-555, contains 1500 atoms. All the systems were tempered by a succession of NVE simulations to introduce a representative thermal disorder before starting the NVT simulations, which were integrated for at least 40 ps in each case.

Results and discussion
3.1. The Pm3m crystal structure as a reference phase The Pm3m cubic phase of CsPbX 3 compounds (Z = 1) is chosen as the reference for the high temperature phase of 3D AIP and HOP with a general formula, AMX 3 . [39][40][41][42][43][44][45][46] The metal (M), the 3 halogens (X) atoms and the cation (A) are located at sites (0,0,0) (M, 1a), (1/2,0,0) (X, 3d) and (1/2,1/2,1/2) (A, 1b), respectively (Fig. 2a). Structural data is available for 3D AIP with A = Cs + or Rb + . In 3D HOP, the orientations of the organic cations are dynamically averaged in the Pm3m cubic phase. The molecular CM is located at the site (A, 1b) of the cubic cell (Fig. 2a). Most of the photovoltaic devices made of perovskites have been developed with the MA cation, to a lesser extent with formamidinium (FA: [HC(NH 2 ) 2 ] + ), but guanidinium (GA: [C(NH 2 ) 3 ] + ), dimethylammonium (DMA), trimethylammonium (TMA) and tetramethylammonium (TeTMA) are also possible. [47][48][49][50] Here we will mainly discuss MA-based 3D HOP, but the methodology can be applied for other organic cations having a different geometry, symmetry or polarity. Noteworthily, the replacement of the molecular CM by a simple cation (Na + or Cs + ) mimics the influence of ionic interactions on the electronic band diagrams close to the band gap. [12][13][14] However, it is necessary to go beyond this simple model in order to describe the crystal packing and structural instabilities, as well as the complex coupling between the disordered orientational dynamics and lattice vibrations. The first step of our approach consists of decoupling translational and rotational degrees of freedom of the cation A. These degrees of freedom can be analyzed separately yielding, for the translational part, a common vibrational analysis for both AIP and HOP bulk solids. The rotational part must be treated carefully to account for the possible disorder effect and will be investigated in a second step. Then rotation-translation coupling allowed by symmetry in the cubic phase will be considered in a third step. Our analysis strongly relies on the use of group theory 39,40 (see for example the pedestrian guide in ref. 19) and reciprocal space. 11,43 The labels of the main high symmetry points related to the first Brillouin Zone (BZ) of the Pm3m space group are shown in Fig. 2c. For the irreducible representation (IR), we will mainly use Miller and Love notations. 19,40 In order to analyze the symmetry properties of lattice vibrations ( phonons) belonging to a Pm3m structure, 5 atoms (among which 1 is a pseudo atom for the molecular CM in HOP) must be considered. 11,19,42  modes, but only the X atoms. These modes related to antiferro-distortive instabilities correspond to the rotations of the metal halide octahedra. The IR decomposition at Γ are useful to analyze electron to optical phonon coupling. For electron to Fig. 2 (a) Simple cubic structure of CH 3 NH 3 PbBr 3 at high temperature. 45 The center of mass (CM) of the CH 3 NH 3 (MA) cation is represented by the large green sphere, and lead atoms by grey spheres. Bromide atoms exhibit anisotropically disordered motions as indicated by the brown ellipsoids. Site labels are also indicated. (b) NaCl-type structure of NaCN at high temperature. The sodium cations are represented by large yellow spheres, the CM of cyanides by small light green spheres. (c) Reciprocal space 3D view of the brillouin zone (BZ) with conventional labels for points of high symmetry (Γ, R, M, X).

Phonon dynamics in AIP
Phonon dynamics has already been thoroughly investigated only in the CsPbCl 3 AIP, especially in its cubic phase that is the most relevant to this work. For instance, an empirical force constant model is in good agreement with the available spectroscopic data and phonon dispersion curves. 51,52 Neutron scattering measurements clearly indicate that the energy of the whole acoustic phonon density of states is remarkably low at 80°C, as compared to other perovskite crystals such as SrTiO 3 or KMnF 3 . According to the direct measurements from Γ to R, M and X points show that the transverse acoustic (TA) branches are limited to about 2 meV over the entire BZ. 51 The optical soft phonon branches, which lie in the same energy range (<2 meV) and are related to the structural instabilities close to the R and M points, are heavily overdamped over half the BZ. It was suggested that anharmonicity plays an important role in the phonon dynamics, as also revealed by the unusually large Debye-Waller factors. Similar observations were reported for the high temperature cubic phase of CsPbI 3 , which exhibits significant anisotropy of the halogen thermal ellipsoids similarly to that of CsPbCl 3 (Fig. 2a). 41,53 The halogen thermal motion is consistent with a small bulging of its distribution that is related to thermal anharmonicity. The central peak observed by Raman scattering in the cubic phase of CsPbCl 3 was originally attributed to order-disorder longrange correlations in the high temperature phases and fitted by a Debye-type relaxor function: I / γ ω 2 þ γ 2 , where γ = 1/τ is the relaxation rate and τ the relaxation time. 54 However, it corresponds to the signature of the strongly overdamped phonon observed by neutron scattering at M and R. 51 In fact, the scattering cross-section for a phonon is given by the phonon frequency and damping. 55,56 For an overdamped phonon, Γ q ≫ ω q , the neutron scattering cross section exhibits a signature similar to that of a relaxor, Sðq; ωÞ / 1=τ q ω 2 þ ð1=τ q Þ 2 , with 1/τ q = ω q 2 Γ q . At high temperature, a phonon signature and a central peak with a true relaxor dynamics may be simultaneously observed by neutron scattering in many materials. [55][56][57][58] The simultaneous observation of the two signatures corresponds to a more complex dynamical susceptibility: , where the phonon (ω q ,Γ q ) is linearly coupled to a true relaxor through the coupling constant δ. When the two signatures lie in different frequency ranges, the neutron scattering cross-section is the sum of two contributions whose intensities are related: and optical phonon ω ∞ 2 = ω q 2 + δ 2 .
Such a refinement was not necessary for CsPbCl 3 , as it was for KMnF 3 or SrTiO 3 , and the neutron scattering data were only described with an overdamped phonon model by the same authors. 51,58 All the phonon peaks are strongly broadened in the high temperature cubic phase. In fact, the high frequency optical phonons (>5 meV) are not observable in CsPbCl 3 by neutron scattering due to broadening. Thus, from the inspection of the available literature on 3D AIP, we may conclude that the phonon spectrum in the pseudocubic phase is split in two parts: (i) heavily damped optical phonons, which are overdamped very close to the M, R and X edges of the BZ, and (ii) acoustic phonons at very low energy (<2 meV). The former contribution corresponds to pretransitional thermal fluctuations. The large damping of phonons observed in 3D AIP can be explained by a strong anharmonicity, which particularly shows off for the thermal vibrational motions of halogen atoms and leads to a relaxor-like signature close to the M and R points of the BZ. However, there is no experimental evidence of an additional quasi-elastic central peak like in KMnF 3 or SrTiO 3 perovskites. Unfortunately, similar detailed phonon spectroscopic data and vibrational models are still lacking for 3D HOP. From the available experimental reports on AIP and HOP, we may only infer that the spectra of optical modes lie in the same energy ranges. 59 Diffraction studies on HOP crystals also report on strongly anisotropic thermal ellipsoids for halide atoms in HOP (Fig. 2a), which can be assigned to a strong anharmonicity. 45,46 The MA cation reorientations 21,23,24 have been investigated recently by inelastic incoherent neutron scattering in MAPbI 3 , 60,61 exhibiting two quasi-elastic responses based on the same technique as used for phonons in AIP. The last neutron scattering investigation 61 (which appeared while this manuscript was at the revision stage) contains a symmetry analysis of the tetragonal phase related to the approach developed in the present work (vide infra). The MA tumbling dynamics is described by discrete jumps between the four equivalent positions in the tetragonal phase, and a fast dynamics around the C-N axis, which is still present in the low temperature orthorhombic phase. 61 A recent neutron scattering study of MAPbBr 3 shows that the low-temperature harmonic ( phonon) fluctuations cross over to fast overdamped relaxational dynamics and coupling to molecular dynamics at higher temperatures. 62 The random reorientations at high temperature thus deserve a specific theoretical framework, which is discussed in the next section. ordered lattice with translational symmetry. 11,20,63,64 The orientations of the molecules undergo at high temperature, random dynamical changes between preferential directions. This random reorientational motion is thus considered as "stochastic", 65 by contrast to the "deterministic" ( predictable) motion of an atom around an equilibrium position in a phonon mode. The PS concept is an approximation, leading to a discretization of the angular variables, keeping only the preferential molecular orientations with maximum occupancies. 20 It is very useful for lattice symmetry analyses, as well as analogies with discrete Ising models. The continuous rotator functions are more general and allow a complete crossover from a PS description to free molecular rotational motion. The discreteness of the angular distribution can be quantified by expanding the symmetric rotator functions in terms of cubic harmonic functions in the case of a cubic lattice. All these concepts will be examined for the high temperature cubic reference phase of hybrid perovskites (vide infra). When plastic crystals are cooled, the rotational motion freezes and an orientational order appears at a certain temperature. Plastic crystals may exhibit two different behaviors at low temperature. In the first case, the rotational motion freezes under careful cooling, keeping quasistatic random molecular orientations. These orientational glasses exhibit most of the properties related to conventional glasses such as specific heat anomalies, ageing and distributions of relaxation times. Prototypes of orientational glasses are cyclohexanol 64 and cyanoadamantane. 66 In the second case, the plastic crystal may undergo a structural transition at low temperature to a partially or completely ordered crystalline phase, if spatial correlations or interactions between the molecular orientations are strong enough. Noteworthily, these interactions can be direct or indirect (vide infra). Prototypes of plastic crystals exhibiting structural phase transitions are the NaCN and KCN compounds. 20 The structural transitions are usually associated to a structural symmetry breaking related to a collective molecular freezing and to an order-disorder character. Displacive and order-disorder dynamical processes are the two limiting cases for the microscopic mechanism driving a continuous structural phase transition. The order parameter dynamics is critical at the transition temperature and for a specific wavevector. For a soft phonon mode and a second order phase transition, the frequency goes to zero at T c , 67 whereas for a critical order-disorder relaxation, a similar behavior is observed for the relaxation rate. 67,68 The high temperature relaxation mechanism in a plastic crystal is not necessarily critical at the structural phase transition corresponding to the molecular freezing. The order parameter of the structural phase transition can be related to a soft phonon mode of the surrounding lattice. The variable (PS or rotator function) describing the stochastic orientational motion is a secondary order-parameter, linearly coupled to the main displacive order parameter, which leads to a long range ordering of the frozen orientations at low temperature. Long-range correlations may be destroyed in alloys such as (NaCN) 1−x (KCN) x , inducing an orientational glass transition at low temperature instead of a structural phase transition. 69 HOP compounds undergo structural transitions to low temperature phases, where the molecular orientations are fixed. The existence of stochastic orientational degrees of freedom at high temperature can be assessed from the experimental results reported in the literature. Onoda and coworkers first demonstrated that the specific heat anomalies and entropy contributions at the structural phase transition temperatures are characteristic of order-disorder processes. 22 Kawamura and coworkers showed that the cubic to tetragonal transition, in MAPbI 3 , is associated to a partial and collective molecular freezing, together with a displacive-like antiferrodistortive rotation of the metal halide octahedra. 70 Moreover, Debye-type functions, which are characteristic of relaxational motions, were used to interpret nicely, both dielectric and millimeter-wave spectroscopy results in the three MAPbX 3 compounds. 21,23 In the prototype MAPbI 3 HOP, the quasielastic neutron scattering signatures 60,61 are also consistent with the simultaneous fast rotations of the methyl and ammonium groups around the C-N axes 24 and the slower tumbling of the C-N axis between the symmetry equivalent positions. 21,23 This is described in detail in ref. 61 with a C 4 × C 3 symmetry-based jump model designed for the tetragonal phase. The same model is used for the cubic phase, although as quoted by the authors, a O × C 3 model expected to be more correct, yields essentially the same result. 61 A reasonable description of the experimental results is also afforded by an isotropic orientational distribution. 61 Both molecular motions are related to entropy contributions at the transition temperatures, 22 but this work focuses on the tumbling of the molecular axis, namely the C-N axis in MA-based HOP, to discuss PS or rotator functions. The existence of a stochastic reorientational motion, leading to a collective freezing at low temperature, and thus a plastic crystal behavior, seems to be consistent with all available experimental data on MAPbX 3 compounds. However the nature of the critical mechanism (displacive or order-disorder) is still not clear from the inspection of the literature. This point will be discussed further in section 3.4.4.
3.3.2. Molecular pseudospins (PS) in hybrid perovskites. This approximate description is able to capture more easily the symmetry properties as well as the translation-rotation coupling or the order-disorder character of phase transitions. In a local cubic environment, the sets of symmetry equivalent molecular orientations can be classified into six scenarios: {1,0,0} (6, A); {1,1,0} (12, B); {1,1,1} (8, C); {a,b,0} (24, D); {a,a,b} (24, E); {a,b,c} (48, F), with corresponding multiplicity (dimension of the point group reducible representation) and scenario labels given in parentheses. In scenarios A, B, and C, the molecular dipole points toward a cubic cell facet, a halogen atom, and the center of an octahedron, respectively. Scenarios D, E and F are mixed, between A and B, B and C, and A, B and C, respectively. In fact, scenarios A, B and C correspond to the simplest sets of unit vectors. They have been suggested earlier to rationalize the influence of molecular disorder in the Pm3m cubic phases of MAPbX 3 . 21,23 Unfortunately, the available experimental data are insufficient to definitely conclude about the prevailing scenario. The {1,0,0} direction (scenario A) appears among the preferred MA orientations: MA cations are considered to be mainly oriented toward the faces of the cube where no atoms are located, but probably also partly toward the halogen anions at the cell edges or the lead atoms in the corner, that is scenario B. 45 The concept of PS is also useful to rationalize properties of other crystals. Table 2 gathers a few examples of plastic crystals showing orientational disorder that can be described by PS in a high temperature cubic phase. 20 NH 4 Br and MAPbX 3 are two examples for a Pm3m lattice, whereas NaCN corresponds to a face-centered Fm3m cubic lattice. NaCN crystallizes in a NaCltype structure (Fig. 2b), where the CN − disordered anion occupies the (1/2,1/2,1/2) (A, 1b) site in the cubic cell. The first, second and third nearest neighbors of CN − are Na + , CN − and Na + , respectively. For this reason preferred CN − orientations correspond to {1,0,0} and, to a less extent, {1,1,1} directions. In addition a CN − head to tail disorder is usually assumed, approximating the molecule to a symmetric dumbell. 20 Finally, in the NH 4 Br crystal, the tetrahedral molecule in the center has eight nearest neighbor Br anions in the corners. The protons point toward the corners of the cube, which define only two preferred orientations. The NH 4 Br example further shows that besides the site symmetry, the molecular symmetry is also important.
In the PS simplified picture, a discrete variable ξ i is defined for each possible molecular orientation i. ξ i is equal to 1 or 0, whether orientation i is occupied or not. Symmetry-adapted combinations of ξ i can be computed for a single crystal cell: (1,0,0), (0,1,0), (0,0,1), (−1,0,0), (0,−1,0) and (0,0,−1) orientations respectively. The symmetry-adapted combinations decompose as Γ 1 + + Γ 3 + + Γ 4 − (Table 1): The totally symmetric combination Ξ Γ 1 + is always equal to 1, whatever the value taken by ξ i . In order to better understand the significance of these symmetrized combinations, we may consider the average values of the ξ i in relation to the orientation occupation numbers. 19,71 If one considers a thermally activated disorder of the molecular orientations, the dynamically averaged value, compatible with lattice symmetry, is <ξ i > = 1/6. It yields 〈Ξ Γ 1 +〉 = 1, and 0 for the other symmetry-adapted combinations. This can be considered as the ground state of the system in scenario A. The Γ 3 + and Γ 4 − IR are relevant to describe dynamically excited configurations corresponding to thermally activated population fluc- symmetrized combination corresponds to δξ 1 = −δξ 4 and δξ 2 = δξ 3 = δξ 5 = δξ 6 = 0, that is, to a net electric dipole along the x axis. In turn, all the other symmetrized excited combinations yield 0. The PS concept is able to account for both dynamic and static orientational disorders. Indeed, a molecule with a frozen (static) random orientation and pointing toward a face (ξ 1 = 1 for example) corresponds to a superposition of where the 4 terms describe respectively, the crystal field potential felt by one molecule in a rigid lattice with the orientations of all the other molecules averaged, the CM vibrational potential energy, the translation-rotation coupling potential and the direct interaction between the orientations of pairs of molecules. Notice that the translation-rotation coupling potential V TR is at the origin of an indirect interaction between the orientations of pairs of molecules. The model proposed by Poglitsch et al. 21 describes the dynamics of individual molecules without lateral coupling, thus implicitly retaining only a few terms in the crystal potential expansion: V R , the single particle susceptibility connected to V RR and a self-energy term derived from the indirect interaction between the molecular orientations (vide infra).
The relaxation rate γ 1 can be measured by experimental techniques related to electric dipoles such as dielectric measurements, 21,23 whereas γ 2 can be probed by experimental techniques relevant to elastic multipoles, like Raman scattering. 72 Both types of relaxation are expected to yield an experimental signature by incoherent neutron scattering (vide infra), 60,61 the best technique being triple axis neutron scattering, but large monocrystals are necessary. 55,56 Noteworthily, more realistic modelling of the interactions between neighboring cells (vide infra section 3.4.2) may change the γ 2 /γ 1 ratio equal to 1.5 in this simplified model. Correlations between the tumbling of cations in neighboring cells also lead to dispersions of relaxation rates. In that case, it is useful to define collective PS variables accounting for the space group symmetry: wave vector in the BZ, n is for a given molecule at the position X n and N is the total number of molecules. The symmetrized decompositions for these collective variables are given in Table 1 as Γ, R and M for scenarios A-D. Finally, we may consider that the orientational disorder in the Pm3m cubic phases of MAPbX 3 is a superposition of the three scenarios A, B and C. Then, the total number of cells is , where N A , N B , and N C are the number of cells undergoing, respectively, A, B, and C type disorder, respectively. But, as it turns out, such a mixing is better described by rotator functions that will be considered in the next section.
In the cubic phase of MAPbI 3 at T = 370 K, the experimental analysis of the tumbling motion is difficult. 61 A fully isotropic random rotation model compatible with the lattice symmetry was used. This continuous model further assumes no correlation between the rotations of the neighbouring molecules. It yields a relaxation rate equal to γ l ¼ l l þ 1 ð Þ 2τ for the rotation with an angular momentum l (τ = 1.62 ps at T = 370 K). According to this model γ 2 /γ 1 is thus equal to 3. By comparison γ 2 /γ 1 is equal to 1.5 and 2, respectively, in the Poglitsch discrete jump models for scenarios A (Γ 3 + /Γ 4 − , see the above derivation) and C (Γ 5 + /Γ 4 − , Table 1). For scenario B, two IR (Γ 3 + and Γ 5 + see Table 1) correspond to an angular momentum l = 2, while only Γ 4 − is relevant for l = 1. γ 2 /γ 1 is equal to 2 and 3 in a Poglitsch discrete jump model, respectively, for Γ 3 + /Γ 4 − and Γ 5 + /Γ 4 − .

Molecular rotator functions in hybrid perovskites.
Rotator functions are symmetry adapted functions that change continuously when varying the orientation. 20 In the case of linear molecules, like the CN − anions in NaCN, they are given by combinations of spherical harmonics. These simple rotator functions are also well adapted to the tumbling motion of the C-N axis in MAPbX 3 materials. Inclusion of other degrees of freedom, like the methyl or ammonium rotations around the C-N axis could be further implemented by means of generalized expansions of the rotator functions in terms of Wigner rotation matrices.
Considering a O h molecular site-symmetry, 19 the probability distribution of the C-N axis orientation as a function of the two polar angles (ϑ,φ), i.e. the unit vector u = (cos(φ)sin(ϑ), sin(φ)sin(ϑ), cos(ϑ)) is given by: where the 〈K 2i 〉 coefficients are the ensemble average of the cubic harmonics. A cubic harmonic K 2i is a specific combination of spherical harmonics of order 2i. 73 The expansion is usually limited to a few terms and may be used to refine diffraction data.  Fig. 3), pure [001] (Fig. 3a) or [111] (Fig. 3b) ordering are qualitatively well described but not the pure [110] one (Fig. 3c), given that the dotted line in Fig. 3c Fig. 3). In real systems, fewer terms are usually needed for cubic harmonics expansions. Fig. 4 reports on the experimental results analyzed in the past for NaCN and KCN with an expansion up to K 6 ( Table 3). 74 For both these materials the [1,1,0] orientation corresponds to an absolute minimum of the probability distribution. NaCN has a preferred orientation along the [001] direction, but significantly different from pure [001] ordering. Indeed, the dashed line in Fig. 4 has a smaller deviation from the horizontal line (corresponding to the isotropic case) than the dashed line in Fig. 3b Fig. 4). The molecular axis exhibits a preferred orientation along [001] but the deviation from the isotropic distribution is less important than in cyanide compounds. This may explain why a reasonable fit of neutron scattering experimental data was obtained in the cubic phase with a purely isotropic rotational model. 61 (Fig. 2b), whereas a repulsive interaction between MA + and Pb 2+ is expected along [111] in MAPbX 3 crystals (Fig. 2a). By analogy to symmetry-adapted combinations that are able to describe PS excited states, it is also possible to design symmetry-adapted rotator functions for rotational excited states. These expansions from spherical harmonics are available for all the IR of the O h group. 73 They are connected to the symmetry-adapted polynomial expansions provided that the unit vector is expressed as u = (x, y, z). The Γ 4 − vectorial (l = 1) IR appears in the expansion of the PS for A, B and C scenarios. The lowest-order symmetry-adapted polynomials with Γ 4 − symmetry appearing as a possible term for a Γ 4 − rotator function are thus (x, y, z). They correspond to the orientation of the CN axis of MA + . The average time dependence of the autocorrelation function of this vector extracted from MD simulation at 400 K is plotted in Fig. 5 for various cubic supercell sizes. Interestingly, the relaxation time depends much on the parity of the number of cells along one direction of the supercell. The CN axis orientation is more stabilized for an even number of cubic cells. The reasons are at least twofold. First, given that the MD simulations are limited to a Γ point sampling of the BZ, hybridization of atomic orbitals related to BZ-folding cannot be adequately addressed with an odd number of cells in one direction. This major effect is exemplified by considering a linear chain in Fig. 5b. 43 Next, the stabilization of the CN axis obtained with an even number most  probably results from enhanced anti-ferroelectric interactions between the molecular dipoles, whereas frustrated interactions are expected with an odd number of cells (Fig. 5c). In the latter case parallel dipole orientations between neighboring supercells may result in unphysical ferroelectric dipole-dipole interactions that are not present in the real bulk crystal. In fact, even though some hybrid perovskites may exhibit phase transitions with a ferroelectric character (vide infra), MAPbI 3 essentially undergoes antiferro-distortive collective phase transitions at low temperatures, where long-range antiparallel arrangements of the molecular dipoles have been evidenced by diffraction experiments. This highlights the need for appropriate design of supercells in MD simulations, which should not break by construction of the overall symmetry, nor such longrange correlations. Fig. 6 reports on the temperature dependence of the average autocorrelation function for Γ 4 − deduced from simulations at three different temperatures. In each case, the relaxation time corresponds to the relaxation rate γ 1 derived within the PS picture for the tumbling of dipoles. As expected, lowering the temperature induces a slowing down of the molecular tumbling. But, the time dependence is strongly non-exponential and can be nicely fitted using stretched exponentials ( Fig. 6): C(t ) = e −(t/τ) β . Due to issues related to affordable time scales, it can hardly be extracted at 300 K. For the other temperatures, 400 and 450 K, the stretching parameter β can be reasonably estimated to 0.56. The relaxation time deduced for the electric dipole (l = 1) dynamics (Γ 4 − ), from the simulation at T = 400 K (Fig. 6a), amounts to 10.5 ps. At high temperature, the available experimental data are not consistent: the dynamics is either described by a single relaxation time with a large value of 12.7 ps at T = 370 K 60 or, a superposition of relaxation dynamics, whose vectorial (l = 1) component has a relaxation time equal to 1.62 ps at the same temperature. 61 The Γ 3 + quadrupolar (l = 2) IR also appears in the expansion of the PS scenarios A and B ( Table 1). The lowest-order symmetry-adapted polynomial with Γ 3 + symmetry appearing as a possible term for an elastic quadrupole is Y 0 2 ∝ (3z 2 − 1). The relaxation time for the autocorrelation function of this rotator extracted from MD simulation at 400 K is equal to 6.6 ps. It corresponds to the γ 2 coefficient of the simple PS analysis reported in section 3.3.2. The actual ratio γ 2 /γ 1 is equal to 1.6, which is closer to PS scenario A than to scenarios B and C (section 3.3.2). This is also consistent with the maximum observed along the [001] direction for the probability distribution function (Fig. 4). The computed relaxation time is smaller than the experimental one (0.54 ps) deduced from a fully isotropic model for the angular momentum l = 2. 61 The analyses of neutron scattering data 60,61 do not consider lateral coupling between the rotational motions of molecules. Individual molecular random reorientations are also described by simple PS relaxation models in section 3.3.2. A similar picture emerges from the MD study presented in section 3.3.3, provided that spurious spatial correlations are removed with supercells containing only an even number of cubic cells (Fig. 5). A random isotropic molecular tumbling is described by a single relaxation time in ref. 60, whereas a more elaborate isotropic model leads to the definition of a relaxation time for each angular momentum in ref. 61. The PS and rotator models   61. This is necessary to account for both electric and elastic multipolar effects, as well as to be able to discuss the possibility of ferroelectric phase transitions. Further analyses of neutron scattering data using symmetry adapted O h rotators should provide a better description of the rotational molecular disorder in HOP.

Linear coupling between collective pseudospins/rotators and lattice vibrations
3.4.1. Simple ferroelectric distortions from the reference pseudocubic phase. The linear coupling between disordered orientational degrees of freedom and lattice vibrations can be predicted using symmetry-adapted PS or rotator functions and group theory. It is useful to separate lattice displacements at the Γ point of the BZ in acoustic and optical phonons (Table 1). 71 The odd Γ 4 − electric dipole IR of the PS is then linearly coupled to the optical phonons. Table 4 gathers some of the basic cubic Pm3m cell distortions relevant to HOP, leading to possible collective ferroelectric phase transitions illustrated in Fig. 7. 10 For MAPbI 3 , the P4mm polar group is a possible alternative to the cubic Pm3m space group for structural refinement of the high temperature diffraction data. 44 Our MD simulations however suggest that local polar symmetry breaking may only occur in relation to frozen configurations of the cations, and not collective fluctuations. The Amm2 appears in the crystallographic analysis of the room temperature phase of FASnI 3 . 44 The R3m space group is relevant for germanium compounds, MAGeX 3 . 43,76 As shown by a systematic study on germanium iodide compounds, CsGeI 3 , MAGeI 3 and FAGeI 3 , 76 the appearance of a ferroelectric trigonal R3m space group and thus a large SHG response was attributed to the high density of states in the valence band due to sp-hybridization of the Ge and I orbitals. It shows that ferroelectricity for this class of HOP and AIP is essentially related to the instability of the Pm3m cubic structure against a phonon mode related to the inorganic lattice, rather than collective reorientational motions of the cations. However, for larger organic cations like GA, DMA or TMA, the influence of the lattice instability is reduced and the size and shape of the cations may induce a variety of crystallographic structures, some of them being non-polar like P2 1 /c for GAGeI 3 . 76 3.4.2. Acoustic phonon anomalies. The acoustic phonon frequencies vanish at the Γ point of the BZ, but the coupling between acoustic phonons and other excitations can be predicted close to the BZ center using the IR of the strain tensor ( Table 1). The Γ 3 + (E g ) and Γ 5 + (T 2g ) excited elastic multipoles are thus expected to be linearly coupled to acoustic phonons through the Γ 3 + (E g ) tetragonal and Γ 5 + (T 2g ) shear strain components. A similar behavior is predicted for cyanide compounds. 75 The elastic free energy is given by: with the symmetry-adapted elastic constants given by C A1g = C 11 + 2C 12 , C Eg ¼ 1 2 ðC 11 À C 12 Þ, C T2g = C 44 and strain tensor components ε A1g = (ε 1 + ε 2 + ε 3 )/√3 ( proportional to the volumetric strain), ε Eg1 = (2ε 3 − ε 1 − ε 2 )/2√3 and ε Eg2 = (ε 1 − ε 2 )/2 (tetragonal strains) and ε T2g1 = ε 4 , ε T2g2 = ε 5 , ε T2g3 = ε 6 (shear  strains). 20 In order to predict the effect of linear coupling between strain and the random orientational degrees of freedom, the total free energy is written as a function of the phonon displacement variables u and the orientational variables where M is the bare k dependent dynamical matrix related to V TT , the single particle susceptibility matrix χ 0 and the bare k dependent coupling matrix J are related to V RR and v the translation-rotation linear coupling matrix to V RT .j is the onsite lattice self-energy connected to the indirect interaction between the molecule orientations and is related to both V RT and V TT . Separating the two degrees of freedom leads to renormalized (or dressed) coupling 20 These renormalizations may lead to ferroelastic structural instabilities 77 that, to the best of our knowledge, are not observed in hybrid perovskites. In more standard cases, the linear coupling leads to renormalized elastic constants and additional temperature depen- where C 0 IR is the elastic constant for a given IR and B is the corresponding constant of the coupling matrix. The additional term may induce an unusual increase of the elastic constant with temperature, 78 which could be tested experimentally in hybrid perovskites to gauge the importance of translation/rotation coupling.
Another possible outcome of linear translation-rotation couplings, even more characteristic of the stochastic molecular dynamics, is an additional acoustic damping related to the tumbling relaxational dynamics. 20 When combining a Lagrangian, which includes the total free energy plus an additional kinetic energy term depending on the density ρ, with a dissipa-tive part accounting for frictional forces, it is possible to predict the relaxation dynamics of the tumbling as well as its effect on acoustic wave propagation. 20,79 Keeping only a damping matrix Γ for the relaxational motion, the set of linearly coupled dynamical equations reads: For a single orientational degree of freedom, the relaxation rate without translation-rotation coupling is: γ = 1/τ = [χ 0 −1 + J + j]/Γ. Separating the two degrees of freedom leads to a where the Laplace transform variable is z = ω + iε + in the limit of vanishing bare phonon damping. This result is similar to the one discussed above (section 3.2) for the coupling between a relaxor and optical modes in KMnF 3 or SrTiO 3 , but the impact of damping renormalization is expected to be very important and close to the Γ point. The coupling will thus lead to two separate signatures, one from relaxor and another from the acoustic phonon, as long as the acoustic phonon energy is much larger than the relaxation energy. However, due to acoustic phonon dispersion, a strong damping affects the whole acoustic phonon branch when k is close to the BZ center. Such a strong damping may affect a large part of the acoustic phonon spectrum in hybrid perovskites, especially the transverse acoustic phonon branches that are limited to about 2 meV over the entire BZ (section 3.2), given that the typical relaxation time of 5 ps (Fig. 6) corresponds to 0.12 meV. Thus, a careful investigation of the temperature dependence of the acoustic phonon spectrum would yield fruitful information about the influence of the stochastic dynamics on the translation-rotation coupling.  3.4.3. Antiferro-distortive phase transitions. Structural phase transitions in HOP can be inspected on the basis of Landau theory, already well developed for perovskite materials such CaTiO 3 , where distortions have been associated with IR of the parent Pm3m space group. Here, we adapt this model to the case of AIP and HOP, with contributions from the disordered organic part in the latter.
Some of the related antiferro-distortive space group changes associated with order parameters for AIP (P4/mbm) and HOP (I4/mcm and Pnma) are summarized in Table 5 and illustrated in Fig. 8 (Table 1). Symmetry analysis confirms that such a vibrational mode, mainly related to the inorganic lattice, exists in the high temperature phase of HOP. As condensation of a vibrational mode is usually related to a displacive mechanism, a phonon mode softening should be observed in HOP: (i) at the BZ boundary (R point) by neutron scattering in the high temperature phase, (ii) at the BZ zone center both by neutron and Raman scattering in the low temperature phase. 62 Symmetry analysis further shows that the phase transition can also be associated to a mode involving the collective tumbling of molecular cations. This is true for the PS scenario B, namely for the C-N axis pointing toward iodine atoms, or the mixed A/B or D scenario, where the C-N axis is pointing toward an intermediate direction between a face and an edge. It shows that an additional order parameter having the same lattice symmetry as the octahedra rotations (R 4 + ) is associated to molecular tumbling. From the Landau theory point of view, these order parameters undergo a linear coupling. But, the mechanism associated to these additional order parameters has an order-disorder nature. In other words, even though a linear coupling is predicted between these order parameters in the Pm3m phase, the characteristic dynamics appear in different frequency ranges. Detailed derivations of the possible Landau free energy expansions are known. 80 A linear-quadratic coupling is also predicted between the two linearly coupled R 4 + order parameters and two components of the strain tensor, ε 1 + ε 2 + ε 3 and 2ε 3 − ε 1 − ε 2 that correspond to the Γ 1 + and Γ 3 + IR, respectively. The cubic to tetragonal phase transition presents an improper ferroelastic character. Such a linear-quadratic coupling is in agreement with the available diffraction data recorded on MAPbI 3 . 70 Below the critical temperature T c , antiferro-distortive octahedra rotation was found to vary as (T c − T ) 0.25 and the tetragonal strain as (T c − T ) 0.42 . 70 The anomalies of elastic constants are predicted as a function of the order-parameter. These anomalies could be tested experimentally in order to quantify the linear-quadratic coupling with strains.    9 Schematic diagram indicating the group-subgroup relationship between some of the space groups observed in HOP and AIP (solid red lines). The tilt systems are given using Glazer's notation, as well as the Pm3m IR related to the lattice instabilities. The dashed red line indicates that a transition from I4/mcm to Pnma space groups must be a 1 st order phase transition.
In oxide perovskites, Pm3m to Pnma and Pm3m to I4/mcm instabilities compete in the same way as in MAPbI 3 and MAPbBr 3 . The Pm3m to Pnma instability is related to the simultaneous condensation of R 4 + and M 3 + rotational modes of the octahedra. Our symmetry analysis confirms that an M 3 + rotational mode, mainly related to the inorganic lattice, exists in the high temperature phase of HOP (Table 1). The Pm3m to Pnma structural instability is indeed directly observed in the [536-602 K] temperature interval for CsPbI 3 , leading to a Pnma phase at room temperature without the intermediate I4/mcm phase. 41 In MAPbX 3 , additional PS order parameters with the M 3 + symmetry exist for disordered collective molecular tumbling in scenarios A and B, as well as in the mixed scenario D (Table 1). They are expected to present an order-disorder character and to be linearly coupled with the displacive order parameter.
Finally, it appears that the Pm3m to P4/mbm instability is only observed in some AIP, namely CsPbCl 3 . 51 This instability is triggered by a heavily over-damped M 3 + rotational displacive mode of the lead halide octahedra in AIP. This result suggests that the potential anharmonicity associated to this normal vibrational mode is large in the cubic phase, and that the phase transition has an intermediate character between displacive and order-disorder. The Pm3m to P4/mbm phase transition is followed by P4/mbm to Cmcm and Cmcm to P2 1 /m phase transitions (Fig. 9). 51

Translational symmetry-breaking and orientational glass behavior
Neutron scattering data have revealed that only a minor fraction of the MA cations is involved in both the rotational motions at high temperature. 60 This is consistent with MD simulations that will be described in this section. Moreover, some authors claim that structural refinements are not easy in pure HOP, even at low temperatures, due to the remaining lattice disorder. 44 It often leads to neglecting a large number of weak Bragg peaks in the structural refinement of diffraction data. This might be indicative of the remaining disorder in the orientations of the cations, thus leading to an incomplete long range freezing of cation orientations. The time dependence of the average autocorrelation function for molecular reorientation in the MD is strongly non-exponential and can be nicely fitted using stretched exponentials (Fig. 6). A detailed inspection of the dynamics of the CH 3 NH 3 + cations at T = 400 K ( Fig. 10b) reveals that it is related to a wide dispersion of reorientation dynamics of individual molecules. A significant part of the cations (roughly 30%, represented in white in Fig. 10a) remain almost frozen with respect to a tumbling motion on the time scale of our MD trajectory. These results are qualitatively consistent with the recent neutron scattering data, which indicate that 60%-77% of the cations have residence times exceeding 200 ps. 60 In pure cyanide compounds, the CN − tumbling motion exhibits long-range correlations, because it is associated to a single orientational degree of freedom. Frustrated interactions and behavior typical of orientational glass are thus observed in mixed cyanide crystals by either anion (KBr x CN 1−x ) or cation (Na x K 1−x CN) alloying. 69,81 We may infer that frustrated cation orientations are more likely in MAPbX 3 compounds for two reasons: the probability distribution functions of molecular orientation exhibits a lower contrast (Fig. 4) and additional disordered rotational degrees of freedom are associated to the rotation of methyl and ammonium groups around the C-N axis. As shown by NMR and neutron scattering, this reorientational motion occurs on a time-scale comparable to the C-N tumbling. 24,60,61 The absence of long range spatial correlations between the individual reorientations could be better assessed from MD simulations using even larger supercells, although it appears to be consistent with the analyses of neutron scattering data (section 3.3.3). 60,61 However, from Fig. 5 we know that one should focus on even numbers of cells, namely 4 × 4 × 4 in the present work. Unfortunately, a 6 × 6 × 6 supercell is already too large, given the available computational resources. Therefore, it is difficult to explore more in detail the long-range critical correlations appearing at the continuous cubic to tetragonal phase transition within a first-principles approach.

Conclusion
The main aim of this work was to show that room temperature properties of hybrid perovskites can be viewed as a merging of properties of all-inorganic perovskites, conventional semiconductors and plastic crystals. Both the stochastic molecular orientational dynamics and the frozen static configurations can be explored rigorously on the basis of simplified and symmetry-based discrete PS as well as using continuous rotator functions. For instance, it has been shown how molecular orientational degrees of freedom can couple with the order parameters, both for ferroelectric and antiferro-distortive phases. Theoretical predictions and experimental guidelines are provided to explore the linear couplings between the various stochastic degrees of freedom and acoustic or optical phonon modes. Thorough Brillouin, Raman and neutron scattering experiments should help gauging the importance of linear translation/rotation coupling in these materials. Finally, as indicated by molecular dynamics simulation, an orientational glass behavior may result from the frustrated molecular interactions and degrees of freedom.