Tensor factorization for elucidating mechanisms of piezoresponse relaxation via dynamic Piezoresponse Force Spectroscopy

Spatially resolved time and voltage-dependent polarization dynamics in PbTiO3 thin films is explored using dynamic piezoresponse force microscopy (D-PFM) in conjunction with interferometric displacement sensing. This approach gives rise to 4D data sets containing information on bias-dependent relaxation dynamics at each spatial location without long-range electrostatic artifacts. To interpret these data sets in the absence of defined physical models, we employ a non-negative tensor factorization method which clearly presents the data as a product of simple behaviors allowing for direct physics interpretation. Correspondingly, we perform phase-field modeling finding the existence of ‘hard’ and ‘soft’ domain wall edges. This approach can be extended to other multidimensional spectroscopies for which even exploratory data analysis leads to unsatisfactory results due to many components in the decomposition.

The advent of Piezoresponse Force Microscopy (PFM) [17][18][19][20][21] in the mid-90s enabled visualization of ferroelectric domains in opensurface films and ferroelectric capacitors, often with sub-10 nm resolution [22][23][24][25][26][27] . These imaging studies have been further extended to map domain evolution under applied fields [28][29][30][31] and as a function of time 32,33 and temperature 34 , providing information on domain wall motion and nucleation mechanisms and frozen disorder 26,[35][36][37][38][39][40][41][42] . In parallel, several spectroscopic imaging techniques in PFM have been developed 43,44 . In Piezoresponse Force Spectroscopy (PFS), a PFM signal is measured as a function of the magnitude of a dc voltage pulse applied to the probe, providing information on the bias dependence of the PFM signal 45 . The latter, in turn, is related to the size of the domain formed below the probe and thus provides information on thermodynamics and kinetics of domain switching [46][47][48] . Nonetheless, several key challenges remain in the interpretation of microscopic mechanisms responsible for ferroelectric behavior.
The proliferation of multidimensional PFM methods brought forth the challenge of efficient and high-veracity data analysis, i.e., reduction of a 3, 4, or 5D data set into 2D maps representing material properties in the most efficient manner. Such analysis is relatively straightforward if the functional model representing materials response to external stimuli is known, allowing for direct least-square functional fits to extract material parameters. Yet this is rarely the case for locally measured responses, necessitating development of strategies based on information theory methods. This complexity and difficulty in interpretation stems from the fact that existing analysis methods address the dimensions of the measurement parameter space individually, as opposed to treating them as linked. Most of the existing methods such as principal component analysis 49,50 , non-negative matrix factorization 51,52 , and independent component analysis 53 are limited to factorizing matrices, whereas for multidimensional data with multiple spectroscopic axes (for instance, time and voltage), the intrinsic form of the data is a tensor. Subsequently, information can be lost during the flattening of the tensor into a matrix that can be factorized using these tools.
A related challenge in multidimensional PFM is signal interpretation 54,55 . It is well known the measured PFM signal is often beset by non-piezoelectric contributions, e.g., long-range electrostatics, as well as artifacts that arise from e.g., contact resonance changes 56,57 . While a number of measurement techniques have been developed to quantify the electrostatic contribution, and mechanisms such as full-band excitation or dual AC resonance tracking can circumvent the resonance shift issue, achieving quantitative measurements in PFM is still a difficult undertaking [58][59][60] .
Here, we use recently developed dynamic piezoresponse force microscopy (D-PFM) 61 , which provides information on the time and voltage dependence of the piezoresponse in each spatial location, to understand the polarization dynamics in a ferroelectric system. As a model material system, we choose tetragonal PbTiO 3 (PTO) as it exhibits a dense a-c domain structure ideal for studying ferroelectric relaxation dynamics (Fig. 1a-c). To isolate the true electromechanical dynamics of the system, D-PFM measurements were employed via interferometric displacement sensing (IDS) 58,62 in a coupled atomic force microscope. Unlike traditional PFM which is an angular measurement of the cantilever motion, interferometric displacement sensing uses a laser focused directly above the tip, which is routed back to an interferometer enabling accurate measurements of tip displacement, as shown in Fig. 2a. Note, IDS eliminates the long-range electrostatics that originate from the interaction between the cantilever shank and sample, which is a key source of spurious background signal and artifacts, and quantifies the tip displacement. D-PFM in conjunction with interferometric displacement sensing provides a new paradigm for measuring local ferroelectric properties and their dynamics under a field.

RESULTS AND DISCUSSION
Understanding polarization dynamics in a ferroelectric system The representative surface topography is illustrated in Fig. 1a, showing a surface roughness of 4.1 nm (RMS) and topographical features as large as 20 nm in height. The corresponding bandexcitation PFM amplitude and phase images are shown in Fig. 1b, c, respectively, exhibiting a well-defined two-level domain structure. For the (001) tetragonal ferroelectric, the domains can be readily classified as a-(in-plane) and c-(out of plane) domains, identified as the low and high amplitude regions in Fig. 1b, respectively.
To understand polarization dynamics in this system, D-PFM measurements with IDS are performed across a spatial grid.
Specifically, the PFM response (see Fig. 2b, blue) at a single position is measured as a function of time after a voltage pulse following a triangular envelope (see Fig. 2b, red), and the measurement is repeated across a grid of points, building up a spatial map of the relaxation dynamics. The size of the generated data set is a 60 × 60 spatial grid with 128 relaxation time steps measured at 6 ms intervals, and 16 distinct voltage pulses. The response averaged over the whole image is illustrated in Fig. 2b, showing clear relaxation processes. Note that the direction and magnitude of relaxation strongly depends on the set voltage amplitude. Importantly, measured with IDS, it is evident that relaxations exist, and given that switching occurs, that domain wall motion is also activated under the applied fields used. We note that long-range electrostatics can also be a dominant cause of relaxation in such spectroscopies measured using traditional optical beam deflection, so the measurement of local relaxations is an improvement on the existing technique 60 .
To understand the polarization dynamics of such a dense domain structure, we choose to perform D-PFM on a smaller area to provide better resolution of a/c domain features (Fig. 3a). Given this is a large multidimensional dataset, for a first visualization we performed k-means clustering on the entire multidimensional dataset (including time) to segment the image into k = 5 clusters; k-means is an unsupervised machine learning algorithm that automatically assigns each pixel to belong to a member of the kth class (or cluster) such that the within-cluster sum of squares is minimized. The k-means map is shown in Fig. 3d, and the associated cluster means (the respective cluster's average piezoresponse versus time) are plotted in Fig. 3b, c, e, f, h. It is evident there are changes in the hysteresis loop depending on the location, i.e., a-domain or c-domain, and it is noteworthy that the relaxation process varies across the a-c domain as shown by the line profile of the loops along the circled region in Fig. 3d, and plotted in Fig. 3g. Here, the relaxation appears to be asymmetric  with respect to the a/c domain boundary. However, it is difficult to extract meaningful quantitative information from such a representation due to high noise and variability. We note also that Fig. 3 shows two distinct types of c-domains: one type embedded within the c/a domain structure, and another (on the right side of image 3(a)) which appears separated and much larger. The reason for low response in the c-domain embedded in the c/a structure 4 is likely due to the geometry: with a large volume fraction of a-domains in the probed volume of the tip, the c-domain fraction is likely much smaller, leading to suppressed response.
To gain quantitative insight into the spatial and voltage variability of the relaxation phenomena, we employed nonnegative tensor factorization (NTF), which is commonly used in a variety of statistics to reduce n-dimensional tensorial data to a more compact representation whilst still reproducing (to some tolerance) the original values 63 . Here, we restricted the study the piezoresponse amplitude, which is always positive, and thus satisfies the non-negativity constraint of NTF. The details of the algorithm are available in ref. 63 , but briefly, we note that we may write the piezoresponse amplitude as a tensor X of order 3 (the respective orders being space, voltage and time), X 2 R M1 M2 M3 . In our case, M1 = 3600 for the vectorized spatial positions, M2 = 128 for the time dimension, and M3 = 16 for the number of voltage steps. Then, NTF aims to determine a decomposition of X such that: where a k 2 R M1 ; b k 2 R M2 ; c k 2 R M3 ; K is the rank, and represents an outer product of vectors. An alternative way to write this problem is that we attempt to find the matrices A,B,C and then find the matrices subject to the constraint that all elements of these matrices remains non-negative, i.e.,: where A; B; C ¼ P K k¼1 a k b k c k ;and F is the Frobenius norm. In general, we wish to choose a low number for K, but still large enough that good approximation for X can still be found. Note, this is the same for traditional multivariate analysis, i.e., nonnegative matrix factorization (NMF) 51,52 . However, unlike in traditional NMF, the NTF algorithm used here results in 3 vectors for each specific k value, i.e., one for the spatial position (interpreted as loadings), one for the voltage vector and one for the time vector. Lastly, it should be noted that unlike NMF, which does not have uniqueness for the approximation, optimum solutions for NTF can be unique under some soft conditions 64 . Figure 4 shows the NTF analysis of the D-PFM data corresponding to the full voltage-time sweep. The larger maps are A matrices, which have been reshaped into 2D matrices and can be considered as loading or abundance maps commonly seen in NMF or principal component analysis (PCA). The smaller insets consist of matrices B and C and are recognized as the voltage and time dependences. It is important to note that the NTF analysis with rank 4 results in a reconstruction error of~10%. Remarkably, the complex time and voltage-dependent dynamics can be represented as a sum of 4 timedependent components with voltage-dependent coefficients. This is especially important because we do not know, a priori, how the voltage pulse will affect the system's relaxations.
Examination of the components reveals component 1 (Fig. 4a) is more prevalent at negative voltages while component 2 (Fig. 4b)  Fig. 3d; the loop labels are colored corresponding to the cluster membership, d k-means clustering on entire multidimensional dataset, g piezoresponse amplitude vs time for line profile which spans a/c-domain boundary to edge of adjacent a-domain (left to right, respectively) in Fig. 3d. Time color bar corresponds to time dependence of piezoresponse displayed in Fig. 3b, c, e, f, h.
is favored at positive voltages with contrasting relaxation behavior, both of which spatially correspond to c-domains as seen in the abundance maps. Figure 4c shows component 3, which is spatially correlated to a-c domain boundaries with a preference towards one a-c domain orientation, in corroboration with clustering seen in Fig. 3d. The strong asymmetry can be attributed to the tilt of the ferroelastic a-c domain walls, with the volume fraction of the c-domain in the probed volume changing asymmetrically as the tip is moved away/towards an a-c domain boundary 65 . Indeed, we observed such an asymmetry in a previous report, where the nonlinearity in the measured piezoresponse as a function of applied AC field was modeled in a similar system, and the asymmetric response across the ferroelastic interface was reproduced 66 . Here, we note that this asymmetry is reproduced from the data with a purely unsupervised analysis technique. The expanded view of the bottom-left corner of the abundance map in Fig. 4c is shown in Fig. 5a, and the asymmetric contrast is clearly visible across the c/a interface in the average   5 Asymmetric relaxation response mechanism. Asymmetric relaxation response from c/a domain boundaries may be ascribed to the existence of a "hard" and "soft" c/a/air junction. a Component 2 from NTF, abundance map, lower left region only. b Average profile along the x-direction of (a), clearly showing the oscillation in the abundance map due to the different relaxation behavior. c Schematic of the "hard" and "soft" junctions. The "soft" junction results in more wall motion during field application, and consequently more relaxation.
profile, plotted in Fig. 5b. Component 4 (Fig. 4d) seems to be isolated to transient effects associated with domain nucleation events. Additionally, all four components exhibit monotonic relaxation behavior and obey dual exponential relaxations (Fig.  4, see fit in insets). Finally, component 2 appears to show some non-physical effect which is not seen in any individual ("pure") pixel; indeed, performing NTF on the c-domains and a-domains separately does not result in such a component. Care must be taken to not ascribe physical mechanisms when the unmixing is not ideal as is the case for this component.
The fits for the dual exponential behavior of all four components are displayed in Table 1 showing distinguishable relaxation times for each component (single exponential fits in Table S1). Component 1, 2, and 4 have orders of magnitude difference in the two relaxation times, where the much larger relaxation time (on the order of minutes) might be attributed to slow polarization dynamics, whether it be electrochemical or purely ferroelectric in origin. Indeed, other studies on electrochemically active materials show the possibility of electrochemical dipoles causing such relaxations in local AFM measurements, and we cannot discount this here 67 .
On the other hand, somewhat faster dynamics can be expected for lateral domain wall motion, although this will still be a creep process 68 . Indeed, if one expects a creep velocity of~100 nm/s, then the domain wall displacement in 5 ms is 5 nm, which should be sufficient to cause detectable changes in the response of a strong ferroelectric for a probing volume dictated by the tip radius (r~10 nm). Thus, we suggest that the fast process is dictated by lateral domain wall motion, whilst the slower process is likely mediated by electrochemical processes.
A possible mechanism for the observed spatial variability in relaxation behavior of component 3 (also Fig. 3g) is the existence of "soft" and "hard" edges at the ferroelastic interface (as described in our previous report 66 ), i.e., when the tip is on the right side of the c/a boundary (Fig. 5c), the domain wall is more labile, and therefore less relaxation is observed when the bias is removed. On the other hand, when the tip is on the left side of the c/a boundary, only the "hard" edge of the ferroelastic wall is activated, resulting in a large relaxation due to the higher energy penalty. To explore the feasibility of the proposed mechanism, we performed phase-field modeling on a PbTiO 3 thin film coherently grown on a KTiO 3 substrate, with a 1 (100) and c(001) tetragonal domains and c/a domain walls in the equilibrium state (Fig. 5c). A constant local tip bias of 2.0 V was applied on the film surface near the "soft" and "hard" edges of the 90 • domain walls with the infield domain structures illustrated in Fig. 5d, e. It is seen that under 2.0 V the switched c-domain (00-1) under the tip near the "soft" edge ( Fig. 5d) grows larger and deeper in size than the one near the hard edge (Fig. 5e), indicating a clear preference for switching on the "soft" side. Furthermore, when the tip biases are removed, the domain structure begins to relax, as seen in Fig. 5f-i. Clearly the relaxation is larger for the "hard" edge (the switched domain shrinks to a much smaller size), and the original domain structure is mostly restored, as compared to the "soft" edge switching. We believe this because the switched domain is less stable on the "hard" side due to the large energy penalty, such that the relaxation is faster. Our phase-field simulation thus clearly shows the difference of domain switching and relaxation near the "soft" and "hard" edges of the 90 • domain wall. The observations of this asymmetry, in combination with the phase-field modeling, suggest that ferroelectric domain wall motion is a strong component of the observed electromechanical relaxation behavior; however as stated earlier, we cannot discount that other mechanisms may also play a role.
Finally, to conclude our investigation of the NTF analysis on D-PFM relaxation data and the associated mechanisms, shown in Fig.  6 is a direct NMF comparison. Figure 6a-d shows the relaxation components plotted alongside the loading maps for −15V pulse (Fig. 6e-h) and +15 V pulse (Fig. 6i-l) for the 4 components. At negative biases, little difference is observed in the abundance maps indicating the relaxation behavior is a convolution of all relaxation components. In contrast, large positive biases show a preference towards component 1 and component 4, as seen in the abundance maps. It is noteworthy the positive bias/ component 4 highlights the interesting relaxation dynamics seen at the a-c domain boundaries. However, it lacks correlation to the larger a-domain/c-domain, i.e., the spatial, temporal, and voltage relationship for this mechanism is not intuitively observed. Nevertheless, NMF produces 64 abundance maps (4 for each voltage) and 4 relaxation components resulting in a hard to interpret representation of the multidimensional dataset and further emphasizing the need for tensor rather than matrix factorization. Furthermore, the technique comparison clearly illustrates NTF's ability to maintain interconnected relationships between multiple variables.
To summarize, spatially resolved time and voltage-dependent polarization dynamics in multidomain PbTiO 3 films are explored using D-PFM in conjunction with interferometric displacement sensing. To extract physically relevant and simple behaviors from the 4D D-PFM data sets, we introduce the non-negative tensor factorization method. We observe distinct relaxation mechanisms and ferroelectric behavior associated with relaxation of domain walls and dipole realignments with an enhanced relaxation behavior at one side of the a-c domain wall. We then perform phase-field modeling to capture this asymmetric relaxation process finding the underlying mechanism to be the presence of "hard" and "soft" domain edges. Finally, to conclude the importance of NTF, we directly compare results derived from the traditional multivariate analysis technique, NFM, finding an underwhelming representation of the multidimensional dataset. This manuscript highlights the advancements in both PFM data acquisition and multidimensional analysis, and the applicability of NTF to broader scientific areas.
METHODS Experimental As a model system, we have chosen a 700 nm-thick PbTiO 3 thin film grown by chemical vapor deposition on (001) KTaO 3 substrates with an SrRuO 3 conducting buffer layer; this film displays a dense polydomain structure as previously reported 69,70 . The PFM was performed using a commercially available atomic force microscopy (Cypher, Oxford Instruments) fitted with an interferometric displacement sensor.

Simulation (phase-field)
In the phase-field simulation we use polarization vector Pi = (P x , P y , P z ) as the order parameter to describe the ferroelectric state in the PbTiO 3 thin film. The temporal evolution of Pi is calculated by minimizing the total free energy F with respect to Pi by numerically solving the time-dependent Landau-Ginzburg-Devonshire (LGD) equations 71 , where x is the spatial position, t is the time, L is the kinetic coefficient related to the domain wall mobility. The total free energy F of the PbTiO 3 thin film includes the Landau, gradient, elastic, and electrostatic energies, i.e., F ¼ V f lan P i ð Þ þ f grad P i;j À Á þ f elas P i ; ε ij À Á þ f elec P i ; where V is the total volume of the system, ε ij and E i denote the components of strain and electric fields, respectively. Detailed expressions of each free energy density can be found in ref. 72 . Equation (3) is numerically solved using a semi-implicit spectral method 73 based on a 3D geometry sampled on a 128Δx × 128Δy × 32Δz system size, with Δx = Δy = Δz = 1.0 nm. The thickness of the film, substrate, and air are 20Δz, 10Δz, and 2Δz, respectively. The temperature is T = 25°C, and an isotropic relative dielectric constant (κ ii ) is chosen to be 50. The gradient energy coefficients are set to be G 11 /G 110 = 0.6, while G 110 = 1.73 × 10 −10 C −2 m 4 N. The Landau coefficients, electrostrictive coefficients, and elasticcompliance constants are collected from ref. 72 .
To model a local contact of a scanning probe tip, we specify the electrical potential distribution on the bottom (ϕ 1 ) and top surface (ϕ 2 ) of the PbTiO 3 film as, ϕ 2 x; y ð Þ ¼ ϕ 0 γ 2 x À x 0 ð Þ 2 þ y À y 0 ð Þ 2 þγ 2 (6) in which ϕ 0 is the tip bias, (x 0 , y 0 ) is the tip position and γ is the half-width half magnitude of the tip. Here we chose γ = 5 nm.

DATA AVAILABILITY
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.