Anharmonic lattice dynamics of superionic lithium nitride

Superionic crystals reach an ionic conductivity comparable to liquid electrolytes following a superionic transition at high temperature. The physical mechanisms that lead to this behaviour remain poorly understood. It has...


Introduction
Superionic conductors are solids that reach a liquid-like ionic conductivity at a temperature T s before reaching their melting point [1,2]. This makes superionic conductors promising candidate solid electrolytes for appli-cations such as fuel cells and all-solid-state batteries [3,4,5]. The fundamental physics of superionic phasetransitions, and the mechanistic details of ion transport in the superionic regime, remain poorly understood, which limits the scope for the rational design and optimisation of superionic solid electrolytes for specific applications. A better understanding of the physical basis of superionic phase-transition behaviour and of the factors that cause some materials to exhibit superionic phases is therefore desirable both as an intellectual endeavour and as the basis for developing improved solid electrolytes for practical applications [6].
Superionic conductors that exhibit a superionic transition can be classified type-I or type-II according to the nature of the transition [1,7]. Type-I superionic conductors, such as AgI, exhibit a discontinuous superionic transition, accompanied by a first-order structural phase transition at T s from a poorly-conducting phase to the superionic phase [1]. In contrast, type-II superionic conductors, such as β-PbF 2 , exhibit a continuous superionic transition, accompanied by a second-order phase transition associated with a constant-pressure heat capacity peak at T s [1,8]. For both type-I and type-II superionic conductors, the superionic transition is associated with a marked increase in ionic conductivity as well as a decrease in ionic conductivity activation energy, which is often interpreted as evidence of a qualitative change in diffusion mechanism of the mobile ionic species across the superionic transition.
Ionic transport in conventional (non-superionic) solids is usually considered as being effected by a sequence of discrete "hops", where individual mobile ions undergo stochastic moves between available crystallographic sites [9,10]. While this model is usually appropriate for ionic solids with low-to-medium ionic conductivities, fast ionic conduction is often associated with highlyconcerted ionic motion, in which groups of ions undergo cooperative near-simultaneous motion [8,11,12,13,14,15]. This observational correlation between fast-ion conduction and cooperative ion-transport mechanisms suggests that ionic motion in superionic conductors might be better described in terms of appropriate collective degrees of freedom, rather than independent single-atom coordinates [16]. One appealing choice of basis for describing collective dynamics in a crystal is the set of normal modes of vibration, which correspond to phonon eigenvectors. Because phonons can be directly probed using established experimental methods, theoretical descriptions of superionic conduction in terms of phonon behaviour can, in principle, also be compared directly to experimental data for real materials.
Previous experimental studies have shown dramatic changes in phonon behaviour for particular solid electrolytes when heated above their superionic transition temperatures. In CuCrSe 2 , AgCrSe 2 , and γ-LiAlO 2 the superionic transition is accompanied by a breakdown of specific phonon modes [17,18,19], characterised by marked softening and extreme broadening of one or more phonon modes. Above the superionic phase-transition temperature, specific features of discrete phonon modes are no longer able to be accurately identified on spectra. This suggests that the associated vibrational degrees of freedom are lost, and that the phonon picture breaks down for these modes in the superionic regime. One proposed interpretation of this behaviour is that the corresponding vibrational degrees of freedom are transferred to translational degrees of freedom [17], with the phonon breakdown suggested as the physical mechanism driving the superionic transition from conventional hopping-mediated diffusion to liquid-like "stochastic" diffusion [17], or correlated ionic motion [18].
Phonon or vibrational density-of-states data have also been used to construct a range of cheap-to-compute descriptors for various properties of fast-ion conductors, including for ionic conductivity, activation energy, and electrochemical stability window [20,21]. Computationally cheap descriptors of specific properties are increasingly used in high-throughput materials discovery, where data from electronic-structure methods are combined with data-mining techniques to screen large numbers of materials to identify candidates with, hopefully, desirable properties [22,23,24,25]. The Debye frequency, obtained by speed-of-sound measurements, and the mobileion phonon band-centre, ω av , calculated from density functional theory (DFT) or from inelastic neutron scattering data, have previously been proposed as descriptors for ionic conductivity and activation energy, respectively, within individual solid electrolyte families [26,27,20]. The host-framework phonon band-centre, meanwhile, has been used as a descriptor for the electrochemical stability window of Li-ion conductors [20]. The utility of cheapto-compute descriptors within high-throughput computational workflows is illustrated by a study by Muy et al. in which screening against a set of descriptors that included the lithium-ion phonon band-centre was used to predict 18 promising lithium-ion conductors, with one, Li 3 ErCl 6 , then shown to display a high room-temperature ionic conductivity of 0.3 mS cm −1 [28]. The phonon band-centre for a given solid electrolyte is defined as the "centre-of-mass" when integrating over the full vibrational density-of-states (DOS), and can be computed for the full set of phonon modes, or for a subset using the atom-projected phonon DOS [20]. While this descriptor is much cheaper to compute than directly calculating ionic conductivities; for example, from molecular dynamics simulations; it neglects the effect of temperature on ion dynamics because the phonon DOS is computed within the harmonic approximation at 0 K [29]. The conventional phonon band-centre is therefore not expected to be able to fully describe a superionic transition-in particular for a type-II transition where the host-framework structure does not change. At nonzero temperatures phonon modes are not equally populated, but instead are populated according to Bose-Einstein statistics. A Bose-Einstein-weighted phonon band-centre, for example, therefore might better describe the effects of changing temperature on ionic conductivity than a conventional unweighted phonon band-centre, and may therefore be able to describe the changes in ionic conductivity that occur across a superionic transition.
Here, we consider the question of how the phonons of a superionic solid electrolyte change across the superionic transition, and what implications this has for various descriptors of ionic conductivity in the context of superionic materials. We consider the archetypal superionic conductor Li 3 N, and compare the phonon behaviour calculated at harmonic, quasi-harmonic, and anharmonic levels of theory.
Li 3 N is a structurally and compositionally-simple longstudied superionic conductor [30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45]. It has three phases: the αand the β-phase (Figure 1), which we consider here, and the γ-phase, which is only stable under high pressure at about 35 GPa [46], and which we therefore will not discuss further. The α-phase consists of alternating Li 2 N layers and pure Li layers. The α-phase conducts ions in its Li 2 N plane at a rate of 1.2 × 10 −3 S cm −1 at 300 K [35]. For several decades, it remained the highest Li ionic conductivity measured in a crystal around room temperature. Commercial applications, however, have been hindered by a low electrochemical stability, 0.45 V versus Li + /Li [47]. The β-phase has a similar structure but exhibits another pure Li layer in between each Li 2 N plane. It conducts ions in its pure Li layer at a rate of 2.085 × 10 −4 S cm −1 [48]. It is clear from the structures of the two phases that they are related by a first order structural phase transition.
Previous X-ray data for α-Li 3 N indicate that this phase exhibits a superionic transition at T s ≈ 678 K [41]. Below this temperature, the thermal mean-squareddisplacements of the Li ions [41] and the ionic conductivity [36] increase both linearly, suggesting Li 3 N exhibits a type-II superionic transition involving only the α-phase.
Based on analysis of X-ray diffraction data, Schulz and Thiemann have proposed that the superionic transition is driven by the increasing amplitude of strongly-directional anharmonic vibrations between Li sites in the Li 2 N plane, and that these same anharmonic vibrations promote diffusion in the superionic phase [41].
In this work, we first compute harmonic phonon spectra using finite difference methods for both α-Li 3 N and β-Li 3 N. Within the harmonic approximation we find that α-Li 3 N is the thermodynamically preferred phase above 300 K, and we therefore discount β-Li 3 N as playing a role in the superionic transition at T s ≈ 678 K. We then show that simple descriptors derived from this harmonic model, such as the phonon band-centre, are not able to describe the superionic transition, even when including a temperature-weighting. We then consider thermal expansion effects within the quasi-harmonic approximation (QHA). We again find that simple descriptors derived from the QHA are unable to describe the superionic transition of Li 3 N and conclude that thermal expansion is not sufficient to explain the superionic transition behaviour of Li 3 N. In contrast, within an anharmonic model, we observe a phonon breakdown above 678 K, which we interpret as a signature of the superionic transition, following previous experimental studies showing phonon breakdown at superionic transitions in CuCrSe 2 , AgCrSe 2 , and γ-LiAlO 2 [17,18,19]. We propose that the origin of this phonon breakdown within our computational model is due to a dramatically decreased residence time for the mobile lithium within individual sites, with this becoming comparable to their characteristic vibrational time period, which prevents the calculation of a well-defined harmonic restoring force for these ions.

Harmonic Lattice Dynamics
First-principles calculations using Density Functional Theory (DFT) were carried out to analyse the vibrational and ionic transport properties of Li 3 N. The crystal potential energy, Φ, can be expressed in terms of the displacement u jl of the j th atom in the l th unit cell. Φ can be expanded as where α, β, and γ are the Cartesian indices, and Φ 0 , Φ α , Φ αβ , and Φ αβγ are zeroth, first, second, and third-order force constants, respectively. In this study, the secondorder, or "harmonic", force constants were calculated using the finite-difference method implemented in the postprocessing software PHONOPY [49]. Within the harmonic approximation, the Helmholtz free energy of a crystal at constant volume can be expressed as where U el (V ) is the electronic internal energy, and F ph (V, T ) is the harmonic Helmholtz phonon free energy.
where the partition function Z can be expressed as .
(4) Here, U is the potential energy of the system and ω qλ are the frequencies of the vibrational modes λ at wave vector q. The force calculations were performed within 2 × 2 × 2 and 2 × 2 × 1 supercells for α-Li 3 N and for β-Li 3 N, respectively.
All of the DFT calculations were performed within the VASP code [50], using the PBEsol generalized gradient approximation exchange-correlation functional [51]. Projector augmented-wave pseudopotentials were used to model the ion cores [52,53], while treating the N 2s and 2p, and the Li 1s and 2s electrons as valence electrons. The electronic wave functions were expanded in a planewave basis with a kinetic-energy cutoff of 950 eV. We used 4 × 4 × 4 and 8 × 8 × 4 Monkhorst-Pack k-point meshes to sample the Brillouin zone for the αand the β-phases, respectively [54]. k-point mesh sizes were reduced for relevant supercell calculations to preserve net k-point density. Geometry optimisations were performed using a quasi-Newton algorithm with a convergence setting of 0.005 eV/Å for the atomic forces. The resulting optimised structures of α-Li 3 N and β-Li 3 N are shown in Figure 1. The calculated lattice parameters for these DFToptimised structures compare well with previous computational and experimental data (full details are given in Table S1 of the Supplementary Information).
The potential energy surface (PES) for a given vibrational mode at a particular wave vector can be obtained by calculating the internal energies of a set of structures along the corresponding normal-mode coordinate. We have used the program MODEMAP [55], which computes the displacement of atoms along a given eigenvector e qλj using N s is the number of atoms in the supercell used to model the displacement, m j are the atomic masses, and r j,l are the the atomic positions.

Phonon Band-Centre
The phonon band-centre, or average vibrational frequency, is a DOS-weighted average frequency defined as [20] This quantity can be calculated for specific species, such as the mobile ions, by using the projected DOS of that species instead of the total DOS. We also define a temperature-dependent Bose-Einstein-weighted phonon band-centre obtained by weighting each DOS(ω) term by the mean phonon occupation number, n(ω):

Quasi-Harmonic Approximation
Within the harmonic approximation, the vibrational DOS is temperature independent, and any temperature depen-dence arises only through considering phonon occupation numbers. The quasi-harmonic approximation (QHA) goes beyond this harmonic model, by considering how phonon frequencies vary in response to thermal changes in equilibrium volume. Within the QHA, phonon frequencies at a given volume are treated within the harmonic approximation. The equilibrium volume at a given temperature can then be found by minimising Equation 2 with respect to volume at each temperature. The mode Grüneisen parameter, γ qλ provides a measure of anharmonicity with respect to changes in volume, and can be obtained from the gradient ∂ω qλ /∂V via where V 0 is the zero-pressure equilibrium volume and ω qλ,0 are the corresponding zero-pressure vibrational frequencies. The heat capacity calculated within the QHA for α-Li 3 N is in good agreement with previous calorimetric measurements [40] (Figure 6(a)).

Ab initio Molecular Dynamics
While the quasi-harmonic approximation does include some temperature dependence, by allowing for thermal expansion, it does not give a complete description of the non-zero temperature phonon behaviour. In particular, phonon modes within the QHA are assumed to be noninteracting and have infinite phonon lifetimes. To go beyond the QHA and to consider a fully anharmonic phonon model, we have performed a series of "ab initio" molecular dynamics simulations (AIMD). We have then used the post-processing software DYNAPHOPY [56] to extract the phonon frequencies and phonon lifetimes as a function of temperature. DYNAPHOPY uses atom velocities, v jl (t), which can be directly derived from AIMD trajectories, to compute the velocity auto-correlation function Taking the Fourier transform of the velocity autocorrelation function gives the power spectrum, which is related to the phonon DOS by a factor of (3N n a k B T ) −1 , where N is the number of lattice points and n a is the number of atoms per lattice point. Anharmonic properties can be derived by fitting the power spectra, G qλ , of individual phonon modes to Lorentzian functions. From this fitting, ω qλ is determined as the peak position, while the linewidth Γ qλ is determined as the fullwidth-at-half-maximum of each peak. Each G qλ is obtained from Equation 10 and the velocity auto-correlation function of v qλ (t) is obtained from Equation 9. Phonon mode velocities v qλ (t) are derived from the normalmode-decomposition technique. Within this technique, v jl (t) are successively projected onto wave vectors q and phonon eigenvectors e qλ . DYNAPHOPY also allows the calculation of renormalised force constants from ω qλ values computed at commensurate q-points to interpolate ω qλ at incommensurate q-points and thus obtain the phonon dispersion.
To characterise the degree of lithium-ion diffusion in our AIMD simulations, we compute mean-squared displacements (MSDs) for the lithium ions as [57] as All MSD analysis was performed using the postprocessing software KINISI. We also compute the timeaveraged probability density of the lithium ions, P (r j ), by projecting the AIMD trajectory onto a uniform 3D grid and calculating the time-average number of atoms at each grid point [58]. Our DFT-AIMD simulations were performed in a canonical (N , V , T ) ensemble. A canonical ensemble was preferred to an isobaric-isothermic (N , p, T ) ensemble. Parinello-Rahman dynamics [59,60] requires complete convergence with respect to the plane wave basis set, which makes it more expensive than simulations performed within a canonical ensemble. Volume effects were also found to be marginal for the properties investigated ( Figure 6). Molecular dynamics simulations with variable cell volume are therefore not necessary. The temperature in the AIMD simulations is kept fluctuating around a given value by using the Nose-Hoover thermostat. A 4 × 4 × 4 supercell was chosen with the Brillouin zone sampled only at the Γ-point. For these AIMD simulations, the plane-wave energy cutoff was reduced to 300 eV to lower the computational cost. We performed simulations at 300 K, 400 K, 500 K, 600 K and 678 K. At each temperature, we performed 10 independent AIMD simulations with different random seeds. Each simulation was run for 8000 steps with a timestep δt = 2 fs. We exclude the first 300 steps of each simulation to allow equilibration, giving a total of 154 ps of independent trajectory data at each temperature. Mean-squared displacements were calculated over a time scale of 15.4 ps, averaged over each set of 10 independent runs to improve our statistics.

Polymorphism of Li 3 N
We first consider the α-Li 3 N and β-Li 3 N polymorphs within the harmonic approximation, to evaluate their relative thermodyamic stability as a function of temperature. Within the harmonic approximation, we calculate the free energy as a function of temperature of both αand β-Li 3 N using Equations 2 and 3. This analysis predicts that the βphase is more stable than the α-phase in the range 0 K to 300 K, with α-Li 3 N becoming the more stable polymorph above 300 K (Figure 2). This phase-transition temperature is broadly consistent with previous neutron powder diffraction data [61], which showed a steady increase in proportion of α-Li 3 N with respect to β-Li 3 N in a α/β mixed phase sample with increasing temperature above 473 K, with only a small fraction, 2.5 wt%, of the β-phase detected at 673 K. Because β-Li 3 N is predicted to be thermodynamically favoured only below 300 K, which is well below the superionic transition temperature T s ≈ 678 K, we conclude that β-Li 3 N is not involved in the superionic transition behaviour of Li 3 N. We therefore solely focus on α-Li 3 N in the remainder of this study. Figure 3 presents calculated Li-ion mean-squared displacements (MSDs) after a simulation time of t = 15.4 ps as a function of temperature and time-averaged Li-ion probability densities in a single ab Li 2 N plane at 600 K and 678 K. The MSD data show the onset of Li-ion diffusion (on our simulation timescale) at 600 K, which can also be seen in a log-log plot of the MSD data (Figure S1 (SI)). These MSD data also show a dramatic increase in lithium-ion displacements over the AIMD simulation timescale between 600 K and 678 K. At 678 K the average Li-ion displacement after 15.4 ps is larger than the Li-site nearest-neighbour distance. The time-averaged lithium-ion probability densities (Figure 3(b)) also show a distinct change in Li-ion behaviour between 600 K and at 678 K. At 600 K, the Li-ion density is localised around the formal lithium-ion sites, which is consistent with a conventional "hopping" ion-diffusion process. At 678 K the Li-ion density is more diffuse with clear "interstitial" contributions, indicative of significantly more mobile lithium ions. These data are consistent with the experimental observations of Schulz and Thiemann [41] who reported a decrease in Li-ion site occupation probabilities and a more diffuse Li-ion density from 600 K to 678 K, which was attributed to the onset of superionic conduction in α-Li 3 N.

Displacive Instability
Figure 4(a) shows the calculated harmonic phonon dispersion of α-Li 3 N. We find an imaginary mode at the Γ-point of the vibrational Brillouin zone that corresponds to out-of-plane vibrations of Li ions in the Li 2 N plane. The associated phonon eigenvector and PES for this mode are shown in Figures 4(b) and (c). The PES shows a double-well potential that is characteristic of displacive phase transitions [62]. The saddle point corresponds to the P 6/mmm space group of the initial crystal structure, while the two local minima have lower P3m1 symmetry. If a displacive phase transition were to happen around T s , this could be the mechanism behind the superionic transition. Out-of-plane motion of lithium ions has previously been identified as an essential part of the diffusion mechanism in α-Li 3 N [63]. Our data, however, show the energy barrier between the two wells is only 1.20 meV, which corresponds to a thermal energy of 14 K. At temperatures above 14 K, therefore, atomic vibrations will behave as for a single harmonic (quadratic) potential [62]. Even below 14 K, the P3m1 phase may not be stable due to zero-point motion and quantum tunnelling; such behaviour has been observed in the quantum paraelectric SrTiO 3 [64,65]. The effects of the double-well PES are therefore expected to be largely negligible within the stability range of α-Li 3 N (> 300 K cf. Figure 2), and to have little effect at the superionic transition temperature T s ≈ 678 K. We conclude from this analysis that the displacive phase transition predicted for α-Li 3 N does not contribute towards the mechanism of the superionic tran-sition observed in Li 3 N. 1 3.3.2 Temperature Dependence of the Phonon Band Centre Figure 5 shows the Li-ion-projected phonon bandcentre, calculated using the conventional temperature independent scheme (Equation 6, and with temperaturedependent phonon occupation weightings (Equation 7)). This analysis gives a temperature-independent ω av = 12.04 THz. Including the temperature-dependent phonon occupation number weightings, we obtain a lower, temperature-dependent value of ω av (298 K) = 8.08 THz, and that ω av (678 K) = 8.31 THz. The difference between the temperature-independent and temperaturedependent ω av values is due to the low occupancy of vibrational modes above 10 THz across the temperature range of interest (see Figure S3). The temperatureindependent phonon band centre calculation is equivalent to computing the temperature-dependent value at T = ∞ K, i.e. all phonons are equally occupied, regardless of their frequency. The behaviour observed here, that the temperature-dependent ω av (T ) < ω av is general; the effective phonon band-centre at 298 K will always be lower than the temperature-independent value. The degree to which ω av (T ) and ω av differ, however, is likely to vary for different materials, and this may cause changes in previously observed trends with respect to room temperature ionic conductivity measurements [27,20,28]. For α-Li 3 N above 300 K, the temperature-dependent phonon band-centre has a nearly constant value. Notably, we do not predict any change in behaviour at 678 K, as might be expected if this descriptor were describing something fundamental about the ionic conduction across the superionic transition. We further expect that using a temperature-dependent fully-anharmonic phonon DOS will not change the overall value of the phonon bandcentre across the superionic transition either [29]. This re-1 An imaginary mode was also observed in the phonon dispersion of β-Li3N at the K-point and is shown in Figure S2(a) (SI). In this case, however, the energy barrier between the two wells of the PES shown in Figure S2(b) (SI) is 40.95 meV. This corresponds to a thermal energy of 475.2 K, which is beyond the stability range of the β-phase according to the calculated phase diagram and neutron powder diffraction data [61]. It is therefore possible that β-Li3N takes its lower energylower symmetry structure with P 6 3 cm symmetry, instead of the accepted P 6 3 /mmc symmetry. This point is discussed further in the SI. sult is perhaps not unexpected because the phonon bandcentre is an average quantity, and it has been suggested that a superionic transition is associated with the breakdown of a single phonon mode [17,18,19]. Our finding that the phonon band-centre fails to capture any signature of the superionic transition in α-Li 3 N, even when temperature-dependent phonon occupations are incorporated, shows that this metric may not be reliable as a descriptor of ionic conductivities in solid electrolytes that exhibit a superionic phase transition. More generally, this result challenges the idea that the phonon band centre provides a useful conceptual framework for thinking about the microscopic dynamics of fast-ion conductors and how these relate to macroscopic ion-transport. Figure 6 presents the constant-pressure heat capacity (panel (a)) and the temperature dependent phonon dispersion of α-Li 3 N (panel (b)) calculated within the QHA from 0 K to 690 K, and the resulting properties are shown in Figure 6. The constant pressure heat capacity, C P , shown in Figure 6(a) is extrapolated from the free energies, while the thermal dependence of the phonon dispersion, shown in Figure 6(b) is calculated using the mode Grüneisen parameters defined in Equation 8. C P is linear around T s ≈ 678 K. There is no peak or discontinuity in the C P close to T s , which reflects the inability of the QHA to capture the α-Li 3 N superionic transition. The calculated values of the weighted average of the mode Grüneisen parameters range from −1.6 to 0.9, which results in mild softening of the phonon dispersion, as shown in Figure 6(b), but, again, gives no distinct mode softening at 678 K that might correspond to the superionic transition. In contrast, a QHA treatment of cubic CeO 2 does predict the emergence of imaginary modes at high temperature [66]. For the case of Li 3 N, thermal expansion effects alone cannot be the physical mechanism driving the superionic transition observed at 678 K. This conclusion stresses that Grüneisen parameters, a common measure of anharmonicity, are not a general descriptor for superionic transitions.

Anharmonic Phonon Breakdown
To go beyond the QHA lattice dynamics description, we consider a higher-order anharmonic model by mapping from the AIMD simulations onto a basis of phonon modes. We calculate velocity auto-correlation functions from phonon mode velocities, v qλ (t), using Equation 9, and power spectra, G qλ , using Equation 10. In contrast to the harmonic and quasi-anharmonic models discussed above, where we found no signature of the superionic transition, here a phonon breakdown occurs at 678 K (Figure 7). Above 600 K there is a drastic increase in the average phonon linewidth (inversely proportional to the mode lifetime), which is correlated with the large MSD of Li ions in the supercell. Such a discontinuous increase in average phonon linewidth cannot solely arise from threephonon scattering contributions. Only a smooth increase in linewidth would then be observed.
This behaviour is consistent with previous experimental studies in other solid electrolytes where the superionic transition has been found to be accompanied by an increase in phonon linewidth [17,18,19]. In contrast to these experimental studies, where the phonon breakdown at the superionic transition was associated with the breakdown of a particular mode, we are not able to assign a single-phonon breakdown at this stage. Instead, as shown in Figure S4 (SI), our model predicts a general phonon breakdown that involves almost all phonon modes at all q. Only acoustic phonon modes near the Γ-point remain coherent, which is expected since the centre of mass is not allowed to drift over the course of an AIMD simulation run with the VASP code. In short, the lattice dynamics description has broken down in α-Li 3 N at the superionic transition.
We interpret the results as a sign that the normal-modedecomposition-technique applied to moderately sized supercells fails in the diffusive regime. It represents an inability to project the velocity auto-correlation functions obtained from the high-temperature AIMD back onto the wave vectors, q, or the eigenvectors, e qλ that represent small harmonic displacements.
Once in the superionic regime, Li ions are not restricted to their regular lattice sites and diffuse, Figure 3(c). The mobile Li ions only vibrate for a short time at each individual lattice site before undergoing a diffusion event; this picture is supported by the large calculated average linewidth at 678 K. In a conventional "hopping" regime, where the average residence time for the mobile ions at individual sites is expected to be significantly greater than the average vibrational period [10], long-range vibrational order is preserved and the average motion of ions can still be described by a set of normal-modes. We also note that ions occupying interstitial positions may disrupt the longranged coherence of any lattice vibrations in a supercell model, which will further contribute to the reduction in mode lifetimes.
Nevertheless, in another study using the normal-modedecomposition-technique to probe the vibrational be-haviour of AgCrSe 2 across its superionic transition, Wang and Chen [67] replicated experimental results [18] and observed the breakdown of Ag-dominated transverse acoustic modes. Other phonon modes with little or no contribution to Ag remained coherent throughout the transition. From a computational point of view, the velocity autocorrelation functions could not be projected back onto the vibrational eigenvectors of the phonon modes where Ag ions were involved because they were diffusing and occupying interstitial positions [18]. On the other hand, where Ag-contributions were absent, the small harmonic displacements remained well-defined and the corresponding phonon modes stayed coherent throughout the superionic transition. Such a selective phonon breakdown is possible because the Ag-projected vibrational DOS is localised at low frequencies [18]. Both CuCrSe 2 and γ-LiAlO 2 also show localised site-projected vibrational DOS, with distinct regions where the mobile species dominates and others where it contributes little or not all to the total vibrational DOS [17,19].
In Li 3 N, however, the site-projected vibrational DOS for the diffusing Li ions is flat (Figure 8). Diffusing Li ions are involved in all phonon modes. In the superionic regime, we calculate a large broadening for all the phonon modes except acoustic modes near the Γ-point because they are all altered to some degree by the diffusion of Li ions. In addition, when Li ions in the Li 2 N plane are diffusing, the structure (Figure 1) may resemble a series of LiN strings. Their vibrations will be different than regular Li 3 N. On the other hand, materials like ternary triangular chalcogenides such as AgCrSe 2 and CuCrSe 2 maintain their well-defined layered structures even in the disordered superionic state [68].
These results suggest an explanation for why the (quasi)harmonic models fail to show any signature of the superionic transition in α-Li 3 N. A quadratic potential, even when it is allowed to soften, or harden, in response to temperature effects misses the physics involved. The superionic state can only be described by including changes in the potential energy surface sampled by the mobile Li ions and thus requires higher-order anharmonic terms from Equation 2 to be included.
Looking beyond lithium nitride, a consequence of this analysis is that descriptors which seek to capture superionic phase transitions must include some features of higher-order anharmonicity. It is doubtful that data from harmonic or QHA models alone will be sufficient for this purpose. Our results do demonstrate that phonon linewidths may be suitable, but that more effort is needed to understand the fundamental physics from a computational lattice dynamics perspective. For example, can the transition be understood in terms of a specific anharmonic term in the standard expansion of the potential energy (Equation 2)? Or can the transition be attributed to specific self-phonon or phonon-phonon interactions? Experimentally, our results also raise the question of whether a general phonon breakdown can be measured in systems where the mobile species is involved in most, or all, phonon modes. We anticipate that the general conclusions drawn here, and these associated open questions, are likely to be relevant for other superionic conductors.

Conclusions
We have computed the lattice dynamics properties of the archetypal superionic conductor Li 3 N using different levels of phonon theory. We show that the physics of superionic conductors is complex, even in a system as compositionally and structurally simple as Li 3 N. Our results confirm that β-Li 3 N is not involved in the superionic transition of Li 3 N, but that the superionic transition is specific to α-Li 3 N. We have demonstrated that both the harmonic and quasi-harmonic approximations fail to provide a signature of the superionic transition at 678 K. By using AIMD simulations to capture the high-order anharmonicity of the system, we identify a breakdown in the computed phonons in the superionic regime as quantified through the computed vibrational linewidths. This observation mirrors previous experimental observations of phonon breakdown at superionic transitions in other systems [17,18,19]. We suggest that the breakdown comes from a combination of short occupation times of the diffusing ions at their formal lattice sites and an increase in the occupation of interstitial positions, with both effects potentially contributing to a radical change in the way the diffusing ions are vibrating, preventing their ion dynamics from being directly mapped onto a set of normal-modes. In particular, diffusing Li ions are involved in all normalmodes, which gives rise to a general phonon breakdown instead of a mode-specific one. We anticipate these conclusions to hold for other superionic conductors, and high-light a number of questions that emerge from our findings regarding the nature of the relationship between anharmonic contributions to ion dynamics and superionic phase transitions.

Author Contributions
G.K. contributed to data curation, formal analysis, investigation, validation, visualisation, writing the original draft, as well as reviewing and editing the manuscript. C.E.K. contributed to formal analysis, as well as reviewing and editing the manuscript. K.T. contributed to validation, as well as reviewing and editing the present manuscript. B.J.M. contributed to supervision, visualisation, as well as reviewing and editing the manuscript. A.W. contributed to conceptualisation, funding acquisition, project administration, resources, supervision, as well as reviewing and editing the manuscript.

Conflict of Interest
There are no conflicts to declare.