We introduce an equivalent-circuit element based on the theory of interface pinning in random systems to analyze the contribution of domain wall motion below the coercive field to the impedance of a ferroelectric, as a function of amplitude $E0$ and frequency $f$ of an applied ac electric field. We demonstrate our model on a bulk $PbZrxTi1\u2212xO3$ (PZT) reference sample and then investigate capacitor stacks, containing ferroelectric $0.5(Ba0.7Ca0.3)TiO3$–$0.5Ba(Zr0.2Ti0.8)O3$ (BCZT) thin films, epitaxially grown by pulsed laser deposition on Nb-doped $SrTiO3$ single crystal substrates and covered with Au electrodes. Impedance spectra from $f=10$ Hz to 1 MHz were collected at different $E0$. Deconvolution of the spectra is achieved by fitting the measured impedance with an equivalent-circuit model of the capacitor stacks, and we extract for $E0=2.5$ kV/cm, a frequency-dependent permittivity of $\epsilon r\u2032(f)=458+7.3ln\u2061(1Hz/2\pi f)$ for the BCZT films from the obtained fit parameters. From an extended Rayleigh analysis, we obtain a coupling strength of 0.187 cm/kV between dielectric nonlinearity and dielectric dispersion in the BCZT films and identify different domain-wall-motion regimes. Finally, we construct a schematic diagram of the different domain-wall-motion regimes and discuss the corresponding domain-wall dynamics. Our approach can be utilized to replace purely phenomenological constant phase elements (CPEs) in modeling the impedance response of ferroelectrics and extracting material properties.

## I. INTRODUCTION

Domain wall motion in ferroelectric materials does not only influence their ferroelectric but also their dielectric and piezoelectric properties.^{1} In general, the complex dynamics of domain wall displacement in ferroelectrics results in a dependence of their dielectric $(\epsilon ij)$ and piezoelectric $(dijk)$ tensor components on the amplitude $E0$ and frequency $f$ of an applied ac electric field.^{1,2}

Regarding perovskite-type ferroelectrics, field and frequency dependent dielectric and piezoelectric response has been mostly investigated in ceramics and thin films based on $PbZrxTi1\u2212xO3$ (PZT).^{1–4} It was found that the field dependence of permittivity and piezoelectric coefficients at subswitching fields can often be described by the linear Rayleigh law, which results in a nonlinear and hysteretic dielectric polarization and piezoelectric strain. Based on this, the linear field dependence of permittivity is also referred to as dielectric nonlinearity.^{4} Moreover, the coupling of dielectric nonlinearity and frequency dispersion has been studied.^{2} However, a unified theory which describes the field and frequency behavior of material properties, e.g., the complex relative permittivity $\epsilon r(E0,f)$ of ferroelectrics, with the same set of equations is not available. Thus, the measurement of $\epsilon r(E0,f)$ is crucial since the permittivity is directly related to the microstructure of the material and essential for device modeling.^{1,5} This is especially true for thin film samples of new and, thus, less investigated lead-free ferroelectrics, such as $K0.5Na0.5NbO3$ (KNN) and $0.5(Ba0.7Ca0.3)TiO3$–$0.5Ba(Zr0.2Ti0.8)O3$ (BCZT).^{6,7} In terms of device applications, the motion of ferroelectric domain walls at subcoercive fields is the key mechanism for the realization of ferroelectric memristors based on ferroelectric tunnel junctions (FTJs).^{8,9}

Recently, a new method for the analysis of the field dependence of complex relative permittivity $\epsilon r(E0)$ due to the motion of ferroelectric domain walls at subcoercive fields (i.e., Rayleigh analysis) was introduced.^{10} This method is based on impedance spectroscopy and equivalent-circuit fitting, which allows one to discriminate different contributions to the permittivity, originating from electronically different components (electrodes, electrode/ferroelectric interfaces, and ferroelectric) of the analyzed thin film devices. This is important, because the interfaces to the electrodes play a more crucial role in thin films as compared to thick ceramic samples,^{10,11} which can lead to extrinsic dielectric relaxation processes,^{12} masking the contributions from the ferroelectric. The superiority of the Rayleigh analysis based on impedance spectroscopy over the conventional approach (i.e., without equivalent-circuit fitting) was demonstrated in Ref. 10 for the measurement of $\epsilon r(E0)$ for a ferroelectric Si:$HfO2$ thin film. However, the frequency dispersion $\epsilon r(f)$ is not accessible with the method described in Ref. 10 because the utilized circuit-element is not based on a physical theory which describes the pinning of ferroelectric domain walls. Note that also the widely utilized constant phase elements (CPEs) are phenomenological circuit elements and not based on a physical model, which makes it difficult to extract physically meaningful parameters from equivalent-circuit fits with CPEs.

In thin film capacitor stacks, with a ferroelectric placed between two metallic electrodes, the dielectric properties can exhibit strong frequency dispersion depending on the resistivity of both, the ferroelectric film and the electrode layers,^{13} thus masking the true dielectric dispersion of the ferroelectric caused by pinning of domain walls on randomly distributed defects occurring in the ferroelectric layer inside the heterostructure.^{14,15} Thus, a circuit-element is needed which model the ferroelectric impedance response and hence allows one to extract material properties from equivalent-circuit fitting.

In this work, we introduce an equivalent-circuit element with impedance $ZDW$ which is based on the theory of interface pinning in random systems^{14} and models the frequency-dependent impedance response of a ferroelectric including the domain-wall-motion induced contribution. Moreover, we report on a method to extract the domain-wall-motion induced amplitude- and frequency-dependent dielectric response $\epsilon r(E0,f)$ in thin film ferroelectrics, embedded in a thin film capacitor stack, from impedance spectroscopy measurements at different $E0$ and equivalent-circuit fits with the introduced element $ZDW$. Finally, we present an extended Rayleigh analysis which allows one to quantify the coupling strength between dielectric nonlinearity and frequency dispersion and the identification of different domain-wall-motion regimes. Our approach will be demonstrated by applying it to thin film capacitor stacks containing the lead-free ferroelectric BCZT and, as a reference sample, to a commercial soft bulk PZT ceramic.

## II. THEORY

### A. Rayleigh behavior and interface pinning in random systems

Domain wall dynamics below the coercive field $EC$ of a ferroelectric can be considered as a superposition of center of mass motion and relaxational motion of internal modes, i.e., domain wall segments.^{15} From theory it is known, that the interaction of domain walls with pinning centers, that form a random energy landscape, gives rise to the so called Rayleigh behavior and the corresponding Rayleigh law.^{16,17} Rayleigh behavior originates from the irreversible center of mass motion of domain walls and results in a linear dependence of permittivity and piezoelectric coefficients on the field amplitude $E0$ above a specific threshold field $ET$ and below the coercive field $EC$, which was also observed experimentally on ceramic and thin film samples.^{4,18–20} As a consequence of Rayleigh behavior, the dielectric polarization and piezoelectric strain responses of ferroelectrics become nonlinear in $E0$ and hysteretic even at subswitching fields $ET<E0<EC$. In the following, we will solely discuss dielectric properties and disregard piezoelectric properties.

The Rayleigh law for the real and imaginary parts of the complex relative permittivity $\epsilon r=\epsilon r\u2032\u2212i\epsilon r\u2033$ is given by Refs. 10 and 20,

Here, $\epsilon r,in\u2032$ and $\epsilon r,in\u2033$ are real and imaginary parts of the reversible (initial) Rayleigh parameter $\epsilon r,in$ due to lattice and reversible domain wall contributions, $\alpha \u2032$ denotes the irreversible Rayleigh parameter (Rayleigh constant), and $\alpha \u2033=43\pi \alpha \u2032$.

In the following, we will denote a general linear field dependence of $\epsilon r\u2032(E0)$ and $\epsilon r\u2033(E0)$ with positive slope as “Rayleigh-like behavior.” In that sense, the Rayleigh law is a special case of Rayleigh-like behavior with the additional requirement $\alpha \u2033=43\pi \alpha \u2032$.

The irreversible domain wall contribution to the real part of the permittivity is given by $\alpha \u2032E0$. It turns out that $\alpha \u2032$ is a direct and quantitative measure for the mobility of domain walls, and it is dependent on the concentration of pinning centers, the film thickness, and the density of domain walls.^{4} However, Rayleigh law Eq. (1) is rate-independent and no frequency dependence of material properties can be derived from it.^{2–4}

Theory predicts for domain wall pinning processes in random systems a logarithmic dependence of permittivity on frequency $f$, which is due to the relaxational contribution of internal modes,^{14,15}

Here, $\Theta $ is related to the roughness of the domain wall and $t0$ denotes a microscopic time constant, which is related to the rate of transitions of domain wall segments between metastable states.^{14,15}

In this pinned regime, the center-of-mass of the domain wall is captured in a potential valley, while domain wall segments can still jump between different metastable states with close energies. The potential valley of a random energy landscape can be described by Ref. 15,

with a constant $Cn$, the transverse spatial coordinate $r\xaf$, a single relaxation time $\tau 1$, and the random pinning force $g$. The first term in Eq. (3) describes the average (parabolic) shape of the well and the second term describes spatial fluctuations due to the residual part of the random pinning force $g$. This term is important because it results in a distribution of relaxation times $\Psi 1(\tau )$ instead of the single relaxation time $\tau 1$.^{15} As a consequence, the dielectric response of a ferroelectric due to domain wall pinning in a random energy landscape is not expected to exhibit a Debye-like dielectric relaxation with a single relaxation time and, hence, one dominant relaxation frequency.

In many ferroelectrics—including relaxors—it was found that $\Theta \u223c1$, which corresponds to a frequency-independent imaginary part of permittivity, and is valid for systems with a broad distribution $\Psi 1(\tau )$ of relaxation times.^{2–4,21–23} In this case, the logarithmic frequency dependence of the complex relative permittivity can be well approximated by the following expression:^{2,3,22,23}

where $\epsilon r,0\u2032$ and $\epsilon r,0\u2033$ denote static contributions to the permittivity and $\Delta \epsilon r\u2032$ is the logarithmic dispersion strength. We note that in Eq. (4), we fixed the value of $t0$ to 1 s, which is also the common choice in the literature.^{4,21–24} We will use this convention throughout the remaining part of this paper.

It has been shown for systems with Rayleigh-like behavior and logarithmic frequency dispersion that the frequency dependence of $\epsilon r\u2032$ in Eq. (4) is due to the logarithmic frequency dependence of both reversible ($\epsilon r,in\u2032$) and irreversible ($\alpha \u2032$) Rayleigh parameters defined in Rayleigh law Eq. (1), which can then be expressed as^{21}

The physical interpretation of the coefficients $\Delta \epsilon r,in\u2032$ and $\Delta \alpha \u2032$ in Eq. (5) will be discussed later in Sec. IV. As a consequence, the Rayleigh relation for the real part of relative permittivity $\epsilon r\u2032$ in Eq. (1) can be written with frequency-dependent Rayleigh parameters given by Eq. (5) as

However, Rayleigh law Eq. (1) comprises also the Rayleigh relation for the imaginary part of relative permittivity $\epsilon r\u2033$ and the relations are coupled by $\alpha \u2033=43\pi \alpha \u2032$. The coupling relation predicts a logarithmic frequency dependence of $\alpha \u2033$ which results in a frequency dispersion of $\epsilon r\u2033$, which is contradictory to the assumption of a frequency independent imaginary part given by Eq. (4), which implies a constant $\alpha \u2033$. Thus, the logarithmic frequency dispersion due to domain wall pinning in random systems and the Rayleigh law are only compatible at a fixed frequency. Taking account of a constant $\alpha \u2033$, we can implicitly define this frequency $f\u2261fR$ by

At $f=fR$, the Rayleigh law is fulfilled for variable field amplitudes $E0$ within the Rayleigh region. The frequency $fR$ will be called the “Rayleigh frequency” in the following.

The more general Rayleigh-like behavior exists in the frequency range where $\alpha \u2032(f)$ is positive. At frequencies $f$ within this Rayleigh-like frequency range, the center of mass of domain walls can move irreversibly via hopping between different minima of the energy landscape since this is a basic assumption in the derivation of the Rayleigh law. Thus, the center of mass of the domain walls is not captured in a potential valley as in the pinned regime and we conclude that the Rayleigh-like frequency range corresponds to a regime with coexisting center of mass motion and relaxational motion of internal modes. In analogy to Ref. 15, we will denote this regime as the stochastic regime of complex relative permittivity.

### B. Derivation of an equivalent-circuit element based on interface pinning in random systems

In an impedance spectroscopy experiment, the complex impedance $Z(f)=Z\u2032(f)\u2212iZ\u2033(f)$ is measured as a function of frequency $f$, typically in the range of $f=10\u2212106$ Hz. The complex impedance is related to the complex relative permittivity by Ref. 25,

with the geometrical capacitance $C0$. For a parallel-plate capacitor, $C0=\epsilon 0Ad$ with the vacuum permittivity $\epsilon 0$, electrode area $A$, and plate distance $d$. Separating real and imaginary parts in Eq. (8) yields

By inserting the logarithmic frequency dependence Eq. (4) into Eq. (9), we obtain the equivalent-circuit element $ZDW$, which models the impedance of a ferroelectric due to pinning of domain walls in a potential valley [cf. Eq. (3)] of a random energy landscape,

with $\epsilon r,0\u2032$, $\epsilon r,0\u2033$ and $\Delta \epsilon r\u2032$ as free parameters. Consequently, $ZDW$ will be called the “domain wall pinning element” in the following.

Impedance spectroscopy allows one to extract the free parameters $\epsilon r,0\u2032$, $\epsilon r,0\u2033$, and $\Delta \epsilon r\u2032$ by fitting the measured impedance over the full frequency range to the impedance of an appropriate equivalent-circuit model of the capacitor stack, including the ferroelectric. Thus, the different contributions (from electrodes, interfaces, and the ferroelectric layer) to the measured impedance of the capacitor stack are deconvoluted by equivalent-circuit fitting with the domain wall pinning element. The deconvolution of impedance data by equivalent-circuit fitting with the domain wall pinning element will be demonstrated below for a bulk PZT reference sample before we apply this to a capacitor stack which includes an epitaxial BCZT film.

Measurements at different excitation field amplitudes $E0$ then allow one to determine the field dependence of the dispersion parameters $\epsilon r,0\u2032(E0)$, $\epsilon r,0\u2033(E0)$, and $\Delta \epsilon r\u2032(E0)$. Subsequently, the field- and frequency-dependent dielectric response $\epsilon r(E0,f)$ of the ferroelectric is obtained by inserting the field-dependent dispersion parameters into Eq. (4). For comparison, the conventional approach for the analysis of complex relative permittivity measurements is described in Appendix A.

Note that for systems in which the assumption of a broad distribution of relaxation times and/or $\Theta \u223c1$ is not fulfilled, similar equivalent-circuit elements can be obtained by inserting the corresponding frequency dispersion into Eq. (9).

## III. EXPERIMENTAL

A ceramic BCZT target was used to grow 200-nm-thick epitaxial single crystalline films on (001)-oriented Nb-doped (0.5 wt%) $SrTiO3$ (Nb:STO) single crystal substrates (SurfaceNet GmbH, Rheine, Germany) by means of pulsed laser deposition (PLD). We used a 248 nm KrF excimer laser with an energy of 110 mJ and 5 Hz repetition rate. For the deposition on Nb:STO substrates, the temperature was fixed at $700\xb0$C with an oxygen partial pressure $pO2=0.1$mbar. After deposition, the films were heated up to $800\xb0$C and kept there for 15 min before cooling down to $700\xb0$C. Subsequently, $pO2$ was increased to 5 mbar and the films were cooled down to room temperature with a rate of 10 K/min.

The microstructure of the films was investigated by x-ray diffraction (XRD) to obtain the crystallographic orientation and to check for possible secondary phases. One sample was further analyzed by scanning transmission electron microscopy (STEM) to investigate possible crystal defects in BCZT, which may act as domain wall pinning-centers.

For electrical measurements, 3 mm diameter disk-shaped Au electrodes ($\u223c$100-nm-thick) were deposited on top of the BCZT films by electron beam evaporation through a shadow mask. Impedance spectra of the capacitor stack from $f=10$ Hz to 1 MHz were collected at different excitation field amplitudes $E0$ using a Solartron 1260 impedance analyzer together with a probe station. For each measurement, $E0$ was successively increased from 2.5 to 25 kV/cm (applied root mean square voltage $VRMS=0.05\u22120.5$ V) with a step-size of 1.25 kV/cm, which results in the measurement of the impedance under 19 different $E0$. Furthermore, impedance spectroscopy measurements were carried out on a 1.6 mm-thick disk-shaped commercial soft PZT ceramic with an electrode diameter of 8.06 mm (with a different impedance analyzer as utilized for the BCZT films).

## IV. RESULTS AND DISCUSSION

In this section, we present XRD data and results from electrical measurements from one representative epitaxial BCZT device. The STEM results were obtained from a different epitaxial BCZT device which was fabricated with the same process parameters (see Sec. III).

### A. XRD and STEM

Figure 1 shows the XRD $\theta $–2$\theta $ diffraction pattern for a BCZT thin film deposited on a (001)-oriented Nb:$SrTiO3$ (Nb:STO) single crystal. No traces of secondary phases are observable within the scanning range. However, we cannot rule out traces of secondary phases at higher angles of $2\theta $ between $60\xb0$ and $100\xb0$. The structural phase diagram of BCZT shows four possible crystallographic phases, i.e., cubic, tetragonal, orthorhombic, and rhombohedral.^{26} Due to the small variations in the lattice parameters of the possible phases, a rigorous specification of the crystallographic phase of the films is not accessible from the XRD data. Thus, we denote the film structure as pseudo-cubic, which represents the four possible phases. For simplicity, we use cubic metrics for the Miller indices and crystal directions. The film grown on a Nb:STO substrate is epitaxial and the full width half maximum (FWHM) of the rocking curve of the (002)-BCZT peak (the inset in Fig. 1) has a value of $\omega =0.08\xb0$, indicating high crystalline quality of the film. The results are in good agreement with BCZT thin films grown by PLD on undoped STO^{27} and $SrRuO3$ coated STO.^{28}

The epitaxial growth of BCZT on Nb:STO was also confirmed by STEM. Moreover, our STEM analysis revealed the occurrence of edge-type misfit dislocations at the epitaxial Nb:STO/BCZT interface, characterized by the Burgers vector $b=a[100]$. Here, $a$ denotes the lattice constant along the [100] direction within the BCZT unit cell (Fig. 2). These misfit dislocations act as possible domain wall pinning centers, which was already observed in similar epitaxial PZT heterostructures.^{29}

### B. Impedance spectroscopy and equivalent-circuit analysis

We first utilized impedance data, collected on a soft bulk PZT ceramic reference sample, to validate the domain wall pinning element $ZDW$ [cf. Eq. (10)]. The PZT reference sample is well-known to exhibit the characteristic logarithmic frequency-dispersion due to domain-wall pinning.^{21} Figure 3 shows the measured impedance response of the PZT reference sample in a Bode plot.^{30} First, the full equivalent-circuit depicted in the inset of Fig. 3 was used to fit the data.^{31} Here, the series resistance $R0$ corresponds to the resistance of the contacts and electrodes. The $R1C1$-element represents the interfaces to the electrodes. The resistance $R2$ and $ZDW$ model the ferroelectric PZT layer. $R2$ is attributed to losses due to possible leakage currents, whereas the losses due to domain wall motion are contained in $ZDW$. Additionally, the dispersion parameters are contained in the domain wall pinning element. The element $L0$ models the inductance of the cables used.^{30} The fitting procedure yields values for $R1$ and $C1$ in the range of $1011\Omega $ and $1015$ F, respectively, and the other fit parameters remain unaffected by setting $R1$ and $C1$ to infinity. This corresponds to a vanishing interface contribution to the impedance response, as expected for a mm-thick bulk sample.^{10} Additionally, no change in the fit parameters could be observed for setting $L0\u22610$. Fitting with the remaining equivalent-circuit of the PZT reference sample yields the following values for the fit parameters: $R0=77\Omega $, $R2=18G\Omega $, $\epsilon r,0\u2032=1400$, $\epsilon r,0\u2033=24.5$, and $\Delta \epsilon r\u2032=15.9$. Setting further $R0\u22610$ and $R2\u2261\u221e$, only changes $\epsilon r,0\u2033$ slightly (25.1 instead of 24.5). Thus, the full equivalent-circuit depicted in the inset of Fig. 3 can be reduced to $ZDW$. The obtained fit with $ZDW$ [cf. Eq. (10)] is depicted as solid lines in Fig. 3. The fit is in good agreement with the measured impedance. However, one might argue that it is not surprising that our model works with bulk PZT, as one can expect that in this kind of sample, the domain wall motion dominates over other extrinsic mechanisms. Hence, to further verify the correct extraction of material properties by domain wall pinning element modeling, we numerically added the impedance of circuit elements $R0=1000\Omega $, $R1=5000\Omega $, $C1=100$ pF, and $L0=1$ mH to the measured impedance of the PZT ceramic as depicted in the inset of Fig. 4(a). The added circuit elements represent extrinsic contributions to the impedance response, which may arise from contacts and interfaces.^{10,11,13} The added circuit elements result in a radically different behavior of the impedance at frequencies >1 kHz, as compared to the impedance measured on the PZT ceramic reference sample [cf. Fig. 4(a)]. Thus, converting the impedance data into complex permittivity data using the standard procedure, Eq. (8) also results in a radically different behavior of the permittivity as can be seen in Fig. 4(b). Here, the real part of the permittivity of the PZT sample with added elements even becomes negative at high frequencies, which is by no means a reflection of the true permittivity of the PZT reference sample. However, the extracted permittivity from the equivalent circuit fit to experimental data plus added elements is in good agreement with the measured permittivity of the PZT ceramic [cf. Fig. 4(b)], and the obtained fit parameters $\epsilon r,0\u2032=1400$, $\epsilon r,0\u2033=24.5$, and $\Delta \epsilon r\u2032=15.9$ are almost identical to the parameters obtained by fitting the impedance of the reference PZT sample with $ZDW$ in Fig. 3. Moreover, the fitting procedure yields fit parameters of $R0=1066$ $\Omega $, $R1=5040\Omega $, $C1=103$ pF, and $L0=1$ mH which are in very good agreement with the numerically added values. Altogether, the data presented in Fig. 4 demonstrate the correct deconvolution of different contributions to the measured impedance response by domain-wall pinning element modeling.

Collected impedance spectroscopy data (for a selection of three different $E0$) on the epitaxial BCZT thin-film device are shown in Fig. 5(a). A validation procedure to rule out artifacts in the measured impedance spectra is given in Appendix B. The equivalent circuit depicted in Fig. 5(b) [i.e., the same circuit as shown in the inset of Fig. 3(a)] was used to fit the collected impedance spectra. Results of the fits can be seen in the Bode plot shown in Fig. 5(a) (for three values of $E0$). In all cases, the fits are in good agreement with the measured data, indicating that the utilized equivalent circuit is a suitable model for the thin film capacitors. It should be noted that the full equivalent-circuit [cf. Fig. 5(b)] is necessary to fit the impedance data of the BCZT device in contrast to the PZT reference sample [cf. Fig. 3], which is due to the fundamental differences in the impedance response of thin-film and bulk samples.^{10} Fit parameters obtained from equivalent-circuit fits at three different excitation field amplitudes are summarized in Table I. The standard errors are calculated from maximum likelihood estimation.^{32} Note that the equivalent-circuit fitting procedure could also be utilized to fit impedance data measured with a fixed (small) excitation field amplitude $E0$ superimposed by a dc-bias voltage $VDC$. This measurement mode can be utilized, e.g., to analyze the properties of possible Schottky barriers at the ferroelectric–metal interfaces, by extracting the interfacial capacitance $C1$ from the equivalent-circuit fits at different dc-bias voltages. Here, the interfacial capacitance due to Schottky barrier formation is expected to be proportional to $(VDC+Vbi\u2032)\u22121/2$, with the apparent built-in potential $Vbi\u2032$.^{33} However, impedance spectroscopy and equivalent-circuit fitting does not reveal the electronic transport mechanism (e.g., thermionic injection, space-charge limited conduction) across metal–ferroelectric interfaces. Thus, a comprehensive analysis of transport properties across metal–ferroelectric interfaces also requires current–voltage (I–V) measurements.

E_{0} (kV/cm)
. | R_{0} (Ω)
. | R_{1} (Ω)
. | C_{1} (nF)
. | R_{2} (MΩ)
. | ɛ_{r,0}′
. | Δɛ_{r}′
. | ɛ_{r,0}″
. | (μH)
. |
---|---|---|---|---|---|---|---|---|

2.5 | 18 ± 1 | 2889 ± 4 | 150.5 ± 0.3 | 29 | 457.9 ± 0.2 | 7.27 ± 0.04 | 10.65 ± 0.03 | 4.8 ± 0.8 |

12.5 | 16 ± 1 | 2880 ± 3 | 151.8 ± 0.2 | 21 | 493.7 ± 0.1 | 8.41 ± 0.03 | 13.23 ± 0.02 | 5.3 ± 0.5 |

22.5 | 14 ± 3 | 2838 ± 7 | 155.1 ± 0.6 | 13 | 526.4 ± 0.4 | 11.41 ± 0.09 | 18.23 ± 0.07 | 5.3 ± 1.4 |

E_{0} (kV/cm)
. | R_{0} (Ω)
. | R_{1} (Ω)
. | C_{1} (nF)
. | R_{2} (MΩ)
. | ɛ_{r,0}′
. | Δɛ_{r}′
. | ɛ_{r,0}″
. | (μH)
. |
---|---|---|---|---|---|---|---|---|

2.5 | 18 ± 1 | 2889 ± 4 | 150.5 ± 0.3 | 29 | 457.9 ± 0.2 | 7.27 ± 0.04 | 10.65 ± 0.03 | 4.8 ± 0.8 |

12.5 | 16 ± 1 | 2880 ± 3 | 151.8 ± 0.2 | 21 | 493.7 ± 0.1 | 8.41 ± 0.03 | 13.23 ± 0.02 | 5.3 ± 0.5 |

22.5 | 14 ± 3 | 2838 ± 7 | 155.1 ± 0.6 | 13 | 526.4 ± 0.4 | 11.41 ± 0.09 | 18.23 ± 0.07 | 5.3 ± 1.4 |

Subsequently, the complex relative permittivity $\epsilon r(f)$ of the BCZT film (i.e., without contributions from the electrodes and ferroelectric/electrode interfaces) was calculated from the obtained fit parameters according to Eq. (4) for each field amplitude $E0$. Due to domain wall pinning element modeling, the extracted frequency-dependent permittivity of the BCZT layer is in agreement with the theory of interface pinning in random systems.^{14} This is different to previous studies on the dielectric dispersion in ferroelectrics, where purely phenomenological CPEs were utilized to extract permittivity data.^{11,12,34,35} The results for the extracted BCZT permittivity for a selection of three different excitation field amplitudes are depicted in Fig. 6. We clearly observe [cf. Fig. 6 and Table I] that $\epsilon r(f)$ and the parameters $\epsilon r,0\u2032$, $\epsilon r,0\u2033$, and $\Delta \epsilon r\u2032$ depend on the field amplitude $E0$, which indicate that dielectric nonlinearity and frequency dispersion are coupled in the BCZT films.

### C. Extended Rayleigh analysis

From the equivalent-circuit fits, the field dependence of the dispersion parameters $\epsilon r,0\u2032(E0)$, $\epsilon r,0\u2033(E0)$, and $\Delta \epsilon r\u2032(E0)$ can be extracted and the results are depicted in Fig. 7.

We observe that all dispersion parameters exhibit a linear dependence on the field amplitude $E0$, i.e., in the guise of the following linear equations:

Comparing the real part $\epsilon r\u2032(E0,f)$ of Eq. (12) with Eqs. (5) and (6) gives rise to the identifications $a0\u2032\u2261\epsilon r,in,0\u2032$, $b0\u2032\u2261\Delta \epsilon r,in\u2032$, $a\u2032\u2261\alpha 0\u2032$, and $b\u2032\u2261\Delta \alpha \u2032$. It should be noted that similar identifications were previously demonstrated to be valid for a PZT-based ferroelectric ceramic.^{21} However, the imaginary part of the material response was not considered in Ref. 21. Here, comparison of the frequency-independent imaginary part in Eq. (12), with the imaginary part of the Rayleigh law Eq. (1), yields the additional identifications $a0\u2033\u2261\epsilon r,in\u2033$ and $a\u2033\u2261\alpha \u2033$.

Altogether, the field dependence of the dispersion parameters Eq. (11) can then be written in the following form:

where the coefficient $\Delta \epsilon r,in\u2032$ is the zero-field contribution to the logarithmic dispersion strength and the field-dependent contribution is given by $\Delta \alpha \u2032E0$.

Subsequently, linear fits according to Eq. (13) were carried out in the range $6.25kV/cm\u2264E0\u226416.25kV/cm$ as indicated by the vertical dashed lines in Fig. 7. The starting point for the linear fits at $E0=6.25kV/cm\u2261ET$ was determined by the criterion that all subsequent data points are required to have a higher value compared to their left neighbor. The end point for the linear fits was then chosen by the longest fitting range for which every linear fit has a correlation factor $R2\u22650.980$, which is fulfilled for $E0=16.25kV/cm\u2261EPS$. At higher fields, all three curves in Fig. 7 show an upward curvature which we attribute to the partial switching of domains.^{20} For the linear fits, the estimated standard errors from the previous equivalent-circuit fitting procedure were taken into account, which is also indicated by the error bars within the fitting range [cf. Fig. 7]. The results of the fits according to Eq. (13) are summarized in Table II.

ɛ_{r,in,0}′
. | α_{0}′ (cm/kV)
. | ɛ_{r,in}″
. | α″ (cm/kV)
. | Δɛ_{r,in}′
. | Δα′ (cm/kV)
. |
---|---|---|---|---|---|

463.9 ± 0.1 | 2.41 ± 0.01 | 9.31 ± 0.02 | 0.318 ± 0.002 | 6.09 ± 0.03 | 0.187 ± 0.003 |

ɛ_{r,in,0}′
. | α_{0}′ (cm/kV)
. | ɛ_{r,in}″
. | α″ (cm/kV)
. | Δɛ_{r,in}′
. | Δα′ (cm/kV)
. |
---|---|---|---|---|---|

463.9 ± 0.1 | 2.41 ± 0.01 | 9.31 ± 0.02 | 0.318 ± 0.002 | 6.09 ± 0.03 | 0.187 ± 0.003 |

From the extracted fit parameters, the functional form of the frequency-dependent reversible and irreversible Rayleigh parameter of the BCZT thin film can be immediately obtained according to Eq. (5) and the results are depicted in Fig. 8. Here, the vertical dotted line at$fR\u223c1$ kHz indicates the Rayleigh frequency for the BCZT thin film, which we calculated according to Eq. (7). As discussed in Sec. II, $fR$ is part of the Rayleigh-like frequency range, where the center of mass of the domain walls can move irreversibly via hopping between potential minima corresponding to the stochastic regime.

The dashed line at $f0\u223c63$ kHz marks the region where the irreversible Rayleigh parameter $\alpha \u2032$ is zero and changes sign to negative values. This might indicate a transition from the stochastic regime to a pinned regime, in which the center of mass motion of the domain walls cannot follow the fast driving field, and, thus, the domain walls are captured in a potential valley.^{15} Similar experimental indications for the freezing of the irreversible center of mass motion of the domain walls at higher frequencies were reported previously.^{24}

By inserting the functional form of the field dependent frequency dispersion parameters [cf. Fig. 7 and Table II] into Eq. (4), we finally obtain the functional form of the field and frequency dependent dielectric response in the epitaxial BCZT thin film, which is given by

Note that $\epsilon r,in\u2033$ in Eq. (14) is not zero, which could be an indication for an internal bias or stress field.^{10}

The coupling between field dependence and frequency dispersion in the BCZT thin film arises explicitly in Eq. (14) as the product of $E0$ and the $f$-dependent term. Thus, this gives rise to the identification of the coefficient $\Delta \alpha \u2032$ as the coupling strength between dielectric nonlinearity and frequency dispersion, which is defined as the change of logarithmic frequency dispersion strength with field amplitude—or equivalent—as the change of the irreversible Rayleigh parameter with the logarithm of frequency, i.e.,

which has a value of 0.187 cm/kV in the epitaxial BCZT thin film. In fact, the negative values of the irreversible Rayleigh parameter are a direct consequence of the coupling between frequency dispersion and dielectric nonlinearity in the BCZT thin film.

Altogether, our analysis yields the schematic diagram depicted in Fig. 9 for the different domain wall motion regimes in the epitaxial BCZT thin film. Note that the diagram constructed in Fig. 9 is not a generic diagram for any epitaxial BCZT film. In fact, it is a schematic diagram of different domain-wall motion regimes, which can now be identified within the framework of the extended Rayleigh analysis.

At low field amplitudes $0<E0<ET$, the center of mass of the domain walls cannot overcome the potential well Eq. (3) of the random energy landscape; however, domain wall segments can jump between metastable states with close energies, which corresponds to the non-coupled pinning regime. This observation has similarities to what has been found theoretically for disordered ferromagnets.^{36} By increasing the field amplitude above the threshold field $ET$, the center of mass of the domain walls can additionally jump between different potential minima, resulting in a coexistence of irreversible center of mass motion and relaxational motion of internal modes in the stochastic regime. Note that this coexistence is not simply additive due to the coupling of dielectric nonlinearity and frequency dispersion. As a consequence of coupling, the irreversible Rayleigh parameter becomes zero at $f0\u223c63$ kHz and changes to negative values in the pinned coupling regime. Note that in this coupled-pinning regime, domain wall segments can still jump between metastable states, which is indicated by the logarithmic frequency dispersion observed in this regime. The domain wall dynamics in the coupled-pinning regime is not fully understood and remains as a task of a unified theory of dielectric nonlinearity and frequency dispersion, which is not yet available.

## V. CONCLUSION

In conclusion, this work demonstrates several key aspects. We introduce the domain wall pinning element $ZDW$ based on the theory of interface pining in random systems, which includes the characteristic logarithmic frequency dispersion due to domain wall pinning in a random energy landscape and hence models the impedance response of a ferroelectric below the coercive field. The domain wall pinning element can be utilized to replace purely phenomenological CPEs in modeling the impedance response of ferroelectrics and extracting material properties. Moreover, we demonstrate the practical application of this new element to extract the field- and frequency-dependent dielectric response in an epitaxial lead-free ferroelectric BCZT thin film embedded in a capacitor stack and we reveal its superiority over the conventional approach. In addition, we perform an extended Rayleigh analysis which results in the quantification of the coupling strength and the explicit functional form $\epsilon r(E0,f)$ of the coupled dielectric response in the BCZT film. Finally, we present a schematic diagram of the different domain wall motion regimes in the BCZT film and discuss the corresponding domain wall dynamics.

The present work is intended to serve as a guideline for future work on ferroelectric materials, which includes to reveal the effect of temperature and phase transitions,^{37} external bias fields,^{4} dopants,^{23} ion-bombardment,^{38,39} and substrate clamping^{24} on domain wall dynamics and related material properties. In terms of ferroelectric devices, our approach may help to characterize and optimize the response of ferroelectric memristors, which is also of practical importance for future interdisciplinary applications of ferroelectrics in bioelectronics, such as ferroelectric microelectrodes^{40} in hybrid neurocomputational systems.

## ACKNOWLEDGMENTS

We acknowledge Dragan Damjanovic for providing data on a ceramic reference sample and gratefully acknowledge fruitful discussions with J. Driscoll (Cambridge), C. Warres (NMI), C. Hofer, and J. C. Meyer and technical support by M. Turad and R. Löffler (LISA$+$). This work was partially funded by the Bundesministerium für Bildung und Forschung (BMBF) under Grant No. 13GW0123E and by the Europäische Fonds für regionale Entwicklung (EFRE) under Grant No. 712303.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Maximilian T. Becker:** Conceptualization (lead); Data curation (lead); Formal analysis (lead); Writing – original draft (lead); Writing – review and editing (equal). **Claus J. Burkhardt:** Funding acquisition (lead); Project administration (lead); Resources (equal). **Reinhold Kleiner:** Formal analysis (supporting); Writing – review and editing (supporting); Resources (equal). **Dieter Koelle:** Formal analysis (supporting); Writing – original draft (supporting); Writing – review and editing (equal); Resources (equal).

## DATA AVAILABILITY

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

### APPENDIX A: CONVENTIONAL APPROACH

For the conventional approach, the measured complex impedance $Z=Z\u2032\u2212iZ\u2033$ of the capacitor stacks is converted into complex permittivity data of the capacitor stacks using the formalism described in Ref. 25, resulting in

The resulting $\epsilon r,stack(f)$ plots are shown in Fig. 10 for a selection of three different $E0$. For comparison, the extracted complex permittivity $\epsilon r(f)$ of the BCZT layer from the equivalent-circuit fits [cf. Fig. 6 and Eq. (4)] is depicted. The capacitor stack exhibits a Debye-like dielectric relaxation at high frequencies ($\u223c105$ Hz) due to resistive losses in the electrodes.^{13} Moreover, the dielectric spectrum of the epitaxial thin film capacitor indicates an additional Debye-like relaxation process at lower frequencies ($\u223c100$Hz). We attribute this second relaxation process to an extrinsic electrode/film interface effect, which was previously reported to occur in other ferroelectric thin films grown on Nb:STO substrates.^{12} This is also consistent with the theoretical discussion of potential well Eq. (3) within a random energy landscape, which leads to a distribution of relaxation times in the BCZT layer and hence the BCZT layer is not expected to exhibit a Debye-like dielectric relaxation which corresponds to a single relaxation time and hence one dominant relaxation frequency $fr$. The relaxation frequency can be calculated from the obtained interface fit parameters $R1$ and $C1$ according to the relation $fr=1/2\pi R1C1$. The calculated relaxation frequencies are in agreement with the peaks in the imaginary part of permittivity of the capacitor stack [cf. Fig. 10] at frequencies ($\u223c$100 Hz), which indicates the correct deconvolution of film and interface contributions by our equivalent-circuit fits. Moreover, for BCZT thin films grown on platinized Si substrates, the second Debye-like dielectric relaxation was absent.^{41} Furthermore, we note that the slightly negative values of $\epsilon r,stack\u2032(f)$ at frequencies around 1 MHz are due to the inductance of the cables used.^{30}

It is clear that the extraction of frequency dispersion of permittivity from [cf. Eq. (A1) and symbols in Fig. 10] does not reflect the true dielectric dispersion of the BCZT thin film since there is no deconvolution of electronically distinct components forming the capacitor heterostructure. The same is true for the extraction of Rayleigh parameters, which has been recently demonstrated.^{10}

### APPENDIX B: VALIDATION OF IMPEDANCE SPECTRA

Impedance spectra might be afflicted with artifacts. To rule out artifacts introduced by experimental setup conditions in the measured impedance spectra, we use the two-pole Hilbert transform (Z-HIT) described in Ref. 42. This approach is beneficial compared to methods based on Kramers–Kronig (KK) relations^{43,44} since the KK relations are strictly defined within the frequency range between zero and infinity and hence, the extrapolation of the measured impedance data is a unavoidable task which may lead to erroneous results.^{42} The Z-HIT avoids the necessity of an extrapolation of the frequency range and allows one to reconstruct the modulus of the impedance $|Z|$ vs $f$ from the measured values of the less artifact-prone^{10,42} phase angle. The results of the Z-HIT for three different excitation field amplitudes $E0$ are depicted in Fig. 11. In all cases, the reconstructed curves match fairly well with the measured data, which indicates that the measured impedance spectra are not afflicted with artifacts and the investigated system is stable, causal, and (sufficiently) linear^{10} and exhibits finite values over the entire frequency range.