---
abstract: 'Ferroelectric relaxors are complex materials with distinct properties. The understanding of their dielectric susceptibility, which strongly depends on both temperature and probing frequency, have interested researchers for many years. Here we report a macroscopic and phenomenological approach based on statistical modeling to investigate and better understand how the dielectric response of a relaxor depends on temperature. Employing the Maxwell-Boltzmann distribution and considering temperature dependent dipolar orientational polarizability, we propose a minimum statistical model and specific equations to understand and fit numerical and experimental dielectric responses versus temperature. We show that the proposed formula can successfully fit the dielectric response of typical relaxors, including Ba(Zr,Ti)O$_{3}$, Pb(Zn$_{1/3}$Nb$_{2/3}$)$_{0.87}$Ti$_{0.13}$O$_{3}$, and Pb(Mg$_{1/3}$Nb$_{2/3}$)O$_{3}$-0.05Pb(Zr$_{0.53}$Ti$_{0.47}$)O$_{3}$, which demonstrates the general applicability of this approach.'
author:
- 'J. Liu'
- 'F. Li'
- 'Y. Zeng'
- 'Z. Jiang'
- 'D. Wang'
- 'Z.-G. Ye'
- 'C.-L. Jia'
title: Insight into dielectric response of ferroelectric relaxors by statistical modeling
---

Introduction
============

Relaxor ferroelectrics are materials that exhibit interesting dielectric responses different from normal ferroelectrics. For instance, they often possess relaxation modes at low frequency ($<1$GHz). Relaxors ferroelectrics have been exploited in many applications such as actuators due to their giant electromechanical couplings [@Uchino1996], and their properties extensively investigated, including structural properties (e.g., polar nanoregions or PNRs) using neutron scattering [@PhononLocalization], dielectric responses [@Bokov2006a; @Nuzhnyy2012; @Petzelt2014; @Kleeman2014; @Wang2014; @Wang2016], the crossover from ferroelectrics to relaxor [@Kleeman2014]. To understand such systems, many theories have been proposed [@Pirc1999; @Bokov2012; @Akbarzadeh2012; @Sherrington2013a; @Uchino2014; @Sherrington2014; @PhononLocalization; @Kleemann2015]. Relaxors are complex systems, to some extent similar to spin glasses [@Sherrington2013a; @Sherrington2014a; @Sherrington2014], in that their compositions are, without exception, made of complex oxides containing different ions and inevitably inhomogeneous. For instance, the B-site ions of typical relaxors Ba(Zr,Ti)O$_{3}$ (BZT) and Pb(Mg$_{1/3}$Nb$_{2/3}$)O$_{3}$(PMN) are randomly distributed.

The dielectric response of ferroelectric relaxors is the defining feature that differentiates them from normal ferroelectrics: (i) large susceptibilities at low frequency (GHz or lower) ; (ii) even more unusual, the characteristic temperature $T_{m}$, at which the susceptibility peaks strongly depends on the frequency of the probing *ac* electric field. In other words, the susceptibility, $\chi$, depends on both temperature, $T$, and the probing frequency, $\nu$. While such phenomena are well known experimentally [@Bokov2007BZT; @Bokov2012; @Gridnev2004; @Bokovo2006; @BZT-susceptibility; @Tagantsev1999; @PMN-BZT], numerical generation of relaxor’s dielectric response from model-based simulations has been a challenging work. For instance, the shift of $T_{m}$ of the lead free relaxor BZT was only numerically achieved recently[@Wang2016]. Since numerous ferroelectric relaxors exist, numerically treating each of them remains a daunting task. One way to mitigate this difficulty is to resort to statistical modeling [@Sethna2010]. For a complex system, a statistical approach can provide intuitive understanding by capturing dominant factors, provide equations to understand experimental results. and help extracting useful information. In the present work, we adopt this approach to treat the dielectric response of relaxors and show that such a statistical model can indeed be applied to understanding how the dielectric constants change with temperatures and probing frequency.

While the susceptibility of relaxors, $\chi\left(T,\nu\right)$, depends on both temperature and frequency, theoretical models are often proposed to treat $\nu$ and $T$ separately [@Uchino2014; @Bokov2012; @Cross1987a; @UniversalRelaxation; @DielRelaxSolids; @Nuzhnyy2012]. For instance, at a given temperature, two processes are employed in the fitting of $\chi\left(\nu\right)$ of Ba (Ti$_{0.675}$Zr$_{0.325}$ )O$_{3}$: the universal relaxor process and the conventional relaxor process, which have different relaxation characteristics employing the Curie-Von Schweidler law at low frequency and the Kohlrausch-Williams-Watts law at higher frequency [@Bokov2007BZT]. Other formula, such as the Cole-Cole and the Havriliak-Negami equations, are also employed to model the dielectric response with respect to frequency at given temperatures. When phonon modes are close or interacting with the relaxation modes, it becomes necessary to use coupled modes to model the dielectric response [@Wang2014; @Wang2016]. On the other hand, there are also many investigations on how the dielectric response, $\chi$, depends on the temperature, $T$, at given frequencies. In addition to the well known Curie law for $\chi\left(T\right)$ at high temperature, most useful equation for fitting around the dielectric peak appears to be the square law. In particular, the formula $$\begin{aligned}
\frac{1}{\varepsilon\left(T\right)}= & \frac{1}{\varepsilon_{A}}+\frac{\left(T-T_{A}\right)^{^{\eta}}}{B}\label{eq:square-law}\end{aligned}$$ was proposed to describe the permittivity at $T>T_{m}$ [@Smolenskii1970; @Kirillov1973]. Initially, $\eta$ was found to be 2, but later was shown to be between 1 and 2 [@Clarke1974; @Uchino1982; @Santos2001; @Bokov2006a; @Bokov2012]. Here, we further the investigation in this direction and attempt to address some important questions regarding relaxor behaviors. We will explain why the dielectric constant has a peak value at $T_{m}$, and what causes the asymmetry around the peak. Moreover, by constructing a statistical model that properly describes how dipoles behave in relaxors, we propose formula to fit experimental results, which further illuminates the physics behind relaxation behavior.

This paper is organized as the follows. In Sec. \[sec:Statistical-modeling\], we introduce the statistical model. In Sec. \[sec:Results\], we apply this model to both lead-free and lead-based relaxors. In Sec. \[subsec:Discussion\], we discuss the implication and limitation of this approach. Finally, in Sec. \[sec:Conclusion\], we present a brief conclusion.

Statistical modeling\[sec:Statistical-modeling\]
================================================

The statistical model starts by considering a critical difference between ferroelectric relaxors and normal ferroelectrics. One crucial observation is that all relaxor ferroelectrics discovered so far are inevitably disordered and inhomogeneous systems. For instance, in BZT Zr and Ti ions are distributed randomly, so are the Mg and Nb ions in PMN, when the samples are treated macroscopically. In addition, PMN possesses the electric field arising from heterovalent Mg and Nb ions, which affects dipole distribution. It is important to further note that well known relaxors can become non-relaxor if their ions are perfectly ordered [@Setter1980; @Bokov1984; @WangJAD].

Individual dipoles
------------------

The randomness of ions and the ensuing lack of long-range correlation has the important consequence that phonon modes may not be the best description to understand relaxor. This fact is evidenced by the effective Hamiltonian that describes the BZT relaxor [@Akbarzadeh2012; @Sherrington2013a; @Sherrington2014] $$\begin{aligned}
E= & \sum_{i}\left(\kappa_{i}\left|u_{i}\right|^{2}+\lambda_{i}\left|u_{i}\right|^{4}\right)+\dots,\label{eq:effHBZT}\end{aligned}$$ where $i$ labels the sites occupied randomly by Zr or Ti, and $\kappa_{i}$ ($\lambda_{i}$) are the second (fourth) order coefficients in the Taylor expansion of energy with respect to $u_{i}$, which is the local dipole on site $i$. For a homogeneous system, where $k_{i}$ and $\lambda_{i}$ are constants, we can usually first consider the harmonic term and construct phonon modes, which are then used to describe the system, especially in low temperature when the system condense to particular phonon modes [@CondenedPhononMode]. In contrast, with the loss of periodicity in relaxors, this approach is no longer profitable. One can insist on using averaged atoms (e.g., replacing Zr and Ti atoms with their average in BZT) to retain the use of phonon modes. However, it is then necessary to consider defect-pinned intrinsic localized modes [@PinnedPhonon] and phonon localization [@PhononLocalization].

The inhomogeneity also has important consequences on ferroelectric phase transitions. In the typical ferroelectric material BaTiO$_{3}$, we may ascribe the temperature-driven phase transition to the condense of phonons to a particular phonon mode [@CondenedPhononMode]. At high temperature, many phonons modes are occupied (occupancy obeying the Bose-Einstein distribution); at low temperatures, due to mode softening, certain mode (often corresponding to the well-known in perovskites) has essentially zero energy, which dominates the system and induce phase transitions. Unlike BaTiO$_{3}$, there is no global phase transition due to the existence of PNRs and/or random electric fields, which eliminates the mode softening phenomenon and renders a global dipole pattern difficult to achieve [@Akbarzadeh2012; @PMNeffHami2015]. In addition, relaxors exhibit strange phonon behavior, such as the “waterfall” effect [@Gehring2000; @Hlinka2003; @PMNWaterFall] and the localization[@PhononLocalization], further showing their difference from normal ferroelectrics. In this work, we focus on individual dipoles and statistically model their dielectric response. This change of view point implies that phonon modes are less important in our analysis. We will show in the following that such change leads to fruitful results, and better understanding of relaxors.

Statistics of individual dipoles
--------------------------------

Individual dipoles can be categorized into different groups based on their dynamics, and each group shall have different contribution to susceptibility. We proceed simplify the interaction between dipoles, assuming that the interaction effectively introduce a potential well of *average* depth, $E_{b}$. We may relate $E_{b}$ to the size of PNRs arising from the clustering of same-type ions [@Akbarzadeh2012] and/or random electric field caused by heterovalent ions [@PMNeffHami2015; @Kleemann2014; @PMNWaterFall].

Since the kinetic energy of individual dipoles obeys the Maxwell-Boltzmann distribution, at temperature $T$, the number of dipoles with kinetic energy $E_{\textrm{kin}}$ is given by $$\begin{aligned}
f\left(E_{\textrm{kin}}\right)= & 2N\sqrt{\frac{E_{\textrm{kin}}}{\pi}}\left(\frac{1}{k_{B}T}\right)^{3/2}\exp\left(-\frac{E_{\textrm{kin}}}{k_{B}T}\right),\label{eq:kinetic-energy}\end{aligned}$$ where $k_{B}$ is the Boltzmann constant, $N$ is the total number of dipoles, and $f\left(E_{\textrm{kin}}\right)dE_{\textrm{kin}}$ is the number of dipoles having a kinetic energy between $E_{\textrm{kin}}$ and $E_{\textrm{kin}}+dE_{\textrm{kin}}$. With this distribution function, we can calculate the number of dipoles with kinetic energy that exceeds the potential well $E_{b}$, which is given by $$\begin{aligned}
N_{1}\left(E_{b},T\right) & =\int_{E_{b}}^{\infty}dE_{\textrm{kin}}f\left(E_{\textrm{kin}}\right)\nonumber \\
= & N\sqrt{\frac{4}{\pi}}\sqrt{\frac{E_{b}}{k_{B}T}}\exp\left(-\frac{E_{b}}{k_{B}T}\right)+N\text{erfc}\left(\sqrt{\frac{E_{b}}{k_{B}T}}\right),\label{eq:particles-above-well}\end{aligned}$$ where erfc is the complementary error function. The number of dipoles confined to the potential well is then given by $$\begin{aligned}
N_{2}\left(E_{b},T\right)= & N-N_{1}\left(E_{b},T\right).\label{eq:particles-inside-well}\end{aligned}$$

The next step is to treat the two sets of dipoles ($N_{1}$ versus $N_{2}$) separately, assigning different susceptibility to them. The total susceptibility is then given by $$\begin{aligned}
\chi(T,\nu)= & \chi_{1}\left(T,\nu\right)P_{1}\left(E_{b},T\right)+\chi_{2}\left(T,\nu\right)P_{2}\left(E_{b},T\right),\label{eq:total-susceptibility}\end{aligned}$$ where $\chi_{1}\left(T,\nu\right)$ and $\chi_{2}\left(T,\nu\right)$ describes the dielectric responses of each dipole group, whose form will be specified later. We also define $$\begin{aligned}
P_{1}\left(E_{b},T\right)\equiv & \frac{1}{N}N_{1}\left(E_{b},T\right),\\
P_{2}\left(E_{b},T\right)\equiv & \frac{1}{N}N_{2}\left(E_{b},T\right),\end{aligned}$$ to normalize the dipoles to unit volume. Equation (\[eq:total-susceptibility\]) is the centerpiece of this work and will be demonstrated to be useful for investigation of various relaxors.

Results \[sec:Results\]
=======================

We now apply Eq. (\[eq:total-susceptibility\]) to fit various $\chi$ versus $T$ obtained experimentally or numerically. The relaxors shown here include both lead-basd (e.g., PMN) and lead-free relaxors (e.g., BZT).

Susceptibility of BZT \[sec:lead-free-relaxor\]
-----------------------------------------------

For the static susceptibility of lead-free relaxor BZT[@BZT-susceptibility; @Prosandeev; @dielectric-susceptibility; @Petzelt2014; @Wang2016], we assume (i) dipoles with kinetic energy that overcomes potential well can be treated as free dipoles, subject only to thermal excitation; (ii) dipoles inside the potential well only contribute a constant susceptibility, $\chi_{2}$. The total susceptibility is thus given by $$\begin{aligned}
\chi(T)= & \chi_{1}\mathcal{L}\left(\frac{\theta}{T}\right)P_{1}\left(E_{b},T\right)+\chi_{2}P_{2}\left(E_{b},T\right)\label{eq:BZT-susceptibility}\end{aligned}$$ where $\mathcal{L}\left(x\right)=\coth\left(x\right)-1/x$ is the Langevin function, known for depicting orientational polarization under thermal fluctuations [@Kasap2006; @CurieLaw]. $E_{b}$, $\chi_{1}$, $\chi_{2}$ and $\theta$ are constants, which will be determined by fitting experimental or numerical data. It can be inferred from equation \[eq:BZT-susceptibility\] that $\chi_{1}$ is the susceptibility of the material at zero Kevin, $\chi_{1}\mathcal{L}\left(\frac{\theta}{T}\right)$ is essentially the Curie law at high temperature, and $\theta$ is proportional to the magnitude of the low-frequency electric field used in experimental measurements.

![Fitting the static susceptibility of Ba(Zr$_{0.5}$,Ti$_{0.5}$)O$_{3}$ obtained from Monte-Carlo simulation using Eq. (\[eq:BZT-susceptibility\]). The blue dots are from Monte-Carlo simulation [@Akbarzadeh2012] and the red solid line is the fitting curve using Eq. (\[eq:BZT-susceptibility\]). \[fig:Fitting-susceptibility-of\]](Fig1){width="8cm"}

We first test Eq. (\[eq:BZT-susceptibility\]) against the static susceptibility versus temperature obtained with Monte-Carlo (MC) simulation in a previous work [@Akbarzadeh2012]. \[fig:Fitting-susceptibility-of\] shows the overall fitting is good enough to reproduce results from MC simulations with parameters shown in and Tab. \[tab:parameters-BZT-static\]. The closeness of $E_{b}$ and $T_{m}$ indicates the average depth of potential wells plays a dominant role in determining $T_{m}$. Close examination of Fig. \[fig:Fitting-susceptibility-of\] also indicates the fitting at the lowest temperature ($\lesssim25$ K) is not as good as the rest. To address this issue, we tried adding a Gaussian distribution to $E_{b}$ and remedied the minor problem. However, the resulting equation is quite complicated and deviates from our original goal of proposing simple analytical formula to fit susceptibility. Therefore this additional step is not adopted here.

            $\chi_{1}$   $\theta$ (K)   $\chi_{2}$   $E_{b}$ (K)
  -------- ------------ -------------- ------------ -------------
   Values     741.6         220.5          64.7         159.1

  : Fitting parameters for the Ba(Zr$_{0.5}$,Ti$_{0.5}$)O$_{3}$ static susceptibility. \[tab:parameters-BZT-static\]

![Maxwell-Boltzmann distribution [\[]{}Panel (a)[\]]{} and the Langevin function, $\mathcal{L}$, [\[]{}Panel (b)[\]]{} versus temperature. Parameters from Tab. \[tab:parameters-BZT-static\] are used in plot each function. \[fig:components\]](Fig_0){width="6cm"}

In order to have a good understanding of BZT’s susceptibility, we show each component of Eq. (\[eq:BZT-susceptibility\]) in Fig. \[fig:components\]. Figure \[fig:components\] (a) shows that $P_{1}\left(E_{b},T\right)$ and $P_{2}\left(E_{b},T\right)$ have opposite trends as temperature increases. The number of dipoles that can overcome the potential confinement ($P_{1}$) steadily increases with temperature, while the number of dipoles inside ($P_{2}$) continuously becomes smaller. Figure \[fig:components\] (b) shows that the Lagevin function is normalized at $T=0$, and decreases with temperature. Such a feature describes the ability of the free dipoles to respond to an external *dc* electric field. Moreover, Fig. \[fig:components\] (b) also shows the product of the Langevin function and $P_{1}$, which already exhibits some resemblance to BZT’s susceptibility [\[]{}Fig. (\[fig:Fitting-susceptibility-of\])[\]]{}.

![Fitting the susceptibility of Ba(Zr$_{0.5}$,Ti$_{0.5}$)O$_{3}$ at $f=1,10,100,100$ GHz, obtained from molecular dynamics simulation using Eq. (\[eq:BZT-frequency\]).\[fig:Fitting-BZT-frequency\]](Fig2){width="8cm"}

Having examined the static susceptibility, we now move to the frequency-dependent dielectric response, which is often taken as a characteristic property of relaxors [@Kleeman2014; @Sommer]. We propose another equation to fit the susceptibility versus temperature:

$$\begin{aligned}
\chi\left(T\right)= & \frac{\chi_{1}}{1+b\exp\left(-\theta/T\right)}P\left(E_{b},T\right)+\chi_{2}\left[1-P\left(E_{b},T\right)\right],\label{eq:BZT-frequency}\end{aligned}$$

where $\chi_{1}$ , $\chi_{2}$, $b$ and $\theta$ are constants at a given frequency (but may change when the frequency changes). The dielectric contribution from dipoles with kinetic energy higher than the potential well is given by $$\begin{aligned}
w_{1}\left(T\right)= & \frac{1}{1+b\exp\left(-\theta/T\right)}.\label{eq:h1T}\end{aligned}$$ which, similar to the Langevin function, monotonically decreases with temperature $T$. The choice of this function reflects two considerations: (i) at very low temperature ($T$ close to 0), such dipoles shall follow the probing $ac$ electric field closely, leading $w_{1}\left(T\right)$ to its maximum; and (ii) at higher temperature, thermal motions of these dipoles hamper their ability to follow the *ac* electric field, leading to smaller $w_{1}\left(T\right)$. We will discuss this equation further in Sec. \[subsec:Discussion\]. With one more parameter ($b$), this function may be taken as an extension to the Langevin function.

                  1 GHz   10 GHz   100 GHz   1000 GHz
  -------------- ------- -------- --------- ----------
   $\theta$ (K)   579.6   762.6    1128.4      2158
    $\chi_{1}$    406.5   312.6     209.3      99.9
    $\chi_{2}$    73.1     63.5     56.2       57.9
       $b$        10.2     9.9       9.5       7.7

  : Fitting parameters of numerically simulated Ba(Zr$_{0.5}$,Ti$_{0.5}$)O$_{3}$ susceptibility at various frequency [@Wang2016] using Eq. (\[eq:BZT-frequency\]). \[tab:Fitting-parameters-of-BZT\]

![Analysis of fitting parameters versus probing frequency. (a) $\chi_{1}$vs $\log\left(f\right)$; (b) $w_{1}\left(T\right)$. \[fig:w1T\]](Fig3_ab){width="8cm"}

We use Eq. (\[eq:BZT-frequency\]) to fit BZT’s susceptibility versus temperature at different frequency and show the results in Fig. \[fig:Fitting-BZT-frequency\]. Since $E_{b}$ is a material parameter, we use the same value ($E_{b}=159.1$ K) obtained by fitting the static susceptibility (cf. Fig. \[fig:Fitting-susceptibility-of\]). In Fig. \[fig:Fitting-BZT-frequency\], the numerical results are obtained from molecular dynamics simulations reported in Ref. [@Wang2016]. As the figure shows, satisfactory fittings are achieved for frequencies between 1 and 1000 GHz. Table \[tab:Fitting-parameters-of-BZT\] shows all the parameters. Among them, $\theta_{1}$ and $\chi_{1}$ have substantial changes over the specified frequency range as shown in Fig \[fig:w1T\](a). Figures \[fig:w1T\] (a) shows that $\theta$ depends on $\log\left(f\right)$ quadratically while $\chi_{1}$ linearly depends on it, and Fig. \[fig:w1T\](b) shows $w_{1}\left(T\right)$. At low frequency ($\lesssim10$ GHz), $w_{1}\left(T\right)$ resembles the Fermi-Dirac function, that is, below$\sim250$ K, its value is close to one but becomes close zero for $T$ above $\sim250$ K. At a higher frequency (e.g., 1000 GHz), however, this function strongly deviates from the Fermi-Dirac function, with a long tail extending to high temperature.

![Experimental susceptibility of Ba(Ti$_{0.675}$Zr$_{0.325}$)O$_{3}$ ceramics versus temperature fitted with Eq. (\[eq:BZT-frequency\]). \[fig:real-bzt\]](Fig4_v2){width="6cm"}

To further verify the suitability of this equation for experimental data, we also fit the result shown in Fig. 1 of Ref. [@Bokov2007BZT], where Ba(Ti$_{0.675}$Zr$_{0.325}$)O$_{3}$ ceramics is measured at $10^{-2}$ and $10^{5}$ Hz. Figure \[fig:real-bzt\] shows that satisfactory fittings are achieved.

Pb-based relaxors \[subsec:Lead-based-relaxors\]
------------------------------------------------

Unlike the lead-free BZT, which possesses PNRs that separate dipole clusters, lead-based ferroelectrics [@Kirillov1973; @Gridnev2004; @PMNeffHami2015; @PNR-PMN] have the Pb-driven dipoles across the system [@Yk], which cause phase transitions in systems such as PbTiO$_{3}$ and Pb(Zr,Ti)O$_{3}$ [@PZTPhase; @PTOeffHam; @JiangPTOFilm-1]. Due to heterovalent ions inside, typical lead-based relaxors are subject to random electric fields, which distort the orientation of dipoles. While the precise consequence of the random field is not all clear [@Sherrington2014; @Kleemann2015], such distracting effect on dipoles appears to lead to a strongly modified phase transition with diffused and smeared peak, in sharp contrast to that of normal ferroelectrics [@Pirc2007; @PMNeffHami2015].

To model such a system and account for the moderate phase transition, we need a function that properly describes the dielectric constant versus temperature. Here, we propose to use the slightly modified well known quadratic relation $\frac{\varepsilon_{A}}{\varepsilon^{\prime}}-1=\frac{\left(T-T_{A}\right)^{2}}{2\delta^{2}}$ proposed by Bokov *et al* [@Bokov2006a] to relate relaxor’s permittivity to temperature [@Smolenskii1970; @Kirillov1973] for dipoles above the average potential well (also see Eq. (\[eq:square-law\]) ). This equation can be rearranged to give the following expression $$\begin{aligned}
w_{2}\left(T\right)= & \frac{1}{1+\left|\frac{T-T_{O}}{\theta}\right|^{\gamma}},\label{eq:w2t}\end{aligned}$$ where $\gamma$ is a critical exponent, $T_{o}$ is associated with the peak position of the moderate phase transition, $\theta$ and $\gamma$ are parameters describing the peak. We note that such choice of $w_{2}\left(T\right)$ also agree with the analysis recently given by Uchino [@Uchino2014]. Combining Eqs. (\[eq:total-susceptibility\]) and (\[eq:w2t\]), we obtain the following equation to fit lead-based relaxors, $$\begin{aligned}
\chi\left(T\right)= & \frac{\chi_{1}}{1+\left|\frac{T-T_{O}}{\theta}\right|^{\gamma}}P\left(E_{b},T\right)+\chi_{2}\left[1-P\left(E_{b},T\right)\right],\label{eq:PZN}\end{aligned}$$ where $E_{b}$, $\chi_{1}$, $\chi_{2}$, $T_{o}$, $\theta$, and $\gamma$ are fitting parameters. The meaning of $\chi_{1}$, $\chi_{2}$, and $E_{b}$ are the same as discussed in Sec. \[sec:lead-free-relaxor\].

![Fitting the relative permittivity of PZN-0.13PT at $f=1,10,100$ kHz using Eq. (\[eq:PZN\]). Note the abnormal decrease in the range below $\sim250$ K. The second panel is $w_{2}\left(T\right)$. \[fig:PZN-relaxor\]](Fig5_v2){width="6cm"}

To verify that Eq. (\[eq:PZN\]) indeed works, we experimentally obtained the permittivity of Pb(Zn$_{1/3}$Nb$_{2/3}$)$_{0.87}$Ti$_{0.13}$O$_{3}$ (PZN-0.13PT) versus temperature at frequencies $f=1,10,100$ kHz. As Fig. \[fig:PZN-relaxor\](a) shows, all three fittings are satisfactory. With Eq. (\[eq:PZN\]) and the fittings, we are able to single out $w_{2}\left(T\right)$, which shows a slight increase of the peak temperature ($T_{o}$, around $220$ K) with frequency [\[]{}Fig. \[fig:PZN-relaxor\](b)[\]]{}. We note that in this fitting there is no need to have two $\gamma$ values above and below $T_{o}$ [@Uchino2014]. The asymmetric peak shown in Fig. \[fig:PZN-relaxor\](a) is naturally caused by the Maxwell-Boltzmann distribution function.

                     1 kHz    10 kHz   100 kHz
  ----------------- -------- -------- ---------
      $\gamma$        2.02     1.82     1.63
     $T_{O}$ (K)     219.6    229.0     240.6
     $E_{b}$ (K)      22.3     22.9     23.5
   $\chi_{1}$ [^1]   56601    55781     53529
     $\chi_{2}$      1320.9   1284.2   1238.0
    $\theta$ (K)     102.4    103.1     104.5

  : Fitting parameters of PZN-0.13PT’s permittivity measured at different frequencies. As the table demonstrates, $\gamma$ and $T_{C}$ are the most important variable that changes a lot with frequency. \[tab:Fitting-parameters-of-PZN\]

Table \[tab:Fitting-parameters-of-PZN\] summarizes fitting parameters of the permittivity measured at various frequencies. Among all the parameters, the critical component $\gamma$ changes most (19.3% from 1 kHz to 100 kHz), and decreases with increasing frequency; similarly, $T_{O}$ also changes by 9.5%. On the other hand, $E_{b}$, $\chi_{1}$, $\chi_{2}$, $\theta$ are relatively constant, which are independent of the frequency, and may be taken as material parameters. Such results hints towards the following formula that describes the dependence of PZN-0.13PT on both temperature and the probing frequency $$\begin{aligned}
\chi\left(T,\nu\right)= & \frac{\chi_{1}}{1+\left|\frac{T-T_{O}\left(\nu\right)}{\theta}\right|^{\gamma\left(\nu\right)}}P\left(E_{b},T\right)+\chi_{2}\left[1-P\left(E_{b},T\right)\right],\label{eq:PZN-1}\end{aligned}$$ where the two functions $T_{O}\left(\nu\right)$ and $\gamma\left(\nu\right)$ are frequency dependent while other parameters are constants for a given material. It is also worth noting that for PZN-0.13PT $\chi_{2}\ll\chi_{1}$, which indicates that dipoles with kinetic energy above the potential well play a more important role, in contrast to the case of BZT (see Tab. \[tab:Fitting-parameters-of-BZT\]).

![Fitting the relative permittivity of Pb(Mg$_{1/3}$Nb$_{2/3}$)O$_{3}$-0.05Pb(Zr$_{0.53}$Ti$_{0.47}$)O$_{3}$ measured at $800$ Hz using Eq. (\[eq:PZN\]). The second panel is $w_{2}\left(T\right)$. \[fig:PMN-PZT\]](Fig6_new_v3){width="6cm"}

To further verify the proposed formula, we also fit the permittivity versus temperature of another lead-based relaxor, Pb(Mg$_{1/3}$Nb$_{2/3}$)O$_{3}$-0.05Pb(Zr$_{0.53}$Ti$_{0.47}$)O$_{3}$ [@Gridnev2004]. It can be seen from Fig. \[fig:PMN-PZT\](a) that the overall fitting is satisfactory. Figure \[fig:PMN-PZT\](b) shows $w_{2}\left(T\right)$ with fitting parameters $\gamma=1.66$ and $T_{O}=256.7$ K. Similar to the PZN-0.13PT case, the results here also shows $\chi_{2}\ll\chi_{1}$.

Discussion\[subsec:Discussion\]
===============================

In the statistical model we divide the dipoles inside ferroelectrics relaxors into two groups, one group being confined in potential wells, while the other having can overcome the potential confinement and exhibiting a more vibrant dynamics. It has been demonstrated that the Maxwell-Boltzmann distribution plays a significant role in determining the profile of $\chi\left(T\right)$. To address a particular type of relaxor, one may only need to adjust the dielectric response function associated with each group of dipoles, while keeping the rest unchanged. In this section, we discuss a few issues of this approach and its limitations.

Characteristic temperature $T_{m}$
----------------------------------

The present analysis helps us to understand why the susceptibility of a relaxor reaches its peak value at some temperature, $T_{m}$. For BZT, the function $\chi_{1}\mathcal{L}\left(\frac{\theta}{T}\right)$ [\[]{}see Eq. (\[eq:BZT-susceptibility\]) and Fig. \[fig:components\][\]]{} or $\chi_{1}/\left[1+b\exp\left(-\theta/T\right)\right]$ [\[]{}see Eq. (\[eq:BZT-frequency\])[\]]{} describes the contribution to susceptibility from dipoles with kinetic energy higher than $E_{b}$. These two functions are both monotonically decreasing with $T$, reflecting the fact that thermal motions prevent dipoles from aligning with the applied electric field. On the other hand, the number of dipoles above potential wells increases with $T$ as governed by the function $P\left(E_{b},T\right)$ (see Fig. \[fig:components\]). The combined effects of the two factors give rise to the susceptibility peak at $T_{m}$. However, the situation for lead-based relaxors is different. The function $\chi_{1}/\left(1+\left|\frac{T-T_{O}}{\theta}\right|^{\gamma}\right)$ [\[]{}see Eq. (\[eq:PZN\])[\]]{}, which largely determines the value of $T_{m}$, manifests the vestige of a true phase transition in normal ferroelectrics, which is torn down by random electric fields and/or PNRs in relaxors.

Rationale for Eq. (\[eq:h1T\])
------------------------------

![The energy well for dipoles is shifted and lowered when an electric field is applied .\[fig:falling\]](Fig7){width="8cm"}

For lead-free BZT, we propose Eq. \[eq:h1T\] to describe the susceptibility of dipoles with kinetic energy higher than $E_{b}$. This choice follows the Debye relaxation, i.e., $\chi\sim1/\left(1+\omega^{2}\tau^{2}\right)$ [@Kasap2006], where $\omega$ is a constant (the probing frequency), and $\tau$ is temperature-dependent relaxation time. For a thermally activated process, the relation between $\tau$ and $T$ is often specified by the Arrhenius law, i.e., $1/\tau=A\exp\left(-E_{a}/T\right)$, where $E_{a}>0$ is the activation energy [@PinnedPhonon; @negativeActivation; @Pirc2007]. In this case, the susceptibility will be $\chi\sim1/\left[1+A^{2}\omega^{2}\exp\left(\frac{2E_{a}}{T}\right)\right]$, which is discussed by Jonscher [@Jonscher1981]. However, the dynamic process considered here describes dipoles falling to a state of lower energy, which temporarily created by the probing electric field (see. Fig. \[fig:falling\]). Therefore the activation energy in this process shall be *negative*, i.e., $\chi\sim1/\left[1+A^{2}\omega^{2}\exp\left(-\frac{2\left|E_{a}\right|}{T}\right)\right]$, which is the form adopted in Eq. (\[eq:h1T\]).

We note that negative activation energy is known in some chemical reactions [@negativeActivation]. Negative activation energy appears here because when an *ac* electric field perturbs dipoles and tilts the relative energy of potential wells, dipoles outside potential wells will move towards the temporary potential minimum. However, the drifting to the potential minimum is hindered by thermal fluctuations of such dipoles. In fact, higher temperature (corresponding to larger kinetic energy) results in a slower relaxation to the energy minimum (corresponding to larger $\tau$), leading to negative activation energy. We also note that since the applied *ac* electric field is responsible for shifting the energy minimum and causing dipoles to drift, the change of its frequency may well alter the negative activation energy, explaining why $\theta$ in Eq. (\[eq:BZT-frequency\]) is dependent on the probing frequency. Similar arguments explains why $T_{O}$ also depends on the probing frequency.

In the above analysis, we focus on dipoles with kinetic energy higher than $E_{b}$. These dipoles are able to drift from one energy minimum to another when an *ac* electric field perturbs the system. It has been proposed that the drifting/hopping of dipoles from one potential well to another causes relaxations . However, without distinguishing dipoles inside and outside the potential well, such proposal seems to have a tendency of confusing the wait time before hopping, $t$ (which reflects the distribution of dipoles at a given temperature) with the relaxation time, $\tau$ (which reflects how fast dipoles drift to the transient energy minimum and relates to the loss peak frequency in the Debye function), leading to some difficulties.

Limitations
-----------

In previous studies [@Wang2014; @Wang2016; @PMNeffHami2015], *ab initio* calculation was used, which prescribes all important interactions between dipoles and other degrees of freedom in relaxors, and MC or MD was used to numerically work out the consequences. In the present work we do not start from *ab initio* calculation, instead, employs statistical and phenomenological arguments. Having shown results and insights obtained with this approach, we now discuss possible limitations to the present approach with respect to treatment of details, accuracy, and prediction power.

First, the proposed equations for lead-free [\[]{}Eq. (\[eq:BZT-frequency\])[\]]{} and lead-based relaxors [\[]{}Eq. (\[eq:PZN-1\])[\]]{} have five and seven parameters, respectively. Ideally, one hopes to be able to reduce this large number and use as few parameters as possible. However, it shall be noted that, among these parameters, many are only material dependent (i.e., they do not depend on frequency or temperature). For instance, for lead-free relaxor, $E_{b}$ is a constant; for lead-based relaxor, $E_{b}$, $\chi_{1}$, $\chi_{2}$, $\theta$ are close to constants (see Tab. \[tab:Fitting-parameters-of-PZN\]). For a given sample, these parameters may only need to be calibrated once. In this way, the number of parameters will be significantly reduced.

Second, in this work we have focused on the temperature dependence of susceptibility. The dependence on frequency needs further investigation. For instance, $T_{C}\left(\nu\right)$ and $\gamma\left(\nu\right)$ in Eq. (\[eq:PZN-1\]) need to be specified explicitly to address this issue. We note that results shown in Tab. \[tab:Fitting-parameters-of-PZN\] will provide clues to $T_{C}$ and $\gamma$’s dependence on $\nu$, and eventually help finding analytical expressions for $\chi\left(T,\nu\right)$. In addition, we generally ignored the long-range correlation of dipoles, which is another limitation to this approach. While such long-range correlation makes relaxor physics so rich, it will bring back Bose-Einstein statistics and make the current formulation more complicated. To what extent the Bose-Einstein and the Maxwell-Boltzmann distribution shall be used for ferroelectric relaxors remains an open question.

Third, at high temperature, Curie law is observed in many ferroelectrics. For the static susceptibility of BZT [\[]{}Fig. \[fig:Fitting-susceptibility-of\][\]]{}, this law can be recovered from the proposed equation, Eq. (\[eq:BZT-susceptibility\]). On the other hand, for Eqs. (\[eq:BZT-frequency\]) and (\[eq:PZN\]), the Curie law cannot be directly recovered. For Eq. (\[eq:PZN\]), we have the asymptotic relation $\text{\ensuremath{\chi}}\sim A/\left(T-T_{C}\right)^{\gamma}+B/T^{3/2}$ at very large $T$. It is unclear how good this relation can fit the Curie law. Therefore, in fitting experimental data at high temperatures, one needs to bear in mind that Eqs. (\[eq:BZT-frequency\]) and (\[eq:PZN\]) should be used with care.

Finally, with Eqs. (\[eq:BZT-susceptibility\]) and (\[eq:PZN\]), in principle we can obtain the relation between $T_{m}$ and $\nu$, which can then be compared to the well-known Vogel-Fulcher law [@Wang2016]. However, we have failed to obtain analytical expressions for $T_{m}\left(\nu\right)$ and believe numerical calculation seems to be the only feasible way to establish the relation between $T_{m}$ and $\nu$.

Conclusion\[sec:Conclusion\]
============================

Instead of working on the atomic level, the present work employs a macroscopic statistical approach to help understanding dielectric properties of relaxors. The effects of disorder, PNRs, and random electric fields are considered statistically by introducing the average potential well, which can trap dipoles of low kinetic energy. An external electric field will mostly increase the magnitude **** of trapped dipoles, but rotate to align dipoles free from such trapping, giving rise to two different types of dielectric responses as shown in Eqs. (\[eq:BZT-susceptibility\]), \[eq:BZT-frequency\], and (\[eq:PZN\]). This approach, by proposing analytical equations, provides insights to experimental and numerical results of relaxors. **** Among other things, it shows that the characteristic temperature, $T_{m}$, is determined by the Maxwell-Boltzmann distribution of dipoles’ kinetic energy, as well as their ability to respond to the applied electric field. We can also conclude that lead-free relaxors (e.g., BZT) are different from lead-based relaxors (e.g., PZN-0.13PT) in that (i) The mechanisms determining $T_{m}$ are different. For lead-based relaxors, it appears $T_{O}$ alone in able to determine $T_{m}$, while for BZT, both $w_{1}\left(T\right)$ and $P\left(E_{b},T\right)$ are important; and (ii) For BZT, $\chi_{1}$ and $\chi_{2}$ are on the same order, in contrast to the fact while $\chi_{1}\gg\chi_{2}$ for the Pb-based relaxors, indicating that dipoles outside the average potential well dominate dielectric response of Pb-based relaxors. With these results, we hope this statistical approach can help better understanding important relaxor systems and the proposed equations be adopted in fitting experimental data.

We thank Drs C.-L. Wang, A. A. Bokov and L. Bellaiche for fruitful discussion. This work is financially supported by the National Natural Science Foundation of China (NSFC), Grant No. 51390472, 11574246, U1537210, and National Basic Research Program of China, Grant No. 2015CB654903. F.L. acknowledges NSFC Grant No. 51572214. Z.J. acknowledges the support from China Scholarship Council (CSC No. 201506280055). We also acknowledge the “111 Project of China (Grant No. B14040), the Natural Science and Engineering Council of Canada (NSERC) and the United States Office of Naval Research (ONR Grants No. N00014-12-1-1045 and N00014-16-1-3106).

urlstyle\[1\][doi: \#1]{}

[Petzelt et al.(2014)Petzelt, Nuzhnyy, Savinov, Bovtun, Kempa, Ostapchuk, Hlinka, Canu, and Buscaglia]{} K. Uchino, . (Springer, 1996).

M. E. Manley, J. W. Lynn, D. L. Abernathy, E. D. Specht, O. Delaire, A. R. Bishop, R. Sahul, and J. D. Budai, Phonon localization drives polar nanoregions in a relaxor ferroelectric. *Nat. Commun.* **5**, 3683 (2014).

A. A. Bokov and Z.-G. Ye, Recent progress in relaxor ferroelectrics with perovskite structure. *J. Mater. Sci.* **41**, 31 (2006).

D. Wang, J. Hlinka, A. A. Bokov, Z. G. Ye, P. Ondrejkovic, J. Petzelt, and L. Bellaiche. Fano resonance and dipolar relaxation in lead-free relaxors. *Nat. Commun.* **5**, 5100 (2014).

D. Wang, A. A. Bokov, Z.-G. Ye, J. Hlinka, and L. Bellaiche, Subterahertz dielectric relaxation in lead-free Ba(Zr,Ti)O$_{3}$ relaxor ferroelectrics. *Nat. Commun.* **7**, 11014 (2016).

D. Nuzhnyy, J. Petzelt, and M. Savinov, Broadband dielectric response of Ba(Zr,Ti)O$_{3}$ ceramics: from incipient via relaxor and diffuse up to classical ferroelectric behavior. *Phys. Rev. B* **86**, 1, (2012).

J. Petzelt, D. Nuzhnyy, M. Savinov, V. Bovtun, M. Kempa, T. Ostapchuk, J. Hlinka, G. Canu, and V. Buscaglia, Broadband dielectric spectroscopy of Ba(Zr,Ti)O$_{3}$: dynamics of relaxors and diffuse ferroelectrics.*Ferroelectrics* **469**, 14 (2014).

W. Kleemann, S. Miga, Z. K. Xu, S. G. Lu, and J. Dec, Non-linear permittivity study of the crossover from ferroelectric to relaxor and cluster glass in BaTi$_{1-x}$Sn$_{x}$O$_{3}$. *Appl. Phys. Lett.* **104**, 182910 (2014).

K. Uchino, Fractal Phenomena in Ferroelectrics. *J. Nanotech. Mater. Sci.* **1**, 1 (2014).

A. A. Bokov and Z.-G. Ye, Dielectric relaxation in relaxor ferroelectrics. *J. Adv. Dielectr.* **2**, 1241010 (2012).

A. Akbarzadeh, S. Prosandeev, E. Walter, A. Al-Barakaty, and L. Bellaiche, Finite-temperature properties of Ba(Zr, Ti)O$_{3}$ relaxors from first principles. *Phys. Rev. Lett.* **108**, 257601 (2012).

D. Sherrington, BZT: a soft pseudospin glass. *Phys. Rev. Lett.* **111**, 227601 (2013).

D. Sherrington, Pb(Mg$_{1/3}$Nb$_{2/3}$)O$_{3}$: A minimal induced-moment soft pseudospin glass perspective. *Phys. Rev. B* **89**, 064105 (2014).

W. Kleemann, J. Dec, and S. Miga, The cluster glass route of relaxor ferroelectrics, *Phase Trans.* **88**, 234 (2015).

R. Pirc, R. Blinc, Spherical random-bond-random-field model of relaxor ferroelectrics. *Phys. Rev. B* **60**, 13470 (1999).

D. Sherrington, A spin glass perspective on ferroic glasses.*Phys. Status Solidi* **251**, 1967 (2014).

A. A. Bokov, M. Maglione, and Z.-G. Ye, Quasi-ferroelectric state in Ba (Ti$_{1-x}$Zr$_{x}$) O3 relaxor: dielectric spectroscopy evidence. *J. Phys.: Condens. Matter* **19**, 092001 (2007).

S. A. Gridnev, A. A. Glazunov, and A. N. Tsotsorin, Nonlinear dielectric response of relaxor PMN-PZT ceramics under dc electric field. *Ferroelectrics*, **307**, 151 (2004).

A. A. Bokov and Z.-G. Ye, ** Double freezing of dielectric response in relaxor Pb (Mg$_{1\slash3}$Nb$_{2\slash3}$)O$_{3}$ crystals. *Phys. Rev. B* **74**, 132102 (2006).

H. N. Tailor, A. A. Bokov, Z. G. Ye, **Freezing of polarization dynamics in relaxor ferroelectric (1x)Pb(Mg$_{1/3}$Nb$_{2/3}$)O$_{3}$xBi(Zn$_{1/2}$Ti$_{1/2}$)O$_{3}$ solid solution.** *Curr. Appl. Phys.* **11**, 175 (2011).

A. Dixit, S. B. Majumder, R. S. Katiyar, A. S. Bhalla, **Studies on the relaxor behavior of sol-gel derived Ba(Zr$_{x}$Ti$_{1-x}$)O$_{3}$ (0.30$\le$x$\le$0.70) thin films*. J. Mater. Sci.* **41**, 87 (2006).

A. K. Tagantsev, A. E. Glazounov, ** **Does freezing in Pb(Mg$_{1/3}$Nb$_{2/3}$)O$_{3}$ relaxor manifest itself in nonlinear dielectric susceptibility?** *Appl. Phys. Lett.* **74**, 1910 (1999).

J. P. Sethna, . (Clarendon Press, Oxford, 2010).

L. Cross. . *Ferroelectrics*, **76**,0241 (1987).

A. K. Josnscher, Universal relaxation law (Chelesea Dielectrics Press, 1996).

A. K. Jonscher, Dielectric relaxation in solids (Chelesea Dielectrics Press, 1983).

G. A. Smolenskii, Physical phenomena in ferroelectrics with diffused phase transition. *J. Phys. Soc. Jpn***.** **28**, 26 (1970).

V. V. Kirillov, V. A. Isupov, **Relaxation polarization of PbMg$_{1/3}$Nb$_{2/3}$O$_{3}$ (PMN)-A ferroelectric with a diffused phase transition*. Ferroelectrics*, **5**, 3 (1973).

K. Uchino, S. Nomura, Critical exponents of the dielectric constants in diffused-phase-transition crystals. *Ferroelectrics* **44**, 55 (1982).

R. Clarke and J. C. Burfoot, The diffuse phase transition in potassium strontium niobate. *Ferroelectrics*, **8**, 505 (1974).

I. A. Santos, J. A. Eiras, Phenomenological description of the diffuse phase transition in ferroelectrics. *J. Phys. Condens. Matter* **13**, 11733 (2001).

D. Wang and Z. Jiang, Dielectric response of BaZrO$_{3}$/BaTiO$_{3}$ superlattice, *J. Adv. Dielectr* **6**, 1650015 (2016).

A. A. Bokov, I. P. Raevskii, and V. G. Smotrakov, Composition, ferroelectric and antiferroelectric ordering in Pb$_{2}$InNbO$_{6}$ crystals, *Sov. Phys. Sol. Stat.* **26**,1708 (1984).

N. Setter and L. Cross, The role of B-site cation disorder in diffuse phase transition behavior of perovskite ferroelectrics. *J. Appl. Phys.* **51**, 4356 (1980).

N. Nakanishi, A. Nagasawa, Y. Murakami, Lattice stability and soft modes.****

A. Bussmann-Holder, A. R. Bishop, T. Egami, Relaxor ferroelectrics and intrinsic inhomogeneity. *Europhys. Lett*. **71**, 249 (2005).

W. Zhong, D. Vanderbilt, K. M. Rabe, Phase Transitions in BaTiO$_{3}$ from First Principles. Phys. Rev. Lett. **73**, 1861 (1994).

****

R. A. Cowley, Structural phase transitions I. Landau theory. *Adv. Phys.* **29**, 1-110 (1980).

A. Al-Barakaty, S. Prosandeev, D. Wang, B. Dhkil, and L. Bellaiche, ** Finite-temperature properties of the relaxor PbMg$_{1/3}$Nb$_{2/3}$O$_{3}$ from atomistic simulations. *Phys. Rev. B* **91**, 214117 (2015).

D. Phelan, C. Stocka, J. A. Rodriguez-Riveraa, S. Chia, J. Leãoa , X. Long, Y. Xie, A. A. Bokov, Z.-G. Ye, P. Ganeshe, and P. M. Gehring, Role of random electric fields in relaxors, *PNAS* **111**, 1754 (2014).

J. Hlinka, S. Kamba, J. Petzelt, J. Kulda, C. A. Randall, S. J. Zhang, Origin of the Waterfall Effect in Phonon Dispersion of Relaxor Perovskites. *Phys. Rev. Lett.* **91**, 107602 (2003).

P. M. Gehring, S.-E. Park, and G. Shirane, Soft phonon anomalies in the relaxor ferroelectric Pb(Zr$_{1/3}$Nb$_{2/3}$)$_{0.92}$Ti$_{0.08}$O$_{3}$. *Phys. Rev. Lett.* **84**, 5216 (2000).

W. Kleemann, Cluster glass ground state via random fields and random bonds. *Phys. Status Solidi* **251**, 1993 (2014). .

S. Prosandeev, D. Wang, A. R. Akbarzadeh, L. Bellaiche, **First-principles-based effective Hamiltonian simulations of bulks and films made of lead-free Ba(Zr,Ti)O$_{3}$ relaxor ferroelectrics*. J. Phys . Condens. Matter*. **27**, 223202 (2015).

S. O. Kasap, Principles of Electronic Materials and Devices, Third Edition (McGraw-Hill, 2006), pp612.

We also verified that replacing the Langevin function with the simpler Curie law ($\sim1/T$) will result in similar fitting results and analysis.

R. Sommer, N. K. Yushin, J. J. Van Der klink. **Dielectric susceptibility of PMN under DC bias*. Ferroelectrics* **127**, 235(1992).

I. K. Jeong, T. W. Darling, J. K. Lee, T. Proffen, R. H. Heffner, Direct observation of the formation of polar nanoregions in Pb(Mg$_{1/3}$Nb$_{2/3}$)O$_{3}$ using neutron pair distribution function analysis. *Phys. Rev. Lett.* **94**, 147602 (2005).

Y. Kuroiwa, S. Aoyagi, A. Sawada, J. Harada, E. Nishibori, M. Takata, M. Sakata, Evidence for Pb-O covalency in tetragonal PbTiO$_{3}$. *Phys. Rev. Lett.* **87**, 217601 (2001).

I. A. Kornev, L. Bellaiche, P. E. Janolin, B. Dkhil, and E. Suard, Phase diagram of Pb(Zr,Ti)O$_{3}$ solid solutions from first principles. *Phys. Rev. Lett.* **97**, 157601 (2006).

D. Sichuga and L. Bellaiche, ** Epitaxial Pb(Zr,Ti)O$_{3}$ ultrathin films under open-circuit electrical boundary conditions. *Phys. Rev. Lett*. **106**, 196102 (2011).

Z. Jiang, R. Zhang, D. Wang, D. Sichuga, C.-L. Jia, and L. Bellaich, Strain-induced control of domain wall morphology in ultrathin PbTiO$_{3}$ films. *Phys. Rev. B* **89**, 214113 (2014).

R. Pirc and R. Blinc, Vogel-Fulcher freezing in relaxor ferroelectrics*. Phys. Rev. B* **76**, 020101 (2007).

J. M. Valverde, On the negative activation energy for limestone calcination at high temperatures nearby equilibrium. *Chem. Eng. Sci.* **132**, 169 (2015).

A. K. Jonscher, A new understanding of the dielectric relaxation of solids. *J. Mater. Sci.* **16**, 2037 (1981).

[^1]: Here, $\chi_{1}$ and $\chi_{2}$ shall be understood as relative permittivity, not susceptibility.
