Local Probe Comparison of Ferroelectric Switching Event Statistics in the Creep and Depinning Regimes in Pb(Zr_{0.2}Ti_{0.8})O_{3} Thin Films.

Ferroelectric materials provide a useful model system to explore the jerky, highly nonlinear dynamics of elastic interfaces in disordered media. The distribution of nanoscale switching event sizes is studied in two Pb(Zr_{0.2}Ti_{0.8})O_{3} thin films with different disorder landscapes using piezoresponse force microscopy. While the switching event statistics show the expected power-law scaling, significant variations in the value of the scaling exponent τ are seen, possibly as a consequence of the different intrinsic disorder landscapes in the samples and of further alterations under high tip bias applied during domain writing. Importantly, higher exponent values (1.98-2.87) are observed when crackling statistics are acquired only for events occurring in the creep regime. The exponents are systematically lowered when all events across both creep and depinning regimes are considered-the first time such a distinction is made in studies of ferroelectric materials. These results show that distinguishing the two regimes is of crucial importance, significantly affecting the exponent value and potentially leading to incorrect assignment of universality class.

Many systems exhibit crackling, responding to a slowly increasing external applied force by discrete jerky events, which span a broad range of sizes and released energies, and where an initial event can trigger much larger ones in an avalanche phenomenon [1][2][3][4].
These systems typically exhibit a power-law distribution PðSÞ ∼ S −τ f S ðS=S c Þ of the size (S) of these jerky events, with characteristic exponent τ, and a sharp cutoff function f S for events above a particular system-dependent threshold. The size refers to the surface spanned by individual events.
Similar power laws describe the statistics of other observables such as the energy released during the events, with characteristic exponent ϵ and cutoff function f E . Theoretical models typically provide predictions for the values of both exponents, allowing independent measurements of different observables to be related [3,5,6]. While the statistics of crackling systems can often be described at small scales with simple system-specific models, the characteristic power-law exponents show scale invariance and universal behavior independent of the microscopic details of the system. Systems with different small-scale physics can share the same exponents, making insight gathered from an experimentally accessible system directly applicable to less tractable ones within the same universality class. For example, magnetization reversal in soft magnets and the seismic activity of earthquakes fall in the mean-field interface depinning universality class [7], sharing the same large scale statistical properties. This independence of the large scale properties with respect to microscopic system details is also technologically relevant, allowing nondestructive materials tests [8]. It is therefore important to determine the conditions under which individual systems fall in particular universality classes.
Elastic models further detailed in the Supplemental Material [9], Sec. I, including Refs. [3,[10][11][12][13][14][15][16][17][18][19][20][21][22][23][24], describe a wide variety of systems in which crackling behavior has been observed, such as fluid contact lines [25,26], slip faults in earthquakes [7,27], and the dynamics of domain walls in soft magnets [28][29][30][31][32][33][34]. In these models, interfaces are described as elastic manifolds pinned by disorder in a random medium [22,23]. Their complex, nonlinear dynamics show a transition from a thermally activated creep regime at low driving force to a depinning regime above a critical force f c , with jerky motion and crackling behavior [35]. In both regimes, the distribution of individual event sizes is expected to follow a power law. However, recent theoretical studies suggest that the two regimes present different characteristic exponents. They predict that in the creep regime, switching events trigger aftershocks and are correlated in space and time, while in the depinning regime, no aftershocks occur and the switching events are spatially uncorrelated [18].
Ferroelectric materials provide a versatile system for studying crackling phenomena, creep, and depinning. These materials are characterized by stable electric polarization states, organized into domains, switchable under an appropriate applied electric field [36]. Various material defects act as pinning and nucleation sites, affecting the motion of ferroelectric domain walls and the size of the jerky switching events [37,38]. The characteristic powerlaw exponents of the energy distribution of these events have been studied in different ferroelectrics using acoustic noise measurements [15], through their jerky switching currents [39], and simultaneously electrically and through optical birefringence imaging [40], showing energy exponents of ϵ ≈ 1.6.
An alternative approach is to use piezoresponse force microscopy (PFM) [41] to access ferroelectric switching dynamics with a resolution of ∼10 nm-at scales inaccessible optically, and maintaining spatial information lost in electrical or acoustic techniques. Such measurements have confirmed creep motion of ferroelectric domain walls in thin films up to a critical field E c [42,43] and have distinguished the creep and depinning regimes [44,45]. Scanning probe tips can also be used to inject or redistribute defects [46] during switching.
In this work, the power-law scaling of the distribution of switching event sizes is studied using PFM in two films of PbðZr 0.2 Ti 0.8 ÞO 3 henceforth labeled PZT-Nuc and PZT-Mot, showing very different disorder landscapes, distinguishing for the first time in ferroelectrics between events occurring in the creep and depinning regimes. In both samples on as-grown and pre-polarized regions, the size distributions of switching events in the creep regime are shown to follow power-law behavior over two decades of event sizes with no visible cutoffs, but with different characteristic size exponent τ ranging from 1.98 AE 0.09 to 2.87 AE 0.12. Moreover, when events occurring above the critical tip bias are also included, the exponent values are systematically lowered to between 1.81 AE 0.05 and 2.56 AE 0.1. The power-law character of the distribution of event sizes is robust against the intrinsic differences in defect landscapes between the samples, as well as to alterations therein triggered by high tip bias during the domain pre-polarizing process.
Two ∼70 nm thick epitaxial films of PbðZr 0.2 Ti 0.8 ÞO 3 with SrRuO 3 back electrodes on (001)-oriented SrTiO 3 substrates were studied. Both samples display an as-grown monodomain polarization oriented normal to the film plane in the direction of growth (up). As detailed in the Supplemental Material [9], Sec. III including Ref. [47], the two samples show high crystalline and surface quality, but varying morphology, relaxation of biaxial strain, and Pb stoichiometry, thus presenting very different defect landscapes. The resulting variations in domain wall pinning can be seen in Fig. 1(a) with noticeably different roughening of domains written under similar conditions in the two films.
Domains of opposite (down) polarization were prepolarized with a high positive tip bias. Stroboscopic measurements gradually switching the polarization in both up and down domains were then performed, in which switching scans with a defined subcoercive dc bias, incremented at each scan by a fixed interval, were alternated with PFM scans imaging the evolution of the domain structure. Each switching scan induced a series of switching events, mapped as the differences between the preceding and successive PFM images. Maps of the tip bias triggering the switching events throughout the measurement series are shown in Fig. 1, and allow the individual event sizes (see Supplemental Material [9], Sec. II) to be extracted.
To analyze crackling statistics in the creep regime, its tip bias range must be determined for the different samples and bias polarities. We extract the effective domain wall displacement by calculating the equivalent disc radius of the total area switched at each tip bias, as shown in Fig. 2. At low bias, the displacement remains below the noise threshold until the driving force is sufficiently high, then increases nonlinearly, defining the lower tip bias cutoff for the creep regime. This threshold is more clearly visible in logarithmic scale, shown in the Supplemental Material [9], Sec. VII.
For the upper tip bias cutoff, since the creep-to-depinning transition is characterized by a sharp upturn of the velocity response to the driving force [20], the intercept of the region of rapid displacement provides a lower bound for the depinning threshold [44], indicated by star markers in Figs. 2(b)-2(d). This is by construction an underestimate. While in the case of Fig. 2(c), only the points at largest bias were used to determine the upper tip bias cutoff in order to include enough statistics for power-law fitting, the resulting estimate is therefore still expected to be lower than its true value.
In the positive tip bias case, switching from the preferred as-grown P up state, the displacement appears to follow a typical creep-to-depinning transition. Under negative bias, however, where the switching occurs in pre-polarized domains with modified defect landscapes, the displacement shows domelike features [in Figs. 2(a), 2(c)], possibly related to different pinning hierarchies. The sudden drop in displacement inside this dome in PZT-Mot is likely a result of surface contaminants temporarily adhering to the tip and changing the effective field applied to the sample.
For PZT-Nuc, where switching is dominated by multiple nucleation sites each giving rise to a separate growing domain, the decreasing displacement after 4.0 V is likely caused by the lack of available unswitched area. Here, the same voltage window is used as for the positive tip bias case.
The distribution of event sizes are then extracted from the switching maps of Figs. 1(b)-1(e), and fitted with a power law following Clauset et al. [57]. In all four scenarios, the event size statistics present the expected power-law scaling, as shown in Fig. 3.
In the case of PZT-Mot under negative tip bias, where the switching occurs exclusively via domain wall motion, successive passes of the biased tip each cause individual polarization reversal events that are connected from line to line. This leads by the end of the switching scan to an overall reversal spanning the entire length of the domain wall, thus recorded as a single event in Fig. 1(d). To extract meaningful event statistics, we divided such large events into smaller units or boxes of a given width.
Power-law fitting was performed for a range of box widths and the final exponent obtained by averaging values for fits passing set quality criteria. The validity of this procedure was verified on the PZT-Nuc V tip > 0 series, with sufficient statistics to compare the size exponents extracted both directly and using boxing. The two techniques yield similar exponent values. The resulting powerlaw exponents are shown in Fig. 3(b). All fitting procedures and fit characterizations are detailed in the Supplemental Material [9], Sec. VIII including Refs. [42,57,58].
In the four scenarios, we find exponent values varying between τ ¼ 1.98 and 2.87. For PZT-Nuc V tip > 0, the size exponent of 1.98 is compatible with field-integrated meanfield models predicting τ ¼ 2 [6]. This model also describes characteristic exponents of other observables such as the energy ϵ released during switching, which was measured using switching currents, parallel-plate PZT capacitors, and in acoustic measurements in BaTiO 3 [15,39,40]. In PZT-Nuc V tip < 0, included for completeness over a bias window equivalent to that of PZT-Nuc V tip > 0, we observe a higher value of τ ¼ 2.85, reflecting a prevalence of smaller-sized switching events. Similarly, in PZT-Mot for V tip > 0 and V tip < 0 we find exponent values of τ ¼ 2.87 and 2.57, respectively.
Moreover, since previous studies of avalanche dynamics in ferroelectric and ferroelastic systems did not distinguish events occurring in different dynamics regimes [15,39,40], we extracted for comparison the scaling exponent values when events in the depinning regime were included in the analysis. As can be seen in Table I, this significantly lowers the exponent values, reflecting the larger size of switching events occurring beyond the upper bias cutoff.
To further compare avalanche statistics in the creep and depinning regimes, we next explored switching in PZT-Nuc under a scanning tip biased at constant voltage, with switching maps for both regimes shown in the Supplemental Material [9], Sec. IX. Since progressive tip contamination by surface adsorbates and defect redistribution during pre-polarization can influence the effective applied field, the tip biases of −3.5 and −5.0 V selected to correspond to these regimes cannot be directly compared with the measurements at incrementally increasing bias.
At −3.5 V, we observe very slow switching dynamics, with few, stochastically occurring nucleation events. Switching events remain small, with displacements of up to a few pixels only. The newly nucleated domains grow logarithmically with time, as shown in Fig. 4(a). At −5.0 V, domain wall motion is linear with time and characterized by much larger displacements, as shown in Fig. 4(b), rapidly leading to merging of the individual domains. The slow logarithmic domain growth at −3.5 V can be seen as a marker of creep dynamics, with a low driving force compared to the energy barriers of the disorder landscape. In this scenario, only small local displacements aided by thermal activation take place, with pinning severely constraining any large scale reconfiguration. At high tip bias, the growing domains expand at a constant rate, suggesting the applied field strength exceeds the characteristic barrier height and is fully in the depinning regime.
From the distribution of event sizes at −5.0 V, we obtain an exponent of τ ¼ 2.10 AE 0.05 over 2.5 decades of event size. At −3.5 V, however, we find much higher exponent values of τ ¼ 3.36 AE 0.16, emphasizing the role of very small switching events, although we note that the fitting is more difficult, and confidence in the exponent value is far lower in this case (see Supplemental Material [9], Sec. VIII). The values of the extracted size exponents are summarized in Table I. Furthermore, we observe evidence of spatial clustering of switching events in the creep regime, that are decorrelated above the critical bias as has been predicted in elastic models [18]. Details can be found in the Supplemental Material [9], Sec. X.
To conclude, we observe markedly different switching dynamics in PbðZr 0.2 Ti 0.8 ÞO 3 films with different defect landscapes established during sample growth, and discriminate between creep and depinning regimes as a function of the applied tip bias. While the switching event statistics all show power-law scaling, we find significant variations in the value of the scaling exponent τ, with higher values ranging from 1.98 AE 0.05 to 2.87 AE 0.12 in the creep regime. Exponents are systematically lower, with values ranging from 1.81 AE 0.05 to 2.56 AE 0.1 when switching events occurring during the entire tip bias window are included. Furthermore, the characteristic size exponents extracted from measurements of switching at constant tip bias in the creep and depinning regimes show higher values in the former. For both regimes, these exponents (3.36 AE 0.16 and 2.1 AE 0.05, respectively) are significantly higher than elastic model predictions. The systematic lowering of the characteristic size exponent during depinning could potentially be used to identify the different regimes, and suggests that studies which do not carefully discriminate between them can lead to significant error in the extracted exponent values and potentially an incorrect assignment of universality class. These conclusions should generalize to all systems described by elastic manifold models.  The switching event sizes cover only a very small surface range and confidence in this particular exponent is low.