Origin of the large entropy change in the molecular caloric and ferroelectric ammonium sulfate

The deceptively simple inorganic salt ammonium sulfate undergoes a ferroelectric phase transition associated with a very large entropy change and both electrocaloric and barocaloric functionality. While the structural origins of the electrical polarisation are now well established, those of the entropy change have been controversial for over fifty years. This question is resolved here using a combination of DFT phonon calculations with inelastic neutron scattering under variable temperature and pressure, supported by complementary total and quasielastic neutron scattering experiments. A simple model of the entropy in which each molecular ion is disordered across the mirror plane in the high symmetry phase, although widely used in the literature, proves to be untenable. Instead, the entropy arises from low-frequency librations of ammonium ions in this phase, with harmonic terms that are very small or even negative. These results suggest that, in the search for molecular materials with functionality derived from large entropy changes, vibrational entropy arising from broad energy minima is likely to be just as important as configurational entropy arising from crystallographic disorder.


Introduction
Phase transitions in orientationally disordered crystals have been known for many years, 1,2 but have recently seen a dramatic resurgence of interest because of their promising caloric and electrical properties. 3 If the entropy associated with the order-disorder transition couples to an external field, such as electric field (the electrocaloric effect) or pressure (the barocaloric effect), then this coupling can be used to produce a solid-state heat pump; these technologies hold great promise to replace environmentally damaging vapour-compression refrigerants. 4 Likewise, if a disordered component has an electrical dipole moment, then dynamically disordered phases are likely to have large dielectric constants and good ionic conductivity, while ordered phases may show macroscopic polarisation and hence pyro-or ferroelectricity. [5][6][7][8]  In the low-symmetry phase, the ammonium ions tilt away from the incipient mirror plane, forming fewer, shorter hydrogen bonds, while in the high-symmetry phase, the ammonium ions sit on this plane (gray lines), forming a more symmetrical arrangement of longer, weaker hydrogen bonds.
Here we revisit one of the oldest known ferroelectrics, the inorganic salt ammonium sulfate, (NH 4 ) 2 SO 4 . This undergoes a first-order phase transition at 223 K and ambient pressure, involving the loss of a mirror plane perpendicular to the crystallographic c axis and, as a result, ferroelectric polarisation along this axis (Fig. 1). 9 This phase transition has been extensively studied for over 50 years. 10 Renewed interest has recently developed in the strong electrocaloric 11 and (inverse) barocaloric effects, 12 due to the unusually high entropy of transition. The caloric behaviour is especially impressive when measured per unit mass, giving an entropy change of ∆S = 60 J K −1 kg −1 , since ammonium sulfate contains no heavy metals and has only a modest density of 1.77 g cm −3 . The corresponding volume change of the unit cell is −4.4Å 3 , corresponding to ∆V = −4.7 cm 3 kg −1 . 12 Since this is rather large, the gradient of the phase boundary |dT /dp| = |∆V /∆S| = 80 K GPa −1 is also impressively steep, giving a large temperature change under modest pressure. There has been substantial discussion of the mechanism of the phase transition, focused mainly on the origin of the electrical polarisation and hence the ferroelectric properties. This question appears to have been resolved by recent single-crystal diffraction measurements, demonstrating that the polarisation is due primarily to distortions of the ammonium ions. The phase transition involves cooperative rotation of the sulfate ions and a concomitant rearrangement of the NH 4 ···SO 4 hydrogen-bonding network. 13 On the other hand, the origin of the caloric properties -the entropy change at the phase transition -remains uncertain. The entropy change was initially explained in terms of an order-disorder model in which the three tetrahedral ions (two ammonium ions and one sulfate ion) in each formula unit are independent "statistical units", showing twofold disorder about the mirror plane above the phase transition but ordering below it. 14 In this case, the Boltzmann formula for configurational entropy predicts a transition entropy of (3 ln 2)R = 17.3 J K −1 mol −1 , in good agreement with the experimental value of 17.6 J K −1 mol −1 . 15 Spectroscopic data were also interpreted as favouring an order-disorder model. 16,17 This model, however, is unsatisfactory for two reasons. First, it is not consistent with the substantial body of accumulated crystallographic data on this material, since there is almost no crystallographic evidence for a disordered structure in the paraelectric phase. One exception is an early X-ray study, which appeared to show some evidence of hydrogen disorder, 18 but this was neither reproduced by the recent X-ray measurements 13 nor shown in neutron experiments, which are far more sensitive to scattering from H atoms. 19 By contrast, crystallographic disorder is clearly visible in materials with similar structures, including (NH 4 ) 2 BeF 4 20 and (N(CH 3 ) 4 ) 2 ZnCl 4 . 21 Similarly, in vibrational spectroscopy data, an order-disorder transition would be expected to produce a sudden change in the width of relevant peaks at the phase transition, whereas Raman spectroscopy shows no such change. 22 Second, the simple order-disorder model makes no attempt to account for dynamics, which will certainly contribute to the entropy change. 23 An alternative explanation is therefore needed for the remarkably large entropy change in ammonium sulfate. This problem was considered very recently by Born-Oppenheimer molecular dynamics simulation 24 but has not yet been tackled experimentally.
To resolve this question, we have studied the local structure and dynamics of ammonium sulfate computationally, by density-functional theory calculations of the phonon spectrum; and experimentally, by inelastic, quasielastic, and total neutron scattering, as functions of temperature and pressure. Inelastic neutron scattering is sensitive to the vibrational dynamics of all atoms in the material; quasielastic neutron scattering reveals the reorientation dynamics of the ammonium ions; and total neutron scattering reflects the correlations between instantaneous atomic positions. This combination of techniques therefore provides access to all sources of our target material's entropy.
In Section 2, we show that the large entropy change in ammonium sulfate cannot meaningfully be described as configurational (i.e., order-disorder). It is better understood as dynamic, resulting from molecular librations that correspond to motion in shallow energy wells in configuration space. In Section 3, we investigate the atomic origins of this energy surface, demonstrating that it results from competition between two different arrangements of the hydrogen bonds between ammonium and sulfate ions.
2 Is the entropy change configurational or vibrational?

Evidence for vibrational entropy
A DFT model of each phase was constructed starting from the experimentally determined structures (Fig. 1). In each case, the two ammonium ions in the formula unit are crystallographically distinct, referred to as N1 and N2. The network of hydrogen bonds around each of these ions is different between the two phases: the low-symmetry phase has fewer, shorter and therefore stronger hydrogen bonds, while the high-symmetry phase has a more symmetrical arrangement with more, longer, weaker hydrogen bonds.
Allowing the unit cell parameters to relax gave reasonable agreement with experiment in the Pna2 1 phase, but an unacceptably large divergence from the experimental structure in the Pnam phase. In particular, there was a 9% expansion along the c axis compared to experimental data at 233 K, which we attribute to neglect of thermal expansion, which will be significant in this material. For our purposes, however, it is sufficient to adopt the quasiharmonic approximation, fixing the cell parameters at their experimental values at temperatures of 5 K (our own work; see below) and 233 K, 25 respectively.
Next, the phonon spectrum was calculated using density-functional perturbation theory. The phonon density of states is shown in Fig. 2a. To characterise individual modes, the GASP algorithm [26][27][28] was used to decompose the DFT eigenvectors at the Γ point into translations, librations, and distortions of the two crystallographically independent ammonium ions and the sulfate ion ( Fig. 2b and c). * A full comparison between these calculations and experimental data will be deferred to Section 3. For now we note that the simulation agrees well with the experimental phonon density of states from inelastic neutron scattering measurements (Fig. 2a). There is excellent agreement up to E = 34 meV, with the dip at 18 meV being accurately reproduced, and only a slight shift in energy between the data and calculation for the next clear dip at 25 meV and 24 meV respectively. Above 34 meV, the calculations show there are no more translational modes. The peaks in the gDOS around 40 meV were identified as ammonium librations as opposed to translations on the basis of their negative Grüneisen parameters (Section 3). The calculations slightly overestimate the energies corresponding to ND 4 libration modes, then slightly underestimate the energies of the internal modes of both the ND 4 and SO 4 ions.
This analysis reveals that the key vibrational difference between the phases is the librational motion of the ammonium ions (red and blue × shapes, Fig. 2b and c). In the Pna2 1 phase, these occur in the 34-50 meV region, where they are well separated from molecular translational modes below about 30 meV. By contrast, in Pnam, the librations mix with the molecular translational modes. Even more strikingly, these decrease in frequency -some modestly, others substantiallysuch that the lowest-frequency of these librations become weakly unstable, within the harmonic approximation, in some or all of the Brillouin zone.
The dispersion curves in the Pna2 1 phase are shown in Fig. 3a, plotted along the path in Fig. 3b. To account approximately for the instabilities in the Pnam phase, the energy of each unstable mode was mapped as a function of the mode coordinate Q at the gamma point. This revealed flat-bottomed, "bathtub"-shaped potentials with very low harmonic terms; the energy is therefore dominated by quartic and higher-order terms (Fig. 3c). Solving the Schrödinger equation for the resulting potentials gave a renormalised frequency that now varies with temperature, 29 falling in the range 11 meV to 16 meV at T = 235 K. For practical reasons, we made the strong assumption that these effective frequencies remain roughly constant across the Brillouin zone. Fig. 3d shows the dispersion curve after this renormalisation; the original, showing harmonic instabilities, is given in Fig. S1.
Finally, the entropy of transition was calculated using the standard relationship 30 where x i =hω i /k B T runs through harmonic vibrational modes with angular frequencies ω i . This in turn allows the Helmholtz free energy A = U − T S to be calculated as a function of temperature for each phase, and the phase transition temperature to be estimated by setting the free energies of each phase to be equal. (The p∆V contribution is negligible at ambient pressure.) This calculation gives an entropy change of 13.2 J K −1 mol −1 across the phase transition, 75% of the experimental value of 17.6 J K −1 mol −1 . The phase transition temperature is predicted to be 1147 K, substantially above the experimental value of 223 K. The cumulative difference in entropy at the phase transition temperature is shown in Fig. 3e. We consider these values to be reasonably consistent with experiment, since the agreement is constrained by the relation ∆U = T ∆S. The DFT model gives ∆U = 144 meV per formula unit while experimentally T ∆S = 40 meV. A 0.1 eV difference, of the order expected for DFT accuracy, requires T and ∆S to differ from their experimental values by a combined factor of 3.6. The accuracy obtained in practice is therefore as good as can be expected from this simple model. Specifically, the internal energy difference between the two phases is reproduced with only modest accuracy, which is perhaps expected given known limitations of DFT and the substantial differences in hydrogen bonding geometry between the phases. On the other hand, the difference in the shapes of the energy surfaces about the respective minima is relatively well reproduced.
To elucidate further the DFT ∆U value, we used GULP 31 to calculate the Coulomb energy of a hypothetical point-charge model in which the full +1 charge of the ammonium ion resides on the N atom and the −2 of the sulfate ion on the S atom. This gives an energy difference of 77.5 meV per formula unit purely from long-range ionic interactions, without considering changes in covalent or hydrogen bonding; this accounts for roughly half of the DFT energy difference.
One possible source of discrepancy between the experimental and calculated ∆S values is the strong approximation that the renormalised frequencies are independent of wavevector. Indeed, only a 17% decrease in the effective frequencies of just the eight renormalised modes would be sufficient to increase the entropy change to its experimental value; this is likely to be within the error of the approximation. However, extending this calculation to consider points away from the zone center and map changes in anharmonicity across the Brillouin zone would be computationally expensive, requiring phonon calculations on large supercells.
Nonetheless, it is clear that the ammonium librational modes differ between the two phases in a way that can account for a large difference in vibrational entropy. In Section 3, we investigate further the atomic features of these modes that are responsible for this behaviour. Before proceeding to this, however, we show that, by contrast, a configurational model cannot account for the observed entropy change.

Evidence against configurational entropy
Considering possible evidence for configurational entropy, we first note that the established structural model of the Pnam phase is fully ordered, 13,19 and thus has a configurational entropy of zero. Our own single-crystal X-ray diffraction data in the high-temperature phase are in good agreement with this published structure. To complement this traditional crystallographic analysis, we therefore measured the pair distribution function, which represents instantaneous local correlations between pairs of atoms and is thus sensitive to local deviations from the crystallographic average.
Total neutron scattering data were collected from a powder (ND 4 ) 2 SO 4 sample using the POLARIS diffractometer (ISIS, U.K.) at several temperatures spanning each phase. First, the crystallographic average structure in each phase was determined by Rietveld refinement of the established, fully ordered, structural model. The refined models gave good fits to the data and agreed well with literature values, without the need to invoke crystallographic disorder.
Next, we used reverse Monte Carlo refinement, as implemented in the RMCprofile code, 32 to derive an ensemble of models, each representing an instantaneous "snapshot" of a 6 × 5 × 8 supercell. In this method, the models are refined simultaneously against the pair distribution function, the scattering function, and the Bragg profile; thus they are consistent with both the local structure and the long-range average structure. Sample fits are shown in Fig. S2; a full analysis of these data is in preparation.
The directions of the N-H bonds were extracted from these models. Plotting these as a spherical histogram in orthographic projection shows isolated peaks in each phase, with no sign of the bimodal shape that would be expected from disorder about the mirror plane (Fig. 4a). However, these distributions are far broader in the high-symmetry than in the low-symmetry phase. To investigate this further, we again used the GASP algorithm to find the rotation that best related each tetrahedron in the configuration to its average crystallographic orientation. This confirmed that the ammonium ions have far more librational freedom in the high-symmetry phase (Fig. 4b). These results are thus consistent both with the established Figure 2: (a) A comparison of the measured (T = 10 K) and calculated generalised phonon density of states of (ND 4 ) 2 SO 4 in Pna2 1 at ambient pressure. The data were measured with four different incident energies E i = 23, 36, 67 and 162 meV. The overlap between these data sets was used to estimate the background in each case and to scale the y-axis by an appropriate factor, allowing us to produce a "continuous" trace for the combined data. The calculated data have been broadened by a gaussian representing the instrumental resolution for the incident energies used. Specific phonon modes have been labelled on the data and simulation. The calculated phonon modes at Γ were separated into librations, bends and stretches of the different molecules using GASP. The energies and displacement magnitudes for the different modes in Pna2 1 and Pnam are shown in (b) and (c). structural model and our conclusion above that ammonium librations are more thermally accessible in the high-symmetry Pnam phase, but inconsistent with the simple model where the entropy of transition is purely configurational. These results are furthermore consistent with our recent measurements of the rotational dynamics of ammonium ions in this material using quasielastic neutron scattering. 33 This technique is sensitive specifically to the ammonium ions because it measures the incoherent scattering, which in this sample is dominated by 1 H; it reflects the "hopping" motion of these ions between energy minima corresponding to different orientations. We refer the interested reader to that paper for full details, but repeat here the two key results relevant to the entropy change.
First, the proportion of elastic scattering as a function of Q, known as the elastic incoherent structure factor (EISF), is determined by the geometry of hopping; that is, by the equilibrium distribution of the moving 1 H nuclei. Hopping between two closely separated sites, on opposite sides of the mirror plane in the high-symmetry phase, would be clearly visible in this signal. However, no such hopping was observed. Again, then, these data are not consistent with the entropy change being caused by an order-disorder transition.
Second, the frequency of reorientation is indicated by the width of the characteristic Lorentzian peak in energy transfer. In the low-symmetry phase at 200 K, this frequency unexpectedly increased with pressure as the phase transition was approached. In a simple model where the ammonium ions move in a static energy landscape, this indicates that the energy barrier to rotation is decreasing. This is consistent with the decrease in this energy barrier with temperature reported on the basis of spin-lattice relaxation times. 14 Of course, in reality the surrounding ions also move, and the observed behaviour may equally be attributed to the barrier collapsing more frequently; this subtlety does not, however, affect our qualitative argument.
An increase in hopping frequency with pressure is expected for single-atom hopping, for instance in asymmetric hydrogen bonds, since forcing two wells closer together will lower the barrier between them. 34 It is, however, very unusual for rotations of whole molecules such as ammonium ions, where hopping typically decreases in frequency with applied pressure. 35 We propose that this highly atypical behaviour reflects a flattening of the energy landscape as the phase transition is approached, both because the existing hydrogen bonds are destabilised by compression beyond their equilibrium length, and because the alternative high-symmetry arrangement with longer bonds is stabilised for the same reason. This result suggests that the ammonium librations may be important in determining the entropy change, an idea we return to in the following section. This change in the energy landscape is also reflected in the volume decrease across the phase transition. Similar behaviour was recently reported in the related material ammonium thiocyanate. 36 Our results so far show that the entropy of transition is likely to be mostly vibrational rather than configurational, due to the librational motion of the ammonium ions. Next, we performed a detailed investigation to characterise why these particular vibrational modes give such a large entropy change.
3 Atomic characterisation of the high-entropy vibrations 3

.1 Ambient-pressure single crystal inelastic neutron scattering
To characterise the vibrational behaviour in more detail, we collected inelastic neutron scattering from a 1.2 g mosaic sample prepared by gluing individual (ND 4 ) 2 SO 4 crystals to aluminium plates (Fig. S6). Inelastic scattering techniques are an excellent probe of coherent excitations such as phonons, while neutrons show a particular benefit when it comes to the study of low-Z elements, as present in ammonium sulfate. In contrast to Raman scattering, inelastic neutron scattering allows the dispersion of phonons to be mapped over a wide range of momentum and energy transfer space. In the particular case of ammonium sulfate, there is no change in systematic absences between the space groups Pnam and Pna2 1 , so we expect the changes in scattering to be more subtle than a mode simply lifting away from the elastic line. For this reason, collecting data over a large region of 4-dimensional Q-E space is particularly important. Indeed, little change in the relative intensities of the different Bragg peaks is observed in the elastic data, as expected given the subtle difference between the two known structures (Fig. S7).
Moving away from the elastic line, acoustic phonons emerge as cones of scattering from the Bragg peaks, appearing as rings of intensity for energy transfers between 4 and 8 meV (Fig. 5a). When the DFT results are weighted by the neutron scattering lengths of the respective nuclei and convolved with the instrumental resolution function, they agree reasonably well with the experimental data (Fig. 5e). In particular, the simulations accurately reproduce the asymmetric shape of the rings of scattering, and hence the dispersion curves of the acoustic modes. For instance, the different orientations of the The distribution is clearly centred on this plane, rather than having separate peaks above and below it. (b) Distribution of ammonium ion rotations φ from the crystallographic average position, from GASP analysis of these configurations. The distributions evolve slightly with temperature, especially in the low-symmetry phase, but differ far more dramatically between phases, again demonstrating the increased freedom of motion in the high-symmetry phase that is responsible for its increased entropy. These data have been normalised by the solid angle differential sin φ dφ dθ . In both parts, there are 9600 ammonium ion geometries, corresponding to five independent RMC runs. The simulation for the cut in Pnam has been horizontally offset by 0.5 units along the x-axis for clarity (h), but no offset has been applied to the equivalent data set (d), which is genuinely more intense than seen for Pna2 1 .
rings around (0, 10, −2), (0, 7, −5), (0, 3, −5), (0, 0, −4), and (0, −3, −5) are well reproduced, as well as the regions of higher scattering intensity at the top of (0, 3, −5) and the top and bottom of (0, −3, −5). Notable differences include first, that many of the rings associated with weaker diffraction peaks are not clearly visible above the noise in the experimental data. Second, the rings from several peaks on the edge of the observable region of reciprocal space, in particular (0, 6, −6) and (0, 10, −4), are predicted to be more intense than seen in the data; we tentatively attribute this to masking of unreliable detectors at the edge of the array. Third, the calculations also predict an intense broad diagonal scattering peak in the region 5 < K < 10, which is not apparent in the experimental data; this may partly be the result of absorption causing a dark ring cutting through the data set at this point. We suggest that all of these discrepancies in intensity at very low energy transfer reflect the intensities of the relevant Bragg peaks more than the dispersion relations. Further comparisons of data and simulation at different energy cuts are given in Fig. S8. Next, we mapped the scattered intensity as a function of energy transfer along the c * axis in each phase. In the low temperature phase data (Fig. 5b), the acoustic modes are clear, emerging from the Bragg peaks on the elastic line at (0 0 L) for L even. The scattering is more intense as L increases, since the 1-phonon scattering intensity includes a factor of |Q| 2 . The bright signal at the base of the plot is bleeding coming from the elastic line due to the colour scale required to make the inelastic signal visible. The agreement between the data and simulation (Fig. 5f) for the acoustic modes is excellent, and it is also reasonable for the optic modes. There is a notable discrepancy at the elastic line, where scattering from (0, 0, L) peaks with L odd is observed despite the fact that these peaks are systematically absent in both Pna2 1 and Pnam. We attribute this to multiple scattering from the large mosaic sample; this effect will be most visible at the elastic line and is not expected to contaminate our data elsewhere.
The agreement can be seen more clearly from a one-dimensional cut through the data at the Z-point (0 0 −3.5) as a function of energy transfer. The blue trace in Fig. 5d shows the tail of the elastic line at low energy transfer, then two broad peaks in the scattering at 12 meV and 22 meV, separated by a dip at 17 meV. Comparing this with the equivalent cut through the simulated scattering (Fig. 5h) including the instrumental resolution, the elastic line tail is absent, since only the inelastic scattering has been calculated, and the trace is more structured, but it too has a broad peak at around 11 meV and a dip at around 18 meV. The calculation has two peaks at 21 and 25 meV as opposed to the single broad peak at 22 meV in the data, but some broadening of the data could mask individual peaks.
The inelastic signal at T = 235 K in the Pnam phase is much more intense (Fig. 5c) due to the increased phonon population associated with the Bose-Einstein distribution, and also appears to be more diffuse, making it harder to identify phonon modes above the low energy acoustic modes. In addition, the multiphonon scattering that contributes to the background increases in intensity with increasing temperature. The equivalent simulation of the scattering for Pnam including the instrumental resolution (Fig. 5g) reveals considerably more structure than observed in the data. The steeply dispersing acoustic modes at the lowest energy transfers are reproduced, but the features above ∼ 10 meV cannot be resolved. One possibility is that, in the higher temperature phase, phonon-phonon interactions reduce the phonon lifetime, which broadens the spectral lineshape, smearing out the signal compared to the calculations performed in the harmonic approximation. However, the fact that well-defined phonons are observed at all is again evidence that there is no orientational disorder, which would lead to very significant broadening from over-damping. Moreover, it is clear that the gap in the scattering at around 18 meV, seen in both the data and simulation in Pna2 1 , is absent in both in Pnam. Again, some signal at systematically absent (0, 0, L) peaks (L odd) is visible at the elastic line, attributable to multiple scattering.
Comparing the same one-dimensional cut as a function of energy transfer at the Z-point (0 0 −3.5) (Fig. 5d), the data in Pnam shows no clear features, just a slowly decreasing count rate with increasing energy transfer. Meanwhile, the calculation shows a series of peaks (Fig. 5h). The calculated peak below 5 meV would be obscured in the data below the elastic line tails, but the calculated peak at 8 meV may be reflected in the data. The broader features seen in the calculation are not visible in the data, but this is likely to be a consequence of the reduced phonon lifetime at higher temperatures.
Our data cover many Brillouin zones, such that the dispersion can be plotted along high symmetry directions other than the z-axis, but this shows the features most clearly. Slices along other high-symmetry directions compared with calculated phonon dispersion curves are shown in Figs. S9 and S10.

Powder inelastic neutron scattering under hydrostatic pressure
Finally, INS measurements were performed on a powder (ND 4 ) 2 SO 4 sample in a clamp pressure cell. Because of powder averaging, this technique gives data only as a function of the modulus of the scattered momentum |Q|. Clear bands of inelastic scattering intensity are seen at base temperature (10 K) and ambient pressure, separated by gaps at 5 and 18 meV (Fig. S11), reminiscent of that seen at Z in the single-crystal data (Fig. 5d). The 18 meV gap is also clearly visible in the generalised phonon density of states (gDOS) calculated according to Eq. 7 (Fig. 2a).
The DFT calculations allow the modes in each phase to be described in detail, consistent with but expanding previous assignments from Raman scattering. 22,37 In the low-symmetry Pna2 1 phase (Table S3), the phonons up to around 30 meV involve whole-body motion, with translations of both ammonium and sulfate coupled with sulfate librations. There are 45 of these modes: 8 ammonium ions with 3 translational degrees of freedom, and 4 sulfate ions with six translational and librational degrees of freedom, minus the 3 zero-frequency modes corresponding to translation of the whole crystal. Because the primary motion is molecular translation, these appear in Fig. 2b as modes with relatively low contributions from libration, bond stretching or angle bending.
The calculated phonons between 18 and 30 meV include several ammonium "rattling" modes. In our powder data, the localised nature corresponds to a vibrational energy that is independent of momentum transfer, so that these appear as flat bands of scattering in S(Q, E), resulting in sharp features in the gDOS (Fig. 2a) between 18 and 30 meV. While the general trend for phonon scattering intensity varies as |Q| 2 , collective motion will result in a departure from this simple relationship. That the measured intensity of these modes varies simply as |Q| 2 (Fig. S11) implies a lack of phase coherence, further evidence for identifying these as "rattling" modes. 38 The calculated modes between 34 and 50 meV are, with one exception, ammonium librations. † The librational modes group together into sets of four adjacent modes, all of which involve the same set of ammoniums rotating about the same axis. For example, the modes at 34.2, 34.3, 34.5 and 34.7 meV all involve the N2 ammonium rotating about the c-axis. The fact that all four have such similar energies means, in effect, that the four motions are almost independent.
Above the ammonium librational modes, our Pna2 1 calculations identify three sets of internal modes for the tetrahedral sulfate ions:  (Fig. S12). The next 36 modes between 120 and 142 meV are SO 4 ν 3 modes coupled with ND 4 ν 4 modes, followed by 16 ND 4 ν 2 modes, 8 ND 4 ν 1 modes and 24 ND 4 ν 3 modes. In these internal ND 4 modes, N1 and N2 are degenerate, unlike in the lower-energy librational modes (Table S1). Neither the ammonium "rattling" modes nor the internal modes of the ammonium and sulfate ions change substantially across the phase transition in our calculation (Table S4).
Due to finite phonon lifetimes, peaks in the density of states become increasingly broad and indistinct as the temperature increases (Fig. S13), again consistent with previous Raman results. 37 To investigate the pressure dependence of the inelastic scattering, we therefore concentrated on low temperature data and hence the Pna2 1 phase. Figure 6a shows the generalised phonon density of states measured at base temperature (T = 10 K) at P = 0.0, 2.6 and 4.8 kbar. This reveals that not all of the phonon modes exhibit the same pressure dependence. The majority of the modes at lower energies stiffen with pressure, as might be expected, since the atoms are brought closer together, increasing the forces between them. However the modes at 35 − 45 meV soften with increasing pressure, indicating a negative mode Grüneisen parameter, while the mode at 34 meV shows no pressure dependence.
The same behavior is observed in our simulations for the Pna2 1 phase (Fig. 6b), with the whole body motion modes mostly stiffening, the polar translation mode at 33 meV unaffected, and the ammonium libration modes having a negative Grüneisen parameter. Extending these calculations to different points in the Brillouin zone (Fig. 3a) gives similar results. The more dispersive, lower-energy translation modes have either zero or positive mode Grüneisen parameters, with the exception of the lowest-frequency branch along the path Γ-X-S: this involves transverse motion (i.e., in the c direction) of the layers of molecules in the ab plane, and has a mildly negative mode Grüneisen parameter. On the other hand, the flatter ammonium librational modes between 38 and 45 meV show negative values of the mode Grüneisen parameter. We † The exception occurs due to longitudinal optic/transverse optic splitting of the translational modes just below this frequency range, in which ammonium and sulfate ions move in opposite directions. The electric polarisation arising from this motion increases the frequency of the respective longitudinal branches in the limit as Γ is approached from X, Y , and Z, so that along each of these branches, that single translational mode has a higher energy than the lowest-energy ammonium libration modes. rationalise the difference in these parameters by observing, analogous to the behaviour of framework solids predicted by the Rigid Unit Mode model, 40 that the ammonium librations involve transverse motion with respect to the N -H···O hydrogen bonds, and a reduced volume makes "buckling" more favourable, lowering the frequency. Negative mode Grüneisen parameters may lead to negative thermal expansion, as is observed in ammonium sulfate just below the phase transition. 41 In Pnam (Fig. 3b), by contrast, almost all of the modes have positive Grüneisen parameter. Although the eight potentially unstable modes for which we calculated effective harmonic frequencies at Γ appear similar, they have rather different mode Grüneisen parameters. For five of the eight modes, these are positive; for two, they are negative only at very low temperature but become positive thereafter; and for a single mode, the mode Grüneisen parameter remains strongly negative (γ = −4.7) up to the experimental decomposition temperature of 523 K. Thus, although modes involving libration of the same molecules about the same axis are similar in energy at ambient pressure, this result demonstrates that these modes couple very differently to strain, depending on the relative phase of the librations.

Discussion and Conclusions
Combining our experimental data and DFT model, we conclude first that there is no central barrier to the ammonium librations in the high-symmetry phase. Thus this phase is far from the classical disordered limit that would be accurately represented by a crystallographic split-site model. In particular, the Boltzmann entropy formula S = R ln n does not apply here, and we suggest that the similarity of the measured entropy to (3 ln 2)R is ultimately a coincidence. This conclusion is in line with the warning given sixty years ago by Guthrie and McCullough, two pioneers of research on orientationally disordered crystals, that "speculations based mostly on thermal data are [best] avoided in public"! 1 Instead, the entropy is best understood as vibrational rather than configurational. Our results demonstrate that a flat-bottomed potential in which there is little resistance to small librations is consistent with the known structure and capable of producing an entropy of the correct magnitude. Such flat potential minima are reminiscent of the well shapes in the crystalline phases of inorganic semiconductor phase-change materials. 42,43 The structural origin of this dynamic entropy lies in the contrast between the relatively rigid hydrogen-bonded network in the low-symmetry phase -as demonstrated by the harmonic librational frequencies in DFT, and by its experimental destabilisation under pressure -and the more loosely held network in the high-symmetry phase that supports low-frequency librations. We note that, although one might loosely attribute these low frequencies to "anharmonicity", it is strictly speaking the small harmonic term in the potential expansion rather than the larger anharmonic terms that are responsible: indeed, a lower frequency, and higher entropy change, still would be possible if the quartic terms were lower.
Second, we have determined in detail the phonon modes of ammonium sulfate.There is excellent agreement between our experimental and simulation results in the low-temperature phase, even in the full four-dimensional Q-E space made available by the single-crystal data, suggesting that this methodology gives a reliable description of both internal and intermolecular contributions to the crystal energy. In the high-temperature phase, low phonon lifetimes obscure the experimental data and comparison becomes more difficult, but the agreement still appears to be reasonable. The good agreement at low temperatures suggests that the DFT simulations are a reliable guide to the phonon dispersion, even at temperatures where this is unclear from the experimental data. The frequency of the ammonium librational modes is overestimated by DFT; these could nonetheless be conclusively identified in experiment by their negative Grüneisen parameters. The frequency of the internal ammonium and sulfate modes is slightly underestimated but, again, identification is clear because these are isolated in frequency.
Our results suggest a structural mechanism for the distinct change in the phonon spectrum between the Pna2 1 and Pnam phases of ammonium sulfate that is responsible for the high entropy of the Pnam phase and hence for the caloric effects in this material. In the low-entropy Pna2 1 phase, the ammonium librational modes from both experiment and simulation form distinct peaks, separate from the other low-frequency collective modes, indicating that hydrogen bonding holds these ions relatively firmly in place. Strong hydrogen bonding also means that the librational ammonium motion is able to drag neighbouring sulfate ions along with it (consistent with a proposed mechanism for the phase transition 25 ), causing the lattice to contract with temperature. This mechanism therefore also explains the negative mode Grüneisen parameters that are, again, observed in both experiment and simulation; similarly, it suggests that these modes are also likely to be responsible for the region of negative thermal expansion just below the phase transition. By contrast, in the high-entropy phase where the ammonium ions sit on the newly created mirror plane, they are less firmly held in place, leading to lower vibrational frequencies that exert less influence over the shape and size of the lattice.
In this work we have not considered the contribution of the phonons to the thermal conductivity, but this too will be important to practical applications, reinforcing the importance of studying the phonon behaviour under working pressure conditions. Ammonium sulfate is a particularly clear example of the perils of ignoring dynamic contributions to entropy, since, as we have shown, the contribution from configurational disorder appears to be negligible. However, these results are also relevant to other materials where configurational contributions are also important. Many materials that undergo high-entropy phase transitions -including the globular organic "plastic" crystals and the molecular perovskites -have an orientationally ordered state held together by weak interactions such as hydrogen bonds and a disordered state in which these bonds are broken. While it is convenient for the purposes of entropy calculations to consider "pure" order-disorder transitions, the distinction between these and displacive phase transitions is in reality a continuum. 44 It thus seems likely that the Boltzmann formula gives a misleading picture of the true origins of entropy in many or even most molecular crystals. The example of ammonium sulfate further suggests a different paradigm for crystal engineers to target in the search for high-entropy and caloric materials: one where competing networks of hydrogen bonds or other weak interactions instead create a complex, high-entropy energy landscape.

Methods Calculation
We used density functional theory to calculate the structures and phonon dispersion relations of both phases of ammonium sulfate, as implemented in the CASTEP software package, v.19.11 academic release. 45,46 We used a plane-wave basis set together with norm-conserving pseudopotentials from the CASTEP standard library. The Perdew-Burke-Ernzerhof functional was used to describe the exchange-correlation energy, 47 with a Tkatchenko-Scheffler empirical dispersion correction to account for the van der Waals interaction between the molecular ions. 48 In the geometry optimization, the energy cutoff was set to 1100 eV, and a 2 × 1 × 2 Monkhorst-Pack k-point grid used. The structures were relaxed until both the forces between atoms became smaller than 0.01 eV/A and the energy change between steps was less than 2 × 10 −5 eV/atom. For the phonon calculation, the force matrix was calculated by density-functional perturbation theory, using a finer 3 × 3 × 3 phonon k-point grid. As discussed in the Results section below, in the Pnam phase using experimental lattice parameters (although not when using relaxed lattice parameters, nor in the Pna2 1 phase) this gave several unstable modes with ω 2 < 0 in the harmonic approximation. To take these modes into account, we mapped the energy V (η) as a function of the normal coordinate η of each mode individually at Γ, then fitted the results to a polynomial potential model. To find the energy eigenstates of this model, we followed the approach of Skelton and co-workers, 29 using the Fourier grid Hamiltonian method 49 to solve the relevant Schrödinger equation 50 numerically: From these eigenstates, we calculated the temperature-dependent partition function and hence an effective harmonic frequency, which is therefore also now a function of temperature 50 p. 226: For simplicity, and because following a mode in this way at an arbitrary point of the Brillouin zone requires a costly supercell calculation, we assumed that these effective frequencies were constant with wavevector. Using a 2 × 1 × 1 supercell to investigate the X point suggested that incorporating more points in reciprocal space has a negligible effect within the error of our model, and does not change our qualitative conclusions.
To compare the DFT calculations and the experimental INS data, we used Euphonic, 51 a Python package that efficiently calculates phonon bandstructures and inelastic neutron scattering intensities from a force constants matrix using the 1-phonon scattering function, such that one can obtain neutron weighted phonon dispersions and the density of states from precalculated phonon frequencies. This can have a stark effect on the range of visible modes in the simulated phonon dispersions, both due to the variation in coherent scattering cross-sections for different elements, and the polarisation factor Q · e in the 1-phonon coherent scattering cross-section, where e is the phonon mode eigenvector. We also used TobyFit within Euphonic to include the experimental resolution function in our simulations.
The mode Grüneisen parameters were calculated according to the standard equation For the phonon calculations performed using experimental rather than optimized cell parameters, we were therefore unable to evaluate this expression under hydrostatic pressure; instead, we applied uniform strains of ±0.2% in each dimension, using the central difference algorithm to numerically evaluate the derivative in (6). Using uniform strain instead of hydrostatic pressure is equivalent to assuming that the material is elastically isotropic, which -especially at the qualitative level of our analysis here -is comparable to the other approximations involved in this calculation.

Inelastic neutron scattering
Ammonium sulfate of natural isotopic abundance was purchased from Sigma Aldrich and deuterated by four successive recrystallisations from D 2 O. Assuming uniform mixing of hydrogen atoms, this should give a 97.7% deuterated sample. In practice, the incoherent background from 1 H was not observed to any significant extent in our neutron scattering data. Single crystals with dimensions ∼ 5 × 30 × 1 mm 3 were obtained by precipitation from a saturated deuterated aqueous solution of (ND 4 ) 2 SO 4 . Two sets of inelastic neutron scattering experiments were performed on fully deuterated ammonium sulfate using the Merlin spectrometer at ISIS. 52 First, single crystal measurements were performed at ambient pressure using an array of crystals of total mass 1.2 g. These were coaligned by crystal habit, confirmed using the ALF alignment facility at ISIS, and attached using cytop to two aluminium plates so that sample covered an area of 30 mm × 40 mm, oriented with a horizontal scattering plane defined by the b and c crystallographic axes and a mosaic spread of less than 5 degrees (Fig. S6). The plates were then stacked and inserted into the CCR.
In the second set of experiments, a 1.3 g sample of polycrystalline (ND 4 ) 2 SO 4 was loaded into a newly commissioned, low-background TAV6 cylindrical clamp pressure cell, 53 using helium as the pressure transmitting medium.This cell has a cylindrical geometry, with walls of 5.3 mm and inner bore with a 7.0 mm diameter. The cell was inserted into a cylindrical shaped radial collimator coated with gadolinium paint, to minimise the background scattering from the pressure cell. This in turn was mounted in a closed cycle refrigerator, allowing data to be collected in the ranges T = 10-280 K and P = 0-4.8 kbar.
Merlin was operated in multi-rep mode, with incident neutron energies of 23, 36, 67 and 162 meV. The raw data were processed using Mantid, 54 following the standard conventions. 55 For the single-crystal measurements, the sample was measured in discrete angular steps (0.5 • ) about the vertical axis [H00] to form a "Horace-scan" over 90 degrees in total, enabling a large region of energy-reciprocal space to be explored. The processed single crystal data was then combined using Horace 56 to create S(Q, ω). Attempts to subtract off the background scattering arising from the aluminium plates proved unsuccessful, but there is little evidence of inelastic scattering from the aluminium in our data, just streaks of elastic scattering which can easily be differentiated from the Bragg peaks from the sample.
For the polycrystalline measurements, background data were collected by measuring the TAV6 cell filled with helium, under the same pressure and temperature conditions. ‡ This was then subtracted from the relevant data sets. The powder data were also corrected for the absorption from the TAV6 cell to give the inelastic scattering function S(Q, ω).
The phonon density of states is calculated in Mantid using the 1-phonon scattering function formula in the incoherent approximation: 57 S inc (Q, E) = exp(−2W (Q)) ‡ It is essential to include the helium in the background measurements, since it has a strong phonon signal in the solid state, and a weaker signal in the liquid state.
where the term in square brackets is the calculated-neutron weighted density of states, and g j (E) is the partial density of states for each element j in the material, and m j is the relative atomic mass of the component. The average Debye-Waller factor exp(−2W (Q)) is calculated using an average mean-square displacement u 2 , using Total scattering and reverse Monte Carlo analysis Neutron total scattering data were measured from a perdeuterated powder sample of ammonium sulfate on the Polaris diffractometer at the ISIS Neutron and Muon Source, U.K. The sample was loaded in a vanadium can of 6 mm diameter within a CCR. A series of 15-minute short runs were measured to extract the lattice parameters, and 6-hour-long runs were measured at six temperatures (150 K, 180 K, 210 K, 240 K, 270 K and 300 K) for total scattering data measurements. The data were focused using Mantid. 54 The Reverse Monte Carlo simulation was performed using the RMCProfile code. 32 Configurations of each phase used a 6 × 5 × 8 supercell containing 14 400 atoms, initially arranged according to the average crystal structure from Rietveld refinement. The RMC simulations were carried out for at least 5 × 10 6 steps in total (at least 1.2 × 10 6 accepted moves) until convergence. To prevent ions from distorting to physically unreasonable shapes, internal molecular bond and angle potentials, taken from the MM3 parametrisation, 58 were used. To avoid unphysical spikes in the partial PDFs, a soft "curvature" restraint was placed on the second derivative. The large entropy change in ammonium sulfate, historically analysed in terms of an order-disorder phase transition, is shown instead to arise from low-frequency librations of ammonium molecules in flat-bottomed energy wells. This suggests that similarly competing hydrogen-bonded networks may be an attractive target for engineering new caloric molecular-ionic materials.