Soft-mode dynamics in the ferroelectric phase transition of GeTe

GeTe that exhibits a strong anharmonicity and a ferroelectric phase transition between the rhombohedral and cubic structures has emerged as one of the leading thermoelectric materials. Herein, combining molecular dynamics simulations and inelastic neutron scattering measurements, the lattice dynamics in GeTe have been investigated to reveal the soft-mode mechanisms across the phase transition. We have constructed a first-principles-based machine-learning interatomic potential, which successfully captures the dynamical ferroelectric phase transition of GeTe by adopting the neural network technique. Although the low-energy acoustic phonons remain relatively unaffected at elevated temperatures, the high-energy optical, and longitudinal acoustic phonons demonstrate strong renormalizations as evidenced from the vibrational phonon spectra, which are attributed to the large anharmonicity accompanying the phase transition. Furthermore, our results reveal a nonmonotonic temperature dependence of the soft-modes beyond the perturbative regime. The insight provided by this work into the soft-modes may pave the way for further phonon engineering of GeTe and the related thermoelectrics.


INTRODUCTION
Thermoelectric materials have spurred a considerable technological interest to alleviate the energy demands for their capability to directly convert thermal gradient to electrical energy and vice versa [1][2][3] . The maximum thermoelectric efficiency is characterized by a dimensionless figure of merit ZT = α 2 T/ρκ, where α is the Seebeck coefficient, T is the absolute temperature, ρ is the electrical resistivity, and κ is the thermal conductivity, constituted by electronic κ e and lattice κ l components. As the strong coupling effects among α, ρ, and κ e , the desired strategy to achieve advancement in thermoelectrics is to minimize the decoupled κ l through the suppression of phonon transport 4 . Although phonon scattering can be enhanced by various defects or nanostructures 5 , single-crystalline semiconductors that possess strong intrinsic lattice anharmonicity, while simultaneously maintaining favorable electronic transport properties have emerged as promising candidates for high-performance thermoelectric applications 6,7 . The IV-VI narrow band-gap semiconductors, Sb-doped GeTe, for example, exhibit a superior thermoelectric figure of merit, ZT >2, in both rhombohedral and cubic phases 5 . The high figure of merit of GeTe is ultimately attributed to band pocket convergence and ultralow κ l , demonstrating a great potential in the substitution of the toxic PbTe-based materials 8 for applications at T < 800 K.
GeTe crystallizes at room temperature with a rhombohedral structure (R3m) and undergoes a ferroelectric phase transition from rhombohedral α-GeTe to a more symmetric rocksalt β-GeTe (Fm3m) at a critical temperature (T c ) of around 650 ± 100 K, depending on the carrier concentration [9][10][11] . Below T c , the cubic rocksalt lattice undergoes a slight distortion with symmetry breaking induced by two distinct soft-modes, which can be represented as angular distortion through shearing along the [111] direction; the low-temperature phase then spontaneously acquires a ferroelectric polarization between the Ge and Te sublattices along the diagonal axis of the original cubic lattice 12 . Although the thermodynamic and vibrational phonon properties of GeTe have been investigated within the perturbation theory, a comprehensive understanding of its lattice dynamics is still hindered by the challenges in handling the lattice instability and high-order anharmonicity across the phase transition [13][14][15] . The validity of perturbation theory is limited to a situation in which the atomic displacement from the equilibrium position is relatively small compared to the interatomic distance 16 . Additionally, the nature of the structural R3m ⇆ Fm3m phase transition remains unclear and controversial. The transition has been categorized as second-order due to the translational group remains unchanged in both phases 17 and extrapolation of the soft-mode frequency to zero based on perturbation and self-consistent phonon theory 13,18 . The lattice discontinuity 19,20 and energy barrier 21 at T c , however, indicate a first-order transition. The vanish of rhombohedral lattice distortion above T c observed in the neutron diffraction experiment 11 , as well as the concurrent softening of the A 1 and E phonon modes with increasing temperature measured for α-GeTe in Raman scattering experiments 12,22 , suggest a displacive nature underlying the GeTe ferroelectric phase transition. The persistence of local structural lattice distortions in β-GeTe from static X-ray absorption fine structure spectroscopy 23 and time-domain terahertz spectroscopy measurements 24 , however, questions the displacive character of the phase transition and suggests an order-disorder type. The microscopic dynamical behavior of the soft-modes driving the lattice instability is a central question. A thorough investigation of temperaturedependent phonon spectra with a nonperturbative treatment of the anharmonic effects; therefore, is of great importance to better understand the vibrational properties of GeTe, as well as the nature of the ferroelectric phase transition.
Ab initio molecular dynamics (AIMD) is the most common nonperturbative approach to study lattice anharmonicity involving large atomic displacements and has demonstrated its capacity in predicting the vibrational phonon properties of crystalline and amorphous GeTe [25][26][27] . The system-size and time-scale limitations due to the cubic algorithmic scaling of density functional theory (DFT) implementations, however, restrict a direct computation of the temperature-dependent and mode-resolved phonon spectra throughout the reciprocal space within the ab initio realm. A feasible solution to circumvent the scaling problem is to develop computationally efficient potential energy surface (PES) functions in place of DFT calculations. Although the physically inspired interatomic potentials have proven to be effective 28,29 , accurate parameterizations can be challenging due to the preimposed constraints of the sophisticated functional forms. Recently, Behler and Parrinello have proposed a high-dimensional neural network (NN) potential 30 , which overcomes the parameterization problem by abandoning all physical constraints and employing a highly flexible functional form that can be accurately fitted to a large set of reference DFT calculations. The validity of the NN approach has been well demonstrated by accurately reproducing the PES of copper 31 , zirconia 32 , and titanium dioxide 33 .
Herein, we have constructed a first-principles-based interatomic potential for GeTe by employing the NN techniques, which successfully reproduces the dynamical nature of phase transition between the rhombohedral and cubic structures. The normal-mode decomposition technique 34,35 , which offers the capability to numerically depict phonon quasiparticles with full anharmonicity from equilibrium molecular dynamic (MD) simulations, has been applied to investigate the vibrational phonon properties of α-GeTe and β-GeTe. The critical mode-resolved anharmonic vibrational phonon spectra, involving structural distortion and thermal expansion contributions, have been computed at elevated temperatures, resolving the temperature dependence of the soft-modes at zone center and zone boundary across the phase transition. Complementing the theoretical study, inelastic neutron scattering (INS) experiments on powder GeTe samples have been performed.

Validation of the GeTe NN potential
The accuracy and transferability of the NN potential have been examined by calculating the correlation of the NN-predicted energies and atomic forces of all structures in the training and validation sets with the corresponding DFT calculations, as shown in Fig. 1. The root-mean-square errors (RMSE) of the total energy are 2.05 and 2.18 meV per atom for the training and validation sets, respectively, suggesting the convergence behavior of a welltrained NN potential. A zero-centered near-Gaussian distribution of the residual energy errors indicates that more than 95% of the randomly distributed validation structures are predicted within an error of less than ± 4.36 meV per atom (2σ), demonstrating the NN potential is capable of achieving high accuracy in interpolating the PES between discrete DFT calculations. The RMSE of the atomic forces for both training (0.074 eV Å −1 ) and validation (0.072 eV Å −1 ) sets have been substantially improved over the previous GeTe NN potential (0.46 and 0.47 eV Å −1 for the training and validation sets, respectively 36 ), which may be attributed to the incorporation of the DFT calculated forces in our training scheme. In addition, no sudden increase of errors is observed for the highenergy structures, which implicates the equal applicability of the NN potential for both α-GeTe and β-GeTe.
Thermal structure evolution of GeTe across the phase transition The thermal structure evolution of the ferroelectric GeTe phase transition is interpreted by analyzing the bond lengths of the average structure as obtained from the equilibrium MD simulations, as shown in Fig. 2a. At 300 K, the short (intra-layer) and long (inter-layer) bond lengths determined from the pair distribution functions are 2.88 and 3.22 Å, respectively, which are in good agreement with 2.83 and 3.18 Å as obtained from neutron diffraction measurements at room temperature 11 . The difference between the two distinct Ge-Te bonds of α-GeTe gradually diminishes, while approaching the phase transition temperature and finally converges at a bond distance of 3.05 Å at T c , in accordance with the single equivalent Ge-Te bond in the centrosymmetric cubic phase observed with the standard diffraction method 11,25 . The gradual disappearance of the Peierls distortion at elevated temperatures in the defect-free GeTe is supportive of the displacive phase transition scenario. However, it should be noted that the phase transition's displacive character remains controversial due to different concentrations of Ge vacancies in real GeTe samples and the averaging effect of conventional diffraction techniques 10,25 . In addition to the static average structure analysis, we have also calculated the probability distribution function (see Supplementary Fig. 3), and no appreciable distorted local structure is observed above T c . The dynamical motions of both Ge and Te atoms are also quantified by calculating the isotropic atomic displacement parameters (ADPs), as illustrated in Fig. 2b. The thermal motion of Ge atoms is found to be more intensive than Te atoms, which agrees well with previous X-ray powder diffraction measurements 21 and atomic onsite force constant calculations 10 . A sudden increase of thermal motion is observed for Ge atoms around T c , whereas it is less obvious for Te atoms, which is attributed to the increased anharmonicity near T c . The thermal expansion coefficients of α-GeTe and β-GeTe are predicted to be 5.36 × 10 −5 and 9.69 × 10 −5 K −1 , respectively, as shown in Fig. 2c, in reasonable agreement with neutron powder diffraction measurements of 4.59 × 10 −5 and 7.67 × 10 −5 K −111 . The increased thermal expansion coefficient above T c manifests the importance of lattice anharmonicity for the high-temperature β-GeTe. As seen in Fig. 2d, a small discontinuity of the rhombohedral angle appears near the phase transition, which is consistent with the first-order characteristic observed by the neutron diffraction experiment 11 . It is also noted that the anomalous volume contraction across the phase transition 11,19 was not observed in our MD simulations, which may be attributed to the differences between our ideal simulation model and real defective materials.
Anharmonic phonon dispersions of GeTe across the phase transition The anharmonic phonon power spectra have been calculated through the Fourier transform of the projected velocity autocorrelation function. With the renormalized phonon eigenvalues extracted from 1000 discrete wave vectors q sampled by a 10 × 10 × 10 supercell, the diagonalized dynamical matrix has been reconstructed, which is then applied to obtain the renormalized phonon dispersions over the whole Brillouin zone, as shown in Fig.  3 (a-c). The low-energy acoustic phonon branches are sharp and well defined at 100 K, while an evident broadening is observed for the optical branches. Upon heating to 400 K, the acoustic branches still demonstrate weak scattering, while the upper A 1 symmetry optical branch exhibits a pronounced softening and broadening at Γ, which is attributed to the structural phase transition. A substantial broadening is also revealed near the top of the longitudinal acoustic (LA) branch and the lower doubly degenerate E symmetry optical branch at elevated temperatures, owing to the increased anharmonic phonon-phonon scatterings. After the transition to the centrosymmetric β-GeTe phase at 700 K, a dramatic change in phonon dispersions is observed, where the optical branches become triply degenerated at Γ. A crossing between the transverse optical (TO) and LA branches is also observed, which is different from the avoided crossing behavior of PbTe [37][38][39] . Furthermore, it is observed that the group velocities of the acoustic phonons slightly increase when the rhombohedral phase transforms to the cubic phase above T c , suggesting the dominance of the phonon lifetime rather than group velocity in suppressing the high-temperature thermal transport, similar to the analogous SnS and SnSe with structural phase transitions 40 . In addition, as shown in Fig. 4 (a-c), the INS contour plots in the transverse acoustic (TA) modes dominated frequency range of α-GeTe remain relatively unaffected upon heating from 100 to 300 K. The imaginary dynamical susceptibility in a momentum transfer range of 0.5 to 4.5 Å −1 , as illustrated in Fig. 4d, demonstrates a nearly temperature-independent feature below 2.0 THz, which corresponds to the weaker anharmonic renormalizations of the TA  modes. The relatively large difference at around 0.60 THz is related to the elastic tail. Furthermore, it is seen from the inset that the corresponding low-energy zone center optical phonon (~2.40 THz at 100 K) softens at elevated temperatures, which indicates relatively large anharmonicity and is consistent with our MD simulations.

Soft-mode dynamics across the phase transition
A group-theoretical analysis of α-GeTe and β-GeTe indicates that the lattice polarization and angular shearing distortion along the diagonal axis across the phase transition correspond to the dynamical atomic motions described by the eigenvectors of the zone center Γ 6 and the zone boundary X 6 modes, where the subscripts denote the number of modes at the high-symmetry point. In order to further elucidate the anomalous temperature dependence of the soft phonon modes across the structural phase transition, the mode-projected power spectra are calculated at q = (0.0, 0.0, 0.0) and (0.5, 0.0, 0.5), as illustrated in Fig. 5. The Γ 6 (~3.80 THz at 100 K) decreases rapidly in energy at elevated temperatures and merges into a single degenerated Ramaninactive F 1u symmetry TO mode with Γ 4 (involving off-diagonal displacements of Ge and Te atoms) at T c . A less pronounced concurrent softening of the Γ 4 mode (~2.40 THz at 100 K) is also observed with increasing temperatures, consistent with previous diffraction measurements on GeTe samples 12,22,41 . Different from the double-peak feature of the TO mode of PbTe at Γ observed around 300 K due to the transverse acoustic (TA) and TO phonon nesting 42 , the zone center TO spectra of β-GeTe are well-defined single peaks. As shown in Fig. 5b, at the zone boundary X, a prominent softening has been revealed for the soft-mode X 6 and the neighboring phonon peak, while the two lowest modes exhibit weaker temperature dependences with less significant frequency shifts.
The mode-projected power spectra are fitted with a Lorentzian function to quantify the temperature dependence of the phonon frequency and linewidth that is inversely proportional to the phonon lifetime, as illustrated in Fig. 6. Below T c , a reasonable agreement of the temperature evolution of phonon frequency is found between MD simulations and Raman measurements 23 . The frequency of Γ 6 decreases continuously up to T c and then hardens with further increasing temperature. The temperature dependence can be fitted to the Cochran's relation in the form of ω(T) ∝ A |T − T c | α with a critical exponent of 0.38 and T c = 635 K. A deviation of the critical exponent from 0.5 for the soft-mode Γ 6 may be attributed to the large atomic fluctuation not incorporated in the linear response within the mean-field theory 43 . It is noted that the frequency of Γ 6 does not drop to 0 THz at T c , which is consistent with the time-domain infrared terahertz spectroscopy    41 , corresponding to a first-order phase transition behavior. Despite the limited softening can be supportive of the first-order scenario, the phase transition's nature has not been quantified thermodynamically in this work. In contrast, the α → β phase transition in the analogous SnTe is a second-order phase transition driven by the collapse of zone center TO phonon at T c ≈ 75 K 44 . In addition, the linewidth of Γ 6 manifests a dramatic increase while approaching T c , with the largest value of 1.43 THz at 500 K, demonstrating a drastically increased scattering rate. This is in reasonable agreement with the Raman measurement of 1.20 THz 22 near T c , at which the anharmonic phonon-phonon scattering mechanism is predominant (see Supplementary Table  I). A divergence from the monotonic temperature dependence within the perturbation theory is observed, suggesting the importance of structural phase transition to the anharmonic lattice dynamics. Above T c , the phonon linewidth of Γ 6 becomes less dependent on temperature. Furthermore, as shown in Fig. 6b, similar but relatively weaker temperature dependences of frequency and linewidth are revealed for the X 6 soft-mode, corresponding to the secondary role of unit cell shearing along the [111] direction in the R3m ⇆ Fm3m phase transition. For the low-frequency phonons at the zone boundary X that are not directly related to the structural phase transition, weaker temperature dependences of the phonon frequencies and much smaller linewidths are observed.

DISCUSSION
In summary, we have explored the lattice dynamics of GeTe across the structural phase transition through MD simulations and INS measurements. A first-principles-based interatomic potential, which is capable of reproducing the ferroelectric phase transition characteristics of GeTe, has been constructed by employing the NN techniques. Vibrational spectra at elevated temperatures have been computed from MD simulations, revealing strong lattice anharmonicity of GeTe related to the structural phase transition. The relatively weak temperature dependence of the low-energy acoustic phonon modes is well supported by our INS measurements on GeTe powder samples. Furthermore, nonmonotonic temperature dependences of phonon frequency and lifetime have been revealed for the soft-modes driving the lattice instability, which supports the first-order nature of the structural transformation. These results provide a detailed microscopic understanding of the relationship between the anharmonic lattice dynamics and phase transition, which may pave the way for further phonon engineering in advancing GeTe thermoelectrics.

Computational details
Collinear DFT calculations of GeTe were performed using the Vienna ab initio simulation package (VASP) 45 , implementing the Perdew-Burke-Ernzerhof generalized gradient approximation (GGA-PBE) exchange-correlation functional 46 . The projector-augmented wave (PAW) pseudopotentials 47 were adopted to represent core electrons, and 4s 2 4p 2 and 5s 2 5p 4 were treated as valence electrons of Ge and Te, respectively. A plane wave energy cutoff of 300 eV and a convergence criterion of 10 −5 eV were used in the electronic self-consistent calculations. A Monkhorst-Pack k-mesh 48 of 3 × 3 × 3 was applied to sample the Brillouin zone of the 3 × 3 × 3 supercell of GeTe (54 atoms). During the structural relaxation, lattice parameters and atomic positions were fully optimized until atomic forces were smaller than 1 meV Å −1 .
Our reference data comprises 6448 DFT calculated configurations, which were generated either by short AIMD simulations at elevated temperatures (up to 900 K) or systematically distorting the ideal structures of α-GeTe and β-GeTe in a range of ± 10% around the ground states. A cutoff radius of 6.8 Å was applied to determine the maximum distance that the NN framework considers atoms to be interactive. The local environment of each atom within this cutoff sphere was represented by thirty radial and angular Gaussian symmetry functions, which are invariant with respect to translations and rotations. In the NN potential, we used the hyperbolic tangent activation function, and it contains two hidden layers with seven nodes per layer. The atomistic machine-learning package (AMP) 49 , implementing the high-dimensional atom-centered scheme proposed by Behler and Parrinello 30 , was applied to train the NN potential. Of the final reference data, 10% of the structures were randomly selected as a validation set to avoid over-fitting and assess the transferability of the NN potential. In addition to the DFT total energies, the atomic forces were also incorporated in the training scheme with a relative weighting of 0.04 to obtain a more predictive fit of the NN potential.
With the first-principles interatomic NN potential of GeTe, equilibrium MD simulations were performed using LAMMPS packages 50 in a temperature range of 100 to 800 K, where a phase transition occurs at around 600 K. We used a 10 × 10 × 10 GeTe supercell (2000 atoms), and the duration of each MD simulation was 3 ns with a 1 fs time step. The atomic velocities and positions were collected every 30 fs under the NVT ensemble for phonon power spectra calculations. A convergence study of supercell size was performed and only insignificant differences were observed for the renormalized phonon dispersions between 10 × 10 × 10 and 12 × 12 × 12 supercells (see Supplementary Fig. 1). Before data collection, all systems were initially relaxed by minimizing the potential energies and then equilibrated for at least 1 ns using the NPT ensemble, where the average structures were applied to account for the thermal expansion. Second-order force constants obtained from the effective harmonic method as implemented in hiPhive 51 , were used to construct the dynamical matrix, which was then diagonalized to produce the temperature-dependent phonon eigenvectors. Power spectra were further decomposed by projecting the atomic velocities onto the temperaturedependent phonon eigenvectors using the normal-mode decomposition technique 52,53 , which allows us to directly extract phonon frequencies and linewidths from MD trajectories.

Experimental details
Polycrystalline GeTe samples (>10 g) were synthesized using high-purity elements (>99.99%). The raw materials were sealed in vacuum quartz ampoules and melted at 1073 K for 6 h, quenching in cold water followed by further annealing at 923 K for 3 days. The obtained ingots were then hand-ground into fine powders. Inelastic neutron scattering (INS) spectra were measured using the FOCUS 54 , a direct geometry cold neutron timeof-flight spectrometer, at SINQ of Paul Scherrer Institut, Switzerland. The sample was encased in a thin-walled cylindrical Al container. The sample temperature was varied from 100 to 300 K using a closed-cycle He refrigerator. We used an incident neutron energy E i = 16.9 meV (with elastic resolution δE 0 ≈ 2 meV). The data were analyzed using the DAVE software 55 .

DATA AVAILABILITY
All data generated or analyzed during this study are included in this published article (and its supplementary information files).