Random Electric Field Instabilities of Relaxor Ferroelectrics

Relaxor ferroelectrics are complex oxide materials which are rather unique to study the effects of compositional disorder on phase transitions. Here, we study the effects of quenched cubic random electric fields on the lattice instabilities that lead to a ferroelectric transition and show that, within a microscopic model and a statistical mechanical solution, even weak compositional disorder can prohibit the development of long-range order and that a random field state with anisotropic and power-law correlations of polarization emerges from the combined effect of their characteristic dipole forces and their inherent charge disorder. We compare and reproduce several key experimental observations in the well- studied relaxor PbMg$_{1/3}$Nb$_{2/3}$O$_3$-PbTiO$_3$.


INTRODUCTION
Relaxors exhibit a myriad of complex phenomena that are both scientifically interesting and technologically important such as diffuse phase transitions where large and frequency dependent dielectric permittivities extend over hundreds of Kelvin degrees [1] without any signature of macroscopic symmetry breaking as well as unltrahigh electromechanical responses [2][3][4][5]. These properties make relaxors attractive material candidates for energy storage and harvesting applications as well as future cooling technologies for integrated microelectronics [6][7][8].
Though relaxors were first synthesised in the 1950s [9,10] and they have been the subject of many theoretical [11][12][13][14][15] and experimental studies [16][17][18][19][20][21][22][23] there is still no consensus on a satisfactory theory of relaxor ferroelectricity [24]. One of the major difficulties in describing relaxors is that they exhibit many characteristic temperatures. From high to low, these are (i) the Burns temperature T B below which its dielectric response deviates from Curie-Weiss law behavior with (ii) a corresponding Curie-Weiss temperature T CW ; (iii) a frequency dependent temperature T max where the susceptibility is maximum but no ferroelectric (FE) transition occurs; and (iv) an induced FE transition temperature T c if sufficiently large electric fields are applied. Crucially, X-ray and neutron scattering studies have found anisotropic quasi-elastic diffuse scattering very near T CW [25][26][27][28][29].
It has been recognized that a central question in the discussion of relaxors is the effect of random electric fields on the FE transition of cubic systems such as the typical perovskite relaxor PbMg 1/3 Nb 2/3 O 3 (PMN) [18,24,30]. The random electric fields originate from charge disorder: cations with different charge valencies are randomly located on the octahedrally coordinated site such as Mg 2 + and Nb 5+ in PMN [30]. These ions do not order with temperature, making the compositional disorder quenched. Unlike the widely studied random fields (RFs) in magnets which linearly couple to an order parameter of the Ising or Heisenberg type [31], the quenched electric RFs of relaxors couple to a cubic order parameter [24]. It is believed that T * is the onset temperature of a RF state in which relaxors exist [24].
In addition to the symmetry of the order parameter, we make the observation that the characteristic dipolar interaction of FEs is equally important. It is well-known that the structural instability that leads to the breaking of lattice inversion symmetry and a spontaneous polarization, is the result of dipolar forces between electric dipole moments induced by the displacements of the ions associated with a zone-center transverse optic (TO) mode [32]. Such dipolar forces are highly anisotropic and long-ranged, which are very much in contrast with the isotropic and short-ranged exchange couplings between the spin degrees of freedom of magnets.
According to the theory of phase transitions [33], FEs and magnets are therefore in different universality classes, rendering the standard models that describe the effects of RFs on magnetic transitions [31] inadequate for relaxors [24].
In a previous paper, [38] we studied the effects on quenched electric RFs in a standard, uniaxial displacive model of the FE transition. Within a statistical mechanical variational solution, we showed that intrinsic polarization fluctuations associated with the dipolar force and RF disorder, result in diffuse phase transitions -a hallmark of relaxor behavior. Typical relaxors such as PMN are cubic, however, and there is no a-priori reason to believe that the results for uniaxial systems will hold in environments with higher symmetries. The purpose of this work is then to study the random electric field problem posed by cubic relaxors within a minimal microscopic model.
We extend the uniaxial model Hamiltonian of Ref. [38] to cubic symmetries by including the usual displacement soft-mode coordinates along each cubic axis, cubic anisotropy, dipole tensor, and cubic RFs. We also extend to cubic symmetries our previously developed variational solution for uniaxial systems. We will show that as a result of the combined effect of dipolar forces and quenched RFs a state with no-long range FE order and anisotropic, long-ranged fluctuations of polarization emerges for any amount of compositional disorder. We identify this disordered state as the RF state of relaxors. We will also show that long-ranged FE order can be induced by application of strong enough electric fields and that such transition ends at a critical point, as it is observed in experiments [34].

RESULTS
We consider a cubic lattice and choose normal mode coordinates that describe local displacements (Q ix , Q iy , Q iz ) in the unit cell i that are associated with the soft TO mode, the condensation of which leads to the FE transition [35]. We consider the model Hamiltonian, with λ, λ = x, y, z. Π iλ is the conjugate momentum of Q iλ ; E 0 λ is an applied electric field; and v λλ ij is the dipolar interaction tensor with Fourier transform v λλ q = 1 3 C 2 − B 2 |q| 2 δ λλ − C 2 q λ q λ |q| 2 , where |q| = q 2 x + q 2 y + q 2 z is the magnitude of the wavevector q; and B and C are constants that depend on the lattice structure [36]. Hereafter, we denote v 0 = C 2 /3 as the component of v λλ q when q → 0 in the direction transverse to λ (the value of v q depends on the direction in which q approaches zero). κ is the lattice stiffness and γ 1,2 are anharmonic coefficients.
For the quenched random fields h iλ , we choose a Gaussian probability distribution of independent random variables with zero mean and variance ∆ 2 . In the absence of compositional disorder, this is a standard minimal model for ferroelectricity in cubic perovskites [35].
To study the statistical mechanics of the Hamiltonian (1), it is necessary to consider thermal and quantum fluctuations at least at the level of the Onsager approximation and random field fluctuations at least at the level of a replica theory [37]. To do so, we generalize a variational method previously developed by one of us [38] to cubic symmetries. Such method allow us to calculate the temperature and disorder dependence of relevant quantities such as the phonon frequencies, the polarization order parameter and the correlation functions in a self-consistent fashion. The details are presented in Methods section.
Our model parameters are κ, γ 1 , γ 2 , v 0 , B and ∆. Throughout this work, we have fixed the values of κ, γ 1 , γ 2 , v 0 , and B to those of typical values of oxide perovskites [40] and to fit the transition temperature of the conventional FE PbTiO 3 (PTO, T 0 c 760 K ) [41] assuming ∆ = 0. The resulting values are given in Table I. Our choice gives a Curie-Weiss constant of C 0 CW 2.4 × 10 5 K and a zone-center TO phonon energy of Ω ⊥ 0 5.2 meV at zero temperature, which are typical of conventional FEs. Depending on the choice of the anharmonic coefficients, the low temperature FE phase predicted by the Hamiltonian (1) in the absence of compositional disorder has tetragonal (γ 1 < γ 2 ) or rhombohedral (γ 1 > γ 2 ) symmetry [35]. In this work, we have chosen γ 1 > γ 2 , as we will study the field-induced FE transition of relaxors, which is typically a cubic-to-rhombohedral structural phase change [20].
We first present our results in the absence of applied electric fields. Figure 1 shows the calculated temperature-disorder phase diagram and the zero temperature free energies where we have identified three regions according to the RF strength. For weak RFs long-range FE order sets in at a transition temperature T c < T 0 c and it is accompanied by a metastable random field disordered state down to T = 0, as it is shown in Fig. 1 (a). For moderate compositional disorder (0.9 ∆ 2 /v 3/2 0 2.2), there is no transition as the the RF state becomes stable at all temperatures and the long-ranged polar state is now metastable, see side of the morphotropic phase boundary are stuck in a metastable disordered random field state below the phase transition line. We will see below that this is also supported by the predicted correlation lengths and static susceptibilities of our model.
The temperature and disorder dependence of the zone-center TO phonon frequency, Ω ⊥ 0 , and the order parameter, A, associated with the RF and FE states are shown in Figure 2 (a) and (b), respectively. While the TO mode softens and condenses at T 0 c for the pure case, as expected, that of the RF state remains finite all the way down to zero temperature for any amount of compositional disorder. The temperature dependence of the metastable states is shown for the sake of completeness. When contrasted to experiments [21], the observed softening of the phonon frequency of the RF state above about T * is in qualitative agreement with our model and we will show that it supports the conclusion that such softening is responsible for the large increase observed in the dielectric constant. Below T * , however, the observed frequencies exhibit a more complex behavior not captured by our model. We believe that some of the discrepancies are due to local spontaneous polarizations in the disordered state, which we do not allow in our model.
We now discuss the correlation functions in the fluctuations of polarization of the RF state.
As it is usually done for conventional cubic FEs [42], we consider mean squared fluctuations on the polarization components Q ⊥ q and Q q that are transverse and longitudinal to a wave-vector q, respectively. For the transverse components we obtain isotropic fluctuations with the following form, where Ω ⊥ q 2 = B 2 ξ −2 + |q| 2 is the doubly degenerate TO mode (see Methods section) and where we have identified ξ = B/Ω ⊥ 0 as the correlation length. . . . denotes thermal and compositional averages taken in that order. In the absence of disorder and in the classical limit (βΩ qλ 1), Eq. (2) reproduces to the fluctuations of pure FEs [42]. In the classical limit, Eq. (2) becomes a Lorentzian plus a Lorentzian squared. While this is analogous to the well-known result of the random field Ising model [43], we will show below that the correlation functions behave very differently in real space due to the anisotropy and long-range nature of the dipole force. The wavevector distribution predicted by Eq.
(2) has been recently observed in diffuse scattering experiments [19] and the quantum fluctuations have been found important to correctly describe the observed static structure factor at low temperatures [37]. For the fluctuations in Q q , we find that they have a similar form to that of Eq. (2) except that the TO frequency is replaced by that of the the longitudinal mode Ω q 2 ∝ ξ −2 + |q| 2 + C 2 . The constant C is related to the depolarizing field [42], which makes these fluctuations significantly smaller than those in Q ⊥ q . Figure 3 shows the calculated temperature dependence of the correlation length ξ of the RF state for several disorder strengths. In the absence of disorder, ξ diverges as expected near the FE transition. In the presence of disorder, the correlation length of the random field state remains finite at all temperatures. At T = 0 it scales with disorder as ξ ∝ 1/∆ 2 , which we identify as the minimum length scale on which domains must appear spontaneously. By a standard procedure [42], it can be shown that the static dielectric susceptibility is given by χ = 3 4π v 0 B 2 ξ 2 . The temperature and disorder dependence of the resulting static dielectric constant = 1 + 4πχ, are shown in the inset of Fig. 3. We find that our model is in fair qualitative and quantitative agreement with the measured correlation length [26] and static dielectric constant [44] in PMN when 0.5 This is consistent with our identification of the Ti-poor region of PMN-PT in our phase diagram (see Fig. 1) and where the RF state is metastable. We now discuss the spatial dependence of the correlation functions of polarization G λλ (r). Figure 4 shows G λλ (r) calculated from Fourier to the correlation length ξ, they fall-off exponentially and then cross over to a power law behavior (∝ r −3 ) for r ξ, see inset in Fig. 4 (a). We verify the large distance behavior by calculating  analytic expressions of the correlation functions of the random field state in the classical limit, where λ = x, y, z, and v BZ is the volume of the Brillouin zone. Note that the corrections to these power laws are exponentially small. The cross-component correlations (λ = λ ) of the random field state are identically zero everywhere, as expected from cubic symmetry.
Note they also increase with decreasing temperature but do not reach long-range order as their correlation length ξ remains finite for all temperatures. This is in stark contrast with the correlations of the pure compound where, while anisotropic, they are strongest near the FE transition and then weaken away from it, as it is shown in Figs. 4 (c) and 4 (d). Previous theoretical work have also found anisotropic correlations [14,15].
We now describe our results in the presence of an applied electric field.

DISCUSSION
We now compare our results to previous theoretical works. When contrasted to uniaxial systems [37,38], we find that they share some similarities at the qualitative level such as the where heuristic arguments are given to conclude that relaxor behavior in heterovalent compounds is mainly due to bond disorder and that RFs only play a secondary role. On the other hand, our results support the view of Takenaka et al. [14] that there is no non-polar matrix in relaxors (our correlations decay as power-laws for r ξ); and that of Al-Barakaty et al. [15] that quenched RF disorder is essential for relaxor behavior. We emphasize, though, that, according to our results, the intrinsic fluctuations associated with the concomitant dipole forces are essential as well.
To summarize, we have studied the effects of cubic random electric fields on the lattice instabilities that lead to a FE transition. We have shown that a RF state emerges from the

METHODS
We now describe our variational solution of our model Hamiltonian of Eq. (1). We consider the trial probability distribution, where H tr is the Hamiltonian of a displaced, cubic harmonic oscillator in a random field h iλ , and Z tr its normalization, where Ω qα (α = 1, 2, 3) are the soft mode frequencies at wave vector q and are given by the squared of the the eigenvalues of the Fourier transform of the dynamical matrix A iλ is the λ-component of the order parameter at site i and it corresponds to the mean displacement averaged over thermal and compositional disorder, We now compute the free energy F = H + T k B ln ρ tr using our probability distribution (3) together with the above equations. The result is the following, where, ψ λλ ij are temperature and disorder induced fluctuations of polarization between local soft mode components Q iλ and Q jλ , with Fourier component, (b q ) λλ is a unitary transformation that takes (D q ) λλ to its diagonal representation.
A standard procedure gives the following dynamical matrix, where ψ αν 0 is given in Eq. (5). The diagonalization of (D q ) αν gives the squared of the soft phonon frequencies (Ω qλ ) 2 .
Minimization of the free energy (4) with respect to A iλ gives the following result, In writing Eq. (7) we have used the property that v αλ ij is translationally invariant so the summation iλ v αλ ij does not depend on the origin i. (D q=0 ) αν depends on the direction in which q → 0 because v αλ q is non-analytic. Eqs. (6) and (7) are the starting point of our analysis.
We first consider the cubic phase. For the cubic phase, there is no long-range order (A x = A y = A z = 0) and, by symmetry, ψ 0 ≡ ψ xx 0 = ψ yy 0 = ψ zz 0 , ψ xy 0 = ψ xz 0 = ψ yz 0 = 0 [35]. Therefore, the dynamical matrix has the form, For an arbitrary direction of q the diagonalization of (D q ) λλ gives a doubly degenerate transverse optic (TO) mode Ω ⊥ q , and a singlet longitudinal optic (LO) mode Ω q , given as follows, where, is the zone-center TO mode frequency and ω 0 ≡ √ v 0 − κ. For cubic symmetry, the transformation matrix b q that diagonalizes the dynamical matrix (8) takes the form, where θ and φ are the usual azimuthal and polar angles in spherical coordinates.
We now calculate the fluctuations ψ 0 , where s φ ≡ sin φ, c φ ≡ cos φ, s θ ≡ sin θ, c θ ≡ cos θ and, By taking the continuum limit over a sphere of wave-vector Q and calculating the angular integrals, we find the result, Thus, Equations (9), (10), (12), (13) determine de temperature and disorder dependence of the TO and LO mode frequencies for the cubic phase.
We now consider the rhombohedral phase. For the rhombohedral phase, we assume a homogenous order parameter along the cube diagonal, Also, by symmetry, [35]. The dynamical matrix is as follows, where α, ν = x, y, z. We first identify the soft mode frequencies. Pure longitudinal and transverse modes are obtained for wavevectors in the (111) direction and the plane transverse to it. For q ⊥ (1, 1, 1) diagonalization of the dynamical matrix gives two distinct TO mode frequencies Ω ⊥

01
and Ω ⊥ 03 and one LO frequency Ω ⊥ 01 + 3(C 2 /3) at the zone-center. For q (1, 1, 1) there is a doubly degenerate TO mode frequency Ω ⊥ 01 and one LO mode frequency Ω ⊥ 03 + 3(C 2 /3) at the zone-center. Ω ⊥ 01 and Ω ⊥ 03 are given as follows, While exact expressions can be derived for the phonon dispersions from the dynamical matrix (14), they are too elaborated and not enlightening. Instead we calculate them from perturbation theory. Our unperturbed basis is that of the cubic phase, therefore making the frequency splitting Ω ⊥ 03 − Ω ⊥ 01 the expansion parameter. It is also convenient to write the wavevector as q = q LqL + q T 1qT 1 + q T 2qT 2 , where {q L ,q T 1 ,q T 2 } is a right-handed coordinate system whereq L is along the where q 2 T = q 2 T 1 + q 2 T 2 and with a transformation matrix given by, For an applied field (E 0 / √ 3)(1, 1, 1), minimization of the free energy with respecto to the order parameter gives the following result, which can be rewritten in terms of the soft mode frequencies as follows, We now calculate ψ xx 0 and ψ xy 0 , To proceed further, we the continuum limit as we did in the previous section and calculate the integrals over φ. The result is the following, Equations (15)- (19) determine de temperature and disorder dependence of the order parameter A and the TO and LO mode frequencies in the rhombohedral phase.