Local and correlated studies of humidity-mediated ferroelectric thin film surface charge dynamics

Electrochemical phenomena in ferroelectrics are of particular interest for catalysis and sensing applications, with recent studies highlighting the combined role of the ferroelectric polarisation, applied surface voltage and overall switching history. Here, we present a systematic Kelvin probe microscopy study of the effect of relative humidity and polarisation switching history on the surface charge dissipation in ferroelectric Pb(Zr$_{0.2}$Ti$_{0.8}$)O$_3$ thin films. We analyze the interaction of surface charges with ferroelectric domains through the framework of physically constrained unsupervised machine learning matrix factorization, Dictionary Learning, and reveal a complex interplay of voltage-mediated physical processes underlying the observed signal decays. Additional insight into the observed behaviours is given by a Fitzhugh-Nagumo reaction-diffusion model, highlighting the lateral spread and charge passivation process contributors within the Dictionary Learning analysis.


Introduction
Surface and bulk electrochemical phenomena are particularly important in ferroelectrics, whose switchable remanent polarisation modulates their electrochemical reactivity, determining fundamental material responses and giving rise to promising catalysis and sensing applications 1,2 . Electrochemistry has been linked to all aspects of ferroelectric switching and screening [3][4][5] , with recent studies highlighting its influence on catalytic effects and functionality-modifying phenomena 6,7 . These effects can be amplified by the additional presence of adsorbates, such as surface water ubiquitous under ambient conditions 8, 9 . As ferroelectric surfaces exhibit varying polarisation orientation and chemical composition, the interplay of electrochemistry with surface water generates a wide range of complex phenomena. Structurally, preferential adsorption can modulate the usual humidity dependent molecular arrangement observed on simple dielectric materials 10 of 1-2 monolayers of ice-like water, followed by less ordered liquid layers. Chemically, varying dissociation rates into OH − and protons vs. molecular water and even modulation of the surface and bulk material composition as a function of polarisation orientation and relative humidity have been reported 11,12 .
The observed effects are further enhanced under the application of an electric field, with radical water-mediated changes and time-dependent processes taking place, as also seen in non-polar materials 13 . In addition, sufficiently high fields will enhance oxidation reactions at the surface and promote charge injection and/or polarisation switching [14][15][16] , leading to a strong dependence of electrochemical response on the electrostatic (switching and charging) history of the material. The combination of all these phenomena gives rise to a highly nontrivial mixed response of ferroelectrics to changes in relative humidity or surface bias, with widely ranging timescales and magnitudes, as demonstrated previously by local probe surface potential imaging 12,17 .
Over the last decade, investigation of complex phenomena with multiple parameters have been made possible by incremental advances in experimental techniques and acquisition hardware, yielding rich datasets of correlated responses at a broad range of timescales. Additionally, tremendous progress has been made in tools designed to analyze the resulting hyperdimensional datasets, based on Big Data approaches such as machine learning and dimensional reduction. These have been successfully applied to studies of ferroelectric materials, enabling a better understanding of the competing and correlated processes involved in switching and surface screening under ambient conditions [18][19][20] .
Here, we investigate the mechanisms of surface charge dissipation on ferroelectric thin film surfaces as a function of time, applied voltage, writing history from as-grown through final polarisation, and relative humidity. Local probe imaging of the surface potential was performed on two Pb(Zr 0.2 Ti 0.8 )O 3 (PZT) thin films of opposite polarisation, after the writing of a predefined structure. The resulting time-dependent datasets were dimensionally-stacked and analyzed by means of Dictionary Learning in order to unravel the underlying physical and chemical behaviours. A simplified reaction-diffusion model complements the experimental work, allowing the results of the machine learning analysis to be interpreted in terms of physical/chemical processes. Our observations unequivocally demonstrate that whilst both positive and negative surface charges decay with time, only negative charges spread laterally on the surface, diffusing far beyond the originally written areas.

Imaging surface charge dissipation
Surface charge dissipation was imaged by Kelvin probe force microscopy (KPFM) on 100 nm thick PZT films (presenting either up-or down-oriented as-grown monodomain polarisation) and tracked as a function of relative humidity, polarisation, applied voltages, and time. Domain structures were patterned by scanning the biased atomic force microscope tip in contact with the surface: first, to define three lateral stripe regions with positive (+8 V), zero, and negative (-8 V) tip bias; then, three 2/18 overlaid vertical stripe regions, again with positive, zero, and negative tip bias, as illustrated in Fig. 1 (a). The sequential writing is performed in order to take into account possible effects of polarisation switching history and varying intensity of charge injection (areas scanned twice with negative or positive voltages). Once the writing was completed, KPFM imaging was carried out until dissipation of the surface potential contrast, approximately 9 to 12 hours. An illustration of resulting KPFM image is shown in Fig. 1 (b), with the writing structure overlaid. We note that whilst the piezoresponse force microscopy image of the ferroelectric domains generated through the writing process follows the writing structure exactly, the KPFM usually extends outwards due to charge diffusion, especially at high humidity. 1 The procedure was repeated at each relative humidity setpoint (ten different values spanning 6% to 80%), for both samples (up-and down-oriented as-grown polarisation). The control of relative humidity was performed with a in-house-built flow-based, low noise humidity controller, which enables fast and precise in-situ control in the atomic force microscope chamber 21,22 . The first five resulting surface potential images (approximately 70 minutes after writing) are shown in Fig. 2 for both samples at all ten humidity levels probed. Several trends can be identified from simple visual inspection. First, in both samples, faster surface charge dissipation and larger spatial spread of the KPFM signal beyond the boundaries of the stripe domains can be observed at higher humidities. Second, the positive charges seem to be much more localized (bright yellow contrast) than the negative charges (dark blue contrast), regardless of the initial sample polarisation. Finally, the charge contrast seems to be more intense in areas where the polarisation was switched from its as-grown state -i.e. where negative tip bias was applied for the down-polarised film and positive tip bias for the up-polarised film -hinting at the role of switching history in the dissipation. A new structure is written at each relative humidity setpoint and then monitored continuously by KPFM until complete surface potential dissipation. At low humidities the charges seem to be more localized and retained longer than at high humidities. Additionally, at higher humidities the horizontal lines corresponding to the first writing step are not visible, suggesting much faster charge dynamics in presence of high water content at the surface. Lastly, the written pattern remains visible in the up-polarised films at much higher humidity levels than in the down-polarised films.

Dictionary Learning analysis
Due to the complex and correlated nature of the full dataset, any analysis beyond simple visual inspection is challenging as it will introduce human error and bias: regions over which observed behaviour is to be averaged may be selected arbitrarily, outliers and/or unexpected data might be ignored, and overall the spatial distribution of information within the images will be lost as most analysis would concentrate on average signals. To address these issues, machine learning-based techniques were leveraged for the correlative analysis of the acquired data.
The KPFM images for all ten humidities and two samples (one up-polarised and one down-polarised) were stacked along the spatial dimension in order to induce a time correlation, and are shown in the Supplementary Materials, Figure S1. Such stacking imposes a constraint of identical behaviour(s) across polarisation, time and humidity values 20 , thus enabling a direct comparison of the evolution of the behaviours as a function of initial sample polarisation and relative humidity. Dictionary Learning was selected as the machine learning dimensional reduction technique for the further analysis of the dataset due to its intrinsically sparse nature 23 , with the expectation of only a few behaviours being spatially co-located. The number of components n is a free parameter and was evaluated from two to ten, allowing us to track which behaviours are associated with the different regions as they successively segregate out, as detailed in the Supplementary Materials, Figure S3.
The resulting decomposition, here presented for n = 10, shows three distinct types of components, each associated with and increasing n simply refines the areas written with negative bias into increasingly fine concentric segments.
To quantify the time evolution of the positive and negative decays, the corresponding components -C2 to C5 -were fitted to a double exponential decay function, CPD(t) = A 0 + A 1 expt/τ 1 + A 2 expt/τ 2 , with CPD the contact potential difference, and A 0 the offset corresponding to the final CPD value, A 1 , A 2 the amplitudes, and τ 1 , τ 2 the time constant of the decays. From the parameters shown in Fig. 3d, two behaviors seem to be present in all decaying components, each with two substantially different time constants. The latter differ by at least an order of magnitude, as reported previously for this material, indicating the coexistence of a material-dependent fast electrochemical process and a slow process mediating the surface charge dissipation 12 .
The fast process, of the order of 1000 seconds, can be related to the re-equilibration of surface screening charges 24 , whereas the slow process, of the order of 10000 seconds, is specifically attributed to the dynamics of tip-injected surface or near-surface charges.
In all cases, however, the overall mechanisms (background, positive and negative decay) are present in both films and share the same parameters, varying in intensity with humidity and initial polarisation, but mostly depending on the polarity of the applied voltage. The positive voltage related component, C2, for instance is found in both films below 50% RH, but seems to be more present at the higher humidity only in the up-polarised films, where a positive voltage application switches the polarisation.
Conversely, negative voltage related components, C3, C4, and C5 seem to be more present within the down-polarised film, where such a voltage polarity switches the polarisation. The coexistence of the same physical dissipation processes (as fitted above) in both samples suggests that the initial polarisation orientation does not affect the dissipation mechanisms which seem to be primarily dominated by voltage and surface electrochemistry. However, the history of polarisation is important, as switching processes during voltage application will result in higher surface charge density and/or longer charge retention, especially at higher humidities.

Reaction-diffusion simulation
To better understand the specific behaviors observed in the experimental data and the resulting Machine Learning analysis, a simplified reaction-diffusion model was constructed, based on the FitzHugh-Nagumo system. This is an excitable system model, ideal for describing the distributions of different species during their relaxation to equilibrium, as for example in neural synapses 25 . As the experimental data suggests the presence of lateral spread acting together with surface and bulk electrochemistry, the parameters of the model were tuned to reproduce the observations qualitatively, with a spatially localized and negligibly diffusive positive charge distribution, and an initially more spread and highly diffusive negative charge distribution.
The resulting time evolution of the combined distributions is shown in Fig. 4a.   The cumulative evidence points to a non-trivial behaviour of surface charges, with the positive and negative charges exhibiting a substantially different diffusion/dissipation mechanism that is modulated by the underlying ferroelectric polarisation. The coexistence of these electrochemical processes opens a pathway towards polarisation-directed catalysis, with a highly localized control of charges.

Materials growth
Films of 100 nm thick Pb(Zr 0.2 Ti 0.8 )O 3 (PZT) have been grown on Nb:SrTiO 3 single crystalline substrates by RF-magnetron sputtering, as described elsewhere 26 . The films used in this study were up-and down-polarised as grown, respectively, with a smooth sub-nanometer roughness topography.

Kelvin Probe Force Microscopy
The scanning probe microscopy imaging was performed on an Asylum Research Cypher microscope equipped with an environmental control chamber, using Bruker SCM-PtSi cantilevers with a resonance frequency of ∼ 75kHz and a spring constant of ∼ 2.8N/m. Sample temperature was kept constant at 30 • C, and relative humidity was varied through the use of a home-build controller 21,22 . The humidity was stabilized for at least an hour prior to each measurement series. Topographic imaging was performed in intermittent contact mode (AC Mode), with a free oscillation amplitude of ∼ 10nm and a setpoint of ∼ 5nm. The contact potential difference between the tip and the sample was acquired over 12x12µm areas at a scan speed of 0.5Hz in single pass Kelvin probe force microscopy mode with the following parameters: 3V excitation voltage, 5kHz frequency, −90 • phase, 5 proportional gain, 10 integral gain. Domain writing was done in contact mode with a 0.5V setpoint over 8x8µm areas at a scan speed of 0.5Hz and applied voltages of ±8V , as detailed in Fig. 1(a).

Dictionary Learning analysis
Time-dependent data was stacked into data cubes, and the data cubes at each humidity and for each sample were concatenated along the spatial dimension x, yielding a correlated dataset. A similar stacking method was performed for the reaction-diffusion modeling dataset. The Dictionary Learning implementation of the scikit-learn library was then used on the flattened data set with the following parameters: 3 − 10 n_components, 1000000 max_iter, 8 n_jobs, lasso_lars transform_algorithm, 15 transform_n_nonzero_coefs, lars fit_algorithm, and 0 random_state. The code and data used for the analysis is freely available on the paper repository [CITE yareta_or_zenodo].

Reaction-diffusion modeling
The reaction-diffusion model used is based on a 2D Fitzhugh-Nagumo excitable system model 25 . The implementation used is adapted from existing Python code, recipe 12.4 in 27 . The specific form of the model is a set of two differential equations: The complete humidity-sample-time dataset is shown in Fig. S1, for the two samples (up and down polarised) and ten humidities.

Spectral decomposition principle
The principle behind spectral decomposition methods such as DictionaryLearning is shown in Fig. S2. The complete dataset is split into a set of one-dimensional components with corresponding two-dimensional weight maps. Dimensional stacking was used to integrate the relative humidity and initial sample polarisation dependence into the maps, thus correlating the data sets and enabling quantitative comparison of time-dependent decays.

Dependence on the number of components
A systematic investigation of the influence of the number of components on the DictionaryLearning analysis is shown in Fig.   S3. As the number of components increases, so does the number of physical behaviours. At first, when n = 2, only a single physical behaviour is visible, alongside a Gaussian noise component with a null weight map. When the number of allowed

12/18
components is increased up to n = 6, the negative voltage process splits from the main component. Increasing the number of components to n = 7, the concentric structure for the negative voltage region starts to form with the split of the latter into two distinct map/component pairs. At n = 8, a background component splits from the positive voltage component, and finally at n = 10, a third concentric negative component makes its appearance.

Reaction-diffusion modeling
The reaction-diffusion modeling, here based on a Fitzhugh-Nagumo excitation-relaxation process, provides insight into the contribution of the various physical processes to the machine learning based analysis. Such a sandbox approach enables a unique way to tailor the interactions present within the two species. In particular, the reaction-diffusion modeling highlights the importance of sparse approaches such as DictionaryLearning when dealing within physical materials presenting very different co-located behaviors.

Initial conditions choice
Initial species distributions were chosen to mimic the positive and negative charges as deposited by biased-tip atomic force microscopy writing. The positive charge species u shown in Fig. S4(a) was initialized with only a reaction interaction with the negative charge species v, and no diffusion component. The negative charge species v, shown in Fig. S4(b) is susceptible to both reaction and diffusion processes. The total distribution u − v is shown in Fig. S4(c), and qualitatively resembles the Kelvin probe force microscopy data in Fig. 2.

Principal component analysis
Whilst principal component analysis (PCA) does not provide a physically interpretable spectral decomposition, it can be used to estimate the complexity of the decomposition. The PCA weight maps for the first 50 components are shown in Fig. S5(a), illustrating the concentric-like pattern generated by the diffusive nature of the negative charge species distribution v. Indeed, the positive charge species distribution u is only present in the first four maps of the decomposition. The corresponding PCA components shown in Fig. S5(b) additionally shows a periodic behaviour with increasing frequency for the first 26 components, as evidenced by their Fourier transform in Fig. S5(c).

Non-negative matrix factorization
A non-negative matrix factorization (NMF) approach was tried concurrent to DictionaryLearning, with results for n = 4 and n = 10 components shown in Fig. S6(b,c) respectively. Although not able to separate the physical processes as cleanly as DictionaryLearning, NMF also highlights the concentric decomposition of the negative voltage species with an increasing number of components.