Theory of relaxor-ferroelectricity

Relaxor-ferroelectrics are fascinating and useful materials, but the mechanism of relaxor-ferroelectricity has been puzzling the scientific community for more than 65 years. Here, a theory of relaxor-ferroelectricity is presented based on 3-dimensional-extended-random-site-Ising-model along with Glauber-dynamics of pseudospins. We propose a new mean-field of pseudospin-strings to solve this kinetic model. The theoretical results show that, with decreasing pseudospin concentration, there are evolutions from normal-ferroelectrics to relaxor-ferroelectrics to paraelectrics, especially indicating by the crossovers from, (a) the sharp to diffuse change at the phase-transition temperature to disappearance in the whole temperature range of order-parameter, and (b) the power-law to Vogel-Fulcher-law to Arrhenius-relation of the average relaxation time. Particularly, the calculated local-order-parameter of the relaxor-ferroelectrics gives the polar-nano-regions appearing far above the diffuse-phase-transition and shows the quasi-fractal characteristic near and below the transition temperature. We also provide a new mechanism of Burns-transformation which stems from not only the polar-nano-regions but also the correlation-function between pseudospins, and put forward a definition of the canonical relaxor-ferroelectrics. The theory accounts for the main facts of relaxor-ferroelectricity, and in addition gives a good quantitative agreement with the experimental results of the order-parameter, specific-heat, high-frequency permittivity, and Burns-transformation of lead magnesium niobate, the canonical relaxor-ferroelectric.

where, in the right side of Eq. 1a, the 1 st -term is the Hamiltonian of the 3-dimensional-random-site-Ising-model (3D-RSIM) [25][26][27][28] , which is a special form of the random interaction Hamiltonian (Eq. 1 of Kleemann et al. 4 ); J the interaction energy constant between the nearest-neighbor PS pairs; σ k the k th -PS on the lattice, and its two states are represented by σ k = ±1; φ the concentration of PS-vacancies or 1 − φ the concentration of PSs on the lattice (φ = 1/3 for PMN and φ = x for BZ x T 1−x ); r k φ the random function that = , and r is a randomly generated number between 0 and 1; {nn} represents the summation of all the nearest-neighbor PS pairs; the 2 nd -and 3 rd -terms are the Hamiltonians related to the RISF and RIEF that originate from the differences in the size and charge of the disorder ions in RFEs, respectively. S k and E k are the random effective Zeeman fields related to the RISF and RIEF 4 . Obviously, E k = 0 in the isovalent RFEs (such as BZ x T 1−x ), but ≠ E 0 k in the heterovalent ones (PMN and Sr x Ba 1−x Nb 2 O 6 ); and N the total number of the unit cells in an RFE. Considering the fact that: (i) Heterovalent and isovalent RFEs have the same characteristics of relaxor-ferroelectricity [31][32][33][34][35][36][39][40][41][48][49][50][51] , so the influence of E k to relaxor-ferroelectricity is the secondary compared with the primary 3D-RSIM as pointed by Kleemann et al. 4 ; and (ii) The similarity of RISF to RIEF 4 , we will use an approximate method to describe their influence (Sec. New mean-field...), and the specific forms of S k and E k are not given here.
Moreover, in order to describe the dynamic parameters of the 3D-ERSIM (such as complex-permittivity), here we use the Glauber-dynamics 29,30 , i.e. the transition probability [w(σ k )] from σ k to −σ k in unit time [Appendix (App.) A of SI] is, where ≡ ∑ σ + + , U B is the energy barrier that PSs stride over during the transition from σ k to −σ k , and v 0 is the orientation vibration frequency of PSs in their local energy valleys; and k B is the Boltzmann constant.
The main reasons for choosing Glauber-dynamics are that: (a) It satisfies the detailed balance condition; (b) The Weiss, i.e. the single-PS, mean-field form of the 3D-RSIM 58,59 along with the Glauber-dynamics (3D-RSIGM) when φ = 0 is the same as the Mason theory that describes the critical-relaxation of 2 nd -order phase-transition 56,57 ; and (c) The corresponding relaxation time predicted by our new mean-field to solve 3D-RSIGM for φ = 0 is consistent with experimental results 57 (Sec. Complex-permittivity…).
New mean-field to solve 3D-ERSIGM. Considering (see SI 3 in details): (a) The existing problems of the approximate theoretical methods for 3D-RSIM [60][61][62] , and especially there is not any feasible method for solving 3D-RSIGM according to the authors' knowledge; (b) The successes of the existing multi-spin mean-field methods for solving Ising-model, and particularly, with increasing the spin number that the mean-fields contain, the corresponding results tend to the exact solution of 2D-Ising-model [63][64][65][66] ; and (c) Inspired by the exact solution of the complex-permittivity of 1D-RSIGM 30 , we propose here a new mean-field of PS-strings (PSSs) that contains more PSs and takes into account the correlation between PSs in 3D-space, referred as PSS-MF, to solve the 3D-ERSIGM. The PSS-MF includes the following four steps: Step-1: PSS construction. There are six kinds of scans, i.e. x-y-z-, y-x-z-, x-z-y-, z-x-y-, y-z-x-, and z-y-x-scans, to construct PSSs for the 3D-RSIM (Fig. 1a-c). For example, the x-y-z-scan is as the follows: (i) Along the x-axis direction of the crystal lattice, connect the nearest-neighbor PSs into short PSSs (Fig. 1b); (ii) Along the y-axis direction, any two nearest-neighbor endpoints of the short PSSs are connected in the x-y plane to form long PSSs (An endpoint already connected to another PSS is no longer reconnected) as indicated in Fig. 1c; and (iii) Along the z-axis direction, continue to connect any two nearest-neighbor endpoints to form longer PSSs in the 3D-lattice. Here, a PSS containing n PSs is expressed as an n-PSS, and the corresponding intra-string Hamiltonian s indicates the i th -PS in the n-PSS, i.e. renumbering the PSs in the 3D-ERSIM Hamiltonian of RFEs (Eq. 1a).
Step-2: Calculation of PSS length distribution. Count the number of n-PSSs in the model (Fig. 1c), obtaining the PSS length distribution function (q n ) versus (vs) n. Due to the spatial isotropy of the 3D-RSIM in Eq. 1a, the calculated q n is the same for the six scans.
Step-3: Calculation of inter-PSS interaction distribution. Count the number of the n-PSSs with the nearest-neighbor number of PSs being g (n-g-PSSs) in the model (Fig. 1d), getting the inter-string interaction (bond) distribution function (ρ n g ) vs n and g.
Step-4: Mean-field of PSSs (App. B, Figs.B1 and B2 of SI for details). The inter-string interaction of an n-g-PSS with its nearest-neighbor PSs, the RISF and RIEF is described by the Weiss-type mean-field 58,59,64 , and the corre- , where b is the factor of the effective interface-effect (due to the interfaces between PS and PS-vacancy groups, the RISF and RIEF), and,  Cyan lines and red solid circles show the crystal lattice and PSs, respectively, while the unlabeled lattice points are PS-vacancies. A blue solid square indicates that its nearest-neighbor two PSs belong to the same string; (c) Connected long PSSs in the x-y plane. They are ring PSSs in the circles of (c); and (d) Three PSSs selected from (c). The red circles, blue squares, and violet diamonds show the PSs, the intra-string and inter-string interaction bonds. PSS-1: = n 11 and = g 11. PSS-2: = n 24 and = g 21. PSS-3: = n 24 and = g 15.
www.nature.com/scientificreports www.nature.com/scientificreports/ Therefore, the total Hamiltonian of n-g-PSSs is,  The normalized condition of q n and ρ n g used here are, respectively, ∑ = = ∞ nq 1 n n 1 and ∑ ρ = = 1 g n n g 0 4 in this paper.
In this article, we use the following computer simulations to show the spatial distribution of PSs and PS-vacancies, and calculate q n and ρ n g . Specifically: (i) Construct an ensemble which includes 10 4 3D-simple-cubic-lattices with 200 × 200 × 200 grid points; and (ii) For any lattice point, a random number r between 0 and 1 is first generated by the Intel-Visual-Fortran 2013 program, and there is a PS-vacancy if r < φ or a PS if r ≥ φ on the point, respectively. Figure 1a shows the distribution of PSs and PS-vacancies in the x-y-plane of a simulated crystal lattice of 3D-RSIM for φ = 1/3. Due to the randomness of the PSs and PS-vacancies on the lattice, there are accumulated regions or clusters of PSs or PS-vacancies with random local structures. As shown in the Sec. Local-order-parameter, it is just these PS clusters that lead to the appearance of the PNRs 15,16,44-47 first proposed by Burns et al. 42,43 , although its definition is quite unclear now as pointed out by Cowley et al. 16 .

PSS length distribution.
Due to the spatial isotropy of the 3D-ERSIM (Eq. 1a), the obtained q n and ρ n g are same by any scan of the six kinds (Sec. New mean-field...). As shown in Fig. 2a, the resulting q n vs n can be described by the following exponential function, n n n 0 / 0 where n 0 is the average length of PSSs, and q 0 is the normalized constant. The probability that a PS belongs to an n-PSS is nq n , and nq n vs n is given in the inset of Fig. 2a. nq n appears as a single peak with n, and the corresponding n value (n p ) of the peak position decreases until n p = 1 as φ increases (Table S4.1 of SI 4).
The calculated n 0 vs φ is illustrated in Fig. 2b, and n 0 decreases with increasing φ. As shown in Fig. 1, n can be used as a size measure of the PS clusters in 3D-RSIM, while n 0 is their characteristic size. In this paper, the ring PSSs in 3D-RSIM are ignored (Fig. 1c). The probability (E R ) that the PSs belong to the ring PSSs is equal to the number of all PSs in the rings divided by the total number of PSs in the model. E R vs φ is shown in the inset of Fig. 2b, and the maximum E R is about 2.7%.
Inter-PSS interaction distribution. The simulated ρ n g vs g of 3D-RSIM for serial φ and n is presented in Fig. 2c-f. There is a single peak of ρ n g vs g for all φ and n, and it could be imagined that ρ n g for φ → 0 is a Dirac δ-function at = 4 g n . By representing the g value corresponding to the maximum of ρ n g as g p n , we could see that: (a) There is a threshold at 0 3 c φ ≈ . , where g p n is nearly irrelevant to n and ≈ . g 2 45 p n (Fig. 2d); and (b) g p n becomes smaller for φ < φ c or larger for φ < φ c as n goes up. In other words, although the sizes of PS clusters in 3D-RSIM of φ ≈ φ c are different, the average coordination PS number per PS in the clusters is nearly same, and this number is equal to + ≈ . . Obviously, this phenomenon is related to the percolation of PS-vacancies if we change the view angle. In fact, φ c is near to the percolation threshold (=0.31) of PS-vacancies in 3D-RSIM 67 , and we consider that they are equal and have the same physical origin in this article. We would like to point out that φ c can be defined as the characteristic concentration of the canonical RFEs (Sec. Definition of canonical RFEs).
Moreover, the ρ n g vs g for n = n p (the magenta diamond symbols in Fig. 2c-f) indicates that there is another threshold at φ ≈ . www.nature.com/scientificreports www.nature.com/scientificreports/ Order-parameter, static-permittivity and specific-heat. When the n-g-PSSs are in thermal equilibrium, the corresponding equilibrium value (η ne g ) of η n g , i.e. the order-parameter of the n-g-PSSs, is, where Z n g is the partition function of n-g-PSSs corresponding to H n g (Eq. 2b), Q n g is an intermediate variable to From Eq. 4a, the static-permittivity (χ s ng ) of the n-g-PSSs in thermal equilibrium is (App. E of SI),  .
., and the inset shows nq n vs n for the corresponding φ values. n p is the peak position of nq n . (b) n 0 vs φ in 3D-RSIM, and the inset illustrates E R vs φ.
(c-f) n n g vs g n / and n in 3D-RSIM with φ = 0.1, 0.3, 0.7, and 0.9. (2020) 10:5060 | https://doi.org/10.1038/s41598-020-61911-5 www.nature.com/scientificreports www.nature.com/scientificreports/ By H n g (Eq. 2b), the average internal energy (u n g ) and specific-heat (c n g ) per PS of the n-g-PSSs in thermal equilibrium are, where ζ ni ge is the equilibrium value of the correlation-function (ζ ni g ) between the i th -and (i + 1) th -PSs, i.e. the expectation value of σ σ + i s i s 1 , in the n-g-PSSs (App. F of SI). Figure 3 shows the η ne g , χ s ng , and c n g vs T for b = 1.5, serial n and g. We obtain that, with decreasing n and for non-zero g, n-g-PSSs have an evolution from 2 nd -order phase-transition to DPT, which is indicated by the diffuse change of η ne g as well as the dispersion peaks of χ s ng and c n g , and the DPT spreads to a wider temperature zone. In this article, we define the temperature corresponding to the maximum value of − η ∂ ∂T ne g as the transition temperature (T p ng ) of the DPT of n-g-PSSs, and T p ng goes up with increasing n and g. ≡ T 0 p ng for n-g-PSSs with g = 0 (n-0-PSSs), giving that they belong to the paraelectric-subsystem in 3D-ERSIM. n-0-PSSs has a small diffuse c n g peak at low-temperature except n = 1 (Fig. 3f). Physically, for finite n and g > 0, the non-zero η ne g value at temperature higher than T p ng originates from that the existence of the PS-vacancy groups, RISF and RIEF makes the probability of the ferroelectric configurations in 3D-ERSIM being higher compared with 3D-Ising-model at high temperature 58,59,64 (App. B of SI).
The order-parameter (η), spontaneous-polarization (P s ), static-permittivity (χ s ps ), as well as the average internal energy (u ps ) and average specific-heat (c ps ) per PS of 3D-ERSIM are, www.nature.com/scientificreports www.nature.com/scientificreports/ In the calculation of macroscopic P s and χ s ps by the microscopic η ne g and χ s ng , this paper uses an approximation similar to the parallel capacitance circuits (Eq. 5b,c).
For serial φ and = . b 1 5, η, χ s ps , and c ps vs T are shown in Fig. 4a-c, and it could be seen that: (i) With varying φ, 3D-ERSIM has the evolutions between the normal-ferroelectrics of (nearly) 2 nd -order phase-transition ↔ RFEs of DPT ↔ paraelectrics of (almost) no PS ordering as indicated by the very small η for φ = .
b 1 5. T p 1 and + T p 2 are the transition temperatures, as well as T d 1 and + T d 2 the diffuse temperatures of n-1-PSSs and n-2 The low-temperature-and high-temperature-DPTs correspond to n-g-PSSs with g = 1 (n-1-PSSs) and n-g-PSSs of g ≥ 2 (n-2 + -PSSs), respectively. The content (R F 1 ), order-parameter (η 1 ), spontaneous-polarization (P s 1 ), static-permittivity (χ s 1 ), and average specific-heat per PS (c ps 1 ) of n-1-PSSs are, , and average specific-heat per PS ( + c ps 2 ) of n-2 + -PSSs are, i.e. the Euler or natural number), as shown in Fig. 4d. 15,16 , see Sec. Burns-transformation) vs φ, which indicates that, as φ increases: (i) + T p 2 first decreases, but it remains almost unchanged after φ ≈ . 0 8; + T d 2 first drops slightly, then rapidly, and keeps as a constant after 0 8 φ ≈ . ; T d 1 first decreases slowly, then rapidly, and increases slightly at the end; as well as + T p 2 is always higher than T p 1 , and T p 1 is almost irrelevant to φ; (ii) R P 0 and + R F 2 , respectively, increases and decreases monotonically, with the maximum growth or drop rate near 0 8 φ = . ; R F 1 shows a diffuse peak with the peak position near φ = .
Local-order-parameter. When the n-g-PSSs are in thermal equilibrium, the corresponding equilibrium value (s ni ge ) of s ni g (App. D of SI) is, shows significant spatial heterogeneity from 0 K to the temperatures far above + T p 2 (Fig. 5b-f) due to the existence of the PS and PS-vacancy clusters (Fig. 5a); (ii) The black regions in Fig. 5b-f that never polarize correspond to = . www.nature.com/scientificreports www.nature.com/scientificreports/ those of PS-vacancies (Fig. 5a); (iii) At very low-temperature ( . + T 0 14 p 2 ), the PS shells adjacent to the PS-vacancy regions have smaller s k e compared with those inside the PS clusters (Fig. 5b); (iv) Even at such high-temperature ( , individual polarized regions of nanoscale about a few a 0 , i.e. PNRs, with small s k e appear in the PS clusters (Fig. 5f), which originates from the effective interface-effect of PS clusters (App. B of SI); and (v) With decreasing T, the size of the PNRs first increases (Fig. 5f,e), then new PNRs appear (Fig. 5e,d), and finally they interconnect to form large polarized regions with quasi-fractal structure (Fig. 5d-b), meanwhile the average of s k e always increases.
Along with the DPT, there are at least three different kinds of interfaces in RFEs: (i) Phase boundary. The interfaces between adjacent paraelectric and ferroelectric regions (Fig. 5c-f); (ii) Sub-phase boundary. The interfaces between adjacent ferroelectric regions of different s k e values (Fig. 5b-d); and (iii) Domain wall. The interfaces in ferroelectric regions with opposite polarization directions of s k e . According to the domain formation theory [68][69][70] , the ferroelectric regions in 3D-ERSIM must become multi-domains due to the influence of the PS-vacancy groups within and adjacent to the PS clusters, the RISF and RIEF.

Complex-permittivity of correlated-relaxation of PSs. The correlated-relaxation of PSs in
normal-ferroelectrics of 2 nd -order phase-transition is also referred to as the critical-relaxation or phase-transition relaxation because it appears in the vicinity of the critical-temperature 56,57 . To date, the most successful theory of the correlated-relaxation is given by Mason 56 , and it is a mean-field of single PS of 3D-RSIGM for 0 φ = . The complex-permittivity (including the linear and higher order) of n-g-PSSs is related to the change of s ni g with time (t). Based on Eq. 2b,d, we obtain that the equation of s ni g vs t is (App. G of SI), . Due to the interaction between PSs, the evolution of s ni g is interrelated to ζ ni The linear complex-permittivity of the correlated-relaxation of n-g-PSSs is directly related to the sufficiently small deviation (δ ni g ) of s ni g from its equilibrium value (  where i c is the imaginary unit, ω angular frequency, and τ nn g the longest relaxation time of the spatial-relaxation-modes of the correlated-relaxation of n-g-PSSs (Figs. H1, H2 and I1 of App. H and I). Figure 6a,b illustrates that, with decreasing n and for nonzero g, ντ nn g shows a crossover from the λ-shape to diffuse peak near T p ng , and from the power-law ντ ∼ − ( ) nn g T T 1 p ng 56,57 to Vogel-Fulcher-law [31][32][33] where T v is Vogel temperature] above T p ng (The relaxation time vs T between the power-law and Arrhenius-relation can be described by the Vogel-Fulcher-law approximately), which indicates that the Vogel-Fulcher-law originates from the effective interface-effect of PS clusters (App. B of SI). ντ nn g with g = 0 always show Arrhenius behavior, and has the same divergent tendency at low-temperature for all g values. Obviously, τ →∞ nn n n 4 corresponds to the case of normal-ferroelectrics with 2 nd -order phase-transition (Fig. 6a), and its power-law behavior is consistent with experimental results 57 . , and = n 10, χ ′ n g and χ ″ n g vs T, serial ω and g calculated by Eq. 11. It could be seen that: (i) χ ′ n g and χ ″ n g have relaxation peaks, and the high-temperature side of χ ′ n g is almost independent of ω; and (ii) The peak temperature (T m ng ) of χ ′ n g vs ω follows the Arrhenius-relation for g = 0 or Vogel-Fulcher-law for ≠ g 0 (Insets of Fig. 6c,d). χ ⁎ n g in 3D-ERSIGM has a distribution with both n and g (Eqs. 11 and 4b). However, there is still no exact method to calculate the complex-permittivity of a heterogeneous system on molecular scale. Here, an extended-Wagner-approximation 71 is used to calculate the complex-permittivity (χ ⁎ ps ) of 3D-ERSIGM, i.e., ps ps c ps n g n n n g n where χ ′ ps and χ ″ ps are the real and imaginary parts of χ ⁎ ps , respectively. Figure 7 shows the calculated χ ′ ps and χ ″ ps of 3D-ERSIGM for = . b 1 5, = U J 20 B , and 0 1, φ = . 0.3, 0.5, 0.7, 0.9 vs T and ω by Eq. 12. The line 5 of Fig. 7a has two peaks of χ ′ ps for φ = . 0 1, which is a typical characteristic of the critical-relaxation 56,57 . Moreover, compared with other φ values, the peak temperature (T m ) of χ ′ ps for 0 1 φ = . changes a little when ω is small, i.e. T m vs ω follows the power-law  (Fig. 4g,h) as the characteristic concentration of the evolution between the normal-ferroelectrics and RFEs.
As indicated obviously in the inset of Fig. 7b, T m vs ω for φ = . . 0 7, 0 9 agrees well with the Arrhenius-relation, and that for φ = .
. 0 3, 0 5 is between the power-law and Arrhenius-relation, which can be described by the (2020) 10:5060 | https://doi.org/10.1038/s41598-020-61911-5 www.nature.com/scientificreports www.nature.com/scientificreports/ approximately, a typical characteristic of RFEs [31][32][33] , which originates from the effective interface-effect of PS clusters (App. B of SI) and the broad interaction distribution between the clusters (Fig. 2d). These results also confirm the rationality that φ p is taken as the characteristic concentration of the evolution between the RFEs and paraelectrics (Fig. 4g,h). Based on Eq. 12 and Fig. 2, it could be expected that the distribution of relaxation time of χ ⁎ ps will first broaden and then narrow with increasing φ, as shown in Fig. 7.

Burns-transformation.
Currently, the interpretation to Burns-transformation of high-temperature thermal-strain (s kl T , = k l , 1, 2, 3) and refractive-index (n kl ) in RFEs is based on the macroscopic quadratic-electro-strictive and Kerr (quadratic-electro-optic) effects 14,42 , and Burns et al. 42 proposed that the transformation stems from the appearance of PNRs at T b being much higher than the DPT 15,16 . However, an unsolved problem is that the calculated spontaneous-polarizations by the above effects are much larger than the data of hysteresis and pyroelectric measurements 14,42 .
Theoretically, the appearance of P s will inevitably lead to the deviations of s kl T and n kl from the high-temperature values. However, this does not rule out the possibility that other factors may also contribute to the transformation. In this paper, according to the coupling between PSs and crystal lattices in 3D-ERSIM, we give a new micro-mechanism of Burns-transformation, which stems from not only the PNRs but also ζ ni ge (App. F of SI) between PSs.       Among them, c kl and d kl are the LI-LD and LI-LE coupling constants; s kl 0 , n kl 0 , α kl and b kl are the thermal-strain, refractive-index, the high-temperature thermal-expansion and thermo-optic coefficients independent of the LI-LD and LI-LE couplings 14,42 , respectively; and T r is a reference temperature.
When the temperature is high enough, → u 0 ps (Eq. 5d and Fig. 4c), . With decreasing T, u ps decreases, resulting in the deviation of s kl T and n kl from their linear behaviors of high-temperatures as shown in Fig. K1 of App. K in SI. Therefore, according to our theory (Eq. 13a,b), Burns-transformation originates from both the PNRs and ζ ni ge between PSs for u ps is related to not only η ne g but also ζ ni ge (Eq. 4c). In addition, the T b increases as c kl and d kl go up (Fig. K1 of App. J-K in SI).
Comparisons with experimental results. At present, the commonly used methods for measuring P s are hysteresis and pyroelectric measurements 14,[34][35][36] . However, for the canonical RFE, PMN, the large external electric field used to polarize the samples will induce the structural phase-transition near 210K 34,36 , which leads to the measured data being not intrinsic. Here, we use the diffuse neutron scattering data (proportional to η) of PMN single crystals 37,38 to compare with the results (Eq. 5a) of 3D-ERSIM, as shown in Fig. 8a. This plot indicates that the model with φ = 1/3, = J 82 K, and = . b 1 5 gives a quite good fit to the data ( = . Moreover, the model results show that PNRs appear well above + T p 2 in PMN (Fig. 8e), and the predicted quasi-fractal characteristic of local-spontaneous-polarization (s k e ) of PMN near and below + T p 2 (Fig. 8f,g) is also consistent with the experimental observations 17,18 .
The experimental data of the specific-heat of PMN, PbMg 1/3 Ta 2/3 O 3 , PbZn 1/3 Nb 2/3 O 3 , and Sr x Ba 1−x NbO 6 show that a specific-heat peak appears in RFEs [39][40][41] , which is also consistent with the model results as shown in Fig. 4c. In Fig. 8b, we give the comparison between the experimental data (c pt ) of the specific-heat of PMN single crystals 39,41 and our theoretical prediction of 3D-ERSIM with 1/3 φ = , = J 102 K, and = . b 1 5 (Eq. 5e). The fit looks reasonable, however the obtained values of J from η and c pt have a relative deviation of ~20%. One possibly direct origination of this deviation is that the preceding c pt data have big measurement error because the data were acquired by subtracting a huge nonlinear background vs T (~100 JK −1 mol −1 at 300 K, and ~25 times larger than the c pt peak height), and this background has some uncertainties, e.g. it was chosen differently by different researchers 39,41 , especially its downturn region with decreasing T overlaps with the c pt peak (Fig. 6 of Tachibana et al. 41 ), which leads to the obtained c pt peak height and temperature of Moriya et al. being ~40% larger and ~7% higher than those of Tachibana et al. 41 , respectively. So, the background selection in PMN is a problem worthy of deep studies. Moreover, we would like to point out that this paper chooses to fit the whole c pt peak of PMN 39 , resulting in the theoretical peak temperature of ~4% higher than that of the experiment, i.e. the corresponding J value may be overestimated by ~4% (Fig. 8b). Another possible origination is the neglection of the mutual Coulomb interaction between PNRs 72 in our theory. E.g., this interaction will induce the adjacent PNRs to tend to the reverse arrangement of spontaneous polarization, leading to the decrease or even disappearance of the local-order-parameter (s k e ) in their adjacent parts (similar to the domain walls in normal-ferroelectrics), i.e. the reduction of η. Moreover, it is worth pointing out that, in addition to s k e , c pt is related to the correlation-function (ζ nk ge ) (Eq. 4c,d). It could be imagined that the above interaction would also result in the decrease of ζ nk ge . As for the quantitative reduction of s k e and ζ nk ge , further researches are needed. Of course, other possible mechanisms may induce the deviation, too.
Due to the spatial distribution of s k e (Fig. 5), the phase boundaries, sub-phase boundaries, and domain walls [68][69][70] appear during the DPT of RFEs, and these movable boundaries and walls will contribute significantly to the complex-permittivity (collectively referred as χ ⁎ b ) below the transition temperature ( + T p 2 ) 15,73-75 . Therefore, the main contributions of the complex-permittivity of RFEs are, respectively, Fig. 8a), and according to the experimental data of complex-permittivity ( c exp ), the frequency (f) is about 1GHz when the peak temperature (T m ) of χ ′ exp is equal to + T p 2 (Bovtun et al. 32 ). Therefore, the χ ⁎ exp of PMN is mainly χ ⁎ ps for f > 1GHz. In Fig. 8c,d, we give the comparison of . The theoretical T 1/ m vs f is shown in the inset of Fig. 8d, and it varies as the Vogel-Fulcher-law. In view of that χ ⁎ exp has large measurement errors when f > 1GHz, especially χ ″ exp 32,76 , the theoretical and experimental results can be considered being consistent with each other. Figure 8h shows the fits of the theoretical predictions (Eq. 13a,b) to experimental data near T b , (i) for n kl : . The good fits illustrate that our theory not only provides a new quantitative micro-mechanism of Burns-transformation, but also avoids the too large spontaneous-polarization data obtained by the quadratic-electro-strictive and Kerr effects (Eqs. J3 and K3 of Apps. J and K in SI that do not consider the contribution of ζ ni ge as shown in Eq. 13a,b) when compared with the hysteresis and pyroelectric measurements 14,42 .
The experimental results show that there is an intermediate temperature (T*) between + T p 2 and T b in RFEs, and according to the current view, T* is the temperature where PNRs change from the high-temperature dynamic to low-temperature static 38,76 . Based on our theory, the relaxation or fluctuation time of the PNRs near T b is much smaller than that near + T d 2 (Insets of Figs. 7b and 8d), and + T d 2 is the characteristic temperature that the PNRs are nearly individual above it (Figs. 5e,f and 8e), while they interconnect to form the quasi-fractal structure of www.nature.com/scientificreports www.nature.com/scientificreports/ local-spontaneous-polarization below it (Figs. 5d-b and 8f,g), so T* just corresponds to + T d 2 as shown in Fig. 8a, and the only difference is their mathematic definitions.
Since x can be continuously adjusted from 0 to 1 49,51 , BZ x T 1−x is an ideal system to verify our theory. However, due to the limitation of crystal growth technology, large-size and high-quality single crystals with x > 0.2 cannot be grown so far 54 , which has affected the measurements of some physical parameters (especially high-frequency permittivity >GHz) to some extent. Nevertheless, with considering the similarity of the low-to high-frequency permittivity of PMN, the analogy of the phase diagram of BZ x T 1−x given by the low frequency (100Hz-500kHz) permittivity 49,51 with our theoretical results has a certain degree of rationality. The obtained phase diagram shows that, with increasing x, BZ x T 1−x evolves from the normal-ferroelectrics to RFEs to paraelectrics (Fig. 11 of Maiti et al. 49 ), while the T m vs ω from the power-law (Fig. 1 of Kleemann et al. 51 ) to Vogel-Fulcher-law to Arrhenius-relation (Fig. 8 of Maiti et al. 49 ), which is consistent with our theoretical predictions as shown in Fig. 4g and the inset of Fig. 7b.
It is worth pointing out that the results of the existing Monte-Carlo-simulations of 2D and 3D-RSIM 25-28 support our PSS-MF method as shown in the SI 6.

Definition of canonical RFEs. PMN is generally viewed as the canonical RFE, but what is a canonical RFE
is not well defined so far 52,53 . The experimental results show that the systems with definite components but adjustable ion distributions, such as PbSc 1/2 Ta 1/2 O 3 , continuously evolve from normal-ferroelectrics to RFEs with the increase of disorder 14 , so it looks reasonable that the random distribution of ions, i.e. the most disordered case, is one necessary condition of canonical RFEs. Clearly, the 3D-ERSIM (Eq. 1a) proposed in this paper satisfies this condition.

Discussion
By the fitting of our theoretical results to the experimental data (Fig. 8), the obtained φ = 1/3 is a reasonable value for PMN. In addition, we can conclude that PbZn 1/3 Nb 2/3 O 3 and PbMg 1/3 Ta 2/3 O 3 also correspond to 3D-ERSIGM of φ = 1/3. Moreover, it could be imagined that, for BZ x T 1−x , J, μ, and U B are almost irrelevant to x when x is small. However, due to the RISF and when x is large enough, J, μ, and U B will change with increasing x, which needs deep studies. For the more complicated and anisotropic RFE, Sr x Ba 1−x Nb 2 O 6 , the relationship of φ, J, μ, and U B with x requires future researches, too.
It could be expected that the domain structure, including domain walls, can be obtained based on the calculated s k e (Figs. 5b-f and 8e-g) and the theory of domain formation [68][69][70] . Moreover, the method for calculating χ ⁎ ps of single-PS flipping in this paper (Sec. Complex-permittivity…) would lay the foundation to obtain χ ⁎ b of the overall movement of PSs, i.e. the multi-PS flipping, in the phase boundaries and domain walls [73][74][75] . Of course, they are the issues that need further studies. We would like to point out that the mutual Coulomb interaction of PNRs at T < T* has recently been evidenced by Kleemann and Dec 72 , to give rise to a super-dipolar glass state below the glass transition temperature (=240 K) in PMN, and the relevant research of χ ⁎ b will undoubtedly deepen the understanding of this state.
It would be speculated that the approximation to calculate P s (Eq. 5b), χ s ps (Eq. 5c), and χ ⁎ ps (Eq. 12) has large errors when φ is high, so it is a potential work to explore more accurate calculation method. Other future studies are: (i) To generalize the PSS-MF for solving isotropic 3D-ERSIGM to anisotropic cases, so that the possible anisotropic RFEs, such as the single-axis tungsten bronze and layered Aurivillius structures 17 , can be described; (ii) Coupling of PSs with crystal lattices and the acoustic properties 16  and (iv) To generalize our theory to the corresponding ferromagnetic systems, in particular spin-glasses, etc. [20][21][22] .
As shown in Figs. 1a, 4, 5, 7 and 8, our theory predicts that RFEs are a special kind of ferroelectrics (observable macro-spontaneous-polarization) with DPT originating from the spatial-dynamical heterogeneity as indicated by the spatial variation of local-spontaneous-polarization [local-order-parameter (s k e )] including the coexistence of para-and ferroelectric regions (Fig. 5c) and corresponding distribution of relaxation time (Figs. 6a,b and 7c-f), which leads to s k e being the key parameter of relaxor-ferroelectricity, instead of the macro-spontaneous-polarization [order-parameter (η)]. In this respect, normal-ferroelectrics are only special ferroelectrics where s k e is equal everywhere. It is worth noting that at low temperature (T → 0), there are macro-spontaneous-polarization regions throughout RFEs (Fig. 5b) according to the percolation theory 67  www.nature.com/scientificreports www.nature.com/scientificreports/ walls 17,18 in RFEs as described in Sec. Local-order-parameter, and by analogy with the experimental results in normal-ferroelectrics 73,74 , the relaxation time of their lateral movement should increase with decreasing temperature and tends to infinity when T → 0, resulting in the freezing or glass transition of domain wall movement at a certain temperature 3,[72][73][74] . In short, based on the present theory, RFEs have spatially-dynamically heterogeneous local-order-parameter and the corresponding mesoscopic defects (domain walls and phase boundaries etc.) of the local-order-parameter show glassy behavior. Of course, the correctness of this relaxor-ferroelectricity picture still needs further theoretical and experimental studies.
The authors would like to point out that, in the preceding theories or models (SI 1), the spherical random-field random-bond model of Pirc et al. 7,8 and the soft-mode theory with random-electric-field of Arce-Gamboa and Guzmán-Verrí 10 give the most quantitative predictions. The theory of Arce-Gamboa and Guzmán-Verrí gives that: (1) For weak random-electric-field, long-range ferroelectric order sets in at a transition temperature where the obtained η changes discontinuously and it is accompanied by metastable paraelectric or random-electric-field state down to T → 0; (2) For moderate one, there is no transition as the paraelectric state becomes stable at all temperatures and the long-ranged polar state is now metastable; and (3) For strong one, only the paraelectric state exists. According this theory, PMN belongs to the case of (1) 10 , so both the theoretically discontinuous change of η with T in the absence of an external electric field and no PNR above the transition are inconsistent with the experimental results 18,37,38 , which is obviously different to the predictions of the present theory (Fig. 8a,e). The physical origination of these differences is that, according to our theory, there are dipole-or PS-vacancies randomly distributing in RFEs (Eq. 1a and Fig. 1a) so that the interfaces between PS and PS-vacancy groups appear, which leads to the dispersion of the phase transition 64 , i.e. interface-effect (App. B), but similar interface-effect term does not exist in the model Hamiltonian as shown in Eq. 1 of Arce-Gamboa and Guzmán-Verrí theory 10 .
The spherical random-field random-bond model predicts that the DPT corresponds a dipolar glass transition, however, in view of the elementary motion and interaction unit being the effective dipole, i.e. PNR, this model is at the mesoscopic or semi-microscopic level and the glass transition might correspond to the super-dipolar glass transition of PNRs 72 , which is an issue of the present theory to be further explored as mentioned in the 2 nd paragraph of this section.