Ferroelectrics are multifunctional smart materials finding applications in sensor technology, micromechanical actuation, digital information storage, etc. Their most fundamental property is the ability of polarization switching under an applied electric field. In particular, understanding of switching kinetics is essential for digital information storage. In this regard, scaling properties of the temporal polarization response are well-known for 180°-switching processes in ferroelectrics characterized by a unique field-dependent local switching time. Unexpectedly, these properties are now observed in multiaxial polycrystalline ferroelectrics, exhibiting a number of parallel and sequential non-180°-switching processes with distinct switching times. This behavior can be explained by a combination of the multistep stochastic mechanism and the inhomogeneous field mechanism models of polarization reversal. Scaling properties are predicted for polycrystalline ferroelectrics of tetragonal, rhombohedral, and orthorhombic symmetries and are exemplarily demonstrated by the measurements of polarization kinetics in (K,Na)NbO_{3}-based ferroelectric ceramic over a timescale of 7 orders of magnitude. Dynamic scaling properties allow insight into the microscopic switching mechanisms, on the one hand, and into statistical material characteristics, on the other hand, thereby providing the description of temporal polarization with high accuracy. The gained deeper insight into the mechanisms of multistep polarization switching is crucial for future ultrafast and multilevel digital information storage.

## I. INTRODUCTION

The most characteristic feature of ferroelectric materials is the appearance of spontaneous polarization below the ferroelectric–paraelectric phase transition temperature. Spontaneous polarization can be switched to another direction by applying an electric field or mechanical stress. The temporal characteristics of the response to an applied electric field are crucial for various applications of these materials, in particular, ferroelectric memories (FeRAM),^{1} future multilevel data storage,^{2–5} and ferroelectric field-effect transistors (FeFET) assembled in neuromorphic computing systems.^{6} The field-driven polarization kinetics in single-crystalline ferroelectrics was successfully described by the stochastic Kolmogorov–Avrami–Ishibashi (KAI) model^{7–9} assuming the formation of 180°-reversed polarization domains characterized by a unique field-dependent switching time. This model, however, failed when applied to thin ferroelectric films and bulk ferroelectric ceramics, both of which exhibit a broad statistical distribution of switching times.^{10–13} To account for the dispersive character of the polarization response, the Nucleation Limited Switching (NLS) model was developed by Tagantsev *et al*.,^{10} and subsequently, applied successfully to lead zirconate titanate (PZT) thin films.^{11} The NLS model, however, did not reveal the physical origin of the wide distribution of switching times. To explain this feature, an idea was suggested by Lupascu *et al*.^{12} that the statistical distribution of switching times may result from the random spatial distribution of an applied field in polycrystalline ferroelectrics, considering the well-known strong field dependence of the switching time.^{14} This idea was first implemented quantitatively by Jo *et al*. using the Lorentzian statistical distribution of local field values in the NLS approach and successfully applied to the PZT thin films,^{15,16} organic ferroelectrics,^{17–19} and composites.^{20,21} Alternatively, the Gaussian field distribution was proposed for bulk PZT ceramics^{22} but also seems to be suitable for HfO_{2}-based thin films.^{6,23–28} It was established, however, that there is no need to *a priori* assume a statistical distribution of local fields because they can be extracted directly from the polarization response, as will be explained below.

^{29}and (Ba

_{0.85}Ca

_{0.15})TiO

_{3}(BCT)

^{30}ceramics of different phase symmetries, 0.94(Bi

_{1/2}Na

_{1/2})TiO

_{3}–0.06BaTiO

_{3}(BNT-BT),

^{31}and (1–x)Ba(Zr

_{0.2}Ti

_{0.8})O

_{3}–x(Ba

_{0.7}Ca

_{0.3})TiO

_{3}(BZT-BCT)

^{32}solid solutions of different phase symmetries, organic ferroelectric copolymers of vinylidenefluoride and trifluoroethylene P(VDF-TrFE),

^{33}and HfO

_{2}-based thin films.

^{23–28}The scaling property is expressed by the fact that the normalized derivative of the polarization variation $\Delta P( E a,t)$ with respect to the applied field, $ E a$, can be represented as a function of the combined variable of $ E a$ and time

*t*,

*t*fall on the same master curve when the field argument is normalized to $ E max(t)$. If relation (1) is experimentally confirmed, that is the case for all the materials mentioned above,

^{23–33}then the polarization response can be described with high accuracy by a formula,

^{34}

Explanation of this behavior was given by the Inhomogeneous Field Mechanism (IFM) model^{34,35} based on the following assumptions:

The field-driven polarization reversal proceeds via statistically independent one-step 180°-switching events (as assumed in the KAI model).

During the polarization reversal, a stable statistical distribution of local fields $F(E/ E a,\theta )$ remains, which depends on the polar angle $\theta $ with respect to the applied field direction and the local field amplitude,

*E*, normalized to the applied field $ E a$.Local polarization switching is characterized by a switching time

*τ*(*E*) dependent on the local electric field magnitude*E.*The KAI time dependence for local polarization switching, $1\u2212 exp[ \u2212 ( t / \tau ) \beta ]$, with the Avrami index

*β*[8,9], is approximated by the Heaviside step function $\u03d1(t\u2212\tau )$.

The assumption of the stable statistical field distribution $F(E/ E a,\theta )$ means neglecting the emerging depolarization fields, a hypothesis supported by the recent self-consistent simulations of polarization reversal in polycrystalline ferroelectrics.^{36} We note that depolarization field variations become locally relevant at the spatial scale of the typical grain size, where they mediate polarization correlations in the neighboring grains. This is in agreement with experimental studies on thin film^{37} and bulk^{38} ferroelectrics. The step-function approximation on the time scale is justified when the local polarization switching step is sharper than the statistical distribution of local fields;^{34} therefore, it works well in polycrystalline systems, including highly textured ceramics,^{30} but fails in single crystals.^{39}

^{34}

*γ*a constant explained later and the mean cosine of

*θ*defined by the relation

A conceptual problem appears when considering polarization reversal in multiaxial ferroelectrics, as was noted in Ref. 40. As was established earlier by x-ray diffraction^{41} and ultrasonic measurements^{42} on tetragonal ferroelectric ceramics, polarization switching occurs at least partially by sequential 90°-reorientations. Distinct characteristic times of two consecutive non-180° switching steps were resolved by *in situ* x-ray diffraction studies.^{43,44} The recently elaborated original Multistep Stochastic Mechanism (MSM) model allowed the identification of fractions, activation fields, and switching times of sequential 90°–90°-switching processes in tetragonal ceramics,^{45} 109°–71°-, 71°–109°-, and 71°–71°–71°-switching processes in rhombohedral ceramics,^{46} and 120°–60°-, 60°–120°-, 90°–90°-, and 60°–60°–60°-switching processes in orthorhombic ceramics.^{47} Considering a large number of polarization switching paths in multiaxial ferroelectrics with switching times of distinct switching events different by orders of the magnitude, the scaling property of the polarization response was not expected in these systems. Paradoxically, the scaling properties^{29} and the multistep switching processes^{48} were clearly observed in the same tetragonal and rhombohedral PZT compositions. In order to resolve this paradox, in the current study, the compatibility of the IFM and the MSM models is investigated. It is shown that, despite multiple parallel and sequential switching processes with different switching times, the scaling property can be preserved under certain conditions. This is proved theoretically for ferroelectrics of tetragonal, rhombohedral, and orthorhombic phase symmetries and experimentally confirmed by an exemplary case of polarization response of an orthorhombic (K,Na)NbO_{3}-based ceramic described in Ref. 47. The theoretical analysis is presented in detail for the latter, the most complicated case, while the derivations for tetragonal and rhombohedral cases are presented in a concise form in Appendixes A and B, respectively.

## II. SCALING THEORY OF THE MULTISTEP POLARIZATION REVERSAL IN ORTHORHOMBIC POLYCRYSTALLINE FERROELECTRICS

In orthorhombic ferroelectrics with a high-temperature cubic parent phase, local spontaneous polarizations may be present in one of the 12 possible directions along the medial plane diagonals of the pseudo-cubic cell ( $ \u27e8 011 \u27e9$ directions). A sufficiently strong electric field opposite to the initial polarization direction causes a reversal of polarization. Let us choose a Cartesian frame with axes collinear with the main axes of a pseudo-cubic cell [Fig. 1(a)]. The initial polarization state is assumed to be $ P s( 0 , \u2212 1 , \u2212 1)/\u221a2$ (polarization of the amplitude $ P s$ pointing to A).

An electric field with an opposite direction drives the spontaneous polarization to the final state $ P s( 0 , 1 , 1)/\u221a2$ (polarization pointing to D). Polarization reversal may occur along different paths exemplarily depicted in Fig. 1(a), namely, by a direct 180° path A–D, by two sequential 120°- and 60°-polarization rotations A–G–D, by two sequential 60°- and 120°-polarization rotations A–B–D, and by a triple sequential 60°-polarization rotation A–B–C–D. We indicate the probabilities for the first switching steps as *η*_{1}, *η*_{2}, *η*_{3}, and *η*_{4} for processes starting with 60°, 120°, 90°, and 180° rotations, respectively, i.e., *η*_{1 }+*η*_{2 }+ *η*_{3 }+ *η*_{4 }= 1. Following the first 60°-polarization rotation, the two following paths are possible, a single 120°-polarization rotation with a weight *η*_{12} (such as the B–D path) and a double 60°-polarization rotation with a weight *η*_{11} (such as the B–C–D path), which satisfy *η*_{11 }+ *η*_{12 }= *η*_{1}. In a stochastic approach, the parameters *η*_{11}, *η*_{12}, *η*_{2}, *η*_{3}, and *η*_{4} are understood as average characteristics of the whole system. The switching probability of a specific path depends on mutual orientations of the domain and the applied electric field. Switching process fractions depend on the energy barriers for different switching paths as was analyzed by means of the Landau–Ginzburg–Devonshire theory in Ref. 46. The energy barriers, in turn, determine the activation fields characterizing different switching events. Expanding to the polycrystalline case is performed later by averaging over a statistical field distribution, which accounts for all possible domain orientations with respect to the electric field. The process fractions are determined, however, not only by the activation barriers for particular processes but also by the microstructure of the material such as grain sizes and shapes, their orientation, properties of the grain boundaries, etc., and, therefore, may be considered as independent fitting parameters.

^{47}At the beginning, we assume a single-domain single crystal with saturation polarization along $[ 0 1 \xaf 1 \xaf]$ [as in Fig. 1(a)]. When an opposite electric field is applied, local polarizations may undergo different switching processes with distinct nucleation rates per unit time and unit volume and develop growing switched domains of different dimensionalities. This results in time-dependent total polarization variation in the field direction,

^{35}the step-function approximation for the local polarization time dependences and apply this to the switching probabilities of particular switching paths derived in Ref. 47. Thus, for the contribution of 180°-switching events,

*β*

_{4}(related to the growing domain dimensionality) and switching time

*τ*

_{4}of this process. For the contribution of two-step 90°–90°-switching events,

*t*are accordingly approximated as

*β*

_{3},

*β*

_{33}and switching times

*τ*

_{3},

*τ*

_{33}of the first and the second switching steps. This results in the total contribution of this switching path,

Using the formulas (10)–(12) and (14) in the following, we will account for the fact that switching times of the sequential switching processes differ from each other by orders of the magnitude,^{47} $ \tau 33<< \tau 3, \tau 21<< \tau 2, \tau 1<< \tau 12, \tau 111<< \tau 11<< \tau 1$, so that the shorter times in additive expressions may be neglected.

^{34,35}we average the total polarization (6) over the statistical distributions of local switching times, which are assumed to follow from the statistical distribution of local electric fields. Thereby, the change from the time to the field scale is required. To this end, the formulas (7), (10)–(12), and (14) are transformed to

^{14}the threshold fields obey the same time dependence with distinct factors (activation fields):

^{47}

The scaling properties (1) and (2) are not expected to apply to each material and are subjected to a range of the mentioned assumptions concerning the width of the statistical distribution of local fields and the field dependence of switching times for different processes. However, they are expected to apply to a wide class of polycrystalline ferroelectrics that should be checked experimentally for each particular material. Indeed, similarly to the above considered ferroelectrics of orthorhombic symmetry, the materials of tetragonal and rhombohedral symmetries also exhibit dynamic scaling properties, though with different master curves $ \Phi (s)$, which are derived in Appendixes A and B, respectively. In the following, the scaling properties of an exemplary KNN-based ceramic^{47} will be investigated experimentally.

## III. EXPERIMENTAL PROOF OF SCALING PROPERTIES OF MULTISTEP SWITCHING

Bulk ceramics with nominal chemical composition (Na_{0.49}K_{0.49}Li_{0.02})(Nb_{0.8}Ta_{0.2})O_{3}, modified with 2 wt. % MnO_{2}, were processed as described in Ref. 47 using the conventional solid-state reaction route.^{49} X-ray powder diffraction revealing an orthorhombic structure and a standard electromechanical characterization of the samples combining large-signal polarization and strain measurements were presented in Ref. 47. The dimensions of the examined ceramic samples were 5 × 5 × 0.5 mm^{3} (length × width × thickness). Based on the compositional similarity between the materials in the current work and previous studies,^{49,50} the domain size and grain size were estimated to be tens to hundreds of nanometers and 1–2 *μ*m, respectively. Due to the polycrystallinity of the studied samples, it can be assumed that they were isotropic, and thus, free of texture before the application of an electric field.

Time-dependent polarization reversal was measured using the pulse switching setup and methodology described in Ref. 44. The sample was poled with an electric field of 4 kV/mm for 10 s, followed by a relaxation time of 30 s. The switching measurement was performed with a rectangular high voltage (HV) pulse, *E*_{a}. Then, the material was poled again at 4 kV/mm for 10 s, followed by a reference switching measurement at *E*_{a} for 10 s to remove the contribution of the sample's dielectric displacement and integrated leakage current from the switching measurement. A broad range of electric field pulses, from *E*_{a} = 0.15 to 1.30 kV/mm, which are approximately in the range of 0.3*E*_{C} < *E*_{a} < 2.6*E*_{C}, were applied to the sample, where *E*_{C} is the coercive field. The coercive field was determined to be 0.54 kV/mm from the macroscopic polarization hysteresis loops measured with a modified Sawyer–Tower circuit and a bipolar electric field of 4 kV/mm at 0.1 Hz.^{47}

Experimental data on time-dependent polarization switching are presented by symbols in Fig. 2(a) for exemplary field values. To study the scaling properties, the polarization variation data are first represented as a function of the applied field at different fixed poling time values in Fig. 2(b). All curves exhibit inflection points related to the maxima of the first polarization derivatives with respect to the field. These derivatives are presented in Fig. 2(c) again using time as a parameter. From Fig. 2(c), a position of maximum, $ E max(t)$, can be derived which is shown in Fig. 3. By calculating the normalized derivatives (1) for different times, using the experimental value ΔP_{max} = 36.1 *μ*C/cm^{2}, and their representation vs the applied field normalized to $ E max(t)$, a master curve $\Phi (u)$ for the studied KNN-ceramic is established in Fig. 2(d), thus confirming the scaling property of polarization response in this material.

The dependence of $ E max(t)$ in Fig. 3 is very well fitted by expression (18) derived from the Merz law,^{13} $ E max(t)= E 0 / [ ln \u2061 ( t / \tau 0 ) ] 1 / \alpha ,$using the parameters $ E 0=1.28 kV / mm$, $ \tau 0=8\xd7 10 \u2212 7 s$, and $\alpha =3.01$ with an overall error of about 0.2%. By implementing this dependence and the master curve of Fig. 2(d) in the numerical form, a theoretical polarization response according to Eq. (2) is evaluated for the respective applied field values and represented in Fig. 2(a) by solid lines. Theoretical curves mostly describe the experiment with high accuracy except for some deviations for the highest field values and shortest measured times. We note that the master curve $ \Phi (u)$ in formula (2) is an experimental one and consists of discrete experimental points of Fig. 2(c) presenting contributions to the normalized field derivative at different times. The shape of this curve is not approximated by splining or assuming a certain distribution type to avoid unnecessary additional hypothesis. Therefore, $ \Phi (u)$ is used as a set of discrete points when performing numerical integration in Eq. (2), which results in irregularities within the experimental accuracy.

Identification of further material characteristics is possible and will be performed in the future. Such microscopic characteristics as fractions of different switching paths $ \eta n$ and activation fields of particular switching events $ E A ( n )$ in the studied KNN material were roughly evaluated in Ref. 47 based on the hypothetic Lorentz distribution of local fields. Having the master curve of this material identified, the microscopic parameters can now be refined on a more rigorous basis. Thus, the process fractions and activation fields, as well as coefficients $ \gamma n$ can be extracted using available dynamic strain measurements.^{47} Identification of the statistical field distributions, however, still presents a challenge because it requires the solution of Eqs. (19), (A6), and (B6) with respect to $f(u)$ that is a nontrivial task.

For ferroelectric ceramics of tetragonal and rhombohedral phase symmetries, the expressions for total polarization, Eqs. (6) and (16), and the master curve $\Phi (u)$, Eq. (19), are different (see in Appendixes A and B), but the principal fact of scaling properties and the validity of Eqs. (1) and (2) remains. This is confirmed by experimental observations in Refs. 29 and 48.

## IV. CONCLUSIONS

The dynamic scaling properties of polarization response in multiaxial polycrystalline ferroelectrics, exhibiting parallel and sequential switching events, are predicted by the combined MSM-IFM model and exemplarily experimentally proven on the orthorhombic KNN-based ceramic. The analysis allows insight into microscopic switching mechanisms, represented by the observed Merz law for local polarization switching, on the one hand, and into macroscopic statistical properties of the polycrystalline media, represented by the master curve related to the weighted statistical distribution of electric field in them, on the other hand. This knowledge substantially improves the description of the temporal polarization response. Thus, the universal properties of polarization reversal, observed before in various polycrystalline ferroelectrics undergoing one-step polarization switching events, are extended to a wider class of multiaxial ferroelectrics undergoing multistep polarization processes.

The experiments performed over years on different polycrystalline ferroelectric materials and interpreted by means of the IFM model^{29–35} did not show an effect of the sample size as long as all samples were of macroscopic dimensions 0.3–3 mm. The domain size effect was not systematically studied. The grain size effect was studied for PZT of different compositions/phase symmetries and grain sizes in Ref. 29. The grain size has an impact on the maximum switchable polarization and activation fields and a less pronounced influence on the statistical field distributions. The accuracy of the model was comparably good for all grain sizes in the range (1–3) *μ*m.

The suggested MSM-IFM approach can be readily extended to further crystalline classes of ferroelectrics, undergoing thermally activated field-driven polarization switching, like hexagonal and monoclinic materials. Furthermore, the insight into multistep polarization switching kinetics is crucial for future ultrafast switching applications, and in particular, for multilevel digital information storage.

## ACKNOWLEDGMENTS

This work was supported by the Deutsche Forschungsgemeinschaft (DFG) Grant No. 405631895 (KO 5100/1-1). K.W. acknowledges the support of the National Nature Science Foundation of China (Grant Nos. 51822206 and 52032005).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Y. A. Genenko:** Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Writing – review & editing (equal). **S. Zhukov:** Data curation (equal); Methodology (equal); Visualization (equal). **M.-H. Zhang:** Data curation (equal); Investigation (equal); Methodology (equal); Visualization (equal). **K. Wang:** Data curation (equal); Funding acquisition (equal); Investigation (equal); Project administration (equal); Supervision (equal). **J. Koruza:** Conceptualization (equal); Data curation (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Supervision (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX A: SCALING PROPERTIES OF POLYCRYSTALLINE TETRAGONAL FERROELECTRICS

^{43}and applying the step-function approximation to them, the formula (A1) for single-crystalline tetragonal ferroelectrics is transformed to

*η*is the fraction of two-step 90°–90° switching events. Considering the fact that $ \tau 1<< \tau 2, \tau 3$ in the second theta-function, an approximation $t\u2212 \tau 1\u2212 \tau 2\u2245t\u2212 \tau 2$ can be adopted.

For establishing the 90°-process fraction *η* and constants $ \gamma n$, additional information is required, such as measurements of the time-dependent strain response.^{42–46} We note that, in the experimentally studied case of a tetragonal PZT ceramic,^{43} the relation $ \tau 2= \tau 3$, implying also $ E t h ( 2 )= E t h ( 3 )$ and $ \gamma 2= \gamma 3$, was observed which further simplifies expression (A6). We again emphasize that, for describing the polarization response using Eq. (2), the knowledge of process fractions and respective constants $ \gamma n$ is not required, as soon as scaling relation (1) is experimentally established as it was done in Ref. 29.

### APPENDIX B: SCALING PROPERTIES OF POLYCRYSTALLINE RHOMBOHEDRAL FERROELECTRICS

^{44}and applying the step-function approximation to them, the contributions to formula (B1) for single-crystalline rhombohedral ferroelectrics can be represented as

where, in notations of Ref. 44, $ \tau 1, \tau 11$, and $ \tau 111$ are the switching times of the first 71°, the second 71°, and the third 71° switching steps, respectively, $ \tau 2$ and $ \tau 21$ are the switching times of the first 109° and the following 71° switching steps, respectively, $ \tau 12$ is the switching time of the second 109° switching step after the first 71° one, and $ \tau 3$ is the switching time of the one-step 180° switching process. The parameters $ \eta 11, \eta 12, \eta 2$, and $ \eta 3$ present the fractions of the three-step 71°–71°–71°, the two-step 71°–109°, the two-step 109°–71°, and the one-step 180° switching paths, respectively. The sequential switching times typically differ by the orders of magnitude. In the case of the studied rhombohedral PZT ceramic,^{44} the following relations were established: $ \tau 1<< \tau 12$, $ \tau 2<< \tau 21$ and $ \tau 1<< \tau 11<< \tau 111$. This allows for the simplification of some arguments in Eq. (B2): $t\u2212 \tau 1\u2212 \tau 12\u2245t\u2212 \tau 12$, $t\u2212 \tau 2\u2212 \tau 21\u2245t\u2212 \tau 21$, $t\u2212 \tau 1\u2212 \tau 11\u2245t\u2212 \tau 11$, and $t\u2212 \tau 1\u2212 \tau 11\u2212 \tau 111\u2245t\u2212 \tau 111$.

For determining the process fractions *η _{n}* and constants $ \gamma n$, additional information is required, such as measurements of the time-dependent strain response.

^{42–46}We note that, in the experimentally studied case of a rhombohedral PZT ceramic,

^{44}the relations $ \tau 111= \tau 21\u2245 \tau 12= \tau 3$, implying also $ E t h ( 111 )= E t h ( 21 )\u2245 E t h ( 12 )= E t h ( 3 )$ and $ \gamma 111= \gamma 21= \gamma 12= \gamma 3$, were established which further simplifies the expression (B6). We again emphasize that, for describing the polarization response using Eq. (2), the knowledge of process fractions

*η*and respective constants $ \gamma n$ is not required, as soon as scaling relation (1) is experimentally established, as it was done in Ref. 29.

_{n}