Activation Energy of Organic Cation Rotation in CH3NH3PbI3 and CD3NH3PbI3: Quasi-Elastic Neutron Scattering Measurements and First-Principles Analysis Including Nuclear Quantum Effects.

The motion of CH3NH3+ cations in the low-temperature phase of the promising photovoltaic material methylammonium lead triiodide (CH3NH3PbI3) is investigated experimentally as well as theoretically, with a particular focus on the activation energy. Inelastic and quasi-elastic neutron scattering measurements reveal an activation energy of ∼48 meV. Through a combination of experiments and first-principles calculations, we attribute this activation energy to the relative rotation of CH3 against an NH3 group that stays bound to the inorganic cage. The inclusion of nuclear quantum effects through path integral molecular dynamics gives an activation energy of ∼42 meV, in good agreement with the neutron scattering experiments. For deuterated samples (CD3NH3PbI3), both theory and experiment observe a higher activation energy for the rotation of CD3 against NH3, which results from the smaller nuclear quantum effects in CD3. The rotation of the NH3 group, which is bound to the inorganic cage via strong hydrogen bonding, is unlikely to occur at low temperatures due to its high energy barrier of ∼120 meV.

Hybrid n class="Chemical">perovskite photovoltaic (HPPV) technology[1,2] has received rapidly increasing interest from the emerging solar-cell community due to its record increase in power-conversion efficiency (PCE) during the last 5 years. The current state-of-the-art HPPV architecture was proposed in 2012, achieving ∼10% PCE at the time.[3,4] With the advancement of materials design and synthesis, the PCE of HPPV cells has already reached 22%, overtaking the best-performing inorganic-based single-junction thin-film cells. The most common photoactive material in HPPV cells is methylammonium (MA ≡ CH3NH3) lead triiodide (CH3NH3PbI3, shortened as MAPbI3 hereafter). This material combines many advantageous features for photovoltaic applications, e.g., optimal band gap for single-junction solar-cell absorbers,[5] excellent absorption strength in the visible part of the solar spectrum,[6] and high mobilities for both electron and hole transport.[7,8] Moreover, hybrid perovskites can be synthesized in solution at low temperature from common materials. Therefore, HPPV cells have become promising candidates for clean, affordable, and sustainable energy generation.

Despite its promises, many funn class="Chemical">damental materials characteristics of MAPbI3 remain poorly understood. Recent experimental and theoretical studies focused on the device architecture and performance,[9−11] different hybrid-perovskite compositions,[5,12−14] the large diffusion-lengths and low recombination rate,[7,8,15,16] the current–voltage hysteresis,[5,17−19] and the stability of hybrid perovskites materials.[20−22] The microstructure of MAPbI3, which is crucial for the initial stages of the light-to-energy conversion process in particular the generation and separation of an electron–hole pair, has received little attention. Although first-principles studies have already revealed that the organic MA+ cations are bound to the inorganic cage via hydrogen bonding between the NH3 group and halide anions,[23,24] information on the motion of the MA+ cations remains controversial. There are two conflicting viewpoints: (1) The MA+ cations can move nearly freely at high temperatures,[25,26] because only a small activation energy (Ea) is required for the cleavage of the hydrogen bonds.[27,28] (2) The MA+ cations are fixed in stable orientations, because Ea is of the order of 100 meV.[29−33]

In this letter, we adn class="Chemical">dress the question of MA+ rotation by combining quasi-elastic neutron scattering (QENS) experiments with first-principles density-functional theory calculations. With QENS measurements we study MA-motion in the low-temperature (<160 K) orthorhombic phase of MAPbI3. QENS is very sensitive to incoherent scattering events involving hydrogen nuclei. It is thus the ideal tool to study the motion of CH3NH3+.

In previous MAPbI3 QEn class="Chemical">NS studies, Leguy et al. measured Ea to be 9.9 meV[27] and Chen et al. 48 meV.[32] Both did not consider the difference between the axial rotation (that is, the rotation around the C–N bond) of only the methyl (Me ≡ CH3) group and that of the whole MA+ cation. In this work, we studied the local dynamics of MA+ cations in MAPbI3 powder samples at temperatures <160 K using incoherent QENS. To better understand the rotation mechanism, we also performed QENS on samples with deuterated methyl groups (CD3NH3PbI3). If deuteration changed Ea, it would indicate that the rotation would mainly involve the methyl group but not the ammonium group. In such a case, CH3 rotates relative to the NH3 group, which in our picture would remain hydrogen-bonded to the inorganic cage.

To obtain an in-depth unn class="Chemical">derstanding of the fundamental atomic processes in MAPbI3, we performed first-principles density-functional theory (DFT) calculations for different MA+ rotations. For each rotation process, we first mapped the potential energy surface onto a 1D hindered-rotor model to extract the barriers. We then performed ab initio path integral molecular dynamics (PIMD) calculations that include nuclear quantum effects beyond the 1D model to better characterize the free energy surface from which an estimate of Ea can be extracted. From the calculated activation energies, we can directly estimate whether a dynamical process is likely to occur at a certain temperature.

Neutron Scattering Experiments

We carried out n class="Chemical">different types of neutron scattering experiments on powder samples of CH3NH3PbI3 and its partly deuterated analogue CD3NH3PbI3 using the backscattering spectrometer IN16B at the Laue-Langevin Institute in Grenoble. The incoherent neutron scattering cross section of 1H is 80.2 b. It provides the major contribution to the total (coherent and incoherent) neutron scattering cross section of both samples. The coherent contribution from 1H is 1.8 b. For comparison, the coherent and incoherent neutron scattering cross section of D ≡ 2H are 5.6 and 2.0 b, respectively. As a result, the total neutron scattering cross section of CH3NH3PbI3 is 523.9 b, much larger than CD3NH3PbI3 (300.7 b). These two samples have quite similar total absorption cross sections (10.8 and 9.8 b for CH3NH3PbI3 and CD3NH3PbI3, respectively) and nearly the same density. Thus, for similar absorption power, the total neutron scattering power of CH3NH3PbI3 is larger than CD3NH3PbI3 by a factor of 1.69.

The chemical synthesis and nucln class="Chemical">ear magnetic resonance (NMR) characterization of samples are described in Section S1 of the Supporting Information. The principles and details of our neutron scattering measurements are given in Section S2 of the Supporting Information. Three types of neutron scattering measurements (elastic, inelastic, and quasi-elastic) were carried out and will be discussed in this letter.

First, we performed elastic fixen class="Chemical">d window scans (EFWS) for a quick overview of the onsets of motions that are faster than the time scale corresponding to the energy resolution. The temperature-resolved results are given in Figure . Panel a shows a ratio of 1.48 for the elastic scattering intensity between the hydrogenated and the deuterated samples at a very low temperature (2 K), where all rotational motions are frozen. This ratio does not deviate much from the aforementioned theoretical ratio of 1.69. 1H NMR spectroscopy obtained after the drying process of the deuterated sample clearly confirmed that the deuteration rate is 99% (see Figure S1 of the Supporting Information).

Figure 1

Elastic fixed window scans (with summing over all selected Q values, see Section S2 of the Supporting Information) of CH3NH3PbI3 (blue) and CD3NH3PbI3 (red): (a) at low temperatures (3, 130) K and (b) within the whole scanned temperature range (3, 350) K with intensities normalized according to that obtained at 3 K.

Elastic fixed winn class="Chemical">dow scans (with summing over all selected Q values, see Section S2 of the Supporting Information) of CH3NH3PbI3 (blue) and CD3NH3PbI3 (red): (a) at low temperatures (3, 130) K and (b) within the whole scanned temperature range (3, 350) K with intensities normalized according to that obtained at 3 K.

EFWS is a recording of the integraten class="Chemical">d intensities of elastic scattering events. Their attenuation with increasing temperature is described by the Debye–Waller factorwhere S is the elastic scattering function, C a constant, the scattering vector, and the displacement of the scatterers. Figure a shows that the mean square displacement ⟨u2⟩ decreases (approximately) linearly for temperatures lower than 58 and 55 K for the all-protonated and deuterated perovskites, respectively. At higher temperatures, the quasi-elastic component of scattering becomes visible within the resolution of the spectrometer. It broadens the scattering function and accordingly the elastic intensity decreases more rapidly, as shown in Figure a. This broadening increases with temperature. Above a certain temperature (80 K for CH3NH3PbI3 and 90 K for CD3NH3PbI3), it becomes too large to be accurately measured (see Figure S3 of the Supporting Information). At this point, the investigated motion becomes so rapid that its quasi-elastic contribution does only weakly contribute to the elastic intensity by its amplitude. Figure a shows that the elastic intensity decreases in the same proportion for the two compounds within the respective temperature ranges where the elastic intensity is most sensitive to the change of the quasi-elastic broadening. This may suggest that the local dynamics involved in these two compounds have the same characteristics.

A scan over a larger temperature range is given in Figure b, in which the scattering intensity of CD3NH3PbI3 is normalizen class="Chemical">d to that of CH3NH3PbI3 at low temperatures. The data illustrates that at temperatures lower than the orthorhombic-tetragonal phase transition at ∼160 K, these two samples exhibit a similar behavior, with the remaining scattering intensity of CD3NH3PbI3 noticeably higher than its all-protonated analogue. It also exhibits that the elastic-intensity decrease becomes slower between 80 and 150 K (before the orthorhombic-tetragonal phase transition at ∼160 K). We have examined this behavior carefully by studying the character of the elastic incoherent structure factor (EISF) and present our analysis in detail in Section S2 (Figures S4–S6) of the Supporting Information. Our analysis supports the assumption that the local dynamics mainly involve the rotation of only the methyl group (while the ammonium rotor remains bound to the inorganic framework via hydrogen bonds). The objective of this work is to couple neutron scattering analysis with first-principles calculations of the atomic structure and potential energy surfaces (in the next section) to verify this assumption.

Next we discuss the inelastic fixen class="Chemical">d window scans (IFWS). Figure shows the results for the normalized intensities of CD3NH3PbI3 scaled by a factor of 1.48. Both data exhibit a quasi-plateau at temperatures lower than 50 K, an increase from about 50 K, and a maximum at 85–90 K, where the peak of the all-protonated system is higher. For a quantitative analysis, we fit the IFWS data with eq following ref (34) and show the results in Figure with solid lines. Specifically, eq describes the relationship between the scattering intensity IIFWS and relaxation time τ at a given temperature T and energy window ωoff. We assume that the dynamics follows a stochastic process so that it can be described by eq (Arrhenius equation), where kB is the Boltzmann constant. The temperature at which IIFWS takes its maximum can be calculated by eq . The fitted activation energy results are listed in Table . We will discuss them together with the quasi-elastic measurement results.

Figure 2

Inelastic fixed window scans of CH3NH3PbI3 (blue) and CD3NH3PbI3 (red). The fits using eq are given in solid lines.

Table 1

Activation Energy (Ea) and Pre-Exponential Factor (τ0) of MA-Motion in Low-Temperature CH3NH3PbI3 and CD3NH3PbI3 Obtained by Fitting Inelastic (IE) and Quasi-Elastic (QE) Neutron Scattering Dataa

		CH3NH3PbI3	CD3NH3PbI3
Tmax [K]	IE	87.9 ± 0.4	90.3 ± 0.3
Ea [meV]	IE	46.9 ± 1.3	57.5 ± 0.7
	QE	48.2 ± 3.9	53.8 ± 7.5
τ0 [fs]	IE	688 ± 122	205 ± 14
	QE	403 ± 270	1540 ± 410

Also listed are the maximal IE-intensity temperature (Tmax).

Inelastic fixed winn class="Chemical">dow scans of CH3NH3PbI3 (blue) and CD3NH3PbI3 (red). The fits using eq are given in solid lines.

To check the consistency of our IFWS Ea results, we also carrien class="Chemical">d out quasi-elastic measurements at different scattering vectors and temperatures. Figure shows that, for both samples in the low-temperature range, the logarithm of the relaxation time depends linearly on the inverse temperature. By fitting these data with the Arrhenius equation, we calculated the activation energies of MA-motion in both cases. The results are indicated in Figure and listed in Table . Note: the signal-to-noise ratio in the spectra recorded for CD3NH3PbI3 is much lower than those for CH3NH3PbI3, as presented in Section S2 of the Supporting Information.

Figure 3

Arrhenius plot (relaxation time vs temperature) of (a) CH3NH3PbI3 and (b) CD3NH3PbI3 from quasi-elastic neutron scattering at low temperatures. Activation energies (Ea) are listed in the plot.

Arrhenius plot (relaxation time vs temperature) of (a) CH3NH3PbI3 ann class="Chemical">d (b) CD3NH3PbI3 from quasi-elastic neutron scattering at low temperatures. Activation energies (Ea) are listed in the plot.

Table illustrates that IFWS and quasi-elastic mn class="Chemical">easurements produce Ea values in good agreement with each other for both all-protonated and partly deuterated MAPbI3. We thus claim that IFWS is an appropriate tool for the study of local dynamics in MAPbI3. Specifically, the fitted Ea for low-temperature CH3NH3PbI3 is 47–48 meV, agreeing well with previous QENS results (48 meV).[32] The activation energy of local dynamics in the deuterated sample is larger than in the all-protonated one by 6 meV (quasi-elastic) or 11 meV (IFWS). The good agreement of this fit with the data indicates that the process is thermally activated. Thus, we can expect that the main contribution to the observed isotope effects comes from zero-point-energy.

Based on our elastic ann class="Chemical">d inelastic measurements, we can already conjecture the atomistic origin of the local dynamics in the low-temperature phase of MAPbI3. As alluded to before, Chen et al. claimed that both CH3 and NH3 rotors participate in a rotation of the MA+ cation.[32] However, considering the discussion of the EISF evolution (see Section S2 of the Supporting Information for details), we instead propose that only CH3 rotors move for the following reasons. First, as the H atoms in the methyl group are not bound to the inorganic cage with hydrogen bonding, we should obtain similar Ea as found in many other systems where CH3 rotors are present. Actually, our Ea of 47–48 meV is close to the QENS results of Me-side-group rotation in polymers (5.3 kJ mol–1 = 55 meV).[35] Second, the activation energy of CD3NH3PbI3 and CH3NH3PbI3 noticeably differ. In the following, we will support our conjecture with first-principles calculations that provide a precise atomistic picture.

First-Principles Density-Functional-Theory Calculations

To model the MA-motion in the orthorhombic phase of MAn class="Chemical">PbI3, we carried out first-principles (DFT) calculations on a series of model systems. The lattice parameters were adopted from single-crystal X-ray-diffraction data a = 8.56 Å, b = 8.84 Å, and c = 12.58 Å.[12] They characterize a supercell consisting of four single (quasi-primitive) cells. In each single cell, Pb2+ ions are located at the corners and I– at the edge-centers. MA+ cations, located close to the cell centers, take different (quasi-)face-to-face directions in different cells. Figure a depicts the DFT-relaxed atomic geometry (for a fixed size of unit cell marked by solid lines), which agrees well with previously reported DFT[33] and ab initio molecular dynamics[36] results. To reduce the coupling between nearest rotating MA+ ions (that result from the periodic boundary conditions), we duplicated the orthorhombic cell into a 2 × 2 × 2 supercell, which contains eight single MAPbI3 unit cells. Its lattice parameters are a = b = 12.30 Å, c = 12.58 Å, and α = β = 90°, γ = 88.15°. Only one of eight MA+ ions (Figure b) is subject to the investigated rotations, which keeps it well separated from other rotating MA+ in the bulk structure.

Figure 4

(a) Optimized structure of the orthorhombic supercell (unit-cell edges marked by black lines) of MAPbI3. Here we duplicated it to a 2 × 2 × 2 right rhombic prism supercell. (b) Starting geometry for the investigated processes. Left: the rotated MA+ cation is highlighted, while the rest part of the supercell is shown in a lighter shade; right: Newman projection of the rotating MA+. (c) Investigated rotation modes of MA+ in the orthorhombic MAPbI3 structure: Me-rotation (upper panel) and MA-rotation (lower panel). C, N, H, Pb, and I atoms are colored in green, yellow, dark gray, blue, and red, respectively. The hydrogen bonds between H atoms in the NH3 group and I– anions are emphasized by gray dashed lines.

(a) Optimized structure of the orthorhombic supercell (unit-cell en class="Chemical">dges marked by black lines) of MAPbI3. Here we duplicated it to a 2 × 2 × 2 right rhombic prism supercell. (b) Starting geometry for the investigated processes. Left: the rotated MA+ cation is highlighted, while the rest part of the supercell is shown in a lighter shade; right: Newman projection of the rotating MA+. (c) Investigated rotation modes of MA+ in the orthorhombic MAPbI3 structure: Me-rotation (upper panel) and MA-rotation (lower panel). C, N, H, Pb, and I atoms are colored in green, yellow, dark gray, blue, and red, respectively. The hydrogen bonds between H atoms in the NH3 group and I– anions are emphasized by gray dashed lines.

Figure c depicts the two rotation mon class="Chemical">des considered in this paper: For the Me-rotation in the upper panel, only the CH3 group was rotated around the C–N bond, whereas for the MA-rotation in the lower panel the whole molecule was rotated around the C–N bond. The rotational angle is well-defined, even when we fully relax the geometry, as we treated the rotating part (CH3 or NH3) as a rigid body. According to the (quasi) C3 symmetry of our model systems, we only scan over the rotational angle range (0,120)°. We also fixed the length and the orientation of the C–N bond. Note: since CH3 is only bonded to NH3 but does not interact with the PbI3– framework via hydrogen bonding, we supposed that it will follow the NH3 rotation. NH3 rotation is therefore equivalent to MA+ rotation, as demonstrated by our results.

We din class="Chemical">d not explicitly consider the flip of the whole MA molecule, i.e., a rotation of the NH3 group around the C–N bond center. Chen et al. had ruled this process out in the low-temperature phase of MAPbI3.[32] Their assessment agrees with our own EFWS results and with test calculations reported in Section S3 of the Supporting Information and in our forthcoming publication.

In our DFT calculations, we usen class="Chemical">d the Perdew–Burke–Ernzerhof (PBE)[37] exchange-correlation functional augmented with long-range van der Waals (vdW) corrections based on the Tkatchenko–Scheffler method.[38] In addition, scalar relativistic effects were included by means of the zero-order regular approximation (ZORA).[39] Our previous study[24] demonstrated that this approach gives lattice constants of hybrid perovskites in good agreement with experiment and can appropriately describe the interaction between the organic cations and the inorganic framework. For a more accurate electronic structure, the total energy of each PBE+vdW relaxed geometry was calculated using the hybrid “PBE0+vdW” functional.[40]

All DFT calculations in this work were performen class="Chemical">d with the all-electron numeric-atom-centered orbital code FHI-aims.[41−44] For all calculations we used tier 2 basis sets. A Γ-centered 6 × 6 × 4 k-point mesh was used for the optimization of the orthorhombic supercell structure, and a 4 × 4 × 4 k-point mesh for the 2 × 2 × 2 supercell models. For the path integral molecular dynamics (PIMD) simulations, we employed the interface of FHI-aims with the i-PI[45] code. The PIMD calculations were run in a cluster model, in which a single PbI3– unit was held fixed in its bulk position and the CH3NH3+ and CD3NH3+ molecules were allowed to move freely. The simulation was run at 150 K with 12 beads and the PIGLET thermostat.[46] The results of all relevant calculations of this work are available from the Novel Materials Discovery (NOMAD) repository.[47]

We start the presentation of our DFT results with the Me-rotation. The n class="Chemical">PBE+vdW relaxed structures at several rotation angles (30, 60, and 90°) of this process are shown in Figure . During the rotation of the CH3 group, the rest of the supercell (shown in light colors) remains largely unaffected. A clearer view is given by the Newman projection in the lower panel. As the CH3 group rotates from the “staggered” (0°, i.e., the fully relaxed structure shown in Figure b) to the “eclipse” (60°, Figure b) configuration, the location of the NH3 group and the surrounding inorganic cage remain nearly unchanged. We can thus conclude that the force on NH3 induced by the torsion around the C–N single bond is much weaker than that related to the hydrogen-bonding interaction with the inorganic framework.

Figure 5

Optimized structures at several Me-rotational angles: (a) 30°, (b) 60°, and (c) 90°. Upper panel: atomic structure of the whole unit cell. Lower panel: Newman projection of the rotating MA+ and the surrounding inorganic framework.

Optimized structures at several Me-rotational angles: (a) 30°, (b) 60°, ann class="Chemical">d (c) 90°. Upper panel: atomic structure of the whole unit cell. Lower panel: Newman projection of the rotating MA+ and the surrounding inorganic framework.

In Figure we plot the potential energy surface (PES) for Me-rotation with both PBE+vn class="Chemical">dW and PBE0+vdW. For clarity, we set the total energy of the 0°-structure to 0 and shifted all other data accordingly. The PBE+vdW PES (green triangles in Figure ) exhibits three minima of equal-height. They are located at 0 (360), 120, and 240°, where the CH3 and NH3 take the staggered configuration in the Newman projection. There are three corresponding equal-height maxima (62 meV) at 60, 180, and 300° that each exhibit eclipsed H atoms. The PBE+vdW total energies of Me-rotation can be perfectly characterized by a cosine function with period 120° (see Section 4 of the Supporting Information). This reflects the nearly perfect C3 symmetry of the Me group, which is not noticeably perturbed by inhomogeneous interactions with the inorganic framework. PBE0+vdW gives an equivalent PES (black squares in Figure ) with slightly higher (by 5 meV) barriers, i.e., torsion energies for the C–N bond.

Figure 6

PES for Me-rotation. The PBE+vdW total-energy data is plotted with green triangles, PBE0+vdW with black squares, and the fitted cosine-potential to the PBE0+vdW results in red.

PES for Me-rotation. The PBE+vn class="Chemical">dW total-energy data is plotted with green triangles, PBE0+vdW with black squares, and the fitted cosine-potential to the PBE0+vdW results in red.

Now we turn to MA-rotation. Figure shows the n class="Chemical">PBE+vdW relaxed structures at several rotation angles of NH3 around the C–N bond. It indicates the major difference between CH3 and NH3 rotation: while Me-rotation does not affect the NH3 group, NH3 rotation forces CH3 to follow the rotation via the torsion potential. This effectively results in the rotation of the whole CH3NH3+ ion, which maintains the staggered configuration of H atoms during the whole process as reflected by the Newman projections (lower panel of Figure a–c). Such a rotation causes the cleavage of H···I hydrogen bonds, which, according to our previous study,[24] would have significant impact on the inorganic-framework deformation (the “chicken-and-egg” interplay). However, since the other seven MA+ ions are held fixed in our calculations, the PbI3– cage cannot deform significantly, as is also seen in Figure .

Figure 7

Optimized structures at several MA-rotational angles: (a) 30°, (b) 60°, and (c) 90°. Upper panel: atomic structure of the whole unit cell. Lower panel: Newman projection of the rotating MA+ and the surrounding inorganic framework.

Optimized structures at several n class="Disease">MA-rotational angles: (a) 30°, (b) 60°, and (c) 90°. Upper panel: atomic structure of the whole unit cell. Lower panel: Newman projection of the rotating MA+ and the surrounding inorganic framework.

Figure depicts the PBE+vdW and PBE0+vdW PESs. They exhibit similar features. The PBE+vdW PES (green triangles in Figure ) has a barrier of 161 meV from 0 to 60°. This corresponds to a breaking of all three hydrogen bonds (Figure b). Unlike for Me-rotation, the inclusion of exact exchange in the density functional has a large effect on this MA-rotation (black squares in Figure ), as it results in a noticeably lower barrier (by 28 meV).

Figure 8

PES for MA-rotation. The PBE+vdW total-energy data is plotted with green triangles, PBE0+vdW with black squares, and the fitted cosine-potential to the PBE0+vdW results in red.

To summarize, the interaction energies between the inorganic PbI3– cage, the n class="Gene">CH3 group, and the NH3+ part of an organic MA+ ion, has the following order:

The major contribution for V(PbI3––n class="Gene">CH3) includes vdW interactions and ionic bonding (the positive charge of MA+ is partially distributed over Me). Both are almost not affected by either Me- or MA-rotation. V(PbI3––NH3+) includes vdW interactions, ionic bonding, and hydrogen bonding. MA-rotation has significant impact on the last component. Note: Vtorsion(CH3–NH3+) is only related to the torsion around the C–N bond but not the covalent bonding itself.

Analysis of Rotational Potential Energy Surfaces Using a One-Dimensional Hindered-Rotor Model

From our DFT (PBE/PBE0+vdW) calculations, we have obtained PESs for both Me- and MA-rotation. They provide us with insight (e.g., energy barriers) into the associated processes. To calculate activation energies, we also need to consider zero-point-energy (ZPE) related to the quantum-mechanical vibrational ground-state energy of the associated PES.

Both Me and MA PESs exhibit a perion class="Chemical">dic and even function. In general, such a function can be expanded in terms of a cosine series in the rotational angle, here θ:Here the period of V(θ) is 360° (2π). V represents the potential energy for a generalized hindered rotor.

By projecting the n class="Chemical">PBE0+vdW PES of Me- and MA-rotation onto the cosine functions cos(nθ) with positive integer numbers n, we only obtained non-negligible coefficients for a few terms (apart from the constant term ): v3 = −32.9 meV,v6 = 2.0 meV for Me-rotation, and v2 = −2.2 meV,v3 = −63.7 meV,v6 = −6.2 meV for MA-rotation. In both cases, we set the constant term for V(0) = 0. The fitted PESs are presented by the red curves in Figures and 8.

The Hamiltonian for the system subject to the generalized-hinn class="Chemical">dered-rotor potential V (eq ) reads

We can calculate the moment of inertia I from the n class="Chemical">DFT-relaxed atomic structure (Figure a). For Me-rotation, we have I = 3.24 u Å2 (note, the NH3 part is not included here as its geometry nearly remains unchanged), which is close to the experimental value of 3.185 u Å2 obtained from rotational spectroscopy experiments of Me-rotation in the isolated CH3NH2 molecule.[48] For MA-rotation, in which all 6 H atoms contribute, we have 6.11 u Å2.

The time-indepenn class="Chemical">dent Schrödinger equationhas the form of a Mathieu equation and can be solved numerically.

By solving the eigenvalue problem eq at 0 K, we calculated the ZPEs of Me- ann class="Chemical">d MA-rotation, respectively. From this we obtain an estimate of the activation energy (note that our activation energy only includes enthalpic contributions and no entropy) by subtracting the ZPE from the energy barrier. Since we model this as a 1D reaction, only the 1D ZPE on the reactant state needs to be included. The details of solving the Mathieu equations eq are given in Section S5 of the Supporting Information. Here we list the results in Table .

Table 2

Results of One-Dimensional Hindered-Rotor-Model Approach to Me- and MA-Rotation (In the Left and Right Sides, Respectively)a

Me group	I	ZPE	Ea	MA group	I	ZPE	Ea
CH3	3.24	10.3	55.5	CH3NH3	6.11	11.4	119.3
CH2D	4.32	9.0	56.8	CD3NH3	8.97	9.4	121.3
CHD2	5.40	8.1	57.7	CH3ND3	9.35	9.3	121.5
CD3	6.48	7.4	58.4	CD3ND3	12.21	8.1	122.6

Listed are the moments of inertia I (in u Å2), the calculated zero-point-energy (ZPE), and activation energy Ea (both in meV). Also different deuterated cases are included for each rotation mode.

Listed are the moments of n class="Disease">inertia I (in u Å2), the calculated zero-point-energy (ZPE), and activation energy Ea (both in meV). Also different deuterated cases are included for each rotation mode.

The ZPEs of CH3- ann class="Chemical">d CH3NH3-rotation are 10.3 and 11.4 meV, respectively. Accordingly, the 1D hindered-rotor-model activation energies are 55.5 meV for CH3-rotation and 119.3 meV for CH3NH3-rotation. The energy required to break the H···I hydrogen-bonds is thus much larger than to rotate the CH3 or CD3 group around the C–N bond.

Our Me-rotation activation energy agrees well with the QENS mn class="Chemical">easurement of the methyl side group rotation in polymers (55 meV).[35] It is also very close to our QENS result (Arrhenius fit) of 48 meV for low-temperature MAPbI3. This activation energy had previously been ascribed to MA-rotation,[32] for which our calculations produce a significantly higher activation energy of 119.3 meV. Therewith the possibility of MA-rotation to occur at <160 K is very low. For example, our simulation indicates that the number of MA-rotation events would be ∼500 times less than the number of Me-rotation events at 150 K and would happen <10 times in 1 ns. Section S6 of the Supporting Information summarizes our analysis.

Table also lists our Ea estimates for the rotation of several n class="Chemical">deuterated Me and MA groups. Deuteration does not change the electronic structure of the system. Only the moment of inertia increases, which subsequently leads to a smaller ZPE (see Section S5 of the Supporting Information). Accordingly, the activation energy increases with the number of D atoms in the rotational group. This trend can be observed for both Me and MA rotation listed in Table . Specifically, for fully deuterated Me, i.e., CD3, our 1D hindered-rotor model approach results in an activation energy of ∼58 meV, agreeing well with our neutron scattering measurement results. Partially deuterated methyl groups have larger ZPEs than the fully deuterated CD3. This will contribute to the lower-energy signals in the neutron scattering spectra of CD3NH3PbI3.

Estimate of Activation Energies from Path Integral Molecular Dynamics

By means of path integral molecular n class="Chemical">dynamics (PIMD), we can estimate the activation barrier of the methyl group directly in three-dimensional space, without restricting ourselves to a 1D model. Due to the large computational cost of PIMD simulations, we here use a model, in which we consider an isolated system consisting of one PbI3– group cut out from the lattice with a methylammonium group attached to it. The PbI3– group is held fixed, while the MA+ can move freely. Due to the absence of the rest of the solid, we expect the barriers to be slightly lower than what they would be in the solid. We performed 30 ps long PIMD simulations with the PBE+vdW functional (tight settings for basis sets and numerical parameters in FHI-aims) for the hydrogenated and the deuterated systems after equilibrating both structures at 150 K, a temperature that is low enough for the orthorhombic structure to be the most stable and high enough for the methyl group rotation to be activated. We then calculated the probability density P projected on the rotation angle θ of the methyl group and calculated free energies as F(θ,T) = −kBT ln[P(θ)]. During our simulations, we never observed the cleavage of an H-bond between the NH3 group and the iodine atoms but instead a slight precession of the C–N bond.

As shown in Figure , from the free energy profiles we can estimate an activation energy of 42 ± 2 meV for CH3NH3 ann class="Chemical">d of 56 ± 4 meV for CD3NH3, assuming a negligible contribution from tunneling effects. Our theoretical results are thus in excellent agreement with our neutron scattering results as alluded to earlier and the experimental results of ref (32). We predict a difference of around 10 meV in activation barriers between the all-protonated and the deuterated systems, which is solely due to the quantum nature of the nuclei. It is worth noting that this difference is more than double the one predicted for the simple 1D model presented in the last section. This difference is probably due to the coupling of this rotation mode to other vibrations of the molecule (most probably the CH/CD stretches), which go beyond the rigid rotor approximation. Regarding these bonds, we find that on average the CH bonds are 0.01 Å longer than the CD bonds in our PIMD simulations.

Figure 9

Free energy profiles for rotation of the CH3 and the CD3 group. The angle was measured as the dihedral angle defined by the H–C–N–H atoms. The simulation was run at 150 K.

Free energy profiles for rotation of the CH3 ann class="Chemical">d the CD3 group. The angle was measured as the dihedral angle defined by the H–C–N–H atoms. The simulation was run at 150 K.

In conclusion, we investigated the motion of n class="Chemical">CH3NH3+ cations in CH3NH3PbI3 at low temperatures both experimentally and theoretically. By fitting our quasi-elastic neutron scattering data, we obtained an activation energy (Ea) of 48.2 meV, which agrees well with previous experiments.[32] Our experiments suggested that this Ea is associated with the rotation of only CH3 rotors and does not involve NH3 rotors. Our findings are corroborated by (a) the temperature and Q dependence of the incoherent elastic intensity and (b) the fact that neutron scattering of samples with deuterated methyl groups results in noticeably larger activation energies.

Theoretically, we modelen class="Chemical">d the rotation of CH3 rotors using density functional theory and path integral molecular dynamics. For CH3NH3PbI3, the prevalent nuclear quantum effects result in Ea = 42 meV for the rotation of the methyl group. This value agrees well with neutron scattering, with a small difference of several meV. For the deuterated (CD3) rotors, we find much weaker nuclear quantum effects. The calculated Ea = 56 meV perfectly matches the value obtained from inelastic neutron scattering. We therefore ascribe the processes observed with neutron scattering at low temperatures to the rotation of only CH3 around the C–N bond. Since the rotation of NH3 breaks hydrogen bonds, its activation energy is much higher (∼120 meV) and the process has only minor contribution to the neutron-scattering intensities at the temperatures of our experiments.

Our study n class="Chemical">demonstrates that nuclear quantum effects are important for the motion of CH3NH3+ cations in MA-based perovskites, inducing significant isotope effects. Our study provides important insight into the microstructure of MAPbI3 that future research and device design can build on.