Perovskite Crystals: Unique Pseudo-Jahn–Teller Origin of Ferroelectricity, Multiferroicity, Permittivity, Flexoelectricity, and Polar Nanoregions

: In a semi-review paper, we show that the local pseudo-Jahn–Teller e ﬀ ect (PJTE) in transition metal B ion center of ABO 3 perovskite crystals, notably BaTiO 3 , is the basis of all their main properties. The vibronic coupling between the ground and excited electronic states of the local BO 6 center results in dipolar distortions, leading to an eight-well adiabatic potential energy surface with local tunneling or over-the-barrier transitions between them. The intercenter interaction between these dipolar dynamic units results in the formation of the temperature-dependent three ferroelectric and one paraelectric phases with order–disorder phase transitions. The local PJTE dipolar distortion is subject to the presence of su ﬃ ciently close in energy local electronic states with opposite parity but the same spin multiplicity, thus limiting the electronic structure and spin of the B( d n ) ions that can trigger ferroelectricity. This allowed us to formulate the necessary conditions for the transition metal perovskites to possess both ferroelectric and magnetic (multiferroic) properties simultaneously. It clariﬁes the role of spin in the spontaneous polarization. We also show that the interaction between the independently rotating dipoles in the paraelectric phase may lead to a self-assembly process resulting in polar nanoregions and relaxor properties. Exploring interactions of PJTE ferroelectrics with external perturbations, we revealed a completely novel property— orientational polarization in solids— a phenomenon ﬁrst noticed by P. Debye in 1912 as a possibility, which was never found till now. The hindered rotation of the local dipole moments and their ordering along an external ﬁeld is qualitatively similar to the behavior of polar molecules in liquids, thus adding a new dimension to the properties of solids—notably, the perovskite ferroelectrics. We estimated the contribution of the orientational polarization to the permittivity and ﬂexoelectricity of perovskite crystals in di ﬀ erent limiting conditions.


Introduction
There are innumerable publications devoted to the study of crystals with perovskite ABO 3 structure, reflecting the rich variety of their properties with wide-ranging applications in physics, chemistry, and materials science. In this paper, we show that there is a fundamental structural feature of such crystals related to the high, cubic-octahedral symmetry of their metallic B center, which defines their main properties. This feature is the pseudo-Jahn-Teller effect (PJTE). Its role in triggering the ferroelectricity of perovskite crystals was revealed first more than half a century ago [1] but gained a more developed form in recent decades [2][3][4][5][6][7]. Similar to the Jahn-Teller effect (JTE), the PJTE, under certain conditions, makes a high-symmetry configuration unstable with respect to

The Pseudo-Jahn-Teller Effect in ABO 3 Perovskite Crystals: B-Center Instability, Mueller's ESR Experiments
As outlined in the Introduction, the high symmetry of perovskite ABO 3 crystals makes them outstanding with respect to a variety of important properties. Among them, properties related to electronic degeneracy and pseudodegeneracy, which are directly controlled by symmetry, seem to be most important. The cases of electronic degeneracy leading to the Jahn-Teller effect (JTE) in the local states and their cooperative interactions in crystals are relatively easily recognized and well-studied [9,10]. Less attention has been paid until recently to the pseudo-JTE (PJTE), which is not seen directly from the initial structural and electronic data of the crystal. Meanwhile, as shown in a series of works [1][2][3][4][5][6][7], the main properties of some important perovskite ABO 3 crystals are controlled by the PJTE.
Consider the local PJTE in the paraelectric phase of BaTiO 3 . Very briefly, in cluster language for the octahedral unit [TiO 6 ] 8− (Figure 1), the local PJTE emerges from the vibronic mixing of the ground state A 1g , formed by the fully occupied six HOMO ((highest occupied (HO) molecular orbitals (MO)) t 1u and t 2u (mostly oxygen 2p orbitals) and three LUMO (lowest unoccupied, LU) t 2g , mostly titanium 3d orbitals (a total of nine MOs with the one-electron energy levels shown in Figure 2), by the polar t 1u type normal coordinates Q x , Q y , and Q z . This PJTE (A 1g + T 1u ) ⊗ t 1u problem [9,10] results in a 9 × 9 secular equation. Its solution yields the following adiabatic potential energy surface (APES) of the TiO 6 8− cluster, obtained already in the first paper on the subject [1]: where Q 2 = Q x 2 + Q y 2 + Q z 2 , 2∆ is the energy gap between the mixing electronic states, K 0 is the primary force constant for the Q displacements (stiffness of the crystal without the vibronic coupling), and F is the vibronic coupling constant (H is the Hamiltonian), In a recent, more rigorous treatment by the Green's functions approach [5], the local PJTE in the cluster unit was appended with the bulk crystal properties by taking into account the interaction of the Ti ion with the whole crystal via its electronic and vibrational bands. It improved the results by yielding appropriately band-averaged parameters instead of the local cluster ones.  The three-dimensional APES (1) has a specific form ( Figure 3). Under the condition the surface (1) has a maximum (meaning instability) when the Ti ion is in the center of the octahedron, eight equivalent minima placed along the four trigonal axes, in which the Ti ion is displaced toward three oxygen ions (away from the other three); higher-in-energy 12 equivalent saddle points along the six C 2v axes, at which the Ti ion is displaced toward two oxygen ions (at the top of the lowest barrier between two near-neighbor minima); and next six higher-in-energy equivalent saddle points, at which the Ti ion is displaced to one of the oxygen ions along the fourfold axes [1, 2,5].
With an additional two experimentally determined structural constants, the band-averaged energy gap 2∆ = 2.8 eV, and the vibrational frequency at the bottom of the trigonal minimumhω Е = 193 cm −1 , all the main parameters of this APES, shown in Table 1, were estimated [5], including K 0 , F, the positions of the minima Q x = Q y = Q z = Q 0 and first saddle points Q x = Q y = q 0 , Q z = 0, their PJTE stabilization energies, and the tunneling splitting δ 0 . The latter is a characteristic measure of the energy barrier between the near-neighbor minima of the APES. Its order of magnitude was first estimated by K. A. Muller and coworkers [27][28][29] in ESR experiments with probing ions by substituting Ti 4+ with the Mn 4+ ion, which produces a very similar APES (Mn 4+ is a "ferroelectric" ion, see below, Section 5), yielding approximate lifetimes in the minima (in sec.): 10 −9 > τ > 10 −10 . It is also seen in the NMR [30,31] and EXAFS [32,33] experiments: in the former, the characteristic "time unit" τ ~10 −15 and the Ti ion are seen in the trigonal minima in all the phases of BaTiO 3 , while in NMR experiments, τ ~10 −8 , and only an averaged picture is revealed (see the discussion in Section 3). Table 1. Numerical values of the PJTE vibronic coupling and APES parameters of the Ti active centers in the BaTiO 3 crystal. E PJT [111] is the PJTE stabilization energy at the trigonal minima, E PJT [110] is the stabilization energy at the maximum of the barrier between two near-neighbor minima [5]. The four phases of BaTiO 3 emerge directly from this APES by gradually populating the states in the minima with temperature and stepwise overcoming its barriers [1,5,6] (see Section 3). Of particular interest is the prediction of disorder in the orthorhombic and tetragonal ferroelectric phases and in the cubic paraelectric phase, and order-disorder phase transitions between them, discussed below in more detail.

Ferroelectricity in ABO 3 Perovskite Crystals-Order-Disorder Phase Transitions
Due to the mentioned above high symmetry, perovskite crystals like BaTiO 3 have an inversion center and no local dipole moments. For such systems, the previous dominant theories of ferroelectricity assumed that the ferroelectric distortion occurs as a result of the compensation of the local repulsions (resisting dipolar displacements) by long-range dipole-dipole attractions. In such "displacive" theories, it is assumed that the local charge-separating displacements are induced by a self-consistent interaction with displacements in the bulk.
On the other hand, it has been proven that the JTE and PJTE are the only source of spontaneous distortions of high-symmetry configurations of polyatomic systems (at T = 0). As shown in the previous section, in perovskites with an inversion center, spontaneous local polar distortions are possible due to the PJTE. This triggered the idea that the spontaneous polarization in ferroelectric crystals is due to the temperature-dependent ordering of considered above PJTE-induced local dipole moments [1], presently developed in the vibronic PJTE theory of ferroelectricity (see the review [2] and references therein). According to this theory, the local dipole moments should be present in the crystal before its spontaneous polarization (meaning in its paraelectric phase), in contrast to the displacive theories, where the dipole moments occur as a result of the spontaneous polarization.
At the time of the first paper on the vibronic theory [1], there were no experiments to test the drastic differences between the two theories. Later on, the number of experimental observations that contradicted the displacive theories of ferroelectricity gradually increased (continuing), including X-ray scattering and diffraction data [34], Raman and optical reflective experiments [35][36][37][38], the mentioned above ESR with probing ions [27][28][29], EXAFS measurements [32,33], NMR experiments [30,31], neutron scattering [39], and later in a variety of experimental studies (see more details and references in [2]; the authors of [34] did not cite the first prediction of the order-disorder phase transitions in ferroelectric perovskites [1] but acknowledged it in a personal letter; see [40]). In these experiments, several "peculiar" (for displacive theories) features in the ferroelectric properties were revealed in some ABO 3 perovskites, notably BaTiO 3 . Among them, we emphasize here the instant trigonal displacement of the Ti ion in BaTiO 3 in all its four phases, ferroelectric and paraelectric [32,33], in strong disagreement with displacive theories, in which the metal off-center displacement occurs as a result of the phase transition to the polarized phase. This fundamental property of perovskite ferroelectrics is confirmed by a variety of experimental data, including the above-mentioned ones, and actually by any observations that implicate the Ti ion. Not only is this ion instantly displaced along one of the [111]-type directions in the paraelectric phase, where the averaged symmetry is cubic, but it is as well displaced in this trigonal direction in the tetragonal phase, where the crystal symmetry and the macroscopic polarization are tetragonal [39]. "The most striking example is the off-center displacement of the Ti atom in barium titanate observed in the paraelectric phase way above the Curie temperature of the tetragonal-to-cubic phase transition. As accurate measurements indicate [32,33], the Ti ion remains instantly displaced closely along [111] directions throughout all the four BaTiO 3 phases, and the magnitude of the off-center displacement decreases monotonically by only 13% when heating from 35 to 590 K, showing no steps at the phase transitions." [32,33].
All these experimental observations and other empirical data (see [2]), as a pattern and in details, confirm the predicted earlier picture of the special, induced by the PJTE local dipolar distortions, which in cooperative interactions produce all observed ferroelectric properties. In a more recent work [6], using the numerical estimates of Table 1, the interaction between the PJTE-induced local dipole moments in BaTiO 3 was taken into account in a mean-field approximation, in which both the local off-center displacement and the mean-field of the environment are interdependent in a self-consistent way (see below). It yields the experimentally observed phase transitions in BaTiO 3 with reasonable values of Curie temperatures [6].
The picture of the B center of the ABO 3 crystal with the PJTE-induced eight-minimum APES in Equation (1) and local dipolar displacement in the minima, and their changes under the influence of the external mean-field E of the environment, is illustrated in Figure 4. White circles correspond to off-center equilibrium positions of the central ion B (the blue ball) in eight trigonal wells. (a) At no polar perturbation, the trigonal wells are symmetry-equivalent and therefore equally populated. The ion B has an equal probability to be off-center shifted (blue arrows) to each well and, on average, its off-center shift is zero. (b) Under a trigonal field E [111], one well is deeper than the others, and at low temperature, it is the only one temperature populated. The vector of average displacement (the blue arrow) is along the field parallel to the axis [111]. (c) In the tetragonal field E [001], the wells 1, 2, 3, and 4 remain symmetry-equivalent to one another; they are the lowest in energy and hence most populated. The mean displacement (red vector) points are along the axis [001]. (d) With the field E along the second-order axis, E [110], the wells 1 and 5 are the lowest and the only ones thermally populated. Therefore, the mean displacement (red vector) is along the electric field, E [110] [6].
Qualitatively, the phenomenon of ferroelectric ordering of the PJTE-induced local dipolar distortions is similar to magnetization in ferromagnetic crystals. The Hamiltonian of the system is H = m H m + mn H mn , where m and n label different unit cells of the lattice. Similar to Section 2, the one-cell term H m = H 0 (m) includes the PJTE in the octahedral unit [BO 6 ]. As for the inter-cell coupling term H mn , we employ the simplified version of dipole-dipole coupling, which follows from the physical nature of ferroelectric ordering. In the mean-field E of all other dipoles in the crystal, the energy of a selected dipole p m equals −p m ·E. At each site, the local dipole is influenced by the mean-field produced by the ordering of all other dipoles. At every other site, the induced electric dipole p n is proportional to the applied electric field E. Therefore, H mn = − 1 2 A mn p m ·p n with m n. Due to the symmetry of the crystal lattice, the inter-cell coupling parameters A mn are symmetric to swapping the indices, A mn = A nm . In what follows, the inter-cell constants A mn are combined into one correlation constant, which is determined by fitting to experimental data.
To decouple H mn , we apply the mean-field approximation. Based on the smallness of the fluctuation ∆p m = p m − p m from the average value p m , it provides reasonable answers far from the phase transition when |∆p m | << p m but does not apply close to the Curie temperature T C . In the ferroelectric phase with a uniform polarization, the average p m is the same at different sites, so p m = p . We insert p m = p + ∆p m into H mn and neglect the terms quadratic in ∆p m . Omitting the constant terms and plugging the resultant approximated expression into H mn yields the decoupled Hamiltonian of the whole crystal as a sum of independent (mean-field) Hamiltonians, H = m H MF (m). In each of its ferroelectric phases, the perovskite crystal has translation symmetry. The one-site Hamiltonian is the same for different unit cells, H MF (m) = H MF = H 0 -p·E. Here, H 0 is the one-site vibronic Hamiltonian yielding the PJTE in Section 2, and E is the mean-field, E = A p , with the inter-cell correlation parameter A = A(m) = n A mn . Due to the translational symmetry, E is the same for all sites, so is A = A(0) = m A m0 .
In the paraelectric phase with zero mean-field, p·E = 0, and there is no intercenter dipolar coupling. The one-site nuclear dipolar displacements are independent of the other sites. It reduces to uncorrelated hindered rotations of the polar distortions Q m at different sites (via tunneling or over-the-barrier transitions between the equivalent minima of the APES), yielding the bulk average Q = 0, meaning a disorder of local dipole moments. Lowering the temperature below T C creates conditions for spontaneous polarization. Even a weak mean-field E can violate the equivalence between the wells and quench the tunneling; by locking the system in the lowest-in-energy potential wells and orienting the local dipoles in the direction E (Figure 4), it polarizes the crystal.
The effect of the mean-field depends on the magnitude of the PJTE induced dipole moment p. The latter is a combination of nuclear and electronic contributions, p = p nucl + p el . It is easy to prove that due to the vibronic coupling, both p nucl and p el are proportional to the polar distortion Q induced by the applied electric field. Therefore, in what follows, we assume p = Z B eQ, where Z B is the so-called Born charge, and e is the elementary charge. The average value p of the local dipole moment is proportional to the off-center nuclear displacement, p = Z B e Q . Induced by the mean-field, it depends on Q in a self-consistent way (see below). Therefore, we present the mean-field Hamiltonian as In the lowest-temperature rhombohedral phase at T = 0 K, the strength of the mean-field is maximal, and the octahedral unit cell [BO 6 ] is locked in one of the trigonal wells at the bottom of the trough with the radius Q 0 . Therefore, at T = 0 K, the magnitude of Q approaches its maximum value Q 0 . At a higher temperature T 0 K, the nuclear motion delocalizes over different wells and the vector Q decreases in magnitude and changes its direction. Therefore, for the sequence of the corresponding ferroelectric phase transitions, we take Q /Q 0 as the order parameter.
In a certain unit cell, the delocalization changes with temperature. Nuclear motion depends on the radius Q 0 of the trough and the height of potential barriers between trigonal wells. Both are due to the one-site PJTE (see Table 1 for BaTiO 3 ), whereas the mean-field constant A is due to the inter-cell interaction, A = m A m0 . Generally, in cubic perovskites, there is a wide range of coupling parameters.
In what follows, we consider two limiting cases: (a) deep trigonal wells where the tunneling model applies, and (b) shallow wells along the warped trough of the APES with the hindered rotation of the nuclei. The two cases correspond to significantly different mechanisms of delocalization and, correspondingly, of the ferroelectric phase transitions. a The tunneling model for deep trigonal wells. Most of the time, the nuclear motion is localized in the trigonal wells with low frequency tunneling through the barriers from one well into another. In a typical perovskite with PJTE, the tunneling energy gaps are of the order of several meV (see the estimates for BaTiO 3 in Table 1). The excited states of hindered rotations along the trough are a few dozen meV higher in energy. In such cases, the tunneling model applies, in which we can (approximately) consider just the eight tunneling states and neglect any quantum entanglement of all other excited states.
At T > T C , in the paraelectric phase with zero mean-field, E = 0, self consistently, the order parameter is zero, Q = 0. The one-site Hamiltonian H MFA = H 0 has cubic symmetry, and all trigonal wells are symmetry-equivalent. With no tunneling (assuming that the potential barriers are infinitely high), the ground state corresponds to eight Born-Oppenheimer states localized in the wells |a , |b , |c , . . . , |h , labeled in Figure 1. Tunneling results in an even spread of the nuclear motion over the potential trough, lifting the eight-fold degeneracy of the ground state. The ground term splits into two triplets and two singlets, A 1g + T 1u + A 2u + T 2g [41]. The tunneling splitting depends on the three overlap integrals, a|b , a|c , and a|g . Evidently, for neighboring wells, a|b dominates. Due to the greater distance, in a good approximation, the other two, a|c and a|g , can be neglected. Measured from the ground-state energy level E 0 in the infinite-deep wells, the tunneling energy levels are equidistant, hω Е being the energy quantum of the transversal mode in a trigonal well, the one perpendicular to the corresponding trigonal axis, U(s) is the potential energy (the APES) in Equation (1), and v and w are the classical turning points where the system hits the barrier wall. The integral is taken over the arc length s (broken curve in Figure 3) along the path of steepest descent from the orthorhombic saddle point to the closest trigonal minimum points of the APES. For BaTiO 3 , the numeric estimate is Г ≈ 35 cm −1 = 4.3 meV [5].
Within the basis set of the eight tunneling states, the one-site Hamiltonian H = H 0 -λQ· Q of the mean-field approximation is represented by the following matrix [41], Here, γ j = λQ 0 Q j /(Г √ 3) and j = x, y, z is the dimensionless mean-field order parameter in terms of the tunneling constant Г. The three order parameters γ j are solutions of the following system of coupled transcendental equations, Z γ Q j γ = Tr(e −βH Q j ) and Z γ = Tr(e −βH ) with j = x, y, and z. Here, β = 1/kT, and the index γ means self-consistent dependence on the components γ j included in the mean-field Hamiltonian (5). In the lowest-temperature rhombohedral phase, all atoms B are coherently shifted into the same trigonal well, say, in the direction [111]. In this case, γ x = γ y = γ z = γ, and there is just one transcendental equation to solve. The mean-field E deepens the respective trigonal well. If E is strong enough, with the ordering energy p·E comparable to Г, it quenches the tunneling and locks the system in this well. Rising temperature weakens the mean-field, reducing its locking effect.
For different values of the correlation parameter j = p 0 A/Г, the calculated temperature dependence of the order parameter γ is shown in Figure 5. Plugging the resultant value of γ into the expression for the Helmholtz free energy, Φ = −∂(ln Z α )/∂β with β = 1/kT, we find the free energy Φ 111 of the rhombohedral phase. Figure 6 shows the temperature dependence of Φ for different ferroelectric phases where we assume j = p 0 A/Г = 5. Figure 5. Temperature dependence of the order parameter γ in the rhombohedral phase when the tunneling mechanism of disorder dominates. The correlation parameter, j = p 0 A/Г, varies from j = 1 for the left-bottom graph to j = 3 for the graph at the right top [6].
In the orthorhombic phase, all atoms B are off-center shifted along the two-fold symmetry axis [110]. As there is no minimum on the APES in this direction, we treat this "displacement" as resulting from an averaged motion evenly spread over two near-neighbor trigonal wells, the well b and the well f in Figure 1, separated by the orthorhombic potential barrier. In this case, γ x = γ y = γ and γ z = 0, and again we have just one equation to solve. Plugging the evaluated values of γ x = γ y = γ and γ z = 0 into the Helmholtz free energy, we find Φ 110 in the orthorhombic phase ( Figure 6). Solving the equation Φ 111 = Φ 110 for temperature, we find T C for the rhombohedral-to-orthorhombic phase transition. At around kT/Γ = 3.8, the graph of Φ 111 intersects with that of Φ 110 . Above T C = 3.8Г/k, the free energy Φ 111 of the rhombohedral phase is above the orthorhombic phase Φ 110 . The latter becomes energy-advantageous, and therefore, at T C = 3.8Г/k, the rhombohedral-to-orthorhombic phase transition takes place. In this model, its dependence on the inter-cell coupling parameter A is close to linear, kT C ≈ 0.8p 0 A − 0.3Г [6]. . Temperature dependence of the Helmholtz free energy Φ (in units of Γ) for the four phases, rhombohedral, orthorhombic, tetragonal, and the paraelectric, at j = 5. At kT < 3.8Г, the free energy of the rhombohedral phase is the lowest. With temperature, at kT > 3.8Г, the free energy of the orthorhombic phase becomes the lowest, causing the rhombohedral-to-orthorhombic phase transition. Above kT ≈ 4.6Г, the tetragonal phase becomes most energy-advantageous [6].
In the tetragonal phase, all atoms B are off-center "shifted" to a tetragonal position, parallel to one of the four-fold symmetry axes, say [001]. Due to tunneling through four potential barriers, similar to the orthorhombic phase, this apparent shift is the average over four close-neighbor trigonal wells, a, b, c, and d in Figure 4, and the respective mean-field order parameter can be set as γ x = γ y = 0 and γ z = γ. Plugging it into the expression of Helmholtz free energy Φ 001 and solving the equation Φ 110 = Φ 001 , we find the temperature of the orthorhombic-to-tetragonal phase transition. For j = 5, it takes place at kT C ≈ 4.7Г. Similar to the previous case, in this phase transition, the dependence of T C on the coupling parameter A is linear, kT C ≈ p 0 A − 0.3Г [6].
For BaTiO 3 , the experimental values of the two low-temperature phase transitions, rhombohedral-to-orthorhombic T C (I) and orthorhombic-to-tetragonal T C (II), are T C (I) = 178 K and T C (II) = 278 K. Repeated with different values of the correlation parameter A, the tunneling model calculated values of T C (I) = 204 K and T C (II) = 256 K are closest to the experimental values when j = 5. The percentage errors are 15% and 8%, respectively. A lower value of the mean-field parameter A better describes the lower-temperature phase transition than the higher-temperature one. This result is quite understandable because of the larger amplitudes of vibration in the wells at higher temperatures and hence the stronger inter-cell interaction of the distortions.
Among the essential results of the tunneling model is the first-order type of ferroelectric phase transition. The equations which we used to evaluate Helmholtz free energy follow from the first-order perturbation approach to the eight-fold degenerate ground state. In this theory, the small parameter is the tunneling overlap Г. In the case of barium titanate, as one can see from Figure 6, the tunneling model fails to provide a reasonable description of the high-temperature phase transition, tetragonal-to-cubic. It could happen by averaging over all eight trigonal wells. Such an averaging implies coherent tunneling through all twelve orthorhombic barriers. Obviously, due to the relatively high entropy factor, the respective probability is too low. An alternative tunneling path is through the high-symmetry point of the APES, which in the oxygen octahedron [BO 6 ] corresponds to the on-center position of the atom B. This potential barrier, of the order of 0.2 eV = 1600 cm −1 , is too high for tunneling.

b.
Shallow trigonal wells, low potential barriers. The warping of the two-dimensional trough on the APES may happen to be relatively weak, which results in moderately low potential barriers between trigonal wells, but the high-symmetry point, like in barium titanate, may be too high in energy, so the first-order tunneling does not hold. In this case, the high-temperature disorder is achieved by second-order processes through a temperature population of the excited states that are close in energy to the top of the corresponding potential barriers between near-neighbor wells. Such over-the-barrier hopping is similar to the Arrhenius-type activation in chemical reactions, and it resembles the Orbach relaxation in magnetic resonance.
The temperature of the ferroelectric phase transition T C is assumed high enough to populate the close in energy excited states, kT C ≥ ∆E JT . Then, the hindered rotation in the warped trough can be treated in terms of classical physics. The temperature averaging contains the exponent with the potential energy U = U 0 -λQ· Q of the mean-field Hamiltonian H = H 0 -λQ· Q . Since it includes Q as a parameter, we come to the following system of coupled transcendental equations: In the tetragonal phase, the atom B is off-center shifted in one of the three tetragonal directions, say, along the axis [001]. In this case, Q x = Q y = 0 and Q z = Q 0. Out of the three equations (6), we have just one to solve. With the same parameter values as used above in the tunneling model, the numerical solution of this equation is shown in Figure 7. With temperature, the order parameter Q /Q 0 decreases smoothly to zero at kT C ≈ 0.1236E JT . In the case of barium titanate with F 0 2 /K 0 ≈ 0.25 eV, this corresponds to kT C ≈ 250 cm −1 = 359 K, reasonably close to the experimental value of 373 K.
Thus, distinguished from the two low-temperature phase transitions, trigonal-to-rhombohedral and rhombohedral-to-tetragonal, the transition to the cubic phase is of second order. A fundamental conclusion, which was outlined already in the first paper of the vibronic PJTE theory [1], is that the ferroelectric phase transitions in such perovskite crystals are of order-disorder nature. This was first confirmed for BaTiO 3 in the experiments with diffuse X-ray scattering [32,33,40]. However, some more recent experimental findings-in particular, the low-frequency (soft) mode at the Curie temperature, observed in neutron scattering [42], IR absorption [43], hyper-Raman scattering [44], etc.-were interpreted as demonstrating displacive features of the phase transition (see also [45]).
As we noted earlier [2,6], there is no contradiction between these experimental observations and the vibronic PJTE theory. To begin with, the latter does not exclude the possibility of a low-frequency mode near the phase transition: it follows directly from the temperature dependence of the Helmholtz' free energy, Φ, with respect to the limiting phonon TO mode, [ω eff (T)] 2 = (∂ 2 Φ/∂q 2 ) 0 . The consideration of the phase transition above in this section is based on the mean-field approximation which does not include any details about phonon dispersion. It follows from the EXAFS experiments [32,33] that in BaTiO 3 , the magnitude of the local instant dipolar distortions does not change significantly as a result of the phase transition, meaning that the polarization process is similar to the magnetic ordering of local spins, which is of pure order-disorder type.
An attempt to include the phonon dispersion in a similar scheme involving the PJTE was undertaken in several papers [46][47][48][49][50], in which the PJTE is applied in a two-band approach, instead of our local two-state (in fact, several state) approach, outlined in Sections 2 and 3. The major difference between the two approaches is in the way the vibronic Hamiltonian is treated. In our paper, the mean-field approximation is used before any Fourier transformation to crystal waves is applied. It addresses the orientation degrees of freedom of the "pseudo spin" at each unit cell. Therefore, obviously, it provides a better description of the properties of local origin and the order-disorder phase transitions. In the two-band theory, the mean-field approximation is applied after the Jahn-Teller Hamiltonian is presented in crystal plane-wave form. It operates with crystal lattice modes that are uniform over the whole crystal, showing that the growing population of the low-energy portion of the electron energy band with lowering temperature strengthens the pseudo Jahn-Teller effect and reduces the respective curvature, (∂ 2 Φ/∂q 2 ) 0 . Evidently, this approach is adjusted to describe displacive transitions, but it fails to explain the variety of the experimental data of local origin [25][26][27][28][29][30][31][32][33][34][35][36][37][38][39]-notably, the instant trigonal displacements of the Ti ion in all the four phases, irrelevant to the phase transition [32,33].
The difference between the two approaches to the solid-state problem with the PJTE, local versus bulk (cooperative), is of a more fundamental nature than it may seem at first sight. Indeed, both the JTE and PJTE are of local origin: they are defined by off-diagonal matrix elements of the vibronic coupling terms in the Hamiltonian, which are significant only for the non-zero overlap of the electronic wavefunctions of near-neighbor atoms. Therefore, if the interaction between the JTE or PJTE centers are not strong enough, meaning before their ordering (before the structural phase transition), their local JT dynamics are not correlated, the translation symmetry of the crystal is obeyed in average, but not in each moment of time, and the traditional presentation of the crystal structure by electron and phonon bands is inadequate. To our knowledge, there are no worked out general methods (or computer programs) to handle such systems beyond our approach, described above based on the cited papers.
In view of this significantly local origin of the PJTE, another important issue emerges with regard to the interpretation of the mentioned above experimental data [42][43][44] as showing elements of displacive phase transitions in BaTiO 3 : it did not take into account the relativity to the means of observation (see [2,6,12]). Indeed, in the eight-minimum APES, the B ion under consideration has a characteristic "lifetime" τ ≈h/Г, where Г is the introduced above tunneling parameter. In cubic perovskites, Г ranges from 0.01 to 50 cm −1 . Respectively, the characteristic time is from τ ≈ 10 −9 s = 1 ns to τ ≈ 10 −13 s = 0.1 ps. Similarly, since (before the ordering) the transitions between the minima at different centers are not correlated, there is a related characteristic dimension l, which is the size of one elementary cell. For barium titanate, it is of the order of l ≈ 4 Å.
On the other hand, the experimental methods of observation have their own characteristic "time of measurement" τ exp , which is directly related to their frequencies. For example, the NMR technique measures the nuclear quadruple transitions and motional average dynamic displacements with a characteristic τ exp = 10 −8 s, whereas in EXAFS experiments, τ exp = 10 −15 . All the other widely used experimental techniques lie in between these limits. Similarly, they have characteristic length limitations l exp mostly determined by the wavelengths. Except for the EXAFS measurements, l exp for all the other widely used experimental methods are several orders of magnitude larger than the unit cell dimension where the disorder begins. Obviously, the experimental results in [42][43][44] reflect the crystal processes averaged over many unit cells, both in space and time, and hence they cannot be used as an indication of displacive versus order-disorder phase transition.
It follows from this discussion that the most appropriate experimental method to reveal the microscopic origin of ferroelectric properties of perovskite crystals is EXAFS, and the already performed experiments with this method fully support the conclusions from the vibronic (PJTE) theory.

PJTE-Induced Orientational Polarization in Solids-Application to Dielectric Susceptibility and Flexoelectricity
One of the most important novel properties of ferroelectric ABO 3 perovskites, revealed by the vibronic (PJTE) theory, is their strongly enhanced interaction with external perturbations (electric fields, pressure, strain, etc.). As outlined in Section 3, in the absence of external influence, the PJTE-induced local dipolar distortions in [BO 6 ] units are of dynamic nature, moving between the eight equivalent minima of its APES by tunneling or over-the-barrier transitions. In fact, the orthorhombic barriers between neighboring trigonal wells are relatively low, so the dynamics of the dipolar distortions may be regarded as some hindered rotations along a three-dimensional trough. In this respect, the crystal acquires properties of polar liquids, with all the consequences for the observable properties, including orientational polarization of crystals, first revealed in the vibronic PJTE theory [51][52][53]. Remarkably, the possible existence of dielectric crystals with randomly oriented dipole moments (similar to ferromagnetics), which may behave like a polar liquid above its freezing point, was first suggested by P. Debye [54] in 1912 (long before the experimental discovery of ferroelectricity), but not discovered till the present work. Similar to polar molecules in liquids, the orientational polarizability of solids with dynamic, PJTE-induced dipolar displacements is expected to be larger than displacive polarizability by orders of magnitude [51][52][53].
Different crystals vary with respect to the barrier height and the radius of the trough Q 0 (see Section 2). At high temperatures, T > T C , in the paraelectric phase with zero mean-field, the trigonal wells are symmetry-equivalent, the tunneling evenly distributes the nuclear motion over the trough, and for the ion B, on average, its off-center shift Q is zero (Figure 4a). The tunneling splitting of localized states in eight trigonal wells is proportional to the tunneling parameter Г. The necessary condition of tunneling is the resonance of the states, localized in the minima. In a typical ferroelectric perovskite with the PJTE, the energy gaps 2Г between tunneling energy levels are of the order of magnitude of several meV. Higher in energy with a few dozen meV are excited states of hindered rotations along the trough (for BaTiO 3 , the estimates are in Table 1). Any polar perturbation W, even not very strong, lowering one or several of the eight wells and/or raising the other ones violates the resonance, quenches the tunneling, and locks the system in the lowest well(s) (Figure 4b-d). The perturbation W causes a non-zero off-center shift Q , which produces anisotropy in the dielectric and elastic properties. A remarkably small magnitude of W~Г, of the order of, or even less than, meV, diminishes the tunneling and locks the system in the deepest minimum. A similar sensitivity, though to a lesser extent, is in the excited states of the hindered rotations. This unique feature of the PJTE in cubic perovskites determines their significant response to polar perturbations resulting in orientational polarization [51][52][53].
There are two types of polar perturbations: (a) electromagnetic fields targeting the electronic subsystem, and (b) applied strain distorting the crystal lattice. In what follows, we discuss both effects in ferroelectric perovskites. In both cases, we consider the paraelectric phase at T > T C with no mean-field and a relatively strong PJTE, when the tunneling approximation applies to the ground state. This allows us to reduce the problem to just eight tunneling states. For all the remaining excited states, we assume that their thermal population is small, and hence their contribution to the corresponding response function is negligible. The effect of locking the PJTE dipolar distortions is shown in Figure 8 [53], where we use an applied electric field E as an example of the polar perturbation, W = − p·E, and the corresponding Hamiltonian is H = H 0 − p·E, with H 0 as the Jahn-Teller Hamiltonian (Section 3). The magnitude of the induced dipole moment p is a combination of nuclear and electronic contributions, p = p n + p e . As shown in [53], in the case of static applied field, the nuclear contribution dominates, |p n | >> |p|. Therefore, as mentioned above, p ≈ p n = Z B eQ, where Z B is the so-called Born charge, and e is the elementary charge constant. The induced off-center nuclear displacement Q is expected to approach its maximum value Q 0 under a trigonal field of infinite strength, when E → ∞ and the octahedral unit cell [BO 6 ] is entirely locked in one of the trigonal wells ( Figure 4b). Accordingly, for the induced dipole moment, its asymptotic value is p 0 ≈ Z B eQ 0. When E [110], the perpendicular component averages out, Q z = 0, and the asymptotic value is Q →(Q 0

Dielectric Susceptibility
In a cubic perovskite, in its paraelectric phase with zero mean-field and no applied field, E = 0, the crystal is not polarized. Under an applied electric field, a non-zero dipole moment p is induced in each unit cell, and the crystal becomes polarized. By definition, the vector of polarization P is the average dipole moment per unit volume, P = p /a 3 . Here, a is the lattice constant, p = Tr(p i e −βH )/Tr(e −βH ), with the averaging over both the temperature and (statistical) over different unit cells, Н = H 0 -p·E is the vibronic Hamiltonian of the octahedral site [BO 6 ] in the applied electric field, and β = 1/kT. Regularly, there is a linear relationship between the generalized "displacement" P and the "force" E, namely Р i = ε 0 Σ j α ij E j . It corresponds to the first terms in the power expansion of P in terms of E. The linear factor α ij is the rank-two tensor of dielectric susceptibility or polarizability. In cubic crystals with no PJTE, the three principal values of the tensor α ij are equivalent, and the polarizability is the same in all directions. The dielectric properties of such a crystal are isotropic, and the principal axes have no special direction. As follows from the numeric results in Figure 8, to be specified below, in cubic perovskites, the PJTE brings about two important properties: (1) the dielectric susceptibility has cubic anisotropy; the principal axes are "tied" to the symmetry axes [100], [010], and [001], and (2) the induced dielectric susceptibility is non-linearly dependent on the applied electric field. The polarizability, and, correspondingly, the dielectric permeability, is not a constant but depends on the strength of the applied field.

a.
Low potential barriers and/or relatively weak electric fields, E << Г/p 0 . This limiting case of E → 0 corresponds to the left side of Figure 8, close to its vertical axis where the field-dependence on Q , and hence on P, is linear. The consideration can be limited to just one localized state in each trigonal well (the potential trough is two-dimensional; even in a very shallow well, there is at least one localized quantum state). Under the applied field, the wells are not exactly symmetry-equivalent. However, since the electric field is relatively weak, the field-induced symmetry breaking is not significant. The localized states in trigonal wells are close to the resonance, with an even distribution of the nuclear motion over the trigonal wells, and the unit-cell octahedron [BO 6 ] carries frequent tunneling transitions or/and hindered rotations close to the bottom of the trough. Accordingly, assuming P i = ε 0 Σ j α ij E j , we have (∂P i /∂E j ) 0 = ε 0 α ij , where the subscript "0" means the derivative at E = 0. Plugging p i = Tr(p i e −βH )/Tr(e −βH ) into the definition P = p /a 3 , we come to α ij = αδ ij with Hence, in a weak electric field, the polarizability is isotropic and the principal axes of the tensor α ij have no specific direction. At relatively low temperatures and/or large Г, when kT << Г (meaning Гβ → ∞), the temperature factor in Equation (7) approaches unity, and α ≈ p 0 2 /(3ε 0 a 3 Г).

b.
Strong electric fields and/or deep potential wells, E ≥ Г/p 0 (this case apparently does not apply to BaTiO 3 , with Г ≈ 35 cm −1 and p 0 ≈ 5 D, E = Г/p 0 = 46 kV/mm, which is above its dielectric breakdown). In this case, as follows from Figure 8, the field-dependence of Q and, hence, of p is non-linear. Therefore, for the dielectric susceptibility, instead of the derivative (∂P i /∂E j ) 0 , we can use the average α ij = (ε 0 ) −1 (∆P i /∆E j ). Since P(0)= 0, we find ∆P i /∆E j = (P i -0)/(E j -0) = P i /E j . For the applied electric field pointing along the trigonal axis, E [111], we have E x = E x = E x = E/ √ 3. According to Figure 4b, it lowers the trigonal well 1 in Figure 1, lifts the opposite well 7, and keeps the remaining potential wells unchanged. At low temperatures, the dipole moment Q x = Q y = Q z ≈ Q 0 / √ 3, so that Q = Q 0 . In this case, the polarization P ≈ p /a 3 = р 0 /a 3 . Similar results can be obtained for the cases when the applied electric field points along the symmetry axes [110] This angular dependence can be approximated by the cubic invariant, where Y 4 (φ, θ) and Y 6 (φ, θ) are the spherical harmonics, Y 4 (θ, ϕ) = (x 4 + y 4 + z 4 − 3 5 r 4 )/r 4 and Y 6 (θ, ϕ) = x 2 y 2 z 2 /r 6 . As follows from Figure 8, if the applied electric field is strong enough, the system is locked at the bottom of the lowest wells of the APES, and the polarization approaches its asymptotic value shown by the dotted line in Figure 8. Therefore, in all three formulas of Equation (8), the polarizability is inversely proportional to E, approaching zero at E → ∞.

c.
High temperature and/or deep potential wells, kT >> Г. This case is of special interest. In most ferroelectric perovskites, the Curie temperature T C is of the order of a few dozen meV, while Г is only a few meV. Therefore, for most cubic perovskites in the paraelectric phase, the condition kT >> Г applies. Tunneling energy gaps are of the order of 2Г. Therefore, in the high-temperature case when β → 0, the Boltzmann exponent ехр(−βH) can be approximated as 1−βH, and the temperature average p = Tr(pe −βH )/Tr(e −βH ) simplifies p ≈ Tr[p(1 − βH)]/Tr(1 − βH) ≈ (−β/8)Tr(pH). At this point, one can take advantage of the invariance of the trace of a matrix to any unitary transformations of the basis set. In particular, we can use the basis of the eight tunneling states when E = 0 with the well-known matrices 8 × 8 for the dipole moment p and the Hamiltonian H. Therefore, under any applied field, Tr(Q) = 0, and Tr(pH) = −8p 0 E/3. Hence, at relatively high temperatures in the applied field of any strength, we have p ≈ p 0 2 /(3kT), resulting in the Langevin-Debye equation α ≈ p 0 2 /(3ε 0 a 3 kT). The dielectric susceptibility is isotropic with the arbitrary direction of the principal axes. The orientational contribution to the dielectric susceptibility α caused by the PJTE-induced dipole moment rotations in the trough in the limits of high temperatures and deep potential wells is thus inversely proportional to temperature. However, the experimentally measured value of α also includes the significant contribution of the PNR (see Equation (16) in Section 6).

Flexoelectricity
Similarly, in cubic perovskites, the PJTE produces an enhanced response to another kind of dipolar perturbation, the strain induced by applied forces, which acts directly upon the crystal lattice. Caused by the stress, the corresponding strain is a rank-two symmetric tensor u with components u mn = 1 ⁄2(∂U m /∂x n + ∂U n /∂x m , where U is the vector of deformation, excluding the rigid-body motions. We follow the traditional notation x n with the index values n = 1, 2, or 3 corresponding to rectangular coordinates x, y, and z. In centrosymmetric crystals, the parity of u is even. Therefore, in the paraelectric phase, a uniform strain does not cause any polarization. A non-uniform strain with a non-zero gradient ∂u mn / ∂x j 0 is odd, and hence it can polarize the crystal. Like in the electric-field case, any non-uniform strain, even a weak one, breaks the resonance between the tunneling states, quenches the tunneling, and locks the system in lowest well(s) (Figure 4). The polarization P induced by the strain gradient is called the flexoelectric effect [56][57][58]. The stronger the perturbation ∂u mn / ∂x j , the greater the induced polarization. As shown in Figure 8, if strong enough, it can induce a non-linear polarization with a cubic anisotropy. In what follows, we limit the consideration with a relatively weak strain gradient ∂u mn / ∂x j and its linear response.
Expanding polarization P in terms of ∂u mn / ∂x j , we can keep just the first non-zero terms, Similar to the influence of an external electric field, due to the break of the local inversion symmetry induced by the PJTE, a relatively small strain gradient can induce a measurable polarization. The rate of change of P i with the strain gradient ∂u mn / ∂x k may be called flexoelectric susceptibility. The strength of the flexoelectric coupling is determined by the four-rank tensor f with components f ijmn , called flexotensor. The elements ∂u mn / ∂x j of the strain gradient can be treated as components of the rank-three tensor u = ∇ ⊗ u.
In the cubic symmetry group O h , the components of the vector operator ∇ form the basis of the irreducible representation T 1u , while the components of the symmetric tensor u transform as the symmetric square [T 1u 2 ] = A 1g + E g + T 2g . Therefore, u = ∇ ⊗ u transforms as the product T 1u The components of the rank-three tensor u form the basis of a reducible representation. To find the irreducible combinations u Λλ (G), we apply the so-called Clebsch-Gordan decomposition of the group theory. Here, Λ is the resultant irreducible representations A 2u , E u , 3T 1u, and 2T 2u , λ being their rows, and G, similar to the quantum number seniority in atomic spectra, labels the original term in the product, T 1u × A 1g , T 1u × E g , or T 1u × T 2g . The left side of Equation (10) is a vector transforming as T 1u . Hence, of the seven representations A 2u , E u , 3T 1u , and 2T 2u , we keep only the three vector representations T 1u , originating from the above products. We present components of the vectors u (A 1g ), u (E g ) and u (Т 2g ) by the expressions u j (Λ)= k,λ (∂u Λλ /∂x k )T 1u kΛλ T 1u j, where Λ 1 γ 1 Λ 2 γ 2 |Λλ are the Clebsch-Gordan coefficients. In these terms, Equation (10) takes the form similar to the one with susceptibility, Instead of just one tensor ε 0 χ ij , it includes three tensors, f ij (A 1g ), f ij (E g ), and f ij (T 2g ). The applied non-uniform strain induces a local electric field E polarizing the crystal, so that the Here, V Λ (with Λ = A 1g , E g, and Т 2g ) are the corresponding coupling constants. Assuming that the strain gradient is relatively weak, and involving the relations ∂P i /∂u j = (∂P i /∂E j )(∂E j /∂u j ) = V Г (∂P i /∂E j ), we get f ij (Λ) = f (Λ)δ ij with Λ of the products T 1u × A 1g , T 1u × E g and T 1u × Т 2g . Quite similar to the procedure that leads to Equation (7), Equation (11) yields the following: Thus, in cubic perovskites, among the 81 components of f ijkl , only three reduced matrix elements f (Λ) with Λ = A 1g , E g , and T 2g remain independent. In addition to the cubic symmetry-related constraints, the flexotensor is invariant to index transpositions [59] f ijkl = f jikl and f ijkl = f ilk j . Therefore, the three reduced matrix elements can be expressed in terms of just two components, √ 2, and f (T 2g ) = 2f 1122 √ 3 [60]. Accordingly, just two measurements are required to determine f (Λ) experimentally. For example, to find the tetragonal component f ( 2)(∂u zz / ∂z) −1 , it is sufficient to measure the polarization that occurs under a non-uniform strain applied along the crystal axis [001], with the strain gradient ∂u zz / ∂z 0. Similarly, if the strain gradient is along the trigonal direction [111], one can find f (T 2g ) = P/u (T 2g ).
The applied strain shifts the equilibrium positions of all ions in the unit-cell octahedron [BO 6 ]. Therefore, in Equation (12), the coupling constants V Г are related to the vibronic coupling constant F in Sections 2 and 3. For example, consider a tetragonal uniaxial strain along [001] with u q = u zz and u e = ( √ 3/2)(u xx -u yy ) = 0 or, equivalently, u xx = u yy = − 1 ⁄3u q and u zz = 2 ⁄3u q .
It shifts the equilibrium position R n j of the jth ion in the nth unit cell to R nj = X (0) nj 1 + 2 3 u nθ , where j = 1, 2, . . . , 7 labels the seven ions in the octahedron [BO 6 ]. We need the site index n in the symmetry-adapted strain u nθ to include its non-uniform nature. In the perovskite crystal lattice, the adjacent octahedrons [BO 6 ] share a bridge oxygen atom. In the next-neighbor octahedron [BO 6 ] n+1 shifted in the direction [001] by one lattice constant a, due to the non-zero strain gradient, the symmetry-adapted coordinates are slightly different. Comparing the tetragonal deformation of the adjacent unit cells, let u θ = ∆u q /∆z = (u n+1,q − u nq )/a where, as above, a is the lattice constant. It follows that the non-uniform tetragonal deformation creates a dipolar distortion Q z → Q z + 0.1a 2 u θ or, in other words, ∆Q z = 0.1a 2 u θ . This gives V E = 0.1Fa 2 = 0.1Z B ea 2 u θ . The corresponding energy increment F∆Q z = V E u θ = 0.1Fa 2 u θ can be treated as being due to a local electric field E induced by the non-zero gradient u θ . Then, 0.1Fa 2 u θ = p 0 E and, therefore, E = 0.1Fa 2 u θ /p 0 . Correspondingly, the induced polarization is P z = ε 0 αE = 0.1ε 0 αFa 2 u θ /p 0 .
The resultant component of the flexotensor, f 1111 = ∂P z /∂u θ = 0.1ε 0 αFa 2 /p 0 , is proportional to the dielectric susceptibility α. For BaTiO 3 , with the experimental value of α of the order of 2 × 10 4 , we find f 1111 = 0.1ε 0 αFa 2 /p 0 ≈ 0.66 µC/m, which is around 2600 times greater than was found in [58] and much nearer to the experimental value. Besides, with the orientational polarizability, as it follows from the outlined here PJTE theory, the flexoelectric factor is positive (in accordance with experimental data), whereas without the latter, it emerges as negative [58]. For other limiting conditions, rough estimates [51][52][53] yield even higher orientational contributions to the flexoelectric coefficients.

Multiferroicity in ABO 3 Perovskites with B(d n ) Configurations
The vibronic theory of ferroelectricity in ABO 3 perovskite crystals based on PJTE-induced local dipolar distortions in the B centers of the perovskite crystal was extended to formulate the necessary condition of coexisting magnetic and ferroelectric (multiferroicity) properties [3,4]. "Multiferroicity implies that the ferroelectric crystal, which is a dielectric, has also a nonzero magnetic moment, meaning unpaired electrons. In the ferroelectric BaTiO 3 the d 0 configuration of the Ti 4+ ion has no unpaired electrons, and attempts to obtain ferroelectricity in perovskites with d n , n>0, configurations of the transition metals B ions were unsuccessful for a long time. This prompted some authors to term the situation as a ferroelectric "d 0 mystery" [61][62][63], accompanied by a conclusion that nonzero spin states are detriment to ferroelectric polarization. However, more recently quite a number of ferroelectrics-multiferroics, mostly perovskites with configurations d 3 -d 7 , were obtained and studied" [3,4].
The origin of these special properties of perovskite ferroelectrics with d n configurations does not follow directly from displacive theories and had no general explanation for a long time. The vibronic (PJTE) theory of ferroelectricity outlined above elucidates the role of spin states in the local polar instability, explains the origin of perovskite multiferroics with proper ferroelectricity, and formulates the necessary conditions that ABO 3 perovskites with a magnetic d n configuration of the B ion may be ferroelectric [2][3][4]. Moreover, because of the spin implication, the multiferroics conditions that emerge from the PJTE for d n ions with n = 3, 4, 5, 6, and 7 are directly influenced by the well-known transition metal high-spin/low-spin crossover, resulting in the coexistence of three phenomena: ferroelectricity (FE), magnetism (M), and spin crossover (SCO). This, in turn, leads to a quite novel phenomenon, magnetic-ferroelectric (multiferroics) crossover (MFCO), creating a rich variety of possible magnetoelectric and related effects [3,4]. The vibronic (PJTE) theory, exclusively, reveals the role of spin in the spontaneous polarization of crystals.
Referring to the typical MO energy scheme for the octahedral cluster BO 6 8− discussed above (Section 2) and shown in Figure 2 for the MO electron population of the d 0 configuration, e.g., when B ≡ Ti, we see that the HOMO in this case is t 1u , which is a three-fold degenerate odd-parity linear combination of mostly oxygen p π orbitals, while the LUMO is t 2g , mostly atomic three d π orbitals of the transition metal ion B, and the next excited MO is the double degenerate one e g (the non-bonding oxygen b 1g MO is not shown as irrelevant). Using the arrows up and down to indicate the two spin states, we find for the d 0 case the HOMO configuration (t 1u ) 6 = (t 1u ↓) 3 (t 1u ↑) 3 , with the energy term 1 A 1g . The excited state with opposite parity is formed by the one-electron, (t 1u ↑) → (t 2g ↑) or (t 1u ↑) → (e g ↑), excitation, resulting in the lowest excited odd-parity term 1 T 1u at the energy gap 2∆. In this case, the PJTE at the B center, under the condition of instability ∆ < 8F 2 /K 0 (Equation (3)), produces a polar displacement of the B atom along [111]-type directions, which triggers the ferroelectric polarization (Section 3). For other d n configurations of the B ion instead of the d 0 one, the electronic structure of the ground and low-lying states changes drastically, and so does the possibility of formation of close in energy ground and excited states of opposite parity but equal multiplicity. In particular, for the B(d 1 ) ions, the HOMO becomes (t 1u ↓) 3 (t 1u ↑) 3 (t 2g ↑) 1 with the term 2 T 2g , while the LUMO (taking into account Hund's rule) is (t 1u ↓) 2 (t 1u ↑) 3 (t 2g ↑) 2 with the lowest excited odd-parity term 4 T 1u . Hence, the two closest terms of different parity possess different spin multiplicity, and hence they do not mix by the vibronic coupling; the latter does not contain spin operators [9,10]. In principle, higher in energy electronic configurations of opposite parity with the same spin as the ground state one are quite possible, but they are at much larger energy gaps ∆ and therefore less appropriate to satisfy the condition of instability (3) (numerical estimates show that the condition (3) may be very restricting). For the next possible d 2 configurations of the B ion, the two lowest terms of opposite parity are 2 T 2g and 5 T 1u , which, again, do not satisfy the condition for the PJTE dipolar instability.
Moving to the case of B(d 3 ) configuration, we find that the HOMO becomes (t 1u ↓) 3 (t 1u ↑) 3 (t 2g ↑) 3 with the ground state term 4 A 1g , and in the low-spin (LS) conditions of the strong ligand fields (sufficiently large t 2g -e g separation in Figure 2), the LUMO is (t 1u ↓) 2 (t 1u ↑) 3 (t 2g ↑) 3 (t 2g ↓) 1 with the lowest odd-parity term 4 T 1u . It follows that in perovskites with B(d 3 ) ions in LS conditions of sufficiently strong ligand fields, the situation becomes again favorable for the PJTE and polar distortions. However, in this case, distinguished from the d 0 case, the system possesses also a magnetic moment created by three unpaired electrons. However, if the ligand field is weak and the separation t 2g -e g is small, the high-spin (HS) arrangement of the excited electronic configuration takes place, and the excitation electron occupies the e g ↑ orbital instead of t 2g ↓; the LUMO configuration under Hund's rule becomes (t 1u ↓) 2 (t 1u ↑) 3 (t 2g ↑) 3 (e g ↑) 1 with the lowest odd-parity state 6 T 1u . Here, again, there is no PJTE on dipolar distortions and no ferroelectric instability. This is one of the examples which shows explicitly how the spin states interfere directly in the possible local polar displacement and ferroelectricity.
On the other hand, in many cases, the two spin states are close in energy, producing the well-known phenomenon of transition metal spin crossover (SCO), in which case the system can be relatively easily transferred from one spin state to another by external perturbations like heat, light, and magnetic fields [3,4]. As shown above, the change in the spin state changes also the possibility of ferroelectric polarization; hence, the SCO in some perovskite crystals is simultaneously a magnetic-ferroelectric (multiferroic) crossover (MFCO). This coexisting magnetic, ferroelectric, and spin-crossover phenomenon opens a variety of new possibilities to manipulate the properties of the system with novel functionalities to electronics and spintronics [3,4]: (1) For d 3 and d 4 (Cr 3+ , Mn 4+ , Mn 3+ , Fe 4+ , etc.) ferroelectrics in the LS state in conditions of MFCO, magnetic fields facilitate the LS→HS transition that destroys the ferroelectricity (and multiferroicity), while an electric field in the HS non-ferroelectric state may transfer the system to the ferroelectric (multiferroic) LS state, thus realizing electric demagnetization; (2) For d 5 ferroelectrics in conditions of MFCO, if the ferroelectricity is (most probably) different in the two spin states, an electric field may change the spin state (electric magnetization or demagnetization); (3) For d 6 and d 7 (Fe 2+ , Co 3+ , Co 2+ , Ni 3+ , etc.) in the non-ferroelectric LS state under conditions of MFCO, magnetic fields facilitate the LS→HS transition that induces ferroelectricity and hence multiferroicity in a strong magnetoelectric effect (the d 6 LS state is nonmagnetic); in the non-ferroelectric LS state in MFCO conditions, an electric field may transfer the system to the multiferroic state (electric magnetization); (4) The SCO phenomenon is well known to be influenced also by stress, heat, light, and cooperative effects in crystals [46], hence these perturbations can be used to manipulate the MFCO and all the consequent properties including those mentioned above. The dependence of the MFCO on pressure adds a ferroelastic order to the magnetic and ferroelectric ones; (5) There is already a long history of attempts to use transition metal SCO systems as units of magnetic bistability. The difficulty is in the fast relaxation (short lifetime) of the higher in energy spin state [64][65][66]. By choosing a system in the MFCO condition, one can increase the lifetime of the excited dipolar (multiferroic) state by applying an external electric field; (6) An important feature of the revealed MFCO is that it is of local origin and hence it does not necessarily require strong cooperative interactions, meaning that, in principle, it may take place as a magnetic-dipolar effect in separate molecular systems, clusters, thin films, etc., provided that the condition of instability (3) is obeyed.
In addition to the mentioned above examples that confirm the (outlined by the PJTE theory) necessary conditions of multiferroicity, several papers [67][68][69] reported observing also the predicted [3,4] magnetic-ferroelectric spin-crossover effect, while in [69], it is shown that in BiCoO 3 , the ferroelectric polarization is greatly enhanced when the Co 3+ ion is in the high-spin state, as compared to the nonmagnetic state with the Co 3+ ion in the low-spin configuration. They demonstrated also the predicted electric magnetization [3] (see point 3 above), when, by means of induced polarization, the spin state changes from low-spin (S = 0) to high-spin (S = 2). The authors concluded that, contrary to the widespread belief, "unpaired electron spins actually drive ferroelectricity, rather than inhibit it, which represents a shift in the understanding of how ferroelectricity and magnetism interact in perovskite oxides" [69].

Origin of Polar Nanoregions and Relaxor Properties
The PJTE theory recently solved also another problem of ferroelectric crystals, the origin of polar nanoregions (PNR) producing relaxor properties [7], which remained unsolved (inconclusive) in spite of over 60 year of studies (see [70][71][72][73][74] and references therein). PNR in the non-polarized, paraelectric phases are observed arguably in all perovskite ferroelectrics. They are formed spontaneously in the paraelectric phases of ferroelectrics above the Curie temperature T C (where the bulk crystal is cubic and no polarization is expected) in the form of small islands (nanoregions), containing a limited number n of unit cells, their size decreasing with increasing temperature T n -T C . Above the so-called Burns temperature T B , the PNR disappear, and the crystal becomes regular, paraelectric. PNR also disappear under sufficiently strong electric fields. Significantly above T C , the paraelectric phase remains ergodic. With cooling, PNR grow in size, and at a temperature T f (closer to T C , T C < T f < T B ), relaxor properties change again to a non-ergodic, glass-like state, which then undergoes the phase transition to the tetragonal polarized state at T = T C . Distinguished from dipole glasses, this non-ergodic state of the crystal can be irreversibly transformed into a regular polarized state by strong external electric fields. These relaxor properties of ferroelectrics influence all their main properties, with direct applications in materials science.
Attempts to explain the formation of PNR as being due to basic structural disorder, or "random fields" produced by the differences in the active centers (especially in mixed perovskites), as well as by other crystal imperfections (see [70][71][72][73][74]), failed because they do not explain the origin of temperature-dependent size effects-in particular, the disappearance of PNR above T B , or under polarizing electric fields. Moreover, PNR are present in ferroelectric perovskites, but they do not show up in very similar non-ferroelectric crystals with structural phase transitions.
The vibronic (PJTE) theory of ferroelectricity explains also the origin of PNR and relaxor properties of ferroelectric crystals [7] as due to the described above local dynamics of the dipolar distortions at the B centers, which are fully disordered (uncorrelated) in the paraelectric phase and partially ordered (correlated) in the ferroelectric phases. At temperatures above T C , the Helmholtz free energy of the cubic phase Φ cub is lower than its value in the tetragonal phase Φ tetr . Accordingly, as Φ = U -TS, where U is the internal energy and S is the entropy, in the temperature interval T C < T < T B , we have (13) This means that at T > T C , the energy gain ∆U = U cub − U tetr in lowering the potential energy from (average) cubic (in the paraelectric phase) to (average) tetragonal (in the PNR) is not enough to compensate the corresponding entropy loss, T(S cub − S tetr ). Therefore, no displacive theory can explain their formation and specific properties.
Based on the basic features of the vibronic PJTE theory of ferroelectricity, outlined above, the local polarized isles (polar nanoregions, PNR, Figure 9) with a limited, relatively small number n of centers are formed in a non-equilibrium self-assembly process of alignment of the local dipolar displacements [7]. As shown below, at T n > T C , the inequality (13) is compensated by the transformation of a part of the tetragonal potential energy ∆U tetr of the PNR into the work of the formation of its surface W n and transfers heat Q to the environment. The gain of energy in the formation of the n-center PNR ∆U n = U n (cub) − U n (tetr) consists of two contributions: ∆U n (in) , the energy of the internal centers, i.e., the energy of ordering the n dipoles inside the PNR (all oriented, on average, along the polarization direction), and ∆U n (surf) , the energy of the centers in the surface layer (the "domain wall"), which are influenced by the neighbor disordered centers of the cubic phase ( Figure 9). Denoting by n and n' the total number of centers in the PNR and in its border surface, respectively, we get the following estimates for the ordering energies based on the mentioned above mean-field approximation (Section 3): ∆U n (in) = (n − n') ∆U 0 and ∆U n (surf) = n'(∆U 0 − g). Here, ∆U 0 is the per-unit energy gain by tetragonal ordering and g is the loss in this energy by the surface centers because of them being in a lower mean-field, E' < E, of the environment; the mean-field contribution from a part of their environment, by the disordered centers of the bulk cubic phase, is zero. A rough estimate yields g (1/6) ∆U 0 [7]. . Polarized nanoregion inside the cubic perovskite with spherical form (center) and a border layer (dark) at temperatures T above the Curie one TC, but below the Burns temperature T B , T C < T < T B . Arrows indicate the direction of the local averaged dipole moments, which are tetragonally ordered inside the PNR and fully disordered in the cubic phase. For the role of the PNR surface layer, see the text (reproduced with permission from [7]. Copyright 2018, American Physical Society).
For the system in a thermostat under consideration, the exchange of energy with the environment may take place without free energy conservation, Φ i Φ f , the total energy balance being preserved by compensation of heat transfer, Q = T(S f − S i ) = Т∆S, with internal energy change ∆U and mechanical work of internal forces, W intern . Similar to other processes with heat (entropy) transfer, the formation of PNR is a non-equilibrium thermodynamic process, described by Gibbs free energy change, ∆G = ∆H − T∆S, where ∆H is the change in enthalpy H. With growing size of PNR, their Gibbs free energy decreases and reaches a minimum value when n satisfies the condition of thermodynamic equilibrium with the environment. At this point, the shape-restoring force, −dG/dr, is close to zero, while according to the fist law of thermodynamics, ∆U = Q − W intern . This determines the size of the PNR with the number of units n at a given temperature T n . Some relatively simple estimates based on the PJTE-induced APES in perovskite crystals briefly outlined above (Sections 2 and 3) yield the following temperature dependence of the size of PNR [7]: T n T C = 1 + 4.8g ∆U 0 n −1/3 , T n = T C 1 + 4.8g ∆U 0 n −1/3 (14) or by introducing the crystal constant A = 36π , we obtain ( Figure 10) the following: Figure 10 shows also the size of PNR by its diameter D, which is obtained under the assumption of its (most probably) approximately spherical form. The most important cubic (T n -T C ) −3 dependence of the PNR size is confirmed experimentally. Subject to the power "−3", the increase in the size of the PNR with decrease in temperature is very fast, especially at smaller increments ∆T n /T C ≈ 0.1-0.2. This explains the formation of the nonergodic "glassy" state of relaxor ferroelectrics [7]. The estimates leading to Equation (14) show also [7] that, albeit small, an initial trigger (fluctuation) is needed to start the growth of a PNR. It clarifies the role of crystal imperfections or structural irregularities (particularly in mixed perovskites) in enhancing the formation of PNR. Another important point is that the formation of the PNR in time is strongly dependent on the speed of cooling [7].
With the PNR included, we can revisit the polarizability of ferroelectric perovskites. Under an external electric field E, the PNR turns as a whole, like a dipolar molecule, aligning with E, its total dipole moment p ≈ np 0 / √ 6. Applying the Debye-Langevin equation to the ensemble of these "molecules," we find the actual PNR-induced polarizability, α = Nn 2 p 2 0 18a 3 ε 0 kT = NA 2 p 2 0 18a 3 ε 0 kT T T C − 1 −6 (16) Here, N is the concentration of PNR, close to the number of precursors (imperfections) in the crystal lattice that trigger their formation. We neglect the contribution of the applied external field to N. Notably, as follows from Equation (16), the PNR contribution to the polarizability of perovskite relaxor ferroelectrics manifests non-Curie-Weiss behavior: due to the strong temperature dependence of the size of PNR, ferroelectric perovskites have a much stronger temperature dependence of α, proportional to (T -T C ) −6 , which may serve as an indication of the PNR contribution (experimental measured polarizabilities contain both PNR and bulk crystal contributions). Note, however, that the employed above Debye-Langevin equation is not applicable for temperatures near the phase transition where the crystal with PNR becomes non-ergodic and glassy-like, and the formation of the latter depends also on the speed of cooling [7]. For BaTiO 3 , at temperatures above the glassy state, assuming (as above) that p 0 ≈ 7 D, a ≈ 4 Å and T ≈ 450 K, with N of the order of 0.1% and A ≈5 [7], we come to a rather realistic estimate, α(PNR) ≈ 5.6 × 10 4 [71][72][73][74].

Conclusions
The main conclusion that emerges from the outlined above studies of perovskite crystals, notably BaTiO 3 , is that its spontaneous polarization and all relevant properties are triggered by the local PJTE. This conclusion is in drastic contrast to the long dominated "displacive" theories, which are based on the assumption that the spontaneous polarization of the crystal occurs due to the compensation of the local repulsions (in the dipolar distortions) with the long-range dipole-dipole attractions. The predicted by the PJTE local dipolar distortions in BaTiO 3 are multiply confirmed experimentally; they take place in all its phases irrelevant to the phase transitions. The latter occur as a result of the temperature-dependent ordering of the local dipolar displacements. Remarkably, this basic finding in perovskite ferroelectrics affects all their main properties and predicts novel features in interaction with external perturbation, important in applications. The latter includes a qualitatively novel property, the orientational polarization of solids, which (similar to the orientational polarization of polar liquids) is orders of magnitude larger than their displacive polarization. The possible presence of ready-made dipoles in solids which, above "freezing temperature", behave similarly to those in polar liquids was first suggested by P. Debye more than a century ago but revealed only by the vibronic theory. Extension of the PJTE theory to other perovskites and other crystal structures seems to be an up-to-date problem.
Another conclusion from these works is related to the general approaches to the study of solid-state problems. The overwhelming way to explore crystal structures and their properties is to start with revealing the cooperative electronic and nuclear motions characterized by electron and phonon bands formed by equivalent atoms or atomic groups, which are uniformly distributed with translation symmetry. The vibronic (PJTE) theory of ferroelectricity shows (convincingly) that this approach to the problem may be inadequate, losing the main properties of the system. The traditional electron and phonon band approach to solid-state problems fails when the translation symmetry groups are equivalent in average, but not equivalent at each moment in time. This takes place when there are additional local nuclear dynamics induced by the JTE or PJTE, which are uncorrelated (before their cooperative ordering). There are no general theories or computer programs to reveal the band structure of such systems. Our works are the first to reveal this problem and to handle it in the particular case of local PJTE that triggers ferroelectricity and related properties in perovskite crystals.
A short note about our meetings with K. A. Muller and his visits to USSR is available on Supplementary Materials