Maximizing Information: A Machine Learning Approach for Analysis of Complex Nanoscale Electromechanical Behavior in Defect‐Rich PZT Films

Scanning Probe Microscopy (SPM) based techniques probe material properties over microscale regions with nanoscale resolution, ultimately resulting in investigation of mesoscale functionalities. Among SPM techniques, piezoresponse force microscopy (PFM) is a highly effective tool in exploring polarization switching in ferroelectric materials. However, its signal is also sensitive to sample‐dependent electrostatic and chemo‐electromechanical changes. Literature reports have often concentrated on the evaluation of the Off‐field piezoresponse, compared to On‐field piezoresponse, based on the latter's increased sensitivity to non‐ferroelectric contributions. Using machine learning approaches incorporating both Off‐ and On‐field piezoresponse response as well as Off‐field resonance frequency to maximize information, switching piezoresponse in a defect‐rich Pb(Zr,Ti)O3 thin film is investigated. As expected, one major contributor to the piezoresponse is mostly ferroelectric, coupled with electrostatic phenomena during On‐field measurements. A second component is electrostatic in nature, while a third component is likely due to a superposition of multiple non‐ferroelectric processes. The proposed approach will enable deeper understanding of switching phenomena in weakly ferroelectric samples and materials with large chemo‐electromechanical response.


Introduction
Scanning Probe Microscopy (SPM) based techniques employ a cantilever to locally probe a materials surface and allow for the acquisition of multiple observable quantities (e.g., The efficacy of PFM is greatly enhanced by resonant modes, which exploit the contact resonant frequency (ω) to increase signal-to-noise ratio and reduce topographic crosstalk. [4,5] Polarization dynamics can be explored with the application of an additional dc field via relaxation or switching experiments. [6,7] Measurements are made both during (On-field, On-F) and immediately after (Off-field, Off-F) the application of the dc bias. The resultant data is subsequently fit to a simple harmonic oscillator model (SHO) to extract information about electromechanical response (amplitude, A), polarization orientation (phase, θ), viscoelastic properties (ω), and energy dissipation (Q).
In prototypical ferroelectrics such as many perovskite oxides, electromechanical deformation is expected to be dominated by true piezoelectric signal. However, as described by Balke, Vasudevan and others, electromechanical deformations can also be induced by non-piezoelectric phenomena such as electrostatic displacement, electrochemical deformation, electrostrictive effects and Joule heating: in fact, non-piezoelectric (i.e., electrostatic or electrochemical) contributions usually dominate the On-F piezoresponse (PR). [8][9][10][11][12] As a result, On-F measurements are typically ignored in the analysis of PFM data, due to the inability to separate or remove the non-piezoelectric contributions. Additionally, instrument-based artifacts can also influence both On-F and Off-F measurements. [13] Multiple different PFM or SPM-based approaches have been reported to separate the piezoelectric and non-piezoelectric contributions to the measured signal. [14][15][16][17][18][19] Some of these involve repeated data collection across multiple specialized experimental techniques (e.g., c-KPFM) or re-scan of the sample with cantilever tips of different stiffness), leveraging differential changes in piezoelectric and non-piezoelectric contributions to distinguish between contributors. An approach for such differentiation without the need for repeated data collection would benefit the analysis of all voltage-modulated SPM variants, including PFM, c-KPFM, ESM, etc.
Machine learning (ML) techniques have recently emerged as a viable method for analysis of PFM and other voltage modulated SPM measurements, addressing some of the above challenges in data interpretation. [1,20] Among others, clustering techniques such as K-means have been used to identify the ergodic and non-ergodic behaviors in ceramic composites, as well as the onset of phase transitions in relaxor-ferroelectric solid solutions. [21,22] A Bayesian optimized support vector ML approach was used to classify switching and stable domains under weak ac electric fields. [23] Neural networks have been leveraged to study polarization dynamics in PZT thin films and competing contributors to the functional response in conductance AFM-based investigation of (Er 0.99 ,Zr 0.01 )MnO 3 . [24,25] Yet, a major limitation of ML techniques is their inherent lack of physical constraints, often leading to debatable physical conclusions. [26] Encoding of physical constraints can be achieved through stacking (concatenation) of data along specific dimensions prior to analysis by ML techniques. [27] In relaxor-ferroelectric solid solutions, composition-and voltage-concatenated piezoresponse relaxation prior to a Non-Negative Matrix Factorization analysis revealed a strong relaxor-like behavior persisting deep into the ferroelectric phase, providing a fingerprint of the relaxor's electromechanical coupling. [27] Similarly, piezoresponse and contact resonance concatenation allowed a simple Dictionary Learning (DL) algorithm to separate ferroelectric and non-ferroelectric contributions to switching in PbZr 0.2 Ti 0.8 O 3 thin films. [28] However, ML-based analysis of PFM data has been often limited to materials with exceptional electromechanical coupling, usually bulk or thin film single crystals. Given the weaker electromechanical properties of emerging ferroelectric materials, it is of paramount importance to establish analysis methods that can be used ubiquitously, across all materials and SPM-based approaches.
As a case study, we investigate and classify the local band excitation piezoelectric spectroscopy (BEPS) switching response in a defect-rich Pb(Zr 0.53 ,Ti 0.47 )O 3 thin film, fabricated via chemical solution deposition and a novel pulsed thermal processing (PTP) crystallization technique on a glass substrate. [29] PTP enables direct processing of ferroelectric films on substrates with limited thermal tolerance, through plasma arc heating delivered in microsecond pulses, and leads to a wide variety of defects (i.e., nanoscale grain sizes, porosity and secondary phases). [29] Such microstructural variations result in non-piezoelectric phenomena coupling with piezoelectricity and dominating the electromechanical response. The resultant decrease in the signal-to-noise ratio of the PFM data results in challenging to impossible extraction of reliable information about the local ferroelectric behavior. However, by carefully removing outliers, incorporating the On-F measurements via dimensional stacking, and by applying clustering and dimensional reduction techniques to our stacked data, we identify different contributors (e.g., ferroelectric switching, electrostatic coupling) to the local switching piezoresponse. We further correlate the contributors to the defect landscape within the sample (as identified by transmission electron microcopy, TEM, and atomic force microscopy, AFM). The correlative approach reveals a superposition of ferroelectric, electrostatic and other complex non-ferroelectric contributors and their probable nanoscale origins in these defect-rich films. Further, the results highlight that analysis algorithms can be substantially strengthened by appropriate data preprocessing, physical constraints imposed through dimensional stacking, and correlative analysis with other characterization methods. This work provides a blueprint for enhanced and robust identification of piezoresponse and hysteretic behavior in defect-rich ferroelectrics and/or materials with weak piezoelectric contributions (i.e., many organic-inorganic perovskites, fluorite and 2D ferroelectrics).

Structural Characterization and Switching
A dark field TEM image (Figure 1a) illustrates the polycrystalline nature of the film. In most cases, the observable bright and dark grains penetrate across the entire 210 nm thickness. The high-angle annular dark-field (HAADF) image ( Figure 1b) reveals a mixture of columnar and triangular/pyramidal shaped grains, as well as several large (≤200 nm) defects (dark pores) in the film, close to the bottom electrode. Such pores could be a result of energy transfer concentration (and localized temperature increase) in proximity of the Pt bottom electrode during crystallization, and hence fast organic removal. [29][30][31] Looking closer at the surface region (highlighted by the red square, Figure 1b), an atomic resolution image ( Figure 1c) shows clear atomic ordering and crystallinity, along with a smaller (≈10 nm) defect within the grain, (highlighted by a red circle). The selected area diffraction pattern (SADP) in Figure 1d shows strong texturing along the <100> zone axis, although the slight streaking of spots implies some disorder with other orientations also present. The determined d-spacing (Table S1, Supporting Information) from the SADP confirms unit cell parameters to be close to the target composition values and the ratio of outof-plane to in-plane d-spacings implies a small tensile strain inplane to the substrate ( Figure S1, Supporting Information). [32] Energy Dispersive X-Ray Spectroscopy mapping ( Figure S2, Supporting Information) also suggests a composition close to target values, with a slight Pb over-stoichiometry. An electron energy loss spectrum (EELS) analysis ( Figure S3, Supporting Information) of the Ti oxidation peak is similarly consistent with Pb excess and/or O vacancies. In general, the film appears to be well-crystallized and ordered, although local defects and regions of disorder are present, as expected due to the large thermal gradients during crystallization. [29] AFM topography (Figure 1e) confirms the polycrystalline structure, with grain sizes varying from tens to hundreds of nano meters, similar to the TEM observations. The specific location studied via AFM techniques has a significantly rougher surface compared to that studied by TEM (Figure 1b), with 37 versus 10 nm max peak-to-peak variation, respectively. The corresponding PFM amplitude and phase images (Figure 1f,g; Figure S5, Supporting Information) further confirm the heterogeneous nature of the film. All PFM images are out of plane images unless otherwise specified. The as-deposited films exhibit a relatively uniform upward polarization, and localized enhancement in the PFM amplitude. After the BEPS experiment, the resultant domain structure varies by location: some locations show an apparent uniform switched polarization, some locations have no indication of switched polarization, and some locations show ringshaped domains (Figure 1i,j; Figure S5, Supporting Information). Such spatial variation in the material's response to the same input signal suggests substantial variability within the sample. However, an analysis of the average Off-F and On-F amplitude or Off-F PR curves does not allow any further insight, given the nature of the averaging process (Figure 1k-n). Similarly, a one-by-one evaluation of the local measured response provides at best circumstantial evidence of any type of behavior. Indeed, only the Off-F ω shows behavior consistent with polarization switching, indicated by the changes in ω in proximity of coercive bias (Figure 1o). [15,33,34] Therefore, in order to gain insight into the origins of the observed variations, we leverage ML approaches enabling correlation of BEPS response with the post-switching domain formation and/or local topographical and microstructural variations.

Data Preprocessing and Machine Learning Analysis
Two different ML approaches were employed. Due to the defective nature of the sample, a large number of outliers were expected. Hence, the removal of outliers while preserving the critical information was the first challenge for both methodologies. As an initial step, both methods identified poor fits to the SHO model. After this step, the two approaches diverged in terms of outlier and noise detection, removal, and replacement for subsequent analysis. After preprocessing, Method 1 utilized Principal Component Analysis (PCA) to reduce noise, before dimensionally stacking the On-F PR and Off-F PR from the second and third switching cycles. [35] K-means clustering was used to identify the possible number of components and classify distinct behaviors within the PFM response. Method 2 dimensionally stacked the On-F PR, Off-F PR, and Off-F ω from all three switching cycles, after initial data preprocessing. K-means clustering and DL were then used to identify localized outliers, and DL was used to identify three principal behaviors within the PFM response. [28]

Method 1: K-Means on PCA-Cleaned Averaged Cycles
Data Cleaning and Preparation: Poor fits to the SHO model accounting for 2.7% of the overall dataset, were identified by unreasonable Q values and removed (see Supplementary Information 2.1 for further details). For these points, the 2nd and 3rd cycle PR were averaged (for each the On-F and Off-F state) before use, while the 1st cycle PR was left out of the analysis. Omission of the first cycle in PFM data analysis is routinely done, based on the fact that by default the switching and nonswitching domains perform differently during the first full period of signal application. Additional noise was removed from the data via PCA analysis ( Figure S7, Supporting Information), with the intention of preserving only the essential information within the data. [35] To encode physical constraints, the On-F and Off-F PR data were dimensionally stacked to provide correlated information about the ferroelectric and non-ferroelectric behaviors. Behavior Classification by K-means Clustering: K-means clustering groups n observations (data points) into k clusters, predetermined by the user. Here, n = 900, corresponding to the spatial locations (i.e., 30 × 30) over which BEPS was conducted. A scree plot (Figure 2a) for the stacked On-F PR and Off-F PR data with inertia (sum-of-squared-distances within a single cluster) as a function of number of clusters was used to determine an appropriate value of k. The sharpest change occured between 3 and 5: hence, k = 4 was selected. We note that an initial manual inspection of the BEPS data also suggested 4 distinct switching behaviors, while k = 5 based analysis exhibited a fifth cluster incompatible with ferroelectric switching ( Figure S8, Supporting Information). The K-means output for k = 4 is presented, including distribution of the various clusters overlaid onto the post-BEPS PFM phase image (Figure 2c; further details in Figure S9, Supporting Information). The corresponding On-F PR (grey) and Off-F PR (black) behaviors for each cluster are shown in the yellow, green, red, and purple boxes (Figure 2d), and are discussed in detail below.
The k 1 cluster (green) exhibits large On-F PR with small hysteresis, small Off-F PR but with larger hysteresis, and no evidence of polarization switching post BEPS (Figure 2c; Figure S9, Supporting Information). The On-F PR exhibits the smallest coercive bias of all identified contributors, while the Off-F PR curves showed relatively shallow slope near coercive values, implying a wide distribution of "switching" energies (voltages) rather than sharp transitions associated with ferroelectric switching. We describe this cluster as nonswitching, due to the absence of a stable, switched domain in the post BEPS phase image. We note that it is possible that the domains are created but fully back-switch after the removal of the dc bias. The k 2 cluster (yellow) has somewhat similar On-F and Off-F PR in absolute value, and a well-saturated Off-F PR response with a counter-clockwise rotation direction. The K 2 cluster has larger coercive biases for both the On-F and Off-F PR compared to k 1 . After BEPS, this cluster often correlates with solid switched domains: hence, it is considered to correspond to ferroelectric switching. The k 3 cluster (red) has a complex On-F PR hysteresis loop and is found to predominantly coincide with ring-shaped domains after BEPS. This domain structure has been previously attributed to partial backswitching of the polarization after the removal of the dc voltage or anomalous switching. [36][37][38] Thus, k 3 is considered to include a "back-switching" component. The k 4 cluster (purple) exhibits "noses" in the Off-F PR, an irregular shaped On-F PR hysteresis loop, and like k 2 , often corresponds with locations exhibiting solid switched domains. Hence, we refer to this component as "irregular switching." While the above observations are generalized, the classified behaviors do not strictly correlate with a single domain type. For example, points classified as "back-switching" could exhibit a small solid switched domain post BEPS ( Figure S9, Supporting Information). The loose correlation to post switching domain type could point to underlying phenomena mediating and leading to irregular local switching behavior, as would be the case for strong pinning sites (defects) interacting with domain walls. [39,40] In general, assigning data that possess features of different behaviors to a single cluster is a major limitation of a K-means approach. [28] The lack of clear separation between clusters ( Figure 2b) indicates a mixture of behaviors, and indeed supports the possible "information reduction" through clustering algorithms. [35] The identification of convoluted contributors is better addressed by dimensional reduction techniques, which allow for the mixing of separate behaviors within a single probed location. [28]

Method 2: Dictionary Learning on Full Cycle Information
Data Cleaning and Preparation: Although often ignored, the first cycle response is inherently most affected by the underlying nanoscale structural or chemical heterogeneities (compared to subsequent cycles). Furthermore, this cycle is the least affected by eventual accumulation of injected charges due to applied external voltages. Thus, all switching cycles were used in the analysis. As mentioned previously, 2.7% of the overall dataset are outliers caused by poor fitting of the raw signal to the SHO model. The percentage of outliers increased to 5.5% after con-sidering packets outside 3 standard deviations of the mean Q values to be out of range. Here, a linear interpolation using the 2 nearest neighbor valid values was performed to fill the missing data values, enabling preservation of the three separate cycles (for each of On-F PR, Off-F PR, and Off-F ω). To incorporate additional physical information, the DL algorithm was applied to the stacked On-F PR, Off-F PR, and Off-F ω data. Polarization switching is accompanied by changes in the viscoelastic properties of the material and hence changes in the tipsample contact resonance frequency, ω. [15,28,33,34] However, the On-F ω data collected was noisy, as is often the case in BEPS measurements ( Figure S10, Supporting Information). Hence, despite containing valuable information, the On-F ω data was not included in the concatenated dataset, in order to maximize the signal-to-noise ratio. The On-F PR, Off-F PR, and Off-F ω data were normalized prior to data stacking. After stacking, further local outliers were removed from the final stacked data set via K-means clustering and DL methods (Figures S11 and S12, Supporting Information). Dimensional Reduction by Dictionary Learning: In the final DL analysis, the model parameters (i.e., number of behaviors, N, and sparsity factor, α) were selected based on minimization of error, as illustrated in Figure 3a,b. [28] The α parameter controls the tolerance for mixing of separate contributors at each probed location and higher values will linearly increase model error. In this analysis, we use the root-mean square (RMS) error plots to enforce the greatest amount of separation without severely comprising model error. Within the range of considered N and α values, the scree plot and error maps show a significant reduction in RMS error between N = 2 and N = 3, and a substantial increase for α values higher than 20. A DL analysis with four components and α = 20, leads to the detection of local outlier responses ( Figure S13, Supporting Information). Therefore, we opt for a choice of three components (N = 3) and α = 20 (Figure 3c). Each map represents the relative spatial intensity of the single component (N 1 , N 2 , or N 3 ) corresponding to the concatenated On-F PR, Off-F PR, and Off-F ω responses. Points in green represent locations where outliers were removed from the dataset. N 1 and N 2 appeared to be somewhat consistently observed across the probed region. Considering first N 1 , the On-F PR response exhibited the weakest hysteresis across all contributions, with substantially smaller coercive voltages compared to the Off-F PR within the same component. Significant softening in ω was observed near the coercive bias, within the first cycle. However, the softening is dampened in the second cycle, transitioning to strong hardening in the third cycle. Such viscoelastic hardening suggests this behavior is related to non-ferroelectric phenomena (accumulated charge injection/electrostatic coupling). Indeed, examination of an individual point with only N 1 behavior (Figure S14a, Supporting Information) clearly indicates a "double V" On-F A (strong linear voltage-dependence of the amplitude beyond the coercive voltage, with no saturation), typical of electrostatic contributions. [8,41] In N 2 , both the Off-F PR and On-F PR showed saturation signs, although the On-F curve had a substantial rotation. The Off-F ω showed softening in proximity of the coercive voltages, enhanced with increasing number of cycles. Such behavior is consistent with a ferroelectric material, with continued unpinning of domain walls through cycling of the applied waveform. [42] The examination of a point dominant in N 2 behavior revealed saturated butterfly loops for the Off-F A, and coercive voltages smaller in the Off-F than On-F ( Figure S14b, Supporting Information), also consistent with ferroelectric switching. [8,11,16] However, the On-F PR showed also characteristics differing from ferroelectric behavior ( Figure S14b, Supporting Information). Specifically, local maxima superimposed on the ferroelectric-like butterfly loops at a low bias, with increasing intensity at increasing number of cycles, suggest the presence of a nonferroelectric component during the On-F measurements.
For N 3 , the On-F PR curves showed continuous PR increase with increasing number of cycles at both positive and negative biases. The Off-F PR loops showed a clockwise rotation, not compatible with simple ferroelectric response. Additionally, while the second half of each cycle showed an increase in (absolute) PR values, the first half of each cycle remained relatively constant after the first cycle. Lastly, ω increased but with smaller increments over time, until it became almost invariant during the last cycle. Looking more closely at a single point with a predominant N 3 behavior ( Figure S14c, Supporting Information), the On-F amplitude and phase (A and θ, respectively) hint at a switching event at low coercive fields, followed by a rise in A at moderate fields and a decrease in the response at higher fields within the first quarter cycle. During the second quarter cycle, there is another strong rise in the On-F A prior to reaching the apparent coercive field. The shape of the On-F PR response might appear to be similar to those associated with ferroelastic switching events or with field-induced phase transitions. [33,34] However, the progression of the N 3 On-F PR response ( Figure S15, Supporting Information) showed the response approaching saturation in the first quarter cycle and then developing the "nose" in the second quarter, unlike the piezoresponse for ferroelastic switching events or field-induced phase transitions. Furthermore, ferroelastic switching or fieldinduced phase transition are also associated with substantial changes in ω near coercive voltages. [33,34] Such changes were clearly absent from N 3 Off-F ω behaviors.
Further insights into the DL identified behaviors were found by the linear addition of a ferroelectric and a non-ferroelectric component (such as electrostatic coupling), shown in Figure S16   Figure S14a,b, Supporting Information). When these two contributors have opposite overall phases (due to a combination of θ and ϕ, phase offset due to instrumentation), the resulting behavior displays "humps" that seem to rotate and expand/contract the overall loop based on the relative amount of non-ferroelectric to ferroelectric contributions ( Figure S16a, Supporting Information). In fact, linear combinations with mostly non-ferroelectric components resulted in very similar shaped curves to those observed in On-F PR for N 2 ( Figure S14b, Supporting Information), illustrating the unavoidable presence of substantial contributions from non-ferroelectric (non-piezoelectric) phenomena during the On-F measurements. However, if the two contributors have the same overall phase, the resultant behavior is closer to a more typical hysteresis loop, with "noses" at saturation ( Figure S16b, Supporting Information). We note that the fifth cluster in the K = 5 analysis in Method 1 ( Figure S8, Supporting Information) also possesses severe "noses." These observations suggest that the contributors identified by DL are all combinations of non-ferroelectric and ferroelectric behaviors. The On-F and Off-F component of N 1 point to strong non-ferroelectric phenomena, possibly electrostatic contributions. However, neither of the non-ferroelectric or the ferroelectric components used here are necessarily characteristic of all possible responses within this sample: changes in coercive voltages, saturation values, and overall hysteresis shape resulting from differing relative contributions from electrostatic, electrochemical, and ferroelectric behaviors are possible. Thus, N 3 (and k 3 ) behavior is presumed to be the superposition of two or more non-ferroelectric pheno mena that cannot be resolved by either of the ML approaches within the current sample.

Variations in Response
To investigate the origins of the local piezoresponse variability within the sample, we examine the relationship between the microstructure of the film and the contributions identified by both analysis methods (Figure 4). For reference, Figure S9 (Supporting Information) (k 1 , k 2 , k 3 , and k 4 ) and 17 (N 1 , N 2 , N 3 ,  and N 4 ) show each of the identified contributor maps overlaid on top of the post BEPS topography, amplitude, and phase images. The spatial distribution of k 1 (Figure 4a) is mostly along the deepest trenches on the surface, while N 1 (Figure 4e) is most intense in the trenches but also contributes to the observed PFM response in the immediate surrounding locations. The non-switching behavior could be due to increased lateral contact of the tip to the film's surface when located in a trench. Specifically, in the PFM-probed area (Figure 4), we observe a maximum trench depth and width of 15 and 250 nm, respectively ( Figure S18, Supporting Information), possibly suggesting that close to the trench edge, the AFM tip may be in lateral as well as normal contact with the sample, driving in-plane switching. The post-BEPS lateral PFM (Figure S19, Supporting Information) offers little to no evidence for stable domain switching. Regardless, lack of normal contact to the sample surface severely limits the intensity of the applied electric field, leading to no or very weak electromechanical contributions from the sample. [43] At the same time, while the apex of the tip is in a surface depression, the cantilever body is more susceptible to non-local electrostatic coupling with the pore side walls. [44] Hence, we ascribe the non-switching behavior to non-ferroelectric contributions from electrostatic coupling of the cantilever body and the sample surface.
The ferroelectric switching is identified by both K-means and DL algorithms (k 2 and N 2 , respectively), however, it occurs in slightly different regions of the sample. In general, these correspond to locations where clear contact to the surface was made by the SPM tip, but that did not show extreme protuberance from the sample. More locations have N 2 identified contributions (Figure 4f), compared to k 2 (Figure 4b). This apparent discrepancy is directly explained by the differences in the clustering algorithms with respect to dimensional reduction approaches such as DL, where the latter allow for partial overlap of different contributions.
Similarly, points with N 3 contributions (Figure 4g) are assigned to k 3 (Figure 4c) or k 4 (Figures 4d) by the K-means analysis. Indeed, the separation between the K 3 and K 4 clusters is not well defined (Figure 2b). Figure S20 (Supporting Information) compares the measured responses for a few locations assigned by K-means to k 4 locations (red circles Figure 4d), with reconstructed responses based on strong N 2 , and different mixing ratios of N 1 and N 3 behaviors. The mixing ratios were determined from relative intensity of each contributor (N 1 , N 2 , and N 3 ) at each location. Specifically, both the nose-like features in the Off-F PR hysteresis loops and the double-loop-like On-F hysteresis curves of k 4 can be reproduced with a mixture of strong N 2 , and weak N 3 contributions. Most notably, a point mostly attributed to N 2 behavior ( Figure S14b, Supporting Information) is classified as k 4 , despite little to no nose-like character in the Off-F PR. The shape of the On-F PR seems to be used as the dominant classification basis.
In general, the N 3 behavior is most intense at locations associated with high PFM amplitude ( Figure S14c, Supporting Information) and local surface perturbances (Figure 4h). The material at these locations is less laterally constrained compared to other locations and hence, can deform in all directions. Additionally, a higher concentration of defects can be expected in these locations due the increased Pb loss at peaks during processing. [45] Indeed, TEM ( Figure S4, Supporting Information) shows regions of enhanced disorder at the surface of the films, possibly resulting in more electrochemically active regions.
Overall, there is a clear correlation between the film's topography and the dominance of the ferroelectric or the non-ferroelectric behavior. Partial contact in pores emphasizes the contribution from electrostatic-like phenomena, while the probable increase in defect concentration at surface protuberances results in contributions from multiple non-ferroelectric, hysteretic behaviors. In between peaks and valleys, ferroelectriclike switching, albeit with an unavoidable On-F non-ferroelectric contribution, dominates the electromechanical response. This ferroelectric-like contributor is also weakly present in close proximity to peaks and valleys. The ability to resolve and differentiate the subtle ferroelectric-like switching behaviors in these microstructural features illustrates how incorporation of On-F PR, as well as Off-F ω measurements via dimensional stacking, can facilitate the differentiation between piezoelectric (ferroelectric) and non-piezoelectric contributions.
Ultimately, the two different approaches yield slightly different results and only Method 2 allows a true mesoscale correlation between topography and obtained behaviors. This correlation is enabled by 1) the preservation and utilization of the maximum amount of information available in the acquired dataset and 2) the use of a dimensional reduction algorithm, which allows for separation of complex (i.e., multiple contributors to) electromechanical response at each location. We recall that in Method 1, similar to many literature reports on machine learning for PFM analysis, the initial data cleaning removes substantial information via exclusion of first cycle information, averaging of the 2nd and 3rd cycles to fill missing data (poor SHO fits) and finally, PCA removal of global noise from the dataset. In an ideal material, the averaging steps may enhance further analysis by removing spurious features unrelated to the piezoelectric response. In a defective material, with substantial variation in microstructure and underlying heterogeneities, such removal of information is instead deleterious to the analysis, ultimately limiting the correlation of classified behaviors with local topography. Method 1 incorporated both On-F and Off-F PR via dimensional stacking, and utilized K-means clustering. K-means clustering forcefully assigned each location to only one behavior, which led to substantial errors in classification and data interpretation challenges, particularly when correlating the classified behaviors with local microstructural features.
Method 2 adopted a different paradigm to noise and outlier removal. In a defective material, outlier responses could still be due to the material's behavior. In the initial data cleaning, Method 2 used linear interpolation to replace poorly fit datapoints, significantly improving information density by preserving all switching cycles. With respect to dimensional stacking, Method 2 correlated On-F PR, Off-F PR and Off-F ω behavior from all cycles. Additionally, this method leveraged K-means clustering and DL for local outlier identification ( Figures S8 and S9, Supporting Information) only after dimensionally stacking. As a result, the identified outliers corresponded to locations with correlated On-F PR, On-F PR and Off-F ω behavior significantly different from the majority, and hence only locations that fully differed from the majority of the sample (e.g., local surface contamination) were removed. Through linear interpolation and local/correlated outlier removal, the final analysis object retained a higher information density compared to averaging and/or global denoising.
For the separation of contributors, Method 2 utilized Dictionary Learning, with multiple possible contributions at each probed location. The relative contributions from different behaviors were user-defined via the sparsity parameter. Maps of relative intensity for each DL-identified contribution were then correlated with local microstructural features. The inclusion of both On-F PR and Off-F ω combined with the correlation to microstructure allowed two of the major contributors to be attributed to known electromechanical phenomena, i.e., electrostatic response and ferroelectricity. Method 2 proved to be successful in differentiating piezoelectric (ferroelectric) and non-piezoelectric contributions through information density maximization (within the dataset): the approach to poor SHO fits, the inclusion of On-F PR, Off-f PR and ω via dimensional stacking and ultimately the analysis were all geared towards preservation of the highest amount of information. In this way, Method 2 built an information bridge between nanoscale materials functionality and microstructural features of the film without the need for pristine materials or repeated data collection. Hence, Method 2 or similar high information density-preserving methods should substantially increase the investigative power of any voltage-modulated SPM technique, where multiple electromechanical contributions are expected.
Here, Method 2 has provided a baseline analysis of a PZT film processed by CSD and PTP. We specifically recommend applying Method 2 to films with precise defect engineering to enable identification of physical and/or chemical electromechanical contributors associated with specific defects and/or microstructures and hence, processing conditions. Furthermore, the proposed stacked On-F PR, Off-F PR and ω can aide in the analysis of voltage-modulated strain response not only in defect-rich materials, but in any material with weak electromechanical contributions such as the novel fluorite-based ferroelectrics, organicinorganic perovskites, and 2D ferroelectric and piezoelectric materials. Finally, a generalized step by step approach consisting of 1) data preprocessing to preserve information, 2) physical constraints imposed via dimensional stacking, 3) localized and correlated outlier removal and 4) behavior identification by ML algorithm most suited to the materials system is by no means limited to SPM techniques, but rather can be applied to any complex, multiparameter dataset.

Conclusion
We employed a combination of electron microscopy, atomic force microscopy as well as band excitation piezoelectric spectroscopy with ML analysis, to characterize the microstructural and functional properties of a defect-rich PZT sample processed by pulsed-thermal-processing. Appropriate data preprocessing, with high information density and dimensional stacking of on-and off-field piezoresponse and off-field contact resonance frequency prior to ML analysis enabled differentiation between ferroelectric switching and non-ferroelectric hysteretic contributions to the nanoscale electromechanical response. K-means clustering identified four apparent contributions to the observed response, while DL identified three contributions. The fourth component identified by K-means was found to be a mix of other contributions (identified by both K-means and DL methods), demonstrating the importance of 1) considerate outlier removal, 2) maximizing information density and 3) ML approaches that can accommodate multiple contributions within a single probed point. Of the three main contributors to the electromechanical response, two were assigned to ferroelectric and non-ferroelectric (electrostatic) phenomena, respectively, while a third behavior was attributed to the local overlap of two or more non-ferroelectric phenomena. The identified contributors showed significant correlation with the surface topography, strongly supporting the electrostatic nature of the primary non-ferroelectric contributor to the response.

Experimental Section
Materials: The 210 nm-thick PZT film examined in this study was prepared via chemical solution deposition and pulsed thermal processing (PTP) on a platinized glass substrate, as described by Yao et al. [29] PTP enables direct processing of ferroelectric films on substrates with limited thermal tolerance, through Xe plasma arc heating delivered in microsecond pulses. The in-house prepared Pb(Zr 0.53 Ti 0.47 )O 3 precursor solution was 2-methoxyethanol-based with 30% Pb excess. The solution was spun coated onto a Pt/Ti coated glass substrate at 3000 rpm for 30 s. Each layer was pyrolyzed at 400 °C for 1 min on a hot plate for partial organic and solvent removal. After deposition and pyrolysis of three subsequent layers under the same conditions, the sample was crystalized in PTP (Novacentrix, PulseForge 3300) through application of 50 pulses at 450 V for 500 µs each.
TEM: Cross-sections for TEM analysis were prepared using an FEI Scios DualBeam instrument. The sample was initially coated with thin layers of evaporated carbon and magnetron-sputtered Au before focused ion beam preparation with milling at 30, 16, and 8 kV, and final polishing at 5 and 2 kV. (S)TEM analyses were performed on a probecorrected FEI Titan Themis operated at 200 kV. Energy dispersive X-ray spectroscopy (EDX) was performed with an FEI SuperX EDX detector and electron energy loss spectroscopy (EELS) with a Gatan Enfinium EEL spectrometer with a fixed 5 mm entrance aperture. For STEM, a convergence angle of 21.2 mrad was used along with high angle annular dark field (HAADF) inner and outer collection angles of 88.4 and 200 mrad. EEL spectra were obtained with a collection angle of 10.3 mrad.
AFM and PFM: AFM-based experiments were carried out using an Asylum Research AFM (Cypher). Conductive tips (Nanosensors, PPP-EFM) with a nominal resonance frequency of 75 kHz and nominal force constant of 2.8 N m −1 were used for all PFM and BE mode measurements. PFM scans were performed at ≈290 kHz, off the contact resonance peak (≈320 kHz), at an excitation voltage of 2 V using a Zurich Instruments lock-in amplifier (HF2LI). A National Instruments module (PXIe-6124) in combination with a chassis (PXIe-1073) and connector block (BNC-2110), controlled via a LabView interface, was used to apply dc bias waveforms and record the resulting piezoresponse in BEPS mode. A train of dc pulses following a symmetric triangular waveform varying between 0 and ±10 V (positive voltages first) was applied in a grid across the film with a setpoint force of ≈60 nN. Three total voltage cycles were applied at each probed point within each waveform. Each cycle consisted of 64 steps lasting 10 ms, with a dc bias applied during the first half (On-field, On-F) and no voltage applied during the second half of the step (Off-field, Off-F). The response was excited at each step with an ac amplitude voltage of 1 V, and a frequency band, 150 kHz in width, centered around the first harmonic contact resonance (≈280 kHz for the PPP-EFM tip used). After fitting to a Simple Harmonic Oscillator (SHO) model, the output data contained A, θ, ω, and Q from each of the On-F and Off-F sets. Additionally, piezoresponse (PR) was also calculated, given by Equation (1), where ϕ represents a phase offset mostly due to instrumentation.

Supporting Information
Supporting Information is available from the Wiley Online Library or from the author.