{'fraction_non_alphanumeric': 0.027663934426229508, 'fraction_numerical': 0.010048865069356873, 'mean_word_length': 4.941699836103957, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dynamic mechanisms of the structural phase transitions in KNbO 3: Molecular dynamics simulations

The question on the dominant driving mechanism (displacive or order-disorder) at each structural phase transition of KNbO3 is investigated by means of molecular dynamics simulations. To this purpose, we first develop a shell model by determining its potential parameters in order to reproduce the ferroelectric instabilities obtained by first-principles total energy calculations. The phase diagram as a function of temperature is obtained through constant-pressure molecular dynamics simulations. The analysis of the dynamical structure factor and the microscopic dynamics of the particles in the different phases allows us to reveal the nature of the dynamics associated with each structural transition. Correlations between local polarizations forming chain-like precursor clusters in the paraelectric phase are examined.


I. INTRODUCTION
Although the ferroelectric perovskite potassium niobate has been extensively studied by a variety of techniques, the dynamical nature of its structural transitions remains an open question. This compound undergoes a sequence of three phase transitions. With decreasing temperature it transforms from cubic paraelectric (C) to tetragonal ferroelectric (T) at 701 K, becomes orthorhombic (O) at 488 K, and finally rhombohedral (R) at 210 K. The question, whether the dominant dynamical character at each transition is displacive or order disorder, is of central importance, and numerous studies have been carried out in order to elucidate this point.
A ferroelectric instability is characterised by a maximum of the potential energy at the average positions of the ions in the high-temperature phase with respect to collective displacements which lead to the low-temperature phase. The average positions in the hightemperature phase correspond to an energy barrier for the collective motion of the ions between equivalent average positions of the low-temperature phase. The dynamical character of the transition is dominantly displacive or order-disorder depending on whether such energy barrier is respectively low oder high compared with the critical thermal energy k B T c . Within a purely displacive picture, there is a mode with the symmetry of the low-temperature phase which softens with decreasing temperature and becomes unstable, i.e. its frequency vanishes at the transition temperature.
Perovskites crystals have long been considered as displacive type ferroelectrics. The main evidence for this type of behavior has been the existence of a Γ-TO soft mode which has been observed in many perovskites [1,2] .
In the case of KNbO 3 , experimental evidences show that the sucessive C-T-O-R phase transitions cannot be driven solely by optical-phonon softening, which was found to be incomplete by means of infrared spectroscopy [3] . Moreover, if hyper-Raman and infrared reflectivity spectra in the cubic and tetragonal phases of KNbO 3 are interpreted in terms of a soft mode, this has to be assigned such a high damping that it can hardly be distinguished from a relaxator [3,4] . In fact, an unified interpretation of the experimental data has been given by assuming the coexistence of a relaxational mode and a soft phonon [3,5] . This picture has been supported by the observation of central components in light spectra [4][5][6] whose linewidths, lineshapes and symmetry properties are consistent with an eight-site model of the sucessive phase transitions. According to this order-disorder model, the potential energy surface has a maximum for the cubic perovskite structure and eight degenerate minima for the [111] displacements of the transition metal ion that correspond to the lowtemperature rhombohedral structure. In the cubic phase, the eight sites are occupied with equal probability, and this symmetry is broken as the temperature is lowered: four sites are occupied in the tetragonal phase, two sites in the orthorhombic phase, and finally, only one site is occupied in the rhombohedral structure. In this way, the relaxation process which displays a critical slowing down when T c is approached in the different phases is interpreted to be the driving mechanism of the phase transitions.
Several features derived from other techniques support the above view. Indeed, the most direct evidence that Nb atoms are displaced along the [111] directions, even in the paraelectric phase, comes from XAFS. Pair distribution functions obtained by this technique in the different phases, show the presence of three short and three long Nb-O bonds which is consistent with a rhombohedral local structure and an order-disorder mechanism for the the sucessive phase transitions [7][8][9] . However, in the orthorhombic phase of KNbO 3 a broad soft mode peak is clearly observed in Raman experiments near the O-R transition, without any quasielastic component [5] , suggesting a more displacive-like dynamics for this transition.
In addition, a more recent femtosecond time-resolved spectroscopy study for this phase rule out relaxational contributions of the same symmetry as the soft mode [10] . Also, other XAFS measurements at room temperature revealed that the direction of the Nb displacement is as close to the polar axis as the order of temperature atomic displacements [11] .
First principles calculations have contributed greatly to understand the origins of the transitions in KNbO 3 [12][13][14][15][16] and related perovskites [17][18][19][20][21] . The local spin-density approximation (LSDA) has provided a useful framework to elucidate important aspects of the underlying physics of these oxides, providing important quantitative information about electronic charge distributions, character of electronic bondings, crystal structure, phonons, structural instabilities, etc., providing considerable insight into the nature of the soft-mode total energy surface.
First-principles techniques are more powerful than any calculation based on empirical models. Nevertheless they are quite computer demanding and simulations of temperature driven structural transitions in perovskites are not computationally feasible at present. A successful approach to study finite temperature properties has been developed on the grounds of effective Hamiltonians whose parameters were fully determined from ab-initio calculations. This scheme has been applied to many ferroelectric perovskites [22][23][24][25] . However, the restricted dynamics considered, due to the lack of an atomistic description, makes this approach inappropiate to investigate the anharmonic lattice dynamics of these compounds.

II. MODEL AND COMPUTATIONAL DETAILS
In spite of the above mencioned experimental evidence, the lattice dynamical approach to ferroelectricity was primarily concerned with the microscopic description of displacive ferroelectrics. In this way, an accurate description of the ferroeletric soft-mode regime in ABO 3 perovskites was given on the grounds of the non-linear oxygen polarizability (NOP) model [26] . This model emphasizes the large anisotropic polarization effects at the oxygens produced by variations of the O-B distance. Such effects are expected in view of the strong environment-dependent oxygen polarizability and its enhancement through hybridization between oxygen p and transition metal d orbitals.
In the NOP model, the A and B ions are considered isotropically polarizable. On the other hand, an anisotropic core-shell interaction is considered for the O −2 ion, reflecting its site symmetry. The above anisotropy is described by two different linear core-shell coupling constants K 2 : one in the directions of the A ions and another in the direction of the B ion.
An additional fourth order contribution to the core-shell interaction along the direction of the B ion is taken into account by a coupling constant K 4 .
The harmonic part of the model is unstable against the ferroelectric mode displacements in the cubic structure. This happens essentially because the strong Coulomb O-B attraction overcomes the repulsive forces that hold the B ion at the cubic cell center. Therefore the stability should be provided by the fourth order term. This has been harmonically approximated by temperature averaging a pair of displacements, which is evaluated selfconsistently. Hence the temperature dependence of the complete phonon dispersion can be calculated. The effective harmonic term stabilizes the ferroelectric mode thus providing its temperature dependence. Applications to KTaO 3 [26,27] , SrTiO 3 [26] , KNbO 3 [28] , BaTiO 3 [29] , lead to results in excellent agreement with the experimental phonon data.
However, the self-consistent treatment of the quartic interaction in the NOP model washes out the potential maximum at the cubic structure in terms of the soft mode coordinate. We have recently shown [30] that an exact numerical treatment of the fourth order interaction at oxygen reproduces quite satisfactorily the ab-initio adiabatic potential for the relevant ionic displacements which lead to the series of ferroelectric transitions in KNbO 3 . A molecular dynamics simulation using this model allowed us to confirm a crossover from a displacive to an order-disorder character for the paraelectric-ferroelectric phase transition in KNbO 3 [31] . However, due to the constant-volume treatment performed, the simulation led to a direct transition from the cubic to the lowest rhombohedral ferroelectric phase.
To obtain a model which describes the phase transition sequence of KNbO 3 , it is necessary to allow for the homogenous strains involved in the phases; a fact that was originally emphasized by Chaves et al. [32] . This can be achieved by replacing harmonic force constants by interatomic potentials. So, we first calculate potential parameters from the force For the ab-initio total energy calculations we used the full-potential linear muffin-tin orbital (LMTO) method, within the LDA and the Ceperly-Alder exchange-correlation potential. Details of the method are given elsewhere [33,34] . The calculations were performed employing the same basis set used by Postnikov et al. [12] and the Brillouin zone integrations have been well converged in the number of k-points.
For the investigation of the temperature driven structural transitions we use constantpressure molecular dynamics (MD) simulations. Shell model molecular dynamics, though it has a long history, is dificult to use due to the treatment of the adiabatic degree of freedoms, i.e. the shells. Most workers have used a steepest descents method to relax the shell positions iteratively to zero-force positions on each step of the molecular dynamics. Although this procedure has been improved based on conjugate gradient relaxation of the shells [35] , in a typical simulation run an average of ten line searches are made within every time step of the simulation, reducing greatly the efficiency of the method in comparisson with rigidion model MD simulations. An alternative approach has been introduced by Mitchell and Fichman [36] in which the shells are given a small mass and their motion, like those of the cores, is found by numerical integration of their equations of motion. They showed that the results of this method are independent of the shell mass, provided it is enough small, and in agreement with those obtained using relaxation of massless shells. Regarding the efficiency of the method, the shortcoming of this approach is that the time step of the simulation must be reduced in order to provide enough accuracy for the integration of the shell coordinates.
We applied the latter approach in the present MD simulation study, which is carried out using the DL-POLY package [37] . The runs were performed employing a Hoover constant-(σ,T) algorithm with external stress set to zero; all cell lengths and cell angles were allowed to fluctuate. Periodic boundary conditions over 4x4x4 primitive cells were considered; the basic molecular dynamics cell therefore contained 320 ions (plus 320 shells which are addi-tional degrees of freedom). The time step was 0.4 fs, which provided enough accuracy for the integration of the shell coordinates. The total time of each simulation, after 5 ps of thermalization, was 45 ps.

A. Potential determination and static properties of the model
We have determined the model potential parameters in order to reproduce LMTO total energy calculations for different ferroelectric distortions, with and without lattice strain.
It is important to remark that one can not afford to reproduce very closely the ab-initio energetics with such simple potentials as the ones employed in our model. In addition, the adjustment of the potentials is not a straightforward and easy procedure, because all pair potentials and core-shell coupling constants contribute to the total energy of a given distorted lattice structure.
To evaluate the adiabatic energy surface for ferroelectric distortions, we have calculated the potential energy of the model for several atomic displacements. To this purpose, the shell coordinates (which represent the electronic degrees of freedom) are evaluated, for a given core configuration, by solving the adiabatic condition. Once the equilibrium solution for the shell coordinates is obtained iteratively by a steepest descent procedure, the potential energy is computed.
As a result of the potential determination, we obtain the curves shown in Figure 1,  Table I.
The bulk moduli, as in ab-initio calculations, can be obtained from the evaluation of the total energy as a function of the uniform volume expansion for the cubic phase. The model calculations yield a lattice parameter of 3.98Å for the static perfect cubic structure. The bulk moduli evaluated at this equilibrium volume is 226 GPa, which agrees fairly well with the LDA value of 208 GPa [12] .

B. Phonon dispersion relations
An interesting test of our modelling concerns the phonon dispersion relations. Recently, Yu and Krakauer [16] have performed first-principles phonon calculations for the ideal cubic  Figure 2. Although a very good agreement is achieved for the stable modes compared with Ref. [16], the imaginary phonon frequency for the ferrolectric mode at Γ turns out to be almost twice as large in our model calculation as compared with the LAPW result. This discrepancy could arise from the fact that the LMTO method produces an overestimation of the ferroelectric instabilities in perovskites compared with LAPW results [21] . Nevertheless, it is remarkable that the model reproduces the wave vector dependence of the instabilities, as obtained by the ab-initio linear response approach (compare Fig.2 with Fig.1 of Ref [16]), indicating chain-like instabilities in real space. As highlighted by Yu and Krakauer [16] , the finite thickness of the slab region of instability corresponds to a minimun correlation length of the displacement required to observe an unstable phonon mode. From the model's phonon dispersion curves, the length of the shortest unstable chain can be estimated to ≈ 4 a.

C. Phase diagram
We have used the above described model to perform constant-pressure MD simulations at several temperatures. In Fig. 3(a) we plot the order parameters (the three components of the mean polarization) as a function of temperature. The corresponding cell parameters are displayed in Fig. 3(b).
At high temperatures, the averaged polarizations p x , p y and p z are all very close to zero and the three lattice constants have almost identical values. As the system is cooled down below 725 K, p x acquires a value clearly different from zero, while p y ≃ p z ≃ 0, and the structure presents a considerable tetragonal strain (see Fig. 3(b)). This indicates the transition from the paraelectric cubic to the ferroelectric tetragonal phase. When the temperature is further reduced, the two lower ferroelectric phases appear: the orthorhombic one below ∼ 475 K, with clearly finite p x ≃ p y and still p z ≃ 0, and finally the rhombohedral phase below ∼ 175 K, with approximately equal values of the three polarization components. Although the model gives slighlty large cell parameters and distortions compared with experimental data [3] , the non-trivial phase transition sequence of KNbO 3 is correctly reproduced.
Regarding the transition temperatures, a good agreement with the experimental values is achieved. However, this can not be taken as a test on the quality of the model because a larger size of the simulation supercell could increase the values of T c . Thus a finite size scaling procedure would be necessary in order to determine the transition temperatures correctly. Such a treatment is beyond the scope of this work.
Finally, we point out that the C-T and T-O transition temperatures obtained from our simulation are considerably larger than the ones obtained using the effective Hamiltonian approach [25] . We believe there are two reasons for such a discrepancy. The first one is the already mentioned overestimation of the ferroelectric instabilities produced by LMTO compared with LAPW. The second one is the thermal expansion which is taken into account in our approach. The thermal increase of volume stabilises the various phases over a wider range of temperatures.

D. Dynamic character of the transitions
To gain insight into the dynamic excitations relevant to the ferroelectric phase transitions, we calculate the dynamical structure factor S(q, ω) at the Γ point. After equilibration, S(q, ω) is calculated with no thermostat and barostat, i.e. using conventional energy and volume conserving dynamics, as the space-time Fourier transform of the core-core displacement correlation function: where Q is an arbitrary wave vector, which can be decomposed in a reciprocal lattice vector G and a wave vector q within the first Brillouin zone, Q = G + q. A gaussian smoothing procedure is applied before the time Fourier transform is performed.
On the left hand side of Figure 4 we show the low-frequency range of S(q = 0, ω) at 800 K (cubic phase), 525 K (tetragonal phase) and 250 K (orthorhombic phase). In both the cubic and tetragonal phases a quasielastic peak is observed, while no peak appears corresponding to the ferroelectric (lowest TO) mode. The peaks observed at ≈ 200 cm −1 and ≈ 280 cm −1 correspond to the second lowest infrared active TO(Γ 15 ) mode and the silent Γ 25 mode of the cubic phase, respectively. The absense the soft-mode peak in our MD-spectra is consistent with the spectroscopic observations mentioned in the Introduction [3,4,6] . Furthermore, in more recent inelastic neutron scattering measurements of cubic KNbO 3 , the TO 1 optic phonon peak was not detected for reduced wave vectors below 0.2 [38] . This experiment also showed the existence of an anomalously flat and low-energy TA mode, strongly coupled with the TO 1 mode at q ∼ 0.2 and extending out to the Brillouin zone boundary. The analogous feature has been observed in KTaO 3 [27] . Krakauer et al. [25] obtained, through a MD simulation with an effective Hamiltonian, that the entire phonon branch softens in the cubic phase. This observation of a single branch softening arises from the fact that the Γ-TO and X-TA eigenvectors have very similar displacement patterns, aside from the wave vector modulation, and they are described as a single effective phonon branch in their approach.
Nevertheless, the dynamical behaviour of the effective degrees of freedom differs from the usual soft-mode picture of a displacive transition, and it is consistent with the previously obtained first-principles LAPW linear response results [16] .
To clarify the nature of the microscopic dynamics leading to the above remarqued features of S(q = 0, ω), we show on the right hand side of Figure 4  simulation study using an interacting polarizable ions model [39] . So, we can conclude that a relaxational slowing-down proccess is therefore mainly responsible for the C-T and T-O phase transitions in KNbO 3 .
While the ferroelectric transition (C-T) is generally discussed in terms of order-disorder models with relaxational-type dynamics, aharmonic damping effects are less severe at room temperature and below, and a simple soft-mode description appears to be adequate for the O-R phase transition [5,10] . The dynamical structure factor for the orthorhombic phase of sheets of the Brillouin zone [41] . These results suggested that the linear correlation of atomic displacements along the [100] directions are dynamic, in agreement with a dynamical model introduced by Hüller [42] a couple of years before. More recent high-resolution diffuse x-ray scattering measurements using synchrotron radiation have supported the dynamical model of correlations [43] .
The existence of correlated atomic displacements forming chains have been also supported for the already mentioned ab-initio linear response calculation of Yu and Krakauer [16] , which showed the existence of Brillouin zone planar instabilities in the cubic phase of KNbO 3 .
Recent MD simulations, using the effective Hamiltonian approach, have revealed preformed dynamic chain-like structures that are related to the softening of a phonon branch over large regions of the Brillouin zone [25] .
As it was showed in Section III-B, our model reproduces the wave vector dependence of the instabilities indicating, as obtained by the linear response approach, chain-like instabili-ties in real space. To further confirm this point we have calculated correlation functions for the three components of the local polarization as a function of the cell-cell distance d along a z-axis chain. This was done by defining the following correlation functions: where P i α (t) and P i+d α (t) are the α component of the instantaneous local polarization of cell i and i + d, respectively, with the condition that both cells belong to the same z-axis chain. The correlation function P αα (d) is defined as the space-time average of the product of those equal time local polarizations. These functions were calculated after equilibration, with no thermostat and barostat, using periodic boundary conditions over 8x8x8 primitive cells. The results at T=800 K are shown in Figure 5. While P xx (d) and P yy (d) decrease strongly when the cell-cell distance (d) increases, P zz (d) shows a slow decrease. This clearly indicates that local polarization components parallel to a given chain are highly correlated along that chain, while the perpendicular components are fully uncorrelated. An estimation of the correlation lenght can be done by fitting the points of Figure 5 to the function P αα (d) = e − d ξαα . However, due to the small size of the simulation supercell, the correlations of relatively more distant cells are affected by the periodic boundary conditions, particularly for longitudinally correlated cells. Therefore we use only the point correponding to the smallest d for the determination of ξ zz . This yields ξ zz ≈ 6a and ξ xx ≈ ξ yy ≈ 0.5a, and the corresponding exponentials are shown in Figure 5.
To visualize chain-like correlations in real space, we show in Figure 6 snapshots of instantaneous local configurations, in an arbitrary slice through the simulation supercell, at three different times. The slice was choosen with the z-axis as the vertical one. Three symbols distinguish the values of the z-coordinate of the local polarizations: ↑ when P i z > 15 µC cm 2 , ↓ when P i z < −15 µC cm 2 and • when |P i z | < 15 µC cm 2 (thus discarding as non z-polarized the cases where |P i z | is small). This plot allows to visualize clearly chain-like correlations of P i z through the appearance of finite lenght chains polarized up or down along the z-axis. It is dificult to estimate quantitatively the lenght of the chains due to their dynamical evolution, but some of them extend over the whole supercell width. It is interesting to point out that cells belonging to a given chain seem to correlate their motion, producing coherent polarization reversals.
To further confirm this point we have plotted in Figure 7 the time evolution of P i z and P i x of four consecutive cells which belong to a given z-axis chain. Chain-like correlations are evident due to the correlated motion (polarization reversals) of P i z and the fully uncorrelated dynamics of P i x .

IV. CONCLUSION
We have developed an atomistic model for KNbO 3 which describes its structural instabilities in good agreement with LMTO total energy calculations. A further molecular dynamics simulation allowed us to evaluate the phase diagram, reproducing correctly the non-trivial phase transition sequence of KNbO 3 .
Regarding the dynamical mechanisms of the phase transitions, were able to identify the dominant character of each one. We find that the C-T and T-O transitions have a predominant order-disorder character, signed by the presence of a central peak in the dynamic response spectra. This excitation appears because of the slow dynamics associated with a relaxational motion of local polarizations, which correlate within chain-like precursor domains in the paraelectric phase. On the other hand, the O-R transition is more displacivelike, showing a more oscillatory dynamics of the non-polar coordinate.
edges support from Consejo de Investigaciones de la Universidad Nacional de Rosario.