The physical significance of imaginary phonon modes in crystals

Atomic motion at small amplitudes is usually well described by the harmonic approximation, where atoms move as if connected by harmonic springs. In this case the total energy can be Taylor expanded with respect to the atomic displacement as:

$$\begin{equation}U(\boldsymbol{u})\approx {U}_{0}+\sum\limits _{j,{j}^{\prime }}\sum\limits _{\mu ,\nu }{{\Phi}}_{j\hspace{0.5pt}{j}^{\prime }}^{\mu \nu }{u}_{j}^{\mu }{u}_{{j}^{\prime }}^{\nu }\end{equation} \tag{ 1 }$$

where u are displacements of the atoms from their equilibrium positions, U₀ is the energy of the static lattice, and Φ are the second-order force constant matrices given by:

$$\begin{equation}{{\Phi}}_{j\hspace{0.5pt}{j}^{\prime }}^{\mu \nu }=\frac{{\delta }^{2}U}{\delta {u}_{j}^{\mu }\delta {u}_{{j}^{\prime }}^{\nu }}.\end{equation} \tag{ 2 }$$

The indices j/j' label the atoms, and μ/ν label the Cartesian directions. It is assumed that the structure is relaxed to a stationary point on the potential-energy surface (PES) and the forces are equal to zero, so that the linear term in u vanishes.

The Φ in equation (2) are used to construct the dynamic matrix D (q) by performing a mass-reduced Fourier transform, which for solids yields:

$$\begin{equation}D{(\mathbf{q})}_{j{j}^{\prime }}^{\mu \nu }=\frac{1}{\sqrt{{m}_{j}{m}_{{j}^{\prime }}}}\sum\limits _{{l}^{\prime }}{{\Phi}}_{j0,{j}^{\prime }{l}^{\prime }}^{\mu \nu }\,\mathrm{exp}\left\{\mathrm{i}\mathbf{q}\cdot [{\boldsymbol{r}}_{{j}^{\prime }{l}^{\prime }}-{\boldsymbol{r}}_{j0}]\right\}\end{equation} \tag{ 3 }$$

where q is the phonon wavevector and we have introduced the indices l/l' to label the crystallographic unit cells of the atoms. Diagonalisation of D (q) is equivalent to solving an eigenvalue equation to determine the set of squared frequencies ω² and associated mass-weighted displacement vectors W (eigenvectors) for the 3n_a phonon modes at q:

$$\begin{equation}\boldsymbol{D}(\mathbf{q})\boldsymbol{W}(\mathbf{q})={\omega }^{2}(q)\boldsymbol{W}(\mathbf{q}).\end{equation} \tag{ 4 }$$

For a given mode, the harmonic energy U(Q) as a function of the normal-mode coordinate (displacement amplitude) Q is given by:

$$\begin{equation}U(Q)=\frac{1}{2}{\omega }^{2}{Q}^{2}\end{equation} \tag{ 5 }$$

Q is related to the Cartesian displacements of the atoms by:

$$\begin{equation}{X}_{j}(Q)=\frac{Q}{\sqrt{{n}_{a}{m}_{j}}}\,\mathrm{R}\mathrm{e}[\mathrm{exp}(\mathrm{i}\phi )\times {\boldsymbol{W}}_{j}\times \mathrm{exp}(\mathbf{q}\cdot {\boldsymbol{r}}_{j})]\end{equation} \tag{ 6 }$$

where ϕ is a phase factor.

There are two common approaches to performing a harmonic phonon calculation. The first is to compute the Φ in real space, typically using numerical differentiation by performing small finite displacements of the atoms. This is also termed the 'direct' approach in the literature, and is similar in principle to the 'frozen phonon' method where atoms are displaced collectively along known eigenvectors (e.g. determined from the crystal symmetry) to compute the corresponding frequencies. The second method is to compute the D (q) in reciprocal space using density-functional perturbation theory (DFPT, also termed 'linear response') [1, 14]. The real-space approach may require a supercell expansion to describe long-wavelength phonons, which can significantly increase the cost of phonon calculations when using methods that scale as a power of the system size. The approach can be implemented using first-principles methods such as DFT, i.e. computing the Φ using the Hellman–Feynman forces [15], but is ultimately agnostic to the underlying method used to calculate the forces. This makes the finite-differences approach compatible with the full range of DFT functionals and also with 'beyond DFT' theories such as dynamical mean field theory. On the other hand, DFPT can be more efficient, but its implementation in DFT codes is a technical challenge and so it tends to be more restricted in its applications, for example to specific families of DFT functionals [16–19].

The truncation of the expansion in equation (1) to second order in u defines the harmonic approximation. This leads to a straightforward method of obtaining the ω and W by diagonalisation of the dynamical matrices, but implicitly yields non-interacting phonon modes with infinite lifetimes and where atoms on average occupy their equilibrium positions regardless of the phonon occupation number (energy level). To understand processes such as thermal expansion and the creation and annihilation processes that give rise to finite lifetimes and thermal conductivity, one needs to consider the anharmonic interactions described by the third- and higher-order terms in u .

Under the harmonic approximation the equilibrium distances between atoms is independent of temperature, so this model cannot formally predict variations in volume with temperature and dependent properties such as thermal expansion. The quasi-harmonic approximation (QHA) extends the harmonic approximation to consider such volume-dependent thermal effects by assuming the harmonic approximation remains valid for a series of volumes about the athermal equilibrium and minimising the combined lattice and harmonic phonon energy as a function of volume and temperature. On a practical level the QHA is a straightforward extension to the harmonic approximation requiring a series of harmonic calculations at a set of compressions and expansions of the unit cell. The QHA provides access to a number of quantities including the Gibbs free energy, thermal-expansion coefficients, and Grüneisen parameters.

Beyond the QHA, higher-order anharmonicity can be incorporated into calculations in a number of ways. A common approach to determining the phonon lifetimes needed to calculate the lattice thermal conductivity is to apply a perturbation derived from the third-order force constants to the harmonic frequencies and eigenvectors [20–22]. Alternatively, a variety of approaches exist that aim to build effective harmonic potentials for a given temperature while incorporating the effects of higher-order terms in the Taylor expansion [7, 23–27]. Some of these approaches are discussed later in section 2.5.

As shown in equation (5), the energy change when atoms are displaced along a phonon mode is quadratic in the amplitude Q with a constant of proportionality ω². A positive ω² indicates a positive local curvature on the PES along the mode, and the energy therefore rises as the atoms are displaced from their equilibrium positions. On the other hand, a negative ω² indicates negative local curvature such that the energy decreases as the atoms are displaced. The harmonic frequencies ω of modes with ω² < 0 are complex numbers and, as such, are termed 'imaginary modes' [28].

Following this, a system for which all the harmonic modes are real is therefore a (local or global) minimum on the structural PES, whereas a systems with one or more imaginary modes is a transition state/saddle point or hilltop. As with real modes, imaginary modes have a well-defined eigenvector W that defines the collective atomic displacements (i.e. the structural distortion) that lowers the energy. By distorting the crystal structure along a particular phonon eigenvector (or group of eigenvectors), it is possible to map out the PES as a function of mode amplitude (figure 1). This 'mode-mapping' procedure is discussed further in section 2.4.

Zoom In Zoom Out Reset image size

For a large class of materials, the imaginary modes take the form of a double-well potential such as that shown in figure 6. This double-well PES is characteristic of materials that undergo the symmetry-lowering displacive phase transitions described by Landau theory [13]. As such, the form of the potential is fundamental to modelling ferroelectric and multiferroic materials. For example, the perovskite PbTiO₃ transitions from a high-temperature non-ferroelectric $Pm\bar{3}m$ cubic phase to a ferroelectric P4mm tetragonal phase at 763 K, driven by a single zone-centre phonon mode with an associated double-well potential [30, 31].

If the imaginary modes result from dynamic instabilities in the structure, further exploration of the PES offers a number of useful insights including predicting phase transitions [1, 6, 32, 33], discovering new crystal polymorphs [34], identifying dynamic Jahn–Teller instabilities [35], and understanding superionic conduction [36, 37]. We will discuss some of these applications as case studies in the second half of this review.

The imaginary modes associated with dynamic instabilities are fundamentally anharmonic in nature and a quantitative description requires moving beyond the harmonic approximation. Both the real-space finite-differences method and DFPT can also be used to compute the force constants for the higher-order terms in the Taylor expansion in equation (1), e.g. the third- and fourth-order force constants ${{\Phi}}_{j{j}^{\prime }{j}^{\prime \prime }}$/${{\Phi}}_{j{j}^{\prime }{j}^{\prime \prime }{j}^{\prime \prime \prime }}$ [20–22, 38–40]. Using the finite-displacement method to calculate these higher-order terms comes at a considerably higher computational cost than their harmonic equivalents, and perturbation theory is not valid for highly anharmonic materials where the anharmonic corrections to the harmonic energies are significant [41]. In these cases, alternative non-perturbative approaches are required to describe the anharmonicity. For example, it is possible to extract force constants and phonon frequencies from molecular-dynamics (MD) trajectories, or to fit force constants from stochastic atomic displacements [25, 26, 42–44]. Although typically more expensive than the perturbation theory-based alternatives, these methods have the advantage of implicitly taking into account higher-order anharmonicity and temperature effects. Non-perturbative methods for treating anharmonic modes are discussed further in sections 2.4 and 2.5.

Within the harmonic approximation, the force in response to an atomic displacement is given by minus the derivative of the harmonic energy with respect to the displacement, i.e. it is a force proportional to −ω². The force is restorative when it takes a negative value. An imaginary frequency corresponds to a non-restorative force leading to a decrease in potential energy as the atoms are displaced away from their equilibrium positions along the mode (figure 1). The presence of such an instability in a crystal structure can, understandably, lead to feelings of discomfort [45]. Being able to distinguish between imaginary modes that are an intrinsic property of the structural dynamics of a system, and those that reflect a technical problem in the calculations, is therefore key. As such, we consider three classes of possible origins for imaginary modes (table 1), viz (i) numerical issues in the calculations; (ii) material-specific physics that is not adequately accounted for when calculating the force constants or dynamical matrices; and (iii) dynamic instabilities that are an intrinsic property of the material.

Table 1. Summary of possible origins of imaginary harmonic modes in calculations, solutions for removing or treating them, and references for examples found in the literature. Each origin is classified as one of three types, viz: (I) technical (implementation) problems in the calculations; (P) material-specific physics that is not adequately accounted for in the calculations; and (D) a dynamic instability intrinsic to the system being studied. Note that the 'mode-mapping' technique proposed as a solution to the latter is discussed in detail in section 2.4.

Origin	Solution	Reference(s)
(I) Inadequately relaxed structure	Tighten force convergence criteria (e.g. to $\leqslant 0.01$ eV Å−1)
(I) Inadequately converged technical parameters	Converge forces with respect to the basis set, k-point sampling density and SCF convergence criteria
(I) Broken translational symmetry	Increase the size of the density grids and/or enforce the acoustic sum rule as a post-hoc correction	[46–49]
(I) Interpolation artefacts	For finite-difference calculations, select supercells commensurate with the wavevector where the imaginary mode is found (consider using non-diagonal supercells if needed). For DFPT calculations, ensure the q-point grid includes the wavevector	[34, 51, 52]
(P) Strongly correlated electrons	Use a suitable exchange–correlation functional for DFT calculations (e.g. DFT + U or a hybrid functional), or a suitable alternative theory (e.g. dynamical mean field theory) to calculate forces	[53, 54]
(P) Spin-phonon coupling	Include the effects of spin and spin–orbit coupling in calculations	[55, 56]
(D) Dynamic instability (Γ wavevector)	Use phonon mode-mapping to locate lower-energy structure along imaginary mode	[29–31, 34, 57]
(D) Dynamic instability (off-Γ wavevector, typically a zone-boundary point)	Ensure point is commensurate with supercell (finite-differences) or included in the q-point sampling grid DFPT, then use mode-mapping as for Γ-point imaginary modes	[1, 29, 34, 57]
(D) Dynamic instability in a defective structure	Use mode mapping or break crystal symmetry to allow for localised distortions during geometry optimisation	[58]

The first class encompasses a range of potential technical problems with calculations. Since the forces on the atoms in the initial structure are assumed to be zero (cf equation (1)), a well-converged geometry optimisation is a prerequisite, and the criteria for force convergence may need to be tightened beyond that typically considered sufficient for relaxations. Accurate forces are essential for computing the Φ needed in finite-displacement calculations, while determining accurate D (q) using DFPT requires tightly-converged ground-state electronic densities. Careful consideration must therefore be given to the electron orbital or plane-wave basis set size and the k-point sampling density used to model the electronic structure. For plane-wave basis sets it is usual to increase the cut-off energy by at least 25% above the default values used for simple electronic-structure calculations, and larger increases of up to 2× may be required in some systems. For numerical atom-centered basis sets, the total energy, density, and force cut-offs may need to be tightened by an order of magnitude.

Auxiliary settings, in particular the sizes of grids used to represent the electron density and to compute quantities such as the exchange–correlation potential, may also need to be considered. The vibrational Hamiltonian is invariant to a rigid translation of the crystal. As a consequence, the sum of the acoustic phonon frequencies at the Γ wavevector (q_Γ = (0, 0, 0)) should sum to zero, which is termed the acoustic-sum rule (ASR). However, in calculations where the electron density is represented on a spatial grid—commonly referred to as FFT grids in plane-wave codes, or integration grids in codes using numerical atom-centered orbitals—the translational invariance can be broken [46]. Deviation from the ASR results in imaginary modes immediately around the Γ wavevector. To remove such imaginary modes, one solution is to increase the size of the grids. A computationally cheaper alternative is to enforce the ASR as a post-hoc correction to the computed forces, although this can in some cases leave 'bowing' artefacts in the dispersion close to the Γ wavevector [47–49].

The choice of supercell in a finite-displacement calculation—or, equivalently, the q-point sampling density in a DFPT calculation—is also important. The sum over the unit-cell index l' in equation (3) means that the supercell used to calculate the force constants is large enough to capture all pairwise interactions for which the corresponding Φ are non-zero. For DFPT calculations, the q-point grid must be dense enough to capture the frequency dispersion of the 3n_a phonon modes in reciprocal space. The point at which these criteria are met is material specific and depends on the nature of the bonding. For example, silicon is a covalent system which requires a relatively long-ranged force-constant matrix (dense q-point grid) for convergence [50], whereas in sodium chloride the bonding is more ionic, the structure is less rigid, and the force constants decay at a faster rate, so a smaller supercell (sparser q-point grid) is sufficient.

A useful way to understand the importance of the supercell size is the concept of 'commensurate' wavevectors. A given supercell expansion will include the set of Φ needed to construct the exact D (q) and determine the ω/ W at q for which the cell contains an integer number of wavelengths (modulation periods). One can evaluate arbitrary dynamical matrices with an incomplete set of Φ using equation (3), but in this case the resulting frequencies and eigenvectors will not be exact. This procedure is known as 'Fourier interpolation'. For a DFPT calculation, a set of D (q) are evaluated exactly by design to obtain phonon frequencies and eigenvectors on a 'coarse' grid of wavevectors. The transformation in equation (3) can then be reversed to obtain a set of Φ to an appropriate real-space range, allowing the ω and W at arbitrary q to be obtained by interpolation in the same way. Therefore, if an imaginary mode is observed at a given wavevector one should check whether the wavevector is commensurate with the chosen supercell expansion (finite-displacement) or q-point grid DFPT—if this is not the case, then the mode may be an interpolation artefact [34, 51, 52]. We note that for finite-displacement calculations the zone-center q = Γ wavevector is always commensurate, so an imaginary mode here is either indicative of issues in the calculations, or indicates a distortion of the atom positions within the unit cell that would lower the energy.

As an example we consider the phonon dispersion of δ-Bi₂Sn₂O₇ (figure 2) [34]. With a 2 × 2 × 1 supercell expansion, an imaginary mode is obtained at the q_Z = (0, 0, 1/2) wavevector. The supercell is however not commensurate with this q point, and when a 2 × 2 × 2 expansion is used instead, which is commensurate with q_Z , the imaginary mode is removed, indicating it to be an artefact of the Fourier interpolation. It is of note that even when using the larger supercell a small imaginary mode remains along the Γ–Z branch of the dispersion—this is also likely an interpolation artefact, but to correct it would require calculations with an even larger suprecell that would be prohibitively expensive. We also note here that the other acoustic and optical modes of this system have converged, indicating that the overall physics of the system remains invariant of this artefact. In cases such as this, it can be useful to turn to non-diagonal supercells, rather than the diagonal (L × M × N) supercells commonly used in finite-differences calculations, which can reduce costs sometimes by orders of magnitude [59].

Zoom In Zoom Out Reset image size

If the theoretical model used in the calculations does not adequately describe the physics of the system, then the resulting force constants will not be correct. We have seen above how chemical bonding determines the supercell size (or q-point grid) required for accurate predictions. Other material-specific properties to consider when calculating the Φ and D (q) include the include the extent of electron correlation [53, 54] and spin-phonon coupling [55, 56]. We address curing of imaginary phonons emerging from these phenomena in table 1.

Finally, imaginary modes can arise when the crystal structure has intrinsic dynamic instabilities. Here the structure corresponds to either a local maximum or saddle point on the PES and may be accessible as an average over nearby minima at elevated temperature. The imaginary modes of SnS shown in figure 1 fall into this category, and it is this type of imaginary mode that forms the focus of this review.

If technical problems in the calculations are ruled out and imaginary modes ascribed to inherent dynamic instabilities in a material, a number of techniques can be used to explore the modes and obtain further physical insight. In this section we discuss the two most basic approaches: 'mapping' the structural PESs along the modes, and using the potentials to determine effective ('renormalised') real frequencies that can be used in harmonic calculations e.g. to determine free energies [6, 32].

To map the PES along a harmonic mode one can generate a sequence of structures with the atomic positions modulated by the mode eigenvector W , with different amplitudes Q using equation (6). The associated energy changes can then be calculated, using DFT or otherwise, to obtain a U(Q) curve. At small amplitudes the energy change is quadratic in Q (cf equation (5)), whereas at larger amplitudes the PES becomes increasingly anharmonic and higher-order polynomials are required to describe it (figure 3). For example, the symmetric double-well potential that is characteristic of many imaginary modes can often be described by the mathematical form U(Q) = aQ² + bQ⁴, as in Landau theory, (where Q is treated as an order parameter) although in principle higher powers of Q may be needed. While the U(Q) associated with harmonic modes in inorganic materials with symmetric bonding environments tend to be symmetric about Q = 0, this need not be the case, and indeed mapping bond-stretching modes in molecular solids, for example, usually yields asymmetric Morse-like potentials [60].

Zoom In Zoom Out Reset image size

For a single imaginary mode, mapping allows the nearest local minima along the mode to be located. In the harmonic approximation, the mode eigenvectors are orthogonal by construction, but this is only valid under the assumption of small displacement amplitudes. Since the minimum on the imaginary-mode PES often occurs at a Q outside the harmonic regime, if the goal of the mapping is to remove the instability it is usually prudent to perform a geometry optimisation on the minimum in order to relax orthogonal degrees of freedom.

One can also map a multidimensional PES U(Q₁, ..., Q_N ) by distorting along multiple harmonic modes simultaneously. This is technically required for degenerate modes because any linear combination of the N eigenvectors is a valid solution to the Schrödinger equation for the harmonic oscillator. In the case of a N-fold degenerate imaginary mode, the minima may be located at an arbitrary combination of displacements along the N eigenvectors W . One can also use this multidimensional mapping procedure to investigate the coupling between specific sets of harmonic modes.

We note in passing that the unit cell used to map the PES along a mode must be commensurate with its wavevector, and for multidimensional maps the unit cell must be commensurate with the q of all the modes. This may require that the mapping is performed in a supercell.

An issue with imaginary modes is that many physical properties derived using the harmonic approximation as a base, such as the phonon free energy, are not defined for imaginary modes, and so these modes are often neglected in such calculations. A convenient way to handle this is to assign the modes a temperature-dependent effective (real) harmonic frequency, a process termed 'renormalisation'.

Consider the typical double-well potential in figure 3. The structure at Q = 0 is a saddle point on the PES between two minima. In many systems, the saddle point corresponds to a symmetric structure observed at elevated temperature, whereas the two minima are equivalent realisations of a symmetry-lowering distortion that occurs on cooling through the phase-transition temperature. At T = 0 K, the system will be confined to the minima and will therefore adopt the low-symmetry structure. As the temperature is raised towards the transition temperature, the system will gain thermal energy and begin to explore the PES in the vicinity of the minima. Above the transition temperature, sufficient thermal energy will be available to traverse the barrier and the system will then be able to 'shuttle' between the equivalent minima. Averaged over time and space, it will appear as though the structure adopts the symmetric Q = 0 structure, but locally there will be large-amplitude distortions towards the lower-symmetry structure. Since the PES is anharmonic, and the region accessible to the system depends strongly on the available thermal energy, the frequency of the mode has a very strong temperature dependence, which can often be used experimentally to track the phase-transition behaviour. These measured frequencies are often relatively low, particularly close to the phase transition. This means the modes are heavily populated, and so it is generally undesirable to neglect them when calculating properties.

Given that we can obtain the potential by mode mapping, a simple scheme for modelling this temperature-dependent frequency is as follows. First, the Schrödinger equation for the PES is solved to give a set of eigenstates and corresponding energies for the nuclear positions in the anharmonic potential. This is shown for the R-point mode in the hybrid perovskite (CH₃NH₃)PbI₃ (MAPbI₃) in figure 4). The occupation of the eigenstates from a Bose–Einstein distribution at a target temperature T is then used calculate the contribution of the anharmonic mode to the thermodynamic partition function. Finally, a harmonic frequency $\tilde{\omega }$ is then assigned to the mode such that it reproduces this contribution when used to compute the phonon free energy. One can also construct a 'density of states' (DoS) showing the probability of the system occupying different Q as a function of temperature. For MAPbI₃, we see that, as expected, the distribution is heavily skewed toward the minima at low temperatures, but gradually becomes broader on warming and becomes close to homogeneous at a very high temperature of 3000 K. The $\tilde{\omega }$ can be used to correct the dynamic matrices D (q) at the q where the imaginary modes are located, which can then be reverse transformed to obtain corrected force constants Φ. This allows the renormalised frequencies to be used in calculations based on the harmonic frequencies/eigenvectors and force constants with minimal effort.

Zoom In Zoom Out Reset image size

This free energy-based renormalisation scheme was developed independently by Skelton et al and Adams and Passerone [6, 32], the latter of whom termed it the 'decoupled anharmonic mode approximation' (DAMA). Both formalisms provide a convenient means to compute temperature-dependent harmonic frequencies for imaginary modes that can be used to estimate phase-transition temperatures and pressures [6, 29]. The only substantial difference between them is that in the scheme proposed in reference [32] only the imaginary modes are renormalised, whereas in DAMA the renormalisation is applied to all modes (i.e. both real and imaginary). DAMA was found to provide a good estimate of the phase-transition temperature in cryolite (Na₃AlF₆), where a high-temperature orthorhombic Immm phases is linked to a lower-temperature monoclinic P2₁/n phase via an imaginary mode with a double-well potential [6]. The approach in reference [32] was also used to explore the impact of imaginary modes on other physical properties including the thermal conductivity of orthorhombic SnSe [32] and the electron–phonon coupling in MAPbI₃ [62]. Some notable studies that make use of mode mapping and renormalisation are explored in detail in the case studies in section 3.

These simple renormalisation approaches offer a number of advantages. Firstly, they are computationally relatively inexpensive—particularly in implementations where they are used only to treat a small number of imaginary modes—as they do not require significant numbers of further calculations. Secondly, since mapping the PES is agnostic to the code or theoretical method being used, the technique is general and easily adapted to e.g. different DFT functionals or post-DFT methods. Finally, the methods do not make any assumptions about the form of the PES and therefore do not require the specification of any model parameters.

On the other hand the main limitation of both approaches is that they retain the orthogonal harmonic eigenvectors, which may be a significant approximation in strongly-anharmonic materials. The approaches in references [6, 32] are defined for single modes, and as such neglect coupling between modes and may not work for degenerate modes where, as noted above, the multidimensional potential should be considered. Generalising these renormalisation schemes to multiple modes is in principle possible, but would significantly increase the cost and complexity of both the mapping and subsequent solution of the DFT. Finally, whereas the method for determining the effective harmonic frequencies is appropriate for thermodynamics, it may not be suitable for other purposes, for example for simulating the outcome of spectroscopic measurements that probe the difference in energy between occupied and virtual eigenstates.

A number of alternative and more sophisticated renormalisation schemes exist, such as the temperature-dependent effective potential (TDEP) [23–25] and the stochastic self-consistent harmonic approximation (SSCHA) [27] methods, which are the subject of the following section.

In systems with very strong intrinsic anharmonicity, or at high temperatures e.g. close to the melting point, the harmonic approximation may break down. Following the previous section, the consequence of this can be thought of as a substantial renormalisation of the harmonic phonon spectrum.

A useful reference point for approaching this is a (MD simulation. In a high-temperature MD simulation, the dynamics in principle incorporate anharmonicity to arbitrary order. A phonon frequency spectrum can be obtained from the Fourier transform of the velocity autocorrelation function and can be thought of as a renormalised harmonic phonon DoS. However, MD simulations, particularly ab initio MD (AIMD) simulations, often incur a significant computational cost overhead compared to harmonic phonon calculations. Furthermore, the detail on individual phonon modes that is readily available from the harmonic approximation is lost, which limits the insight available from the simulations, and also the scope for calculating properties that cannot be formulated as integrals over the DoS.

There have therefore been a number of approaches developed that aim to account for anharmonicity within an effective harmonic potential that can then be used in standard harmonic calculations. While such methods perform a renormalisation of the full phonon spectrum by design, and are thus not specifically targeted at treating imaginary modes, we review them here because intrinsic dynamic instabilities often lead to substantial deviations from the crystallographic average structure and can and do therefore impact upon other modes. Some full-spectrum renormalisation methods also yield higher-order (e.g. third- and fourth-order) force constants as a useful by-product. Full-spectrum methods include the self-consistent ab initio lattice dynamics (SCAILD) [7], TDEP [23–25], and SSCHA methods [27], the phonon quasiparticle approach [63], and the Anharmonic LAttice MODEl (ALAMODE) method [42].

The SCAILD method is based on a self-consistent approach to determining effective harmonic frequencies [7]. An initial estimate of the phonon frequencies is made using the harmonic approximation as outlined in section 2.1. Atomic displacements with random amplitudes are then performed according to equation (6) and the resulting forces computed. These forces are then used to determine new phonon frequencies, and the procedure iterated until the frequencies converge. SCAILD has been shown accurately to reproduce the phonon spectra of the group IVB elements and the group IIIB elements Sc and Y, and is expected to be applicable to the body-centered cubic (bcc) phases of the lanthanides and actinides, numerous ferroelectric materials, and some transition metals. However, as in the methods outlined in section 2.4 above, SCAILD uses the eigenvectors from the initial harmonic calculation, and, as a consequence, may not give accurate results for highly-anharmonic systems where the displacements deviate from this. However, SCAILD is less computationally intensive than alternative renormalisation methods that require MD simulations [7].

The TDEP method [23–25] uses MD simulations to sample the PES of the material at a finite temperature, and the atomic positions and calculated forces are then mapped onto a model Hamiltonian that describes the lattice dynamics and allows extraction of force constants [23–25, 64]. These can be used within the harmonic approximation to compute properties at finite temperature, and has been shown to yield accurate results for a number of systems including bcc Zr [25], PbTe [65], MgO [66] and SnSe [67]. The model Hamiltonian also allows for higher-order (e.g. third-order) force constants to be extracted. In a typical TDEP calculation, the unit cell volume would be equilibrated at the target renormalisation temperatures. In contrast to SCAILD, and also to the mode-mapping approach, the method therefore does not assume that the PES is independent of temperature and as such will provide more accurate predictions of thermal properties such as free energies. The main limitations of TDEP are that the use of MD—typically AIMD with a first-principles method such as DFT—comes at a significant computational cost, as noted above, and also that as a result of using Newtonian dynamics in MD calculations the method does not account for nuclear quantum effects, which can become significant at low temperatures [23–25, 64]. The issue of cost can be partially mitigated with the stochastic TDEP approach, where force-displacement data to fit the force constants is generated by stochastic displacements of the atoms drawn from a Boltzmann distribution corresponding to the temperature of interest, rather than by performing AIMD simulations [43].

In the SSCHA, a variational technique is used to minimise the free energy with respect to the density matrix of an auxiliary quadratic Hamiltonian, which allows it to account for anharmonicity at all orders and to incorporate quantum and thermal fluctuations [27]. The term 'stochastic' in the name comes from the use of Monte-Carlo to obtain averages of the minimum free energy and anharmonic stress tensor. The anharmonic force constants are evaluated by minimizing the residual sum of squares of the difference between the calculated forces and those from the anharmonic lattice model. In principle the SSCHA avoids the limitation of the SCAILD and (S)TDEP models and enables accurate prediction of phonon properties even for highly-anharmonic systems. For example, it has been shown to effectively predict second-order phase transitions and to accurately reproduce phase diagrams [68]. It is also a reliable method for structure searching, as it can locate saddle points in complex systems with large numbers of atoms [69, 70]. The scheme also has the capability to determine accurate phonon linewidths e.g. for thermal-conductivity calculations. The biggest drawback of this method is however its computational cost [27, 64].

It is worth noting that if applied to all the modes in a calculation, as in DAMA, the mapping and renormalisation approaches discussed in the previous section could also be used for full-spectrum renormalisation [6]. However, to do this the mapping and subsequent post-processing would invariably represent a substantial workload, at which point the cost may be competitive with some of the other, more sophisticated, methods discussed in this section.

As can be seen, there are a growing number of methods in the literature for going beyond the (quasi-)harmonic approximation by describing anharmonicity in lattice dynamics, including imaginary modes, but most of these come at a substantially higher computational cost over the simpler mode-mapping and renormalisation approaches described in the previous section. When choosing an appropriate treatment, there is generally therefore a trade-off between computational cost and accuracy, and the balance will be determined by both the system under study, the properties to be calculated, and the required level of accuracy. For systems with strongly-coupled and highly-anharmonic modes, techniques such as SCAILD, TDEP or the SSCHA will give accurate results but at an increased computational cost. On the other hand, in systems where the anharmonicity is largely restricted to a small number of imaginary modes, or for studies where cost is a primary concern (e.g. high-throughput modelling), the significantly cheaper mode-mapping approaches may suffice.

It may be surprising that the harmonic phonon dispersion of the high-temperature cubic $(Pm\bar{3}m)$ phase of the halide perovskite CsPbI₃ has several imaginary phonon modes (figure 5) when the structure is confirmed by x-ray diffraction (XRD) [71, 72]. This discrepancy however highlights the disconnect between experimental methods that probe a system at finite temperature and average over space and time, and the harmonic model of structures with an intrinsic dynamic instability.

Zoom In Zoom Out Reset image size

CsPbI₃ is in fact dynamically disordered at room temperature, with large anharmonic displacements of the atoms away from the ideal cubic sites as the crystal explores the anharmonic PES associated with the imaginary modes [72]. However, standard elastic-scattering techniques—including x-ray and neutron diffraction—probe the mean atomic positions and interatomic distances, averaged over the entire crystal and over the duration of data collection, and will thus identify the average positions albeit with large thermal displacements. On the other hand, the harmonic phonon calculation does not take account of the higher-order anharmonic terms that would stabilise the cubic structure and produce renormalised (real) frequencies at higher temperatures. These factors lead to a discrepancy between experiment and harmonic phonon calculations.

Mapping the PES along the imaginary modes in cubic CsPbI₃ demonstrates that this high-symmetry, high-temperature structure is not a local minimum at low temperatures, but a 'hilltop' along multiple double-well modes [57]. Displacements along the combination of ${M}_{3}^{+}$ and ${R}_{4}^{+}$ imaginary modes lead to a lower-symmetry and dynamically-stable orthorhombic phase. The calculated height of the energy barrier between this and the cubic structure is comparable to k_B T, indicating the cubic perovskite structure measured in XRD [71] is actually an average over the lower-symmetry minima.

Dynamic instabilities in the cubic perovskite structure are present across the inorganic halide perovskite series [35], and are also observed in their hybrid analogues such as methylammonium lead iodide ((CH₃NH₃)PbI₃, MAPbI₃). To probe the dynamics directly, inelastic-scattering techniques, such as inelastic neutron scattering [72] or high energy resolution inelastic x-ray scattering (HERIX) [33] can be used, in conjunction with atomistic modelling, to identify dynamic instabilities. Another alternative is to use methods sensitive to the local atomic structure such as extended x-ray absorption fine structure or x-ray absorption near edge structure [72, 73], again in conjunction with atomistic modelling.

Experimental measurements on MAPbI₃ identified a displacive phase transition driven by large-amplitude anharmonic rotations of the PbI₆ octahedra [33]. HERIX measurements yielded phonon dispersion curves that showed strong softening at the zone-edge points M and R, appointed to the in-phase (${M}_{{\mathrm{T}\mathrm{A}}_{1}}$) and out-of-phase (R_TA) octahedral tilting away from the high-symmetry cubic orientation. X-ray atomic pair distribution function measurements further demonstrated that the local structure was best described by the low-temperature orthorhombic phase whereas the longer-range structure was best described by the cubic phase. Both these outcomes were confirmed by theoretical calculations, where mapping of the imaginary ${M}_{{\mathrm{T}\mathrm{A}}_{1}}$ and R_TA modes produced two shallow double-well potentials leading to local metastable symmetry-broken minima. The relative well depths supported the experimentally-observed phase transition sequence of a cubic-to-tetragonal, transition resulting from condensation of the R_TA mode, followed by a tetragonal-to-orthorhombic transition from the condensation of the ${M}_{{\mathrm{T}\mathrm{A}}_{1}}$ mode [33].

In the context of applications of (hybrid) halide perovskites to opto-electronic devices, dynamic disorder has several repercussions at the device level. The connection between dynamic disorder and device performance is mediated by several physical phenomena including band gap broadening [62], spin-splitting [74], and carrier recombination [75]. For example, Niesner et al used resonant excitations of the photocurrent to provide direct experimental evidence for a dynamic Rashba effect (a spin splitting caused by dynamic disorder) in MAPbI₃ [74]. Munson et al used temperature-dependent and time-resolved infrared spectroscopy to provide evidence for increased dynamic disorder leading to large polaron formation in this system. The Rashba effect and large polaron formation are significant as they are thought to enhance the lifetimes and diffusion lengths of photocarriers in hybrid perovskite solar cells.

The mode-mapping technique can also be used to investigate how the imaginary phonon modes associated with dynamic disorder affect the electronic structure of a material. For example, mapping was used to calculate the change in band gap, ΔE_g, with respect to the imaginary-mode amplitude Q in MAPbI₃ [62]. Summing ΔE_g(Q) multiplied by the probability DoS along Q (figure 4) then yielded a thermally-averaged electron–phonon coupling for each imaginary mode and estimated positive band gap shifts of 35.5 and 27.9 meV for the R and M modes, respectively, at room temperature, which are comparable in magnitude to the photoluminescence broadening measured in experiments [76–78].

The tin monochalcogenides SnS and SnSe have received significant attention as emerging energy materials made from earth-abundant elements. SnS has long been studied as a potential thin-film absorber for photovoltaics applications due to its near-ideal bandgap E_g and large absorption coefficient α [79–81]. More recently, SnSe has been intensively studied for its potential as a thermoelectric material, as it has the ideal combination of physical properties needed for high heat-to-electricity conversion efficiencies, viz a high Seebeck coefficient S, low thermal conductivity κ, and high electrical conductivity σ [82–84]. SnS and SnSe have been proposed to adopt one of five different crystal structures, viz orthorhombic Pnma and Cmcm, cubic rocksalt (RS), zincblende (ZB) and 'π' phases. The Pnma structure is the energetic ground-state of both chalcogenides and transitions reversibly to the Cmcm phase at elevated temperature [85]. Epitaxial growth of RS SnS and SnSe on NaCl was first reported as far back as the 1960s [86, 87]. Much more recently, a nanoparticulate form of SnS was assigned to a ZB structure based on x-ray diffraction and the particle morphology, despite theoretical calculations indicating this to be unstable [88, 89]. This debate was eventually settled with the discovery of the π-cubic phase in the P2₁3 spacegroup [90–92]. It is well known that the Cmcm phase has intrinsic dynamic instabilities and, like the cubic perovskite structure, is an average over equivalent Pnma structures separated by a small energetic barrier [29, 32]. A comprehensive study of the five structures of SnS analysed the phonon dispersions and further determined that the π SnS is dynamically stable, RS SnS has a dynamic instability at the X wavevector, and ZB SnS is highly unstable with ∼50% of the modes being imaginary [93]. It was further found that the imaginary modes in RS SnS harden under compression and eventually become real, providing an explanation for how this phase is stabilised under epitaxial growth. More recent calculations found that both the π and RS phases of SnSe are dynamically stable at equilibrium, a contrast which can be explained by the stereochemical activity of the Sn lone pair [29].

The relationship between the low-temperature Pnma and high-temperature Cmcm phases of SnSe was established using the mode-mapping approach outlined in section 2.4 [32]. The phonon dispersion of Cmcm SnSe has imaginary modes at the Y and Γ wavevectors, each of which correspond to two double-well potentials (figure 6). Evaluation of the 2D PES spanned by both modes showed that the Cmcm phase is a hilltop on the PES. The global minima lie along the Y-point mode and correspond to the equivalent distorted Pnma phases, whereas the Γ-point mode leads to a further saddle point. Interestingly, this is despite the steeper curvature of the PES along this mode at the Cmcm phase. Again, as for cubic perovskites, at elevated temperature sufficient thermal energy is available for the system to access the Cmcm phase as a hilltop between Pnma minima, and the Cmcm structure is thus observed as a crystallographic average [32]. This finding is consistent with experimental inelastic neutron scattering measurements showing a strong softening of the low-energy optic modes in the Pnma phase towards the phase-transition temperature [94].

Zoom In Zoom Out Reset image size

A more recent theoretical study using the QHA [29] further examined the effect of pressure on the Pnma → Cmcm phase transition (figure 6). The barrier between the Pnma and Cmcm phase decreases under pressure and eventually disappears, resulting in the Cmcm phase becoming dynamically stable. This was predicted to occur at pressures around 11.5 and 8 GPa for SnS and SnSe, respectively, and the phase transition temperature, predicted using renormalised harmonic frequencies for the imaginary modes, was also found to decrease with pressure. The study also examined the impact of soft modes on the Helmholtz and Gibbs free energies, and it was shown that in this case the free energies calculated using renormalised frequencies for the soft modes were only very slightly different to those calculated with them being omitted from the partition functions.

Finally, the study in reference [29] also examined the structural relationship between the RS and π phases, which are based on a distorted 2 × 2 × 2 supercell of the eight-atom rocksalt conventional cell [91, 92]. Despite the imaginary mode in RS SnS, an energetic barrier was found between the RS and π phases of both monochalcogenides.

These two case studies—dynamic disorder in perovskites and phase transitions in the tin monochalcogenides—highlight the advantages of combining theory and experiment for materials characterisation. Returning briefly to CsPbI₃, as noted by Bertolotti et al the structural similarity of the various CsPbI₃ polymorphs results in very small differences in the XRD patterns, and they are almost indistinguishable by other widely-used characterisation techniques such as transmission electron microscopy [73]. In addition to this, single-crystal studies are challenging as the orthorhombic-to-cubic phase transition in CsPbI₃ is accompanied by a large change in volume, which is typical of materials with significant anharmonicity, leading to fracturing [71]. In these cases, imaginary modes in the harmonic phonon dispersions calculated for known experimental structures can provide important insight into the complex physics underlying phase transitions and dynamic disorder. The extension of this to crystal-structure prediction is explored in our third case study below.

The drive to find new materials with improved properties for specific target applications (e.g. efficient absorbers for solar cells) has led to an increasing interest in materials discovery [96]. The scale of this challenge is enormous: the chemical space of multi-component materials is vast, and only a fraction of potential ternary and higher compounds have thus far been made and characterised [97]. While an initial set of candidate materials can be identified using metrics based on the chemical composition alone, later-stage screening requires knowledge of the crystal structure [98]. Efficient methods for crystal-structure prediction are therefore an important challenge in contemporary materials modelling.

Given a composition and possibly a target unit-cell size, the goal of a structure-prediction algorithm is to traverse the structural PES and locate the global energy minimum, a very general problem to which a number of potential approaches have been examined [99, 100]. The simplest conceptually is the so-called ab initio random structure searching method [101], which aims at an unbiased sampling of the PES by generating random structures within some universal chemical constraints such as reasonable minimum bond lengths. Building on this are a variety of 'accelerated sampling' methods that attempt to accelerate the search over the PES, with three popular approaches being basin hopping, particle-swarm optimisation and genetic/evolutionary algorithms [102–104]. More recently, lower-cost alternatives have been developed based on data-mining of known crystal structures, such as using atomic substitution probabilities to identify known structures that could be adopted by a target composition [105] and building new structures from 'modules' of existing ones [106].

A relatively common type of structure-prediction sub-problem is that where a high-temperature/high-symmetry structure is known from experiments (or can be guessed) but shows dynamic instabilities, and we wish to find the corresponding stable minima. Returning again to the CsPbI₃ example discussed above, we might assume that an ABX₃ compound would adopt an analogous, dynamically-unstable cubic structure, and then attempt to find the low-temperature/low-symmetry phases. A powerful approach to this problem is the following. We compute the phonon spectrum of the starting structure to identify the imaginary modes at the high-symmetry wavevectors (i.e. the zone-centre and/or the zone-boundary points). We then map each of the modes sequentially to locate the nearest local minima, perform full geometry optimisations to relax the orthogonal degrees of freedom, and hence obtain a second 'generation' of structures (one per imaginary mode). This procedure is then applied to each of the new structures, and iterated until each branch terminates in a dynamically-stable phase with no imaginary modes.

There are a handful of examples where this approach has been applied successfully. As a proof-of-concept, Togo and Tanaka studied the elemental metals Na, Mg, Al, Ti, Cu, Zr and Hf [95]. In each case the search was started from a sc parent structure 0.35–1 eV atom⁻¹ higher in energy than the known stable phases and with multiple phonon instabilities.

For Cu, the instabilities in the SC phase led to stable fcc/hcp minima via a set of lower-energy unstable maxima including simple hexagonal (SH), body-centered tetragonal (BCT), body-centred cubic (bcc) and ω phases (figure 7). The resulting 'tree-like diagram' also enables phase-transition pathways between the stable phases to be identified. The ω structure is a 'junction' between the ω → bcc → hcp and ω → fcc pathways and thus serves as a previously unreported low-energy barrier in the fcc ← → bcc transition pathway. The bcc ← → hcp transition was previously thought to be complex [107], but was found from the PES mapping to proceed via an intermediate with a Cmcm space group.

Zoom In Zoom Out Reset image size

A similar analysis of Mg found hcp, fcc and two 9R and 18R long-period stacking (LPS) structures as minima together with a metastable $P6\bar{2}m$ structure that has yet to be reported experimentally but is formed by Mg-rich alloys. Comparison of  Ti and Hf reveals an interesting contrast wherein the low-energy ω structure is linked to a SH phase in Hf but not in Ti, implying the presence of an additional energy barrier in the latter that explains why this phase is not accessible under ambient pressure, and also why the transition under pressure occurs with a large hysteresis [108]. Searches performed under applied pressure found that pressure stabilises the bcc phase in both Cu and Mg and introduces several LPS minima in Cu, and therefore has a significant effect on the PES.

Across the seven elements examined the fcc and hcp structures were identified as common stable phases together with five LPS structures, but, notably, the ω phases of  Ti and Zr was found to be more stable than the fcc phases despite not being a close-packed structure. This highlights the level of insight that can potentially be obtained from this structure-prediction approach, even for comparatively simple and well-studied systems.

Rahim et al applied the same approach to the ternary oxide Bi₂Sn₂O₇ in conjunction with an exhaustive analysis of variable-temperature neutron-diffraction measurements [34]. Experiments indicate that a cubic pyrochlore structure (γ-Bi₂Sn₂O₇; figure 8(a)) is formed at high temperature, which on cooling converts to an intermediate β phase and a low-temperature α phase with two proposed structures (here denoted α and α_new). The dispersion of γ-Bi₂Sn₂O₇ has imaginary modes at four different wavevectors (figure 8(b)). The PES map obtained from the structural search is reproduced in figure 8(c). The structure search recovers the β phase and the α_new structure, providing strong evidence that the latter is the correct model. Imaginary modes were also found in the β dispersion and followed to a new ζ phase, which proved to be a better match to the experimental diffraction data. The mapping also identified two new metastable phases, termed δ and , that have yet to be identified experimentally. The parent γ phase has 88 atoms in the unit cell, while the intermediate and stable phases have between 44 and 176 atoms, which would make finding these structures with conventional structure-prediction algorithms extremely challenging.

Zoom In Zoom Out Reset image size

The low-temperature α_new structure was subsequently predicted to be a high-performance thermoelectric material [109], due to the unusually-low lattice thermal conductivity arising from the complex structure and heavy metals, which further highlights the potential of this approach for materials discovery.

Kayastha et al used a similar mode-mapping approach as part of a high-throughput workflow to predict dynamically stable quasi one-dimensional (Q1D) polymers [110]. The study considered hypothetical materials produced by bonding one of eleven monovalent metals to one of 109 substituted cyclopentadienyl anions in a vertical (stacked) packing arrangement. The phonon dispersion curves were calculated for the resulting 1199 Q1D chain polymers along the Γ–X direction.

Based on the form of the dispersion the materials were divided into three classes, viz: (i) dynamically stable systems with no imaginary harmonic modes; (ii) materials with a single or degenerate instability; and (iii) materials with multiple independent phonon instabilities. Approximately half of the materials (587) were identified as having maximum instabilities at the Brillouin-zone boundary with $q=\frac{1}{2}$, of which 33 were selected for further analysis to examine trends in the effect of the distortions on the electronic structure. A further 28 materials were predicted to have charge-density wave like distortions with a instabilities at $q=\frac{1}{3}$. Figure 9 shows an example where the maximum instability occurs at $q=\frac{1}{3}$. To remove the imaginary mode and obtain a dynamically-stable structure, the unit cell was expanded threefold in the stacking direction and optimised with the atomic positions modulated by the imaginary mode wavevector and geometry relaxed. Subsequent phonon calculations confirmed the resulting structure to be dynamically stable.

Zoom In Zoom Out Reset image size

To conclude our discussion of using imaginary phonon modes for crystal structure prediction it is interesting to consider how the mode-following technique compares to other structure-prediction methods. In general, complex materials are challenging both because structure prediction requires searching over more independent variables, and because the number of minima on the PES grows exponentially with the system size [111]. In the approach adopted in these studies, it can be seen that the search over 3n_a atomic positions is reduced to a series of sequential searches over a much smaller number of imaginary modes. An additional subtlety is that in the mode-mapping approach any changes required to the cell volume are indicated by the wavevectors of the imaginary modes, as was exploited in the example of the Q1D polymer discussed above [110]. This is in contrast to many prediction algorithms, where the number of formula units is treated as a fixed variable or restricted to a range of values. As highlighted in the study on metallic elements, the imaginary modes linking structures in each generation can also provide valuable information on the connectivity of the PES (figure 7) [95]. Finally, the successive freezing in ('condensation') of imaginary modes mimics how a material naturally explores its PES on cooling, so if the mode-following procedure is started from a high-temperature experimental structure it is not unreasonable to expect that the matching low-temperature structure(s) will be found, as in the case of Bi₂Sn₂O₇ [34].

On the other hand, the requirement for an initial structure means the mode-mapping approach cannot be used to predict the structures of compositions for which there are no experimental structures or suitable known analogues. Also, since the stable structures found using this technique must be accessible from the starting structure by following imaginary modes downhill in energy, if the global minimum and/or other local minima are not connected to the starting structure in this way they will not be found. This is the reason why for example the study in reference [95] does not locate the stable ω phase of Ti. If the starting structure is taken from experiment, i.e. it is the phase obtained from the synthesis, this may not be an issue. However, if the initial structure is guessed then an appropriate exploration of the PES is not guaranteed. Finally, this technique requires a full phonon calculation on each structure, which can easily impose an overhead of 1–2 orders of magnitude on the geometry-optimisation steps.

If there are other systems like Bi₂Sn₂O₇, where a high-temperature structure has been solved but the low-temperature phase(s) are unknown, exploring these with mode-mapping has the potential to yield new materials that should be synthetically accessible. Alternatively, there may be some scope for using this approach in combination with other structure-prediction algorithms to refine dynamically-unstable structures obtained from a coarse initial sampling of the structural PES.