Atomic motion in solids with dimpled potentials

Polymorphic solids of the same chemical composition can have different atomic structures; in each polymorph atoms vibrate around a local potential energy minimum (LPEM). If transformations to other structures have sufficiently high enthalpy barriers, then each polymorph is either stable or metastable; it is stationary and does not spontaneously change with time. But what happens, if those barriers are low? As examples, we consider NiTi shape memory alloy exhibiting a large elastocaloric effect, and selected elemental solids. We suggest a model for dynamically polymorphic solids, where multiple LPEMs are visited by ergodic motion of a single atom. We predict that upon cooling a dynamically polymorphic phase should undergo a symmetry-breaking first-order phase transition, accompanied by a finite change of the lattice entropy. We discuss 3 methods used to calculate phonons in solids with non-harmonic dimpled atomic potentials, and compare theoretical predictions to experiment.


I. INTRODUCTION
A non-harmonic atomic potential presents a challenge for those who use a harmonic or quasiharmonic approximation for addressing solids. A dimpled potential is not harmonic, can cause a lattice instability, is not straightforward to deal with, but is very common in practical materials [1] with more than one local potential energy minimum (LPEM), covered by ergodic atomic motion. There are well-developed methods for a harmonic potential, where parabola has a single minimum. However, a dimpled potential dramatically differs from a harmonic one. Quite often, the high-symmetry atomic position, which used to be an energy minimum in a harmonic potential, happen to be a local energy maximum or a saddle in a dimpled potential, with several LPEMs surrounding this position. How to calculate phonons and describe atomic motion in such non-harmonic potentials? We will consider 3 algorithms, obtain ab initio results, and compare them to experiment.
Crystal is a solid, in which atoms are arranged in a definite pattern and whose surface regularity reflects its internal symmetry [2]. A crystal can be described as a Bravais lattice -an infinite periodic array of discrete points [3]. The crystalline periodic arrangement of atoms is manifested by an x-ray [4] and neutron diffraction crystallography. According to Laue [5], the Bragg spots [4] are observed due to constructive interference, which occurs if the change in the wave vector belongs to the reciprocal lattice. Nevertheless, in some crystals the ideal atomic positions on a Bravais lattice are unstable [1].
Comparing an experimentally observed diffraction pattern with the one predicted for a particular Bravais lattice, crystallographers suggest a crystal structure [4]. However, a correspondence between crystal structures and diffraction patterns is many-to-one: more than one crystal structure can produce the same pattern, while each particular fully ordered crystal produces a unique (one and only one) diffraction pattern. To add confusion, a partially disordered crystal can produce a similar pattern. In addition to thermal atomic motion, materials can have athermal atomic disorder, which can be chemical or displacive, local (due to lattice defects) or nonlocal. An ordered crystal and a solid with a substantial displacive disorder can produce similar diffraction patterns with the same positions but different broadening of peaks. One example is NiTi austenite with the assumed unstable B2 structure and multiple stable representative structures [6,7]. Lattice instabilities were also found in antiferromagnetic (AFM) phase of B2 FeRh [8] and in body-centered cubic (bcc) phases of Ti, Zr, Hf, and Li [1,9,10]. A crystal-like diffraction pattern can be produced by solids, which are not periodic (for example, quasicrystals) [11].
Atoms in solids are trapped in deep potential energy (PE) basins ( Fig. 1), separated by PE barriers E B , which are high compared to kT , where T is temperature and k is the Boltzmann constant. We define N L to be the number of the local potential energy minima per PE basin. We refer to a crystal with N L = 1 as conventional; a 1-dimensional (1D) example is in Fig. 1(a).
A solid with a dimpled potential can have multiple LPEMs and hence N L > 1 [e.g., N L = 2 in Fig. 1(b)]. In general, potential energy is a functional of atomic positions, and a path from one LPEM to another can be a collective atomic motion.
Each LPEM is a stable or metastable arrangement of atoms. Multiple stable chemical structures are known as isomers in molecules and polymorphs in solids. If the enthalpy barriers between polymorphs are high compared to kT , then each polymorph is stationary and does not spontaneously transform to other polymorphs. However, if the barriers between LPEMs are low compared to kT , then thermal atomic motion covers several LPEMs; this state of matter can be called "dynamically polymorphic", and this phenomenon -a "dynamic polymorphism." (A similar semantics is used in computer programming to refer to the runtime polymorphism).
Below we discuss atomic motion in a dimpled potential (section II) and compare phonon methods (section III). As a model, we consider a periodic potential with deep PE basins; only one atom occupies each basin (we assume that a strong interatomic repuslion at short distances makes presence of another atom in the same basin energetically unfavorable). In our examples (section IV) we consider one [Figs. 1(a) and 2] or more [from two in Fig. 1(b) to 48 in Fig. 7(c,d)] LPEMs per basin. Properties of the simplified models and real materials are discussed in sections II, IV, V and summarized in section VI. Computational details are provided in Appendix A.

A. Model
Let us consider a solid with several (N L > 1) LPEMs per PE basin. In a single PE basin, a set of LPEMs connected by minimal enthalpy paths (MEPs) forms a network, which might include (Figs. 1(b) and 3) or exclude (Figs. 4, 6 and 7) the high-symmetry crystallographic position at the basin center (x = 0, with energy E L above LPEM). Let the enthalpy barriers along MEPs be E l (relative to LPEM); for simplicity we assume that enthalpies of all LPEMs are the same and the barriers have the same height (one can index individual enthalpies and generalize our consideration to less symmetric cases). At a high enough temperature (kT ≥ E l ), an atom can move from one LPEM to another within the same basin.
We assume that the enthalpy barrier E l is not higher than the PE E L at the high-symmetry point x = 0, and both enthalpies are low compared to the barrier E B be-tween the PE basins, see Fig. 1 In our examples, E l = E L = 0 for the harmonic potential with N L = 1 in Figs. 1(a) and 2; E l ≡ E L > 0 for the double-well potential with N L = 2 in Figs. 1(b) and 3; 0 = E l < E L for the muffin-tin sombrero potential with N L = ∞ in Figs. 4 and 5; 0 < E l E L for the 2D and 3D corrugated sombrero potentials with a finite countable N L > 1 in Figs. 6 and 7.
In general, one of the LPEMs can be at x = 0 (for example, the only LPEM in a harmonic potential in Fig. 2).
Here we focus on a less trivial atomic motion with an instability at x = 0; the MEP between adjacent stable LPEMs can either include or bypass the high-symmetry point x = 0. The double-well basin in Fig. 3 is a 1D example of a potential with N L = 2 [ Fig. 1 Assuming that each LPEM is a point (and not a line, like in the muffin-tin sombrero potential in Fig. 4), one can define a harmonic limit E h < E l around it, such that atomic vibration around a single LPEM is harmonic at sufficiently small displacements and potential is approximately parabolic at E ≤ E h .

B. Thermal atomic motion
Atomic motion in a solid with N L > 1 depends on T and can cover vicinity of one or several LPEMs.
• kT ≤ E h E l : harmonic vibration around a single local potential energy minimum. The small atomic displacement method can be used to calculate phonons at a LPEM, see Section III.
• E h < kT < E l : anharmonic vibration around a single LPEM.
• kT ∼ E l : a phase transition happens.
• E l ≤ kT < E L : atomic motion covers several LPEMs in the same PE basin. If such LPEMs are distributed symmetrically around x = 0, then the time-averaged atomic position is zero: x = 0.
• E L ≤ kT < E B : atomic motion covers a significant part of the PE basin, including the center (x = 0) and multiple LPEMs. If LPEMs can be interpreted as a negligible roughness (E L kT ) at the bottom of a nearly harmonic potential, then a finite atomic displacement method can be used to calculate phonons around x = 0, see Section III. Snapshots of thermal atomic motion in molecular dynamics at fixed T are finite collective atomic displacements in a solid phase.
• kT ≥ E B : atomic motion is no longer restricted by a PE basin; the solid has melted or sublimated.
A dynamically polymorphic solid phase exists at temperatures E l ≤ kT < E B . Upon cooling, it transforms to a lower-symmetry phase below We expect this transformation to be of the first order, because it is accompanied by a discontinuous change of the lattice entropy ∆S L . Under certain conditions, a change of the total entropy is responsible for an isentropic temperature change (caloric effect). In alloys like NiTi, this effect is quite large.

C. Representative structures
At kT E l , atomic arrangement is stationary: each atom vibrates around a single LPEM. If at each PE basin an atom randomly chooses one of the LPEMs [e.g., either right or left LPEM in each double-well basin in Fig. 1(b)], then this atomic structure is aperiodic, but its diffraction pattern coincides with that produced by a crystal with partially occupied LPEMs in periodic PE basins. Even if there is no atomic periodicity, it is possible to consider a "large enough" representative periodic unit cell, which correctly represents energy, phonon spectrum, average LPEM occupations, and other physical properties.
Let us consider a 1D example in Fig. 1(b). A periodic unit cell with even number of atoms, where each (right or left) LPEM is occupied by half of the atoms, provides the same energy per atom (at T = 0), the same occupancy (c i = 1/2) of each (right or left) LPEM, and a vibrational spectrum similar to that of the whole aperiodic solid. Thus, for the purpose of computing its physical properties, a solid can be approximated by a periodic atomic configuration with a representative structure. Any particular representative unit cell is not unique, there are many others. All representative structures have similar properties, which approximate those of the solid.

D. Displacive pattern
Displacements of atoms from the high-symmetry positions (PE basin centers at x = 0) create a pattern, which can be either stationary (at kT < E l ) or dynamic (at E l ≤ kT < E B ), and has a characteristic length [12], imported from cosmology [13] to condensed matter [14].
At kT ≥ E l , motion of each atom covers several LPEMs and the displacive pattern is dynamic. A broadened distribution of the interatomic distances differs from that in a conventional crystal (e.g., NiTi austenite [6] has such broadening due to athermal atomic displacements). Although there is no periodicity of the instantly occupied LPEMs, a diffraction pattern produced by such a solid reminds that of a conventional crystal (with peaks at the same positions, but not of the same width).

E. Effect of interatomic interactions
Interactions between atoms can change relative energies and positions of LPEMs. A shift of energies can force atoms to choose one particular LPEM at each basin, thus forming a fully ordered crystal (an example is the γ-Se, hP3, A8, P 3 1 21 structure of Te and Se-Te alloys, which consists of 3 triangular lattices stacked along z [15], and can be considered as a distortion of a simple cubic lattice [16]). A shift of the equilibrium atomic positions affects a distribution of the interatomic distances and a diffraction pattern. Nevertheless, atomic displacements and interatomic distances in a dynamically polymorphic phase will differ from thermal ones in a harmonic crystal.

A. Small atomic displacement method
The quasi-harmonic approximation (QHA) is often used to calculate phonons in conventional crystals [17]. Atoms are assumed to be at stable equilibrium at 0K. Small atomic displacements u from this stable equilibrium result in the increase in the potential energy E, which is quantified in QHA using a Taylor expansion: Here D αβ ij is the force-constant matrix: The latin indices i, j enumerate atoms, while the greek letters α, β label directions. The higher-order terms beyond the second order can be neglected, if the amplitudes of all atomic displacements are small: |u| ≤ u h , where u h is the harmonic limit. The instant atomic forces for the near-equilibrium atomic configuration τ n are Given a sufficient number N of independent atomic configurations τ n (n = 1...N ) with known atomic displacements u j and forces F i , one can solve a system of linear equations (5) and find matrix D ij , which can be used to find the phonon spectrum and density of states (DOS). The minimal number N of independent atomic configurations [and linearly independent equations (5)] is equal to the number N D of independent components of D ij . The system (5) might be over-determined, if N > N D . For either well-determined or over-determined system (5), the effective force-constant matrix D e ij can be found by minimizing the sum of the differences between the actual and predicted forces [18]: Within the small atomic displacement method, expansion is around a LPEM, and displacement is assumed to be within a harmonic limit (e.g., |u| ≤ x H in Fig. 2). Any displacement |u| = |x − 0| ≤ x H (including infinitesimal) provides the same vibrational frequency for the harmonic potential in Fig. 2 To emphasize that harmonicity is just an approximation, we exaggerated the harmonic region in Fig. 2, where the force F (x) = Kx is precisely linear and the energy The result of the QHA is correct, if each displacement u = x − x L from a stable equilibrium at x L is indeed within the harmonic limit at |u| < u h , and the potential To get a correct result, one must know the stable equilibrium coordinates x L at each occupied LPEM (Figs. 1(b) and 3).

B. Finite atomic displacement method
Alternatively, the harmonic potential in Fig. 2 can be comprehended as a limit of the double-well potential in Fig. 3 with x L → 0, so that x 0 = x L = 0 and E(x L ) = E(0) = 0. If multiple LPEMs can be interpreted as a negligibly small roughness at the bottom of a nearly harmonic potential (approximately parabolic at −x H < x < x H , with x L → 0, i.e., x L x H ), then one can approximately calculate phonons using the finite displacement method (example is bcc Li, section IV F), which avoids unstable phonons, if F (x)/x > 0; this happens when atomic displacements u = x − 0 are larger than the distance between unstable (at zero) and stable (at x L ) atomic positions, see Fig. 3].
In the finite atomic displacement method, expansion (3) can be around the high-symmetry crystallographic position at x = 0, which might (Fig. 2) or might not (Fig. 3) be a LPEM. If the finite displacements u ≡ x − 0 are sufficiently large, i.e., x L < |u| ≤ x H and none of them is in the region −x L < x < x L , where a destabilizing force pushes an atom away from unstable equilibrium at x = 0 (see Fig. 3), then there are no imaginary frequencies in the calculated phonon spectrum. Because this method deliberately avoids the region −x L < x < x L , unstable phonons are overlooked [19].  Fig. 1(a)], harmonic at E < EH . 3. (color online) Potential energy E [same as in Fig.  1(b)] and force F = dE/dx versus displacement x for a potential with several (two) LPEMs per potential energy basin. Highlighted are forces with F/x < 0 at −xL < x < xL, which push an atom away from the unstable equilibrium at x = 0. At displacement x0 (filled circle), effective harmonic linear force and parabolic potential are shown (dashed orange lines). finite displacements. MD sampling gives a set of atomic positions u j (t n ) and forces F i (t n ) for a large number N of time steps t n . ThermoPhonon code [18] solves the over-determined set of equations (5) and finds the effective force constant matrix D e ij (6), which is then used to construct the phonon spectrum.

A. Harmonic potential
Any smooth curve is approximately parabolic near its local minimum; it can be expanded in a Taylor series around it: . Hence, a vibration with a small enough amplitude around a LPEM is expected to be harmonic (see Figs. 1(a) and 2 and discussion in section III A). Average amplitude of atomic vibrations depends on T, and so does the harmonicity.
For the potential in Fig. 2 [same as Fig. 1(a)], one can categorize temperature dependence of atomic motion.
kT ≤ E H : harmonic vibration; E H < kT < E B : anharmonic vibration; T melt ∼ E B /k: phase transformation [20] from a solid to a liquid or a gas at kT ≥ E B .
In one dimension, a 1×1 matrix D αβ (4) has only one element D 11 = K. For the harmonic potential the force obeys the Hooke's law it is linear vs. x at any "small" displacement |x| < x H . Motion of a mass M in such potential is harmonic and has a frequency which does not depend on x at |x| < x H . The Thermo-Phonon method [18] at kT < E H and a displacement method at |x| < x H provide the same frequency (9).

B. Double-well potential
A dimpled potential in Figs. 1(b) and 3 is a less trivial example with two LPEMs per basin. It has a crystallographic high-symmetry point at x = 0, but E(0) = 0 is a local energy maximum, with two minima nearby at E(±x L ) = −E L . The point (x, E) = (0, 0) is the energy barrier E l = E L between those two LPEMs.
At kT E l ≤ E L , one can apply the smalldisplacement method (section III A) to one of the representative stable structures (section II C) to find phonons. Alternatively, at E L kT < E B , if multiple LPEMs could be interpreted as a negligibly small roughness at the bottom of an otherwise nearly harmonic potential (similar to that in Fig. 2), then one could use a finite displacement method to approximate phonons at the unstable atomic position at x = 0, see section III (B,C).
The 1D matrix D e in eq. 6 has dimension 1 × 1 (i.e., α = β = 1), and the value of its single element is K e . For the double-well potential in Fig. 3, the "effective" spring stiffness K e = F (x)/x depends on displacement x. A negative K e for a small displacement x < x L results in imaginary frequency ω e = K e /M , which characterizes an unstable phonon mode. The displacement x L with F (x L ) = 0 gives K e = 0 and ω e = 0. A large displacement x > x L leads to a positive K e and a stable effective phonon frequency ω e > 0, which depends on displacement x. In other words, the frequency ω e of vibrations around an unstable equilibrium at x = 0 is not welldefined. However, a choice of a finite x T ≥ x 0 , related to thermal motion of atoms at temperature T , such that can result in a phonon spectrum, which compares well with experiment at the same T . In eq. 10 we assumed E L kT and neglected the dimples in Fig. 3. The ThermoPhonon method [18] at small kT < E L returns K e < 0 and an unstable (imaginary) phonon frequency for the expansion (3) around unstable equilibrium at x = 0. With increasing T , the amplitudes of the negative K e and of the unstable imaginary phonon frequency ω e = K e /M become smaller, until at a sufficiently large kT ≥ E L the effective K e becomes positive. If K e ≥ 0, than ω e ≥ 0 is real and the effective phonon mode looks stable. Temperature T uniquely determines the phonon frequency ω e (T ) within this method [18].

C. A muffin-tin sombrero potential
To illustrate that a MEP between LPEMs can bypass a high-symmetry crystallographic point, we show a 2D sombrero potential in Fig. 4. This muffin-tin type model potential is symmetric around x = 0. In the radial coordinates (r, φ), it has a local energy maximum at r = 0 and an extended energy minimum at the radius r = x L (at any angle φ); the MEP between LPEMs is a circle r = x L . The energy barrier E l is zero for the muffintin sombrero potential (Fig. 4), but not for a corrugated sombrero (Fig. 6). A muffin-tin type potential can form a lattice, see Fig. 5. At a low temperature kT < E L , atoms (black dots in Fig. 5) vibrate near the PE minima at r = x L (blue circles in Fig. 5), while averaged over the angles φ atomic positions happen to be at r = 0. By construction, the central cross-section of the muffin-tin sombrero potential in Fig. 4 is identical to the double-well potential in Fig. 3, and the cross-section of this periodically repeated potential in Fig. 5 is identical to Fig. 1(b).

D. Corrugated sombrero potentials
Due to interatomic interactions (discussed in section II E), an actual atomic potential differs from an idealized muffin-tin sombrero model potential, which approximates a corrugated sombrero with E l E L in the limit E l → 0. A corrugated sombrero potential has multiple local LPEMs linked by a MEP forming a closed path, with small but finite barriers 0 ≤ E l ≤ E L , see Fig. 6.
At kT ≥ E l , atomic motion covers multiple LPEMs in a corrugated sombrero potential. However, at kT < E l an atom is trapped in the neighborhood of a single LPEM. Thus, there is a symmetry-breaking first-oder phase transition at T c , estimated by eq. 2.
In Fig. 7, MEPs (lines) between LPEMs (dots) form a loop in 2D and a network in 3D. Similarity between the muffin-tin and corrugated sombreros improves with increasing number of LPEMs, see Fig. 6. A cubic austenite with an unstable high-symmetry atomic position, such as NiTi, can have 48 symmetry-equivalent stable collective atomic displacements (LPEMs). The model sombrero potentials and a qualitative distribution of LPEMs linked by chains of MEPs in Fig. 7 help to understand atomic motion with large athermal displacements in real materials, such as NiTi B2-type austenite. Many solids with lattice instabilities have atomic potentials, which remind a corrugated sombrero. Examples include the austenitic phase of NiTi and the bcc phases of Ti, Zr, Hf, and Li.

E. NiTi austenite
We find that the austenitic phase of NiTi shape memory alloy above T c = 313 K has multiple LPEMs, separated by low energy barriers E l kT ∼ E L . While B2 structure is unstable [7], a representative stable structure [6] can be used for an approximate description of this solid.
To obtain a representative structure, we tried several unit cells of various shapes, with various number of atoms (Fig. 8). Computational details are in [6,7] and in Appendix A. Each supercell was heated to 800 K for 100 fs, cooled to 0 K in 800 fs (using ab initio molecular dynamics), and than fully relaxed to the nearest local potential energy minimum (using the conjugate-gradient algorithm). Among the results, we selected the structure with the lowest energy; the smallest one had 54 atoms per representative supercell. We checked that repeating the whole procedure with a twice larger, 108-atom unit cell (doubled along one of the 3 lattice constants) results in a different structure with the same energy per atom, see Fig. 8(d). We verified that this hexagonal 54-atom structure is indeed a local potential energy minimum with a stable phonon spectrum (Fig. 2b in [6]), and we found that its phonon DOS compares well to that obtained from the neutron scattering experiment [21]. Again, this representative hexagonal 54-atom structure is just an approximate representation of NiTi austenite.
Interatomic interactions can affect not only energy, but also atomic positions. Atomic displacements from the unstable ideal B2 positions are shown in Fig. 6 in [6] and in Fig. 8 (c,d). The largest one (0.66Å) is 25% of the nearest neighbor (NN) distance and 22% of the B2 lattice constant (3Å); it is above the Lindemann criterion for melting [22]. Nickel atomic radius r(Ni) = 1.24Å is smaller than r(Ti) = 1.47Å. An average displacements of Ni from B2 is larger than that of Ti ( Fig. 6 in [6]). The athermal NN-pair distribution function (Fig. 5 in [6]) is skewed and has a substantial width (it spreads from 2.42 to 2.88Å, and from 2.48 to 2.65Å at halfmaximum), in contrast to a δ-function at 2.60Å for an ideal B2 single crystal. An average next-nearest neighbor (NNN) distance for Ni-Ni is smaller than that for Ti-Ti, while NNN Ni-Ni and Ti-Ti distances are the same in ideal B2. There are more short Ni-Ni bonds than short Ti-Ti bonds, see Fig. 8 (c,d). The smallest NNN distance (for both Ni-Ni and Ti-Ti) in the austenite is shorter than the NN (Ni-Ti) distance in ideal B2 crystal. The NNN bonds form chains, which are linear for Ti-Ti and branching for Ni-Ni, see Fig. 8 (c,d). A representative LPEM has a smaller energy than B2 due to optimization of interatomic distances and bond angles. The cross section of the potential energy E vs. collective atomic displacement x from B2 [ Fig. 8(a,b)] at x = 0 to a representative 54-atom NiTi structure [ Fig. 8(c)] at x = ±x L and beyond is a double-well E(x) curve [see the shaded part of Fig. 9], which reminds Fig. 3. B2 also transforms without a barriers to the BCO ground state (see Fig. 4 in [7]). The barrier (B19 in Fig. 9) for the BCO-to-BCO transformation is well below B2, see Fig. 3 in [7] and Table I. In general, atoms are displaced from B2 along directions, which are not high-symmetry ones, and a cubic structure has 48 isometries, which form the O h octahedral symmetry group, isomorphic to S 4 × C 2 . Thus, for each atom there are at least 48 LPEMs around the energy maximum at ideal B2 [ Fig. 7(c)]. The barriers E l ∼ 1 meV/atom between those LPEMs are quite low: they are comparable to the barriers between the NiTi austenite (middle of thin black line in Fig. 9) and the BCO ground state, which vary from only 1 to 5 meV/atom, depending on the transformation path [23].
Stoichiometric NiTi austenite exists at temperatures between 313 K and 1293 K: it transforms to the low-T B19' martensite below T c = 313 K (40 C, kT c = 27 meV) [24,25], segregates above T s = 1293 K (1020 C, kT s = 111 meV), and melts at T melt = 1586 K (1314 C, kT melt = 137 meV). The LPEM representing the NiTi austenitic structure is ∼ 30 meV/atom above the ground-state BCO and ∼ 20 meV/atom below unstable B2. At T c < T < T s , E l kT , while E L ∼ kT . Thus, atoms in the NiTi austenite not only vibrate around a particular LPEM, but also move from one LPEM to another, forming a  [18]), and in a stable representative austenitic NiTi structure (LPEM) [6], compared to that from neutron scattering experiment [21].
dynamically changing pattern.
In spite of atomic motion across multiple LPEMs, the small displacement method applied to a stable representative austenitic NiTi structure [6] provides phonon DOS, which resembles experimental one [21], see Fig. 10. From the other hand, increasingly large atomic displacements from B2 structure suppress the relative weight of unstable phonons. One can compare results of 3 phonon methods (small displacement at B2, finite displacement at B2, and small displacement at a representative LPEM) with the  [26]) vary from 0 to 8 meV/atom [7,27,28], and increase to 15.4 meV/atom for B19 at ϑ = 90 • . B19 is the energy barrier for BCO-to-BCO MEP [7]. The barrier E l for LPEM-to-BCO [23] and LPEM-to-LPEM MEP is ∼ 1 meV/atom (10 K) above LPEM. Measured temperatures of martensitic transformation Tc, segregation Ts, and melting T melt [24,25]. assessment based on the neutron diffraction experiment [21] in Fig. 10.

F. Lithium
The bcc Li transforms to the ground-state close-packed rhombohedral 9R structure below T c = 70 K [9] and melts at T melt = 453.65 K. Its bcc structure is unstable. At T c T < T melt , including room T, atomic motion covers the whole central part of a PE basin, including the bcc high-symmetry point at x = 0 and the nearby LPEMs. At such conditions one can use a finite atomic displacement method to calculate phonons.
We compare various methods from section III in Fig. 11. Because the bcc Li structure is unstable, the small displacement method returns unstable phonon modes around N. As atomic displacements become larger, the relative weight of those unstable phonons reduces, until they completely disappear for a sufficiently large displacement (0.1Å).
The molecular dynamics (MD) at a sufficiently high temperature T T c provides predominantly large collective atomic displacements and a stable phonon spectrum [29]. The periodic boundary conditions in the 4×4×4 supercell [29], which is incommensurate with the 9R ground state, might also suppress a phonon instability. Using MD in a smaller 3 × 3 × 3 supercell, we find a very minor phonon instability at the ΓN branch, which reminds the result of a medium atomic displacement (0.05Å).
Overall, all 3 methods (based on the small displacements, finite displacements, and MD at a finite T) give phonon spectra, which are comparable everywhere, except for the ΓN branch. A phonon instability from the small displacement method is not a mistake, but a property of the unstable bcc Li structure. Suppression of  [29]. Phonon frequencies from small and medium displacement methods were rescaled to match the large displacement method at P. Unstable phonon modes are shown as negative frequencies.
those unstable phonons by using finite atomic displacements is just a computational trick.

G. Group 4 bcc metals: Ti, Zr, and Hf
Materials with a lattice instability are quite common. In particular, bcc Ti, Zr, and Hf are unstable [19,30,31]. They transform from the high-T bcc β-phase to the low-T hcp α-phase at cooling, and to the ω-phase at pressure. Using the SS-NEB method, we find that both β − α (see Fig. 12) and β − ω transformations are barrierless for all 3 metals, in agreement with the previous calculations for Ti [30]. Thus, ideal bcc is either a local energy maximum or a saddle point in these metals.
The calculated interatomic force in bcc Zr for displacements δ from bcc at 0 to ω at x ω was found to be negative (i.e., unstable) for small δ < 0.25x ω , but not for larger δ [38]. Our Fig. 3 explains how the finite displacement method provides stable phonons in such cases.
Using MD above T c , we find that an average PE of atoms at T of experimental bcc existence is well above that of an ideal bcc structure in Ti (Table II), Zr, and Hf. Thus, atomic motion covers the whole central part of a PE basin, including its center and all nearby LPEMs, which look like dimples, responsible for the lattice in-  stability. These PE dimples are shallow compared to kT , and this justifies the use of the finite displacement method for these bcc metals.

H. Stable bcc metals: Fe, Nb
Our story would be incomplete without mentioning a few conventional crystals. Some bcc crystals are stable. For example, magnetic iron has a bcc ground state [39] with a stable phonon spectrum ( Fig. 1 in [40], Fig. 4 in [41]); its transformation to hcp has a barrier [42]. The calculated phonons in bcc Nb (Fig. 2 in [43], Fig. 3 in [41]) are also stable and reasonably agree with experiment [44].

I. Tellurium
The Te A8, γ-Se, hP3, P 3 1 21 trigonal structure [45] can be interpreted as a hexagonal deformation of the simple cubic (SC) structure [16]. The SC structure is unstable and transforms to A8 without a PE barrier. The barriers E l between various orientations of A8 are so high, that atoms vibrate around a single LPEM, which is a stable A8 crystal, observed in experiment [46]. Tellurium is a provocative example, which can be interpreted either as a stable Te A8 structure with a single LPEM (see section IV A), or as a highly unstable SC structure with multiple LPEMs (section IV B), separated by high barriers E l (kT < E l < E L , section II).

V. DISCUSSION
In general, a dynamically polymorphic solid has a higher lattice entropy than a conventional crystal. Entropy is proportional to the logarithm of the number of states, and vibration around a single LPEM has fewer states than atomic motion across multiple LPEMs (see Figs. 1(b) and 7). Many solid phases are stabilized by entropy at a finite T, and dynamic polymorphism is very common in those high-temperature solid phases.
Examples of dynamically polymorphic phases include solids with lattice instabilities (such as bcc Li, B2 AFM FeRh, and NiTi austenite), crystals with a mobile interstitial dopant (some of metal hydrides and boron steels), polymers and organic molecules with rotating molecular units, and numerous solid phases, dynamically stabilized by entropy at a finite T.

VI. SUMMARY
Solids with dimpled atomic potentials are quite ubiquitous among natural and industrial materials. Some of them have multiple local minima of the atomic potential energy. They include many high-temperature solid phases, which have a higher lattice entropy than conventional crystals. Examples include many anharmonic crystals with lattice instabilities, such as bcc Li, Ti, Zr, and Hf elemental solids, B2-type antiferromagnetic FeRh and NiTi austenite. To understand properties of these materials, we have constructed simplified models. We compared three methods for calculating phonons, applied them to conventional harmonic crystals and solids with a lattice instability, and discussed various types of atomic motion. We predicted a first-order phase transition in cooled dynamically polymorphic phases, and provided an estimate of the transition temperature.