Multivariate temperature-series analysis of stress-induced ferroelectricity in SrTiO₃: a machine learning approach with K-shape clustering and hierarchical Bayesian estimation

ABSTRACT

A new machine learning approach that transforms time-series analysis into temperature-series analysis is introduced to analyze stress-induced ferroelectricity in SrTiO₃ at 231 MPa using birefringence images observed at successive temperatures. The spatial distribution of the temperature-series data for each of the 42,280 pixels was clustered using the multivariate $\mathit{K}$-shape clustering method based on the shape similarity of the temperature dependence. In addition, to obtain the structural and ferroelectric phase transition temperatures, $T_{\mathrm{c}}$ and $T_{\mathrm{F}}$, hierarchical Bayesian temperature-series estimation was performed at each pixel (as a lower level) constrained over the entire cluster (as a higher level) considering the measurement error. Consequently, the K-shape clustering method revealed four clusters corresponding to elongated ferroelectric domains, explained by slight differences in retardance and fast-axis direction. Statistical analysis of the Bayesian posterior probability distribution showed a uniform distribution of $T_{\mathrm{c}}$ over the sample, but an inhomogeneous distribution of $T_{\mathrm{F}}$. The higher $T_{\mathrm{F}}$ regions exhibited a concentration of stress and/or strain. The Pearson correlation coefficient calculations suggested a strong to moderate relationship between the distribution of T_F and the ferroelectric state, while the correlation between T_c and the ferroelectric state was weak or nonexistent. The combination of machine learning and statistics provides a more reliable and less arbitrary approach to analyzing temperature-series data. These multilevel analyses are particularly useful in studying critical phenomena near phase transition temperatures in condensed matter physics.

GRAPHICAL ABSTRACT

IMPACT STATEMENT

Temperature-series analysis, combined with statistics, provides deeper insights into numerous imaging data observed at successive temperatures. This method will bring innovation to study of critical phenomena in condensed matter physics.

KEYWORDS:

• Temperature-series analysis
• $\mathit{K}$-shape clustering
• hierarchical Bayesian estimation
• birefringence image
• stress-induced ferroelectricity
• SrTiO${ }_3$

1. Introduction

In condensed matter physics, one of the most important research topics is investigation of the complex relationship between temperature ($\mathit{T}$) and physical properties, particularly critical phenomena occurring near phase transition temperatures [1]. These phenomena gain an additional layer of complexity when intertwined with changes in external parameters, such as magnetic field, electric field, and pressure. An increase in the number of external parameters generates complex, multivariate data structures that escalate in volume and complexity [2–4]. Traditional analysis methods, which are typically limited to examining variables in isolation, struggle to unravel these complex interactions. To bridge this gap, we have proposed an innovative machine learning approach for analyzing $\mathit{T}$-dependent data, which is optimized for the study of critical phenomena in condensed matter physics.

Since the 2000s, remarkable progress has been made in imaging techniques for the detection of spatial information regarding physical properties [5–10]. These developments are essential for the study of functional materials and focus on the formation and dynamics of macroscopic structures, such as domains or domain walls that affect magnetic, electrical, and transport properties. Therefore, advanced imaging techniques have become a fundamental part of materials research; however, such techniques have triggered an era of data explosion [6–8]. For example, spectral imaging collects the wavelength ($\mathit{\lambda }$) or energy dependence at each pixel to obtain a single image [2,4,5]. Extracting meaningful insights from a large dataset remains challenging despite advances in measurement techniques. Traditional analysis methods, based on averaging data from selective regions of interest, are limited in scope and prone to arbitrary bias. They often fail to analyze the overall trends of complex distributions and may miss local but important details. This generates an opportunity for machine learning, i.e. by applying advanced algorithms, we can extract deeper insights from a large dataset that previously remained inaccessible [11]. Pioneering research in condensed matter physics has involved the use of multivariate analysis methods based on statistics, i.e. principal component analysis [12,13], a multivariate curve resolution technique [14], and nonnegative matrix factorization [15], to analyze spectral imaging data. Although these methods are adept at extracting features by reducing multivariate data to a lower dimension, they prove insufficient at learning from data to make predictions or decisions for new data. Our machine learning approach is designed to analyze $\mathit{T}$-dependent imaging data, thus providing more accurate and flexible analysis than before.

We developed a novel machine learning approach for analyzing $\mathit{T}$-dependent birefringence imaging data, focusing on stress-induced ferroelectricity in SrTiO${ }_3$ [16–18]. Although such phase transition phenomena have been known for a long time [19–22], the high electrical conductivity caused by strain-induced ferroelectricity at the interface between SrTiO${ }_3$ and LaAlO${ }_3$ has attracted renewed attention [23–26]. Birefringence is a technique that uses the polarization properties of light to reveal both the internal structure and dielectric properties of materials. In this study, to effectively manage and interpret numerous birefringence images at different temperatures, we employ machine learning techniques that focus on unsupervised learning methods such as clustering algorithms [27]. These techniques are important for identifying underlying patterns and correlations not readily apparent in a large dataset. Pioneering research has been conducted on the application of machine learning to the analysis of reflection high-energy electron diffraction images [28–31]. Our previous study reported the application of a $\mathit{K}$-means clustering technique to the analysis of a birefringence image [18]. However, the analysis of $\mathit{T}$ dependence over numerous images remains a substantial challenge. In this study, we focus on the similarity between time-series and $\mathit{T}$-series data for efficient analysis of successive images.

Because $\mathit{T}$-series analysis methods have not yet been developed, our approach attempts to simply convert the time variables into the $\mathit{T}$ variables. Time-series analysis is a major topic in machine learning with enormous related knowledge available [32,33]; therefore, the application of these methods to $\mathit{T}$-series analysis can lead to considerable development in condensed matter physics. First, we recognize that series data have the following unique properties: (1) a distinct sequential order, which is essential for accurate interpretation; (2) a relational dynamic, where each data point depends on its neighbors; (3) a contextual essence, where the true value of the series is unlocked when viewed as a cohesive whole. With these insights, the $\mathit{T}$-dependent birefringence images of the stress-induced ferroelectricity in SrTiO${ }_3$ are analyzed. For spatial distribution analysis, we use the $\mathit{K}$-shape clustering method to identify the ferroelectric domain structures [34]. In traditional time-series analysis, Euclidean distance and dynamic time warping are commonly used although they are sensitive to variations in the data scale and can overlook similarities in shapes [32,33]. The $\mathit{K}$-shape clustering method is a specialized technique for clustering based on the shape of time-series data. Shape similarity is one of the most important features that forms the scaling law in condensed matter physics [1]. Furthermore, for the pixel-by-pixel $\mathit{T}$-series analysis, we use Bayesian $\mathit{T}$-series estimation to identify the structural and ferroelectric phase transition temperatures ($T_{\mathrm{c}}$ and $T_{\mathrm{F}}$) [35,36]. The Bayesian approach has an important advantage of directly incorporating prior knowledge and uncertainties such as measurement error directly into the modeling process, allowing us to perform comprehensive and probabilistic parameter inference based on the posterior probability distribution. Our combined machine learning approach toward spatial distributions and $\mathit{T}$ variations provides a deeper understanding of complex data structures in earlier impossible ways. This approach overcomes the limitations of traditional analysis methods in condensed matter physics and opens new avenues for exploring various critical phenomena, potentially revolutionizing the way researchers approach a large dataset.

2. Dataset for analysis

2.1. Birefringence images in SrTiO₃

Machine learning was applied to $\mathit{T}$-series birefringence imaging data for the stress-induced ferroelectricity of SrTiO${ }_3$. Under stress-free conditions, SrTiO${ }_3$ undergoes a structural phase transition at $T_{\mathrm{c}}=105$ K from the cubic paraelectric state to the tetragonal quantum paraelectric state, where the ferroelectric state is suppressed by quantum fluctuations [37–39]. The quantum paraelectric state undergoes ferroelectric phase transition upon the application of stress or strain [19–22,24–26,40–42]. For a SrTiO${ }_3$ thin film, strain-induced ferroelectricity appears at room temperature because of the lattice mismatch between the substrate and the thin film [24–26,40–42]. Similarly, for a SrTiO${ }_3$ bulk sample, when the external force ($\mathit{\sigma }$) is applied along $[001{]}_{\mathrm{c}}$, whose direction corresponds to the cubic phase, $T_{\mathrm{c}}$ increases from 105 K with increasing $\mathit{\sigma }$ [20–22]. At $\mathit{T}\le 30$ K, ferroelectric phase transition occurs, resulting in spontaneous polarization along $[010{]}_{\mathrm{c}}$; moreover, $T_{\mathrm{F}}$ increases with increasing $\mathit{\sigma }$ [20–22]. In the previous experimental studies of the stress-induced ferroelectricity in SrTiO${ }_3$, the value of $\mathit{\sigma }$ has been considered as stress because the external and internal forces ($=$ stress) are balanced. Unlike hydrostatic pressure experiments wherein $\mathit{\sigma }$ is applied uniformly in all directions, the stress caused by $\mathit{\sigma }$ is unevenly distributed throughout the sample. This is because $\mathit{\sigma }$ is applied from only one axial direction, and stress tends to be distributed perpendicularly to $\mathit{\sigma }$. Therefore, the stress-induced ferroelectric state may also be inhomogeneous because of the stress distribution. According to a study conducted by Müller et al. in Ref [19], when $\mathit{\sigma }\parallel [111{]}_{\mathrm{c}}$, the relationship between the converted stress (${\sigma }_{\mathrm{Cv}}$) and $T_{\mathrm{c}}$ is given by

(1) ${\sigma }_{\mathrm{Cv}}=30\times \left(T_{\mathrm{c}}-105\right)\left(\mathrm{MPa}\right).$(1)

Reportedly, this relationship also holds when $\mathit{\sigma }\parallel [001{]}_{\mathrm{c}}$ [16]. In this study, we used imaging techniques to evaluate the local stress and ferroelectric domain structure from the distributions of $T_{\mathrm{c}}$ and $T_{\mathrm{F}}$ on a pixel-by-pixel basis.

We obtained $\mathit{T}$-dependent birefringence images of the SrTiO${ }_3(110)$ substrate at $\mathit{\sigma }=231$ MPa applied along $[001{]}_{\mathrm{c}}$ when $\mathit{T}$ decreases from 300.0 to 14.1 K, as already reported in Ref [16]. To fully understand the stress-induced ferroelectric phase transition, our analysis focuses on capturing the inhomogeneous distribution of critical phenomena in the substrate. Because the SrTiO${ }_3(110)$ substrate at $\mathit{\sigma }=231$ MPa belongs to the plastic region, the linear relationship between stress and strain disappears [17]. According to theoretical studies on SrTiO${ }_3$ thin films [43–47], the lattice distortion functions as a driving force for ferroelectricity and $T_{\mathrm{F}}$ is expected to increase with increasing strain. However, because our experimental method cannot help determine whether the observed birefringence is more influenced by stress or strain, we do not distinguish between stress and strain.

The measurement principles required for machine learning with respect to analyzing the $\mathit{T}$-series birefringence imaging data are described as follows. Incident white light with circular polarization expressed as the Stokes parameters ($S_1^{\mathrm{I}\mathrm{n}},S_2^{\mathrm{I}\mathrm{n}},S_3^{\mathrm{I}\mathrm{n}})\simeq (0,0,-1$) is passed through the SrTiO${ }_3$(110) substrate along $[110{]}_{\mathrm{c}}$, and the polarization state of the transmitted light is obtained as ($S_1^{\mathrm{O}\mathrm{u}\mathrm{t}},S_2^{\mathrm{O}\mathrm{u}\mathrm{t}},S_3^{\mathrm{O}\mathrm{u}\mathrm{t}}$) [48,49]. The polarization state at each point of $384\times 288$ pixels can be obtained through a single measurement using a CCD camera with polarizing plates. When the refractive indices along the orthogonal principal optical axes are defined as $n_1$ and $n_2$ for the slow and fast axes with $n_1\ge n_2\ge 1$, respectively, the birefringence ($\mathrm{\Delta }\mathit{n}$) and retardance ($\mathit{\delta }$) are expressed as follows:

(2) $\mathrm{\Delta }\mathit{n}\equiv n_1-n_2,$(2)

(3) $\mathit{\delta }\equiv \mathrm{\Delta }\mathit{n}\times \mathit{t},$(3)

where $\mathit{t}$ represents the sample thickness ($\mathit{t}=0.333$ mm), and the fast-axis direction ($\mathit{\psi }$) corresponds to the direction for $n_2$. To elucidate the large value of $\mathit{\delta }$ exceeding $\mathit{\lambda }/2$, we observed the transmitted light separated into three different $\mathit{\lambda }$s at 523, 543, and 575 nm using bandpass filters [50]. In case of the SrTiO${ }_3(110)$ substrate at $\mathit{\sigma }=231$ MPa, the value of $\mathit{\delta }$ helps evaluate the amount of spontaneous polarization due to the Pockels effect [51,52] and stress and/or strain due to the photoelastic effect [16,17]. In addition, the $\mathit{\psi }$ direction reflects the crystal orientation; therefore, its distribution image can be used to infer the distribution of the spontaneous polarization direction [53–56]. As reported in Ref [16], Figures 1(a,b) show images of $\mathit{\delta }$ and $\mathit{\psi }$ at 130.9 K, corresponding to the cubic paraelectric state. Notably, these images are converted at $\mathit{\lambda }=543$ nm. Stripe patterns due to the slip planes appear even at 300.0 K, and their directions coincide with the directions of the projection of $⟨111{⟩}_{\mathrm{c}}$ onto $(110)$. Figure 1(c) shows an image of the $⟨111{⟩}_{\mathrm{c}}$ stripe pattern; the origin of this image can be attributed to the fact that the fast-axis directions indicated by the arrows differ by a few degrees between the elongated domains. These images show considerably similar distributions to those observed at 30.0 (paraelectric tetragonal state) and 14.1 K (ferroelectric state) shown in Figure 1 of Ref [18]. Because the $\mathit{\delta }$ and $\mathit{\psi }$ values correspond to a certain cross-section of the Poincaré sphere comprising the Stokes parameters [18,48,49], the analysis results from different approaches may appear different. Furthermore, it is challenging to quantitatively assess the amount and direction of spontaneous polarization by simultaneously examining the $\mathit{\delta }$ and $\mathit{\psi }$ distribution images. Therefore, a multivariate analysis using machine learning must be developed to mine the hidden information in $\mathit{T}$-series birefringence imaging data.

Figure 1. In SrTiO${ }_3$(110) substrate at an external force of $\mathit{\sigma }=231$ MPa applied along $[001{]}_{\mathrm{c}}$, images at 130.9 K of (a) retardance $\mathit{\delta }$ and (b) fast-axis direction $\mathit{\psi }$. The large open arrows schematically indicate $\mathit{\sigma }$. In (a), the cross marks indicate the analysis position, as shown in Figure 4. The order is defined from the left end to the right end as a1 to a7. (c) Schematic illustration of the origin of the stripe pattern along $⟨111{⟩}_{\mathrm{c}}$. The arrows are the fast-axis directions and differ by a few degrees between the elongated domains.

In our previous research, reported in Ref [18], we applied the $\mathit{K}$-means clustering method to the birefringence imaging data to observe the ferroelectric domains in SrTiO${ }_3$ at 231 MPa. This approach involved subtraction between the 30.0-K and 14.1-K datasets on a pixel-by-pixel basis, and it effectively simplified the $\mathit{T}$-series data for spatial clustering. Our results identified three distinct clusters (D1–D3, see Supplementary materials), each with unique characteristics. Cluster D1 exhibited large spontaneous polarization and disordered crystal orientation, indicating stress and/or strain concentration. Conversely, cluster D3 showed small spontaneous polarization but a more uniform distribution and less disordered crystal orientation. Furthermore, cluster D2 showed intermediate properties. This difference between clusters highlights the complex state of the stress-induced ferroelectricity in SrTiO${ }_3$. However, our method could not adequately visualize the domain structures along $⟨111{⟩}_{\mathrm{c}}$; this limitation is attributed to the reduction of $\mathit{T}$-series information in the clustering process. Therefore, there is an urgent need for advanced analysis techniques that comprehensively integrate spatial information and $\mathit{T}$-dependent data. The development of such methods could considerably improve the accuracy of visualization of ferroelectric domain structures and provide deeper insights into the complex $\mathit{T}$ dependence of birefringence in case of SrTiO${ }_3$. This approach will provide a new perspective for the analysis of critical phenomena in condensed matter physics.

2.2. Preprocessing of birefringence imaging data

We preprocessed our data sets to efficiently perform machine learning without sacrificing accuracy. Focusing on the SrTiO${ }_3(110)$ substrate, birefringence imaging data at each $\mathit{T}$ were extracted at 42,280 pixels (302$\times$140 pixels). As previously described in Ref [18], the Stokes parameters ($S_1^{\mathrm{I}\mathrm{n}},S_2^{\mathrm{I}\mathrm{n}},S_3^{\mathrm{I}\mathrm{n}}$) and ($S_1^{\mathrm{O}\mathrm{u}\mathrm{t}},S_2^{\mathrm{O}\mathrm{u}\mathrm{t}},S_3^{\mathrm{O}\mathrm{u}\mathrm{t}}$) at each $\mathit{T}$ were averaged over a 3$\times$3 pixel area using the moving average method. Because the sample position shifts with decreasing $\mathit{T}$, the analysis locations must be adjusted by one pixel unit and averaging is most effective in suppressing the effect of deviations of less than one pixel. Using these averaged Stokes parameters, we subsequently defined the important variables $\mathit{\theta }$ and $\mathit{\varphi }$, setting the stage for our subsequent analysis on the Poincaré sphere [18]. The angle $\mathit{\theta }$, corresponding to $\mathit{\delta }$, was defined within the range $[0^{\circ },{180}^{\circ }]$ with a periodicity of ${360}^{\circ }$. The axial direction $\mathit{\varphi }$, reflecting the direction of the fast or slow axis, was restricted to $[0^{\circ },{180}^{\circ })$, also with a periodicity of ${180}^{\circ }$. With these definitions, we generally expressed the retardance as $\mathit{\delta }=|\mathit{l}\pm \mathit{\theta }÷{360}^{\circ }|\times \mathit{\lambda }$, where $\mathit{l}=0,1,2,\cdots$, and the fast-axis direction as $\mathit{\psi }=(\mathit{\varphi }\pm {90}^{\circ })/2$. Whenever $\mathit{\delta }$ exceeds $\mathit{\lambda }/2,\mathit{\lambda },3\mathit{\lambda }/2,\cdots$, the angle $\mathit{\theta }$ folds back at 0${ }^{\circ }$ or 180${ }^{\circ }$; therefore, the variables for $\mathit{l}$ and $\pm$ change concurrently. The accuracy of machine learning predictions can be improved using the $\mathit{T}$-series data for the three $\mathit{\lambda }$s without corrections for these variables on a pixel-by-pixel basis [18].

Because the variables $\mathit{\theta }$ and $\mathit{\varphi }$ are angular measures, we adopted circular statistics for our analysis [57,58]. This approach allows us to consider angles as vectors, which simplifies their computation and interpretation. From Ref [18], the mean (or maximum likelihood estimate) of $\mathit{\theta }$ and $\mathit{\varphi }$ is expressed as follows:

(4) $cos\bar{\theta }=\sum _{\mathit{i}=1}^m\left(\frac{cos{\theta }_i}{{\bar{R}}_{\mathit{\theta }}}\right)\times \frac{1}{m},$(4)

(5) $sin\bar{\theta }=\sum _{\mathit{i}=1}^m\left(\frac{sin{\theta }_i}{{\bar{R}}_{\mathit{\theta }}}\right)\times \frac{1}{m},$(5)

(6) $cos2\bar{\varphi }=\sum _{\mathit{i}=1}^m\left(\frac{cos2{\varphi }_i}{{\bar{R}}_{\mathit{\varphi }}}\right)\times \frac{1}{m},$(6)

(7) $sin2\bar{\varphi }=\sum _{\mathit{i}=1}^m\left(\frac{sin2{\varphi }_i}{{\bar{R}}_{\mathit{\varphi }}}\right)\times \frac{1}{m},$(7)

(8) ${\bar{R}}_{\mathit{\theta }}=\frac{\sqrt{{\left({\sum }_{i=1}^{m}cos{\theta }_i\right)}^2+{\left({\sum }_{i=1}^{m}sin{\theta }_i\right)}^2}}{\mathit{m}},$(8)

(9) ${\bar{R}}_{\mathit{\varphi }}=\frac{\sqrt{{\left({\sum }_{i=1}^{m}cos2{\varphi }_i\right)}^2+{\left({\sum }_{i=1}^{m}sin2{\varphi }_i\right)}^2}}{\mathit{m}},$(9)

where $\mathit{m}$ denotes the number of data points and $\bar{R}$ represents the mean resultant length. The value of $\bar{R}$ varies between 0 and 1 and indicates the alignment of the distribution, i.e. values closer to 1 signify a well-aligned distribution. The maximum likelihood estimation of $\mathit{\theta }$ and $\mathit{\varphi }$ can be reformulated as follows:

(10) $tan\bar{\theta }=\frac{{\sum }_{i=1}^{m}sin{\theta }_i}{{\sum }_{i=1}^{m}cos{\theta }_i},$(10)

(11) $tan2\bar{\varphi }=\frac{{\sum }_{i=1}^{m}sin2{\varphi }_i}{{\sum }_{i=1}^{m}cos2{\varphi }_i}.$(11)

This allows angles to be treated as linear statistics [57,58]. Consequently, the four variables ($cos{\theta }_i/{\bar{R}}_{\mathit{\theta }}$, $sin{\theta }_i/{\bar{R}}_{\mathit{\theta }}$, $cos2{\varphi }_i/{\bar{R}}_{\mathit{\varphi }}$, and $sin2{\varphi }_i/{\bar{R}}_{\mathit{\varphi }}$) are derived from ${\theta }_i$ and ${\varphi }_i$ at the $\mathit{i}$-th measurement point. After successfully applying circular statistics, we extended this procedure to our datasets for three different $\mathit{\lambda }$s, creating a 42,280$\times$12 matrix for a birefringence image. In this study, it is necessary to further develop a method for analyzing the numerous images observed at successive temperatures. Because all $\mathit{T}$-series data are standardized for each variable before machine learning, as discussed later, ${\bar{R}}_{\mathit{\theta }}$ and ${\bar{R}}_{\mathit{\varphi }}$ are not necessary for analyzing the $\mathit{T}$ dependence only at a pixel. However, the $\mathit{K}$-shape clustering would be more reliable by preserving information regarding the distribution within the substrate. Therefore, as a preprocessing of the 12-variable $\mathit{T}$-series data for 42,280 pixels, we calculated ${\bar{R}}_{\mathit{\theta }}$ and ${\bar{R}}_{\mathit{\varphi }}$ for each birefringence image. These careful preprocessing steps help accurately analyze our $\mathit{T}$-series birefringence imaging data and set the stage for insightful machine learning.

The stress apparatus made of Cu thermally shrinks with decreasing $\mathit{T}$, causing a shift in the sample position [59]. With 300.0 K as reference, the birefringence image at 14.1 K showed shifts of 23 pixels to the left and 4 pixels to the bottom. Therefore, we manually corrected the sample position by one pixel unit in all images. Previously, we selected as homogeneous a region of the sample as possible and averaged the data within a 10–30$\times$10–30-pixel area to ensure that the shift of less than one pixel was negligible [16]. However, since the $\mathit{T}$-series data is analyzed on a pixel-by-pixel basis, the 3$\times$3 pixel averaging preprocessing proved insufficient in regions where the amount of birefringence changes abruptly. In these regions, because small jumps appear in $\mathit{\delta }(\mathit{T})$ and $\mathit{\psi }(\mathit{T})$ curves, conventional time-series analysis methods, which mainly focus on the temporal differences of the changes, would cause serious problems [32,33]. In this study, we adopted the multivariate $\mathit{K}$-shape clustering method [34] because it is robust to changes in data scale and time deviations, effectively classifying data into the same cluster based on shape similarity. It also maximizes the cross-correlation within a cluster, increasing the likelihood of grouping similar data. Furthermore, for the $\mathit{T}$-series analysis using the $\mathit{\delta }(\mathit{T})$ and $\mathit{\psi }(\mathit{T})$ curves, such jumps will inevitably have a negative impact on the estimation of $T_{\mathrm{c}}$ and $T_{\mathrm{F}}$ values. Because it is impossible to check the $\mathit{T}$-series data for all 42,280 pixels, the Bayesian $\mathit{T}$-series estimation should be used, considering error factors based on certain criteria [35,36]. Using the Bayesian posterior probability distribution, the confidence and accuracy of the estimates can be quantitatively understood, allowing discussion of the relationship between the stress and/or strain concentration and the ferroelectric state.

3. Data analysis procedure

Our computer was equipped with an Intel(R) Core(TM) i9-13900K CPU running up to 5.80 GHz and implemented 128 GB of memory. To meet the demands of datasets larger than 20 GB, we avoided data input/output bottlenecks on hard disk drives and increased data throughput many times by installing a solid-state drive. For the $\mathit{K}$-shape clustering analysis, we used the tslearn package (version 0.5.3.2), a common tool for time-series analysis in the Python environment (version 3.11.1). To determine the optimal number of clusters $\mathit{k}$, a combination of methods was used, including the elbow method, silhouette method, and gap statistic [60], all of which are integrated in the tslearn package. The computation time, inherently dependent on $\mathit{k}$, was approximately half a day for computing nine different cluster types ranging from $\mathit{k}=2$ to 10. For Bayesian $\mathit{T}$-series estimation, we used the computational power of the Hamiltonian Monte Carlo (HMC) algorithm, a sophisticated technique for Markov chain Monte Carlo (MCMC) simulations, run through the rstan package (version 2.26.23) in the R environment (version 4.1.2). To efficiently process the Bayesian $\mathit{T}$-series estimation for each of the 42,280 pixels, we employed a strategy that used the parallel package (version 4.1.2), foreach package (version 1.5.2), and doParallel package (version 1.0.17) for parallel computations. This methodology helped us take advantage of parallel processing on up to 32 cores. The computational time depended not only on the number of iterations and chains in the MCMC simulations, but also on the settings of the prior distribution. After extensive optimization, we could estimate the fast-axis directions in a few minutes, as well as $T_{\mathrm{F}}$ and $T_{\mathrm{c}}$ in approximately a day and a week, respectively. By applying $\mathit{T}$-series analysis using hierarchical Bayesian $\mathit{T}$-series estimation, we aim to complement our spatial data analysis using $\mathit{K}$-shape clustering and provide a more holistic understanding of stress-induced ferroelectricity.

Our approach to $\mathit{T}$-series analysis, which relates changes in sample $\mathit{T}$ to changes in time, presents unique challenges and opportunities for deeper insight. To address these challenges, we focus on statistical significance in time-series analysis. In time-series analysis, direct tests of statistical significance often require complex procedures tailored to the unique characteristics of the data. Nevertheless, in our approach, we demonstrate the utility of the exploratory values obtained from the MCMC simulations within the Bayesian framework. These values provide an indirect but an effective approach to assess differences between groups by providing insights into posterior distributions. Because of the large sample size in the MCMC simulations, the $\mathit{t}$ test and the Wilcoxon rank sum test are prone to type 1 errors, but Cohen’s $\mathit{d}$ remains robust to them [61,62]. Cohen’s $\mathit{d}$ with Hedges’ correction was calculated using the effsize package (version 0.8.1) in the R environment. Cohen classified effect sizes as small ($|\mathit{d}|=0.2$), medium ($|\mathit{d}|=0.5$), and large ($|\mathit{d}|\ge 0.8$). According to Cohen’s guidelines, $0.2\le |\mathit{d}|<0.5$ indicates a subtle but noticeable difference, $0.5\le |\mathit{d}|<0.8$ indicates a moderate and observable difference, and $|\mathit{d}|\ge 0.8$ indicates a substantial and practically significant difference between two groups [61,63]. Statistical treatment using Cohen’s $\mathit{d}$ with Hedges’ correction is important for drawing robust conclusions from $\mathit{T}$-dependent patterns. The normal-based 95% confidence interval calculated using the normal distribution method is based on the standard error and the assumption that our estimates follow a normal distribution. However, because the data may not follow a parametric distribution, we used the bootstrap 95% confidence interval, which is derived by bootstrap resampling and does not assume a specific distribution. This analysis was performed in the R environment using the boot package (version 1.3.28).

4. Results and discussion

Figure 2 shows the time and measuring point evolution of sample $\mathit{T}$. As described in Ref [16], our measurements continuously recorded 3,362 sheets $\times$ 3 $\mathit{\lambda }$s of the birefringence imaging data from 300.0 to 14.1 K while carefully controlling the decrease of $\mathit{T}$ at a fixed rate of $-0.25$ K/min. Because of the consistent relationship between $\mathit{T}$ and time, which accurately reflects the order and interval characteristics of the series data, it is feasible to apply the $\mathit{T}$-series birefringence imaging data to machine learning using $\mathit{K}$-shape clustering and hierarchical Bayesian $\mathit{T}$-series estimation.

Figure 2. (a) Sample temperature variation with time and measurement points. Regions of insufficient temperature control, (b) affecting the determination of the structural phase transition temperature, and (c) affecting the determination of ferroelectric phase transition temperature.

4.1. K-shape clustering

To comprehensively characterize the spatial information, the $\mathit{T}$-series birefringence imaging data were analyzed using the $\mathit{K}$-shape clustering method. Figure 3 shows the $\mathit{T}$ dependence of ${\bar{R}}_{\mathit{\theta }}$ and ${\bar{R}}_{\mathit{\varphi }}$ calculated at each $\mathit{T}$. The ${\bar{R}}_{\mathit{\theta }}(\mathit{T})$ curves indicate that the distribution of $\mathit{\delta }$ becomes large below 30 K because ferroelectric phase transition occurs and spontaneous polarization is inhomogeneously distributed over the substrate, as reported in Ref [18]. The ${\bar{R}}_{\mathit{\varphi }}(\mathit{T})$ curves for different $\mathit{\lambda }$s show a remarkable difference near $T_{\mathrm{c}}$ and $T_{\mathrm{F}}$. Furthermore, the ${\bar{R}}_{\mathit{\theta }}(\mathit{T})$ curves also seem to slightly differ in the same $\mathit{T}$ regions. This is attributed to the $\mathit{\theta }$ folding back at 180${ }^{\circ }$ when $\mathit{\delta }$ exceeds $\mathit{\lambda }/2$ at 523 and 543 nm. Figure 4 shows the $\mathit{T}$ dependence of the 4-variable $\times$ 3-$\mathit{\lambda }$ for each pixel, specifically, the seven pixels shown in Figure 1(a), selected from a total of 42,280 pixels. These curves show numerous small jumps due to the continuous shift of the sample position. Based on a comparison of Figures 3 with 4, owing to the considerably small variation in ${\bar{R}}_{\mathit{\theta }}$ and ${\bar{R}}_{\mathit{\varphi }}$ below 30 K, the shape of the $\mathit{T}$-dependent curves of 12 variables divided by their coefficients remains almost unchanged. Furthermore, around $T_{\mathrm{c}}$, ${\bar{R}}_{\mathit{\varphi }}$ changes considerably because of the folding back of $\mathit{\theta }$. These coefficients tend to effectively expand the difference between the pixels in the $\mathit{T}$-series data. The ${\bar{R}}_{\mathit{\theta }}$ and ${\bar{R}}_{\mathit{\varphi }}$ estimates are expected to further improve the accuracy of the $\mathit{K}$-shape clustering process. However, as shown in Figures 2(b,c), it is difficult to completely control sample $\mathit{T}$ owing to the limitations of experimental techniques, resulting in small fluctuations within certain $\mathit{T}$ ranges. Such fluctuations temporarily impose a nonintrinsic external force on the substrate, resulting in small birefringence anomalies that can be observed across all pixels. This is a considerable difference between the $\mathit{T}$-series data and the time-series data. In conventional time-series analysis, the order can easily affect the results; therefore, the $\mathit{K}$-shape clustering method, which focuses on shape similarity, is preferred in condensed matter physics [32–34]. Although there is no intrinsic thermal hysteresis in SrTiO${ }_3$, we entered machine learning in its measured order in this study.

Figure 3. Temperature dependence of the mean resultant length (a) ${\bar{R}}_{\mathit{\theta }}$ and (b) ${\bar{R}}_{\mathit{\varphi }}$ for $\mathit{\lambda }=523$, 543, and 575 nm.

Figure 4. Four variables ($cos{\theta }_i/{\bar{R}}_{\mathit{\theta }}$, $sin{\theta }_i/{\bar{R}}_{\mathit{\theta }}$, $cos2{\varphi }_i/{\bar{R}}_{\mathit{\varphi }}$, and $sin2{\varphi }_i/{\bar{R}}_{\mathit{\varphi }}$) as a function of temperature at (a–d) $\mathit{\lambda }=523$ nm, (e–h) 543 nm, and (i–l) 575 nm. Analysis positions a1 to a7 and the corresponding cluster numbers to which they belong.

To optimize computational efficiency, we preprocessed the 42,280-pixel $\times$ 12-variable $\mathit{T}$-series data between 130.9 and 14.1 K by averaging every four $\mathit{T}$ points, resulting in each $\mathit{T}$-series containing 346 points with an approximate step of 0.34 K. For the multivariate $\mathit{K}$-shape clustering method, determining the number of clusters $\mathit{k}$ in advance is imperative. Figure 5 shows the results of identifying the optimal values of $\mathit{k}$ using the elbow method, silhouette method, and gap statistic [60]. Typically, these methods rarely yield consistent results because each emphasizes different aspects of cluster cohesion and separation. In other words, the elbow method focuses on the variance of the data in the clusters, the silhouette method focuses on the high cohesion in the clusters and the large separation between the clusters, and the gap statistic technique focuses on the difference between the random and real datasets. In this study, because the elbow and silhouette methods indicate $\mathit{k}=4$, as shown in Figure 5, we selected this value. Although these indicators shown in Figure 5 serve as guidelines [60], re-evaluating the validity of $\mathit{k}=4$ after clustering is essential. Figure 6 shows the division of the 12-variable $\mathit{T}$-series data into four groups (E1–E4), revealing the elongated clusters along $⟨111{⟩}_{\mathrm{c}}$ corresponding to the dislocation. Comparing the multivariate $\mathit{K}$-shape clustering shown in Figure 6 with the multivariate $\mathit{K}$-means clustering (D1-D3, see Supplementary materials) [18], we found that clusters E1 and E2 coincide with clusters D1 and D2, indicating large spontaneous polarization in these regions due to stress and/or strain concentration. Conversely, clusters E3 and E4, which correspond to cluster D3, could be regions of low spontaneous polarization and less disordered crystal orientation. Table 1 shows the number of pixels for each cluster obtained from the $\mathit{K}$-means and $\mathit{K}$-shape clustering analyses. This result also corroborates that E1 and E2 correspond to D1 and D2, respectively, E3 and E4 correspond to D3.

Figure 5. Evaluation results for the optimal value of $\mathit{k}$ using the (a) elbow method, (b) silhouette method, and (c) gap statistic.

Figure 6. Result of the K-shape clustering obtained from 12 variables of temperature-series data between 130.9 and 14.1 K.

Table 1. Number of pixels in respective clusters determined by K-means and K-shape clustering methods. The K-means clustering data are referenced in Ref. [18].

K-means		K-shape
Cluster	Number	Cluster	Number
Total	42,280	Total	42,280
D1	794	E1	5,503
D2	10,110	E2	5,068
D3	31,376	E3	15,669
		E4	16,040

4.2. Analysis for each cluster

In Figure 6, it appears that domains E1 and E3 along $⟨111{⟩}_{\mathrm{c}}$ form one group, while domains E2 and E4 form another. The $\mathit{T}$ dependence of $\mathit{\delta }$ and $\mathit{\psi }$ was summed of each cluster to confirm the characteristics for each cluster. The averaged values of the retardance, ${\delta }^{\mathrm{Av}}$, and the fast-axis direction, ${\psi }^{\mathrm{A}v}$, for each cluster were obtained by averaging the difference between the incident and transmitted light of the Stokes parameters at each $\mathit{T}$, defined as $(S_1^{\mathrm{A}\mathrm{v}},S_2^{\mathrm{A}\mathrm{v}},S_3^{\mathrm{A}\mathrm{v}})$, as follows:

$S_1^{\mathrm{A}\mathrm{v}}\equiv \frac{{\sum }_{i=1}^m{\left(S_1^{\mathrm{O}\mathrm{u}\mathrm{t}}-S_1^{\mathrm{I}\mathrm{n}}\right)}_{\mathit{i}}}{\mathit{m}},$

(12) $S_2^{\mathrm{A}\mathrm{v}}\equiv \frac{{\sum }_{i=1}^m{\left(S_2^{\mathrm{O}\mathrm{u}\mathrm{t}}-S_2^{\mathrm{I}\mathrm{n}}\right)}_{\mathit{i}}}{\mathit{m}},$(12)

$S_3^{\mathrm{A}\mathrm{v}}\equiv \frac{{\sum }_{i=1}^m{\left(S_3^{\mathrm{O}\mathrm{u}\mathrm{t}}-S_3^{\mathrm{I}\mathrm{n}}\right)}_{\mathit{i}}}{\mathit{m}},$

where $(\cdots {)}_i$ represents the difference between the Stokes parameters at the $\mathit{i}$-th pixel. Figure 7 shows the $\mathit{T}$ dependence of ${\delta }^{\mathrm{Av}}$ and ${\psi }^{\mathrm{Av}}$ at 575 nm when $\mathit{\theta }$ and $\mathit{\varphi }$ are rewritten as follows:

(13) $\mathit{\delta }=\left(1-\frac{\theta }{{360}^{\circ }}\right)\times \mathit{\lambda },$(13)

(14) $\mathit{\psi }=\frac{\mathit{\varphi }+{90}^{\circ }}{2}.$(14)

Figure 7. Temperature dependence of the (a) averaged retardance, ${\delta }^{\mathrm{Av}}$, and (b) averaged fast-axis direction, ${\psi }^{\mathrm{Av}}$, for $\mathit{\lambda }=575$ nm. In (a), the ${\delta }^{\mathrm{Av}}(\mathit{T})$ curves for E1 and E2 overlap, and the curves for E3 and E4 overlap.

These curves are almost identical to those for 523 and 543 nm, except for the $\sim$$T_{\mathrm{c}}$ and $\sim$$T_{\mathrm{F}}$ regions. In these regions, it is difficult to correctly estimate the value of ${\delta }^{\mathrm{Av}}$ from $(S_1^{\mathrm{A}\mathrm{v}},S_2^{\mathrm{A}\mathrm{v}},S_3^{\mathrm{A}\mathrm{v}})$ because we cannot separate the two components of $\mathit{\theta }$ that fold back at 180${ }^{\circ }$ and do not fold back on a pixel-by-pixel basis. In this study, we focused our analysis on $\mathit{T}$-series data at 575 nm.

In Figure 7(a), the ${\delta }^{\mathrm{Av}}(\mathit{T})$ curves exhibit a similar shape regardless of the clusters, with two inflection points appearing at $\sim$$T_{\mathrm{c}}$ and $\sim$$T_{\mathrm{F}}$. In addition, the values of ${\delta }^{\mathrm{Av}}$ are overall higher for E1 and E2 than for E3 and E4 because of stress and/or strain concentration. However, for E1 and E2, the direction of spontaneous polarization may be disturbed because the crystal orientations are disturbed, as discussed in Ref [18]. The slopes of the ${\delta }^{\mathrm{Av}}(\mathit{T})$ curves at $\mathit{T}<T_{\mathrm{F}}$ (the ferroelectric state) appear to be steeper for E3 and E4 than for E1 and E2 because the increase in $\mathit{\delta }$ due to the Pockels effect may be proportional to the superposition of spontaneous polarizations resulting from disturbed crystal orientations. In Figure 7(b), the ${\psi }^{\mathrm{Av}}(\mathit{T})$ curves show numerous small jumps due to shifts in the sample position. In particular, the jumps of the ${\psi }^{\mathrm{Av}}(\mathit{T})$ curves for E1 and E2 exhibit opposite trends with almost identical magnitudes in the same $\mathit{T}$ range. As shown in Figure 6, the effect of less than one pixel misalignment is pronounced because the vertically elongated domains for E1 and E2 alternate along the horizontal direction where the displacement is large. From a crystallographic viewpoint, it is unlikely that the optical principal axis rotates by a few degrees in the cubic or tetragonal phase with decreasing $\mathit{T}$. However, if the crystal orientations are disturbed by stress and/or strain concentration, the averaged direction may rotate slightly with decreasing $\mathit{T}$. Consequently, it remains inconclusive whether the ${\psi }^{\mathrm{Av}}(\mathit{T})$ curves vary with decreasing $\mathit{T}$. In this study, assuming that the values of ${\psi }^{\mathrm{Av}}$ remain constant with respect to $\mathit{T}$ between 130.9 and 14.1 K, we performed a Bayesian $\mathit{T}$-series estimation of the $\mathit{T}$-independent value of ${\psi }^{\mathrm{Av}}$ for each cluster considering the measurement error.

The HMC algorithm for the MCMC simulations was used to derive the $\mathit{T}$-independent value of ${\psi }^{\mathrm{Av}}$ for each cluster. To ensure a comprehensive exploration of the parameter space, 16 independent MCMC chains were run, each with different random initial values. A burn-in period of 1,000 iterations was used, followed by the collection of 100,000 samples for analysis. The prior probability and model Hamiltonian for ${\psi }^{\mathrm{Av}}$ are expressed as follows:

(15) $tan2{\psi }_n^{\mathrm{A}\mathrm{v}}\sim \mathit{N}({\hat{t}}_{2\mathit{\psi }},{{\mathit{\sigma }}_{\mathrm{obs}}}^2),$(15)

${\hat{t}}_{2\mathit{\psi }}\sim \mathit{N}(0,1^2),$

${\sigma }_{\mathrm{obs}}\sim \left|\mathit{N}{(0,1}^2)\right|,$

where the normal distribution $\mathit{N}(\bar{x},{\sigma }^2)$ in linear statistics is a mean $\bar{x}$ and a standard deviation $\mathit{\sigma }\ge 0$. The variable ${\psi }_n^{\mathrm{A}\mathrm{v}}$ denotes the experimental data of ${\psi }^{\mathrm{Av}}$ at the $\mathit{n}$-th $\mathit{T}$ point, while ${\hat{t}}_{2\mathit{\psi }}$ serves as a parameter in the MCMC simulations. Using circular statistics, the mean (maximum likelihood estimate) of the angles is expressed as shown in Equations (10)(10) $tan\bar{\theta }=\frac{{\sum }_{i=1}^{m}sin{\theta }_i}{{\sum }_{i=1}^{m}cos{\theta }_i},$(10) and (11(11) $tan2\bar{\varphi }=\frac{{\sum }_{i=1}^{m}sin2{\varphi }_i}{{\sum }_{i=1}^{m}cos2{\varphi }_i}.$(11) ). Because of the summation of trigonometric functions, it is not easy for the Bayesian approach to successively evaluate the error of each measurement.

To incorporate the error in each measurement into the Bayesian $\mathit{T}$-series estimation, the derivation process of $\mathit{\varphi }$ and $\mathit{\psi }$ should be considered. As explained in Ref [18], the angle $\mathit{\varphi }$ is represented by a periodicity of 180${ }^{\circ }$ in the $S_1{S}_2$-plane. Because all data between 130.9 and 14.1 K, regardless of the clusters, are concentrated around $\mathit{\varphi }\simeq {90}^{\circ }$ in the third and fourth quadrants, as shown schematically in Figure 8, the angles of ${\varphi }_{+}$ and ${\varphi }_{-}$ can be calculated as follows: if $S_1\ge 0$ and $S_2<0$,

(16) ${\varphi }_{+}={180}^{\circ }+arctan\left(\frac{S_2^{\mathrm{O}\mathrm{u}\mathrm{t}}-S_2^{\mathrm{I}\mathrm{n}}}{S_1^{\mathrm{O}\mathrm{u}\mathrm{t}}-S_1^{\mathrm{I}\mathrm{n}}}\right),$(16)

Figure 8. Schematic drawing of the Stokes parameters $S_1^{\mathrm{O}\mathrm{u}\mathrm{t}}$ and $S_2^{\mathrm{O}\mathrm{u}\mathrm{t}}$ in the $S_1{S}_2$-plane. All our experimental data for $\mathit{\lambda }=575$ nm between 130.9 and 14.1 K show $\mathit{\varphi }\simeq {90}^{\circ }$ in the third and fourth quadrants. The definitions of ${\varphi }_{+}$ and ${\varphi }_{-}$ differ depending on the sign of $S_1^{\mathrm{O}\mathrm{u}\mathrm{t}}$ as in Equations (16) and (17).

while if $S_1<0$ and $S_2<0$,

(17) ${\varphi }_{-}=arctan\left(\frac{S_2^{\mathrm{O}\mathrm{u}\mathrm{t}}-S_2^{\mathrm{I}\mathrm{n}}}{S_1^{\mathrm{O}\mathrm{u}\mathrm{t}}-S_1^{\mathrm{I}\mathrm{n}}}\right).$(17)

However, since “$\mathit{arctan}$’’ is a multivalued function, it is necessary to adjust $\pm {180}^{\circ }$ to match the angle as shown in Figure 8. Using Equation (14)(14) $\mathit{\psi }=\frac{\mathit{\varphi }+{90}^{\circ }}{2}.$(14) , Equations (16)(16) ${\varphi }_{+}={180}^{\circ }+arctan\left(\frac{S_2^{\mathrm{O}\mathrm{u}\mathrm{t}}-S_2^{\mathrm{I}\mathrm{n}}}{S_1^{\mathrm{O}\mathrm{u}\mathrm{t}}-S_1^{\mathrm{I}\mathrm{n}}}\right),$(16) and (17(17) ${\varphi }_{-}=arctan\left(\frac{S_2^{\mathrm{O}\mathrm{u}\mathrm{t}}-S_2^{\mathrm{I}\mathrm{n}}}{S_1^{\mathrm{O}\mathrm{u}\mathrm{t}}-S_1^{\mathrm{I}\mathrm{n}}}\right).$(17) ) can be rewritten as follows:

(18) $tan{\varphi }_{\pm }=-\frac{1}{tan2\mathit{\psi }}.$(18)

Because the angle $\mathit{\varphi }$ has the ${180}^{\circ }$ periodicity, to avoid “$\mathit{tangent}$’’ singularities, the relationship between the Stokes parameters and $\mathit{\psi }$ expressed as

(19) $tan2{\psi }_n^{\mathrm{A}\mathrm{v}}=-\frac{{\left[{\sum }_i{\left(S_1^{\mathrm{O}\mathrm{u}\mathrm{t}}-S_1^{\mathrm{I}\mathrm{n}}\right)}_{\mathit{i}}\right]}_{\mathit{n}}}{{\left[{\sum }_i{\left(S_2^{\mathrm{O}\mathrm{u}\mathrm{t}}-S_2^{\mathrm{I}\mathrm{n}}\right)}_{\mathit{i}}\right]}_{\mathit{n}}}\simeq 0,$(19)

where $[\cdots {]}_n$ is the sum of the differences in the Stokes parameters at the $\mathit{n}$-th $\mathit{T}$ point. In our birefringence imaging measurements, the variables $(S_1^{\mathrm{O}\mathrm{u}\mathrm{t}},S_2^{\mathrm{O}\mathrm{u}\mathrm{t}},S_3^{\mathrm{O}\mathrm{u}\mathrm{t}})$ at all pixels are identified simultaneously in a single measurement, assuming the condition $(S_1^{\mathrm{O}\mathrm{u}\mathrm{t}}{)}^2+(S_2^{\mathrm{O}\mathrm{u}\mathrm{t}}{)}^2+(S_3^{\mathrm{O}\mathrm{u}\mathrm{t}}{)}^2=1$ for the respective pixels. Furthermore, as mentioned in $\mathcal{\S }$ 2.1, the incident light was set to $(S_1^{\mathrm{I}\mathrm{n}},S_2^{\mathrm{I}\mathrm{n}},S_2^{\mathrm{I}\mathrm{n}})\simeq (0,0,-1)$ [16]. Therefore, the measurement error of $(S_1^{\mathrm{O}\mathrm{u}\mathrm{t}}/S_2^{\mathrm{O}\mathrm{u}\mathrm{t}})$ in a single measurement at a certain $\mathit{T}$ is important, and the Bayesian $\mathit{T}$-series estimation based on linear statistics was performed in a simplified manner using Equation (19)(19) $tan2{\psi }_n^{\mathrm{A}\mathrm{v}}=-\frac{{\left[{\sum }_i{\left(S_1^{\mathrm{O}\mathrm{u}\mathrm{t}}-S_1^{\mathrm{I}\mathrm{n}}\right)}_{\mathit{i}}\right]}_{\mathit{n}}}{{\left[{\sum }_i{\left(S_2^{\mathrm{O}\mathrm{u}\mathrm{t}}-S_2^{\mathrm{I}\mathrm{n}}\right)}_{\mathit{i}}\right]}_{\mathit{n}}}\simeq 0,$(19) . The convergence of the MCMC simulations was assessed using trace plots and Gelman – Rubin statistics [64]. The trace plots showed effective mixing and overlap between the chains. As detailed in Supplementary materials, the effective sample size ($n_{\mathrm{eff}}$) derived from the MCMC simulations exceeded $\sim$40% of the total sample size, which consisted of 1,584,000 iterations. An $n_{\mathrm{eff}}$ lower than 10% of the total iterations typically indicates a high degree of autocorrelation, suggesting that the MCMC chain may not be exploring the parameter space efficiently. Conversely, our results with $n_{\mathrm{eff}}$ well above this threshold indicate that the MCMC simulations were sufficiently efficient, with lower autocorrelation, ensuring a robust exploration of the parameter space. The Gelman – Rubin diagnostic value ($\hat{R}$) was 1.00 for all three parameters including the log probability of the unconstrained parameters, indicating satisfactory convergence between the chains (see Supplementary materials).

Figure 9 shows the results of the analysis for ${\psi }^{\mathrm{Av}}$ to visually represent the posterior distributions for the MCMC simulations of ${\hat{t}}_{2\mathit{\psi }}$. Figure 9(a) shows a box and whisker plot of the exploration values in the fast-axis direction, expressed as ${\hat{\psi }}^{\mathrm{Av}}$ $(=(arctan{\hat{t}}_{2\mathit{\psi }})÷2\simeq {90}^{\circ }$) [65]. Figure 9(b) shows the probability density histogram for ${\hat{\psi }}^{\mathrm{Av}}$ for E3 and E4. The scales of the angular measures in these figures are very narrow, ranging from ${89.4}^{\circ }$ to ${91.7}^{\circ }$; thus, they can be treated as linear statistics. Therefore, these figures allow a direct comparison to detect any significant differences between the four clusters. Table 2 shows the five-number summary, i.e. minimum, first quartile $(\equiv Q_1)$, median, third quartile $(\equiv Q_3)$, and maximum, and mean for ${\hat{\psi }}^{\mathrm{Av}}$ of each cluster needed to construct the plot in Figure 9(a). In our analysis, outliers in the box and whisker plots are identified based on a common statistical criterion involving the interquartile range (IQR) defined as $(Q_3-Q_1)$. Exploration values outside the range from $(Q_1-1.5\times \mathrm{I}\mathit{QR})$ to $(Q_3+1.5\times \mathrm{I}\mathit{QR})$ are considered outliers [65]. This approach is widely accepted for outlier detection because it provides a balance between sensitivity to extreme values and robustness to the natural variability of the data. Although the exploration of the parameter space for ${\hat{\psi }}^{\mathrm{Av}}$ was set wide, i.e. $2{\hat{\psi }}^{\mathrm{Av}}\sim (180\pm arctan1.96)=$ 117–243${ }^{\circ }$, both the minimum and maximum values in Table 2 fell within a narrower range than the variation in the ${\psi }^{\mathrm{Av}}(\mathit{T})$ curves shown in Figure 7(b). In Bayesian estimation, the prior probability distribution is used to reduce the computational cost by narrowing the range of variables. However, if this range is extremely restrictive, it can limit the exploration of the parameter space, potentially leading to incomplete or biased results. Thus, it is crucial to balance the reduction in computational cost with the need for sufficient exploration. The relationship between exploration values and range limits is a key consideration in the Bayesian approach (see Supplementary materials). The difference in the mean of ${\hat{\psi }}^{\mathrm{Av}}$ between E1 and E2 is about 1.9${ }^{\circ }$, and this becomes the origin of the stripe pattern along $⟨111{⟩}_{\mathrm{c}}$, as shown in Figure 1(c). To further assess the significant difference in ${\hat{t}}_{2\mathit{\psi }}$ between E3 and E4, Cohen’s $\mathit{d}$ with Hedges’ correction was calculated as $|\mathit{d}|=5.603$ [61,62]. In addition, after 100,000 resamplings of the $|\mathit{d}|$ value, the bootstrap 95% confidence interval is $[5.549,5.659]$. This indicates a significant difference between E3 and E4, and agrees well with the probability density histogram shown in Figure 9(b). From ${\hat{\psi }}^{\mathrm{Av}}$, it can be quantitatively seen that E1 and E3 form one pair and E2 and E4 form another pair in the elongated domains along $⟨111{⟩}_{\mathrm{c}}$. These results confirm the successful separation of $\mathit{K}$-shape clustering with $\mathit{k}=4$. This Bayesian approach, although indirect, simplifies the assessment of significant differences between clusters. In addition, the stripe pattern along $⟨111{⟩}_{\mathrm{c}}$, which remained ambiguous in the $\mathit{K}$-means clustering, was separated by the $\mathit{K}$-shape clustering. This may be because the $\mathit{K}$-means clustering used the absolute value of the difference between the 30.0 and 14.1 K datasets to reduce the computational cost, possibly eliminating some information in $\mathit{\psi }$ [18].

Figure 9. (a) A box and whisker plot of the exploration values of the averaged fast-axis direction, expressed as ${\hat{\psi }}^{\mathrm{Av}}$. (b) Probability density histogram for ${\hat{\psi }}^{\mathrm{Av}}$ for E3 and E4. The scales of the angular measures in these figures are highly narrow; thus, they can be treated as linear statistics.

Table 2. Five-number summary (minimum, first quartile, median, third quartile, and maximum) and mean of exploration values of fast-axis direction for respective clusters expressed as ${\hat{\psi }}^{\mathrm{Av}}$.

Cluster	Mean (deg.)	Min (deg.)	1st (deg.)	Median (deg.)	3rd (deg.)	Max (deg.)
E1	89.59	89.49	89.58	89.59	89.61	89.70
E2	91.49	91.44	91.48	91.49	91.49	91.53
E3	90.40	90.33	90.39	90.40	90.41	90.46
E4	90.46	90.42	90.45	90.46	90.46	90.49

4.3. Hierarchical Bayesian temperature-series estimation

From previous studies of stress-induced ferroelectricity in SrTiO${ }_3$, it was found that both $T_{\mathrm{c}}$ and $T_{\mathrm{F}}$ increase with increasing stress [19–22]. As reported in Ref [16], the value of $\mathit{\delta }$ also increases due to increased stress and/or strain due to the photoelastic effect. The purpose of this study is to clarify the distribution of stress and/or strain from $T_{\mathrm{c}}$ and the distribution of the ferroelectric state from $T_{\mathrm{F}}$ on a pixel-by-pixel basis. The previous method typically estimates the phase transition temperature from the intersection of two straight lines on the $\mathit{\delta }(\mathit{T})$ curve using the least squares method [16–18]. However, this method is sensitive to outliers and cannot account for measurement errors due to its point estimation methodology. To ensure consistent quality of analysis across all 42,280 different $\mathit{T}$-series data, we used hierarchical Bayesian $\mathit{T}$-series estimation. This fits individual models at lower levels, guided by constraints from the higher-level model. In other words, as shown in Figure 10(a), the phase transition temperatures were estimated for each pixel while simultaneously estimating the overall phase transition temperature within the cluster.

Figure 10. (a) Representation of the hierarchical Bayesian $\mathit{T}$-series estimation for phase transition temperatures at each pixel ($T_{\mathrm{c}}^i$ and $T_{\mathrm{F}}^i$ for $\mathit{i}=1$–42,280) and for the entire cluster ($T_{\mathrm{c}}^{\mathrm{E}j}$ and $T_{\mathrm{F}}^{\mathrm{E}j}$ for $\mathit{j}=1$–4). Schematic representation of the model Hamiltonian used for the hierarchical Bayesian $\mathit{T}$-series estimation of (b) $T_{\mathrm{F}}$ and (c) $T_{\mathrm{c}}$ from the $\mathit{\delta }(\mathit{T})$ curve at position a3.

Because the estimation of $T_{\mathrm{F}}$ is a more intuitive treatment than that of $T_{\mathrm{c}}$, the hierarchical Bayesian $\mathit{T}$-series estimation method for $T_{\mathrm{F}}$ is explained first. For the $T_{\mathrm{F}}$ estimation, we used the intersection of two straight lines, as shown schematically in Figure 10(b). The notations $T_{\mathrm{F}}^i$ for the pixels ($\mathit{i}=1\sim 42,280$) and $T_{\mathrm{F}}^{\mathrm{E}j}$ for the clusters ($\mathit{j}=1\sim 4$) denote the ferroelectric transition temperature at the $\mathit{i}$-th pixel and as a whole in the E$\mathit{j}$ cluster, respectively. However, because the number $\mathit{i}$ is a sequential number unrelated to the cluster number $\mathit{j}$, it is necessary to count only the skipped number $\mathit{i}$ belonging to each cluster based on the $\mathit{K}$-shape clustering result shown in Figure 6. As shown in Figure 10(b), a positive slope is assumed for $\mathit{T}\ge T_{\mathrm{F}}^i$ (the paraelectric region), and a negative slope is assumed for $\mathit{T}<T_{\mathrm{F}}^i$ (the ferroelectric region). For this analysis, the $\mathit{\delta }(\mathit{T})$ data at 575 nm from 40.0 to 14.1 K were entered into the machine learning in their measured order. To reduce the computational cost, the $\mathit{\delta }(\mathit{T})$ curves were smoothed by averaging every two $\mathit{T}$ points for the $\mathit{\delta }$ values, resulting in 159 points with an interval of approximately 0.17 K between 40.0 and 14.1 K. In addition, the smoothed $\mathit{\delta }(\mathit{T})$ curve at the $\mathit{i}$-th pixel was standardized, denoted as $Y^i(\mathit{n})$, where $\mathit{n}$ is the $\mathit{n}$-th $\mathit{T}$ point. Standardization makes it easier to define a more appropriate prior probability distribution because all data can be treated on the same scale without affecting the estimation of $T_{\mathrm{F}}^{\mathrm{E}j}$ and $T_{\mathrm{F}}^i$. After extensive optimization, the prior probability distribution was defined as follows:

$T_{\mathrm{F}}^{\mathrm{E}j}\sim \mathit{f}(T_{\mathrm{F}}^{\mathrm{E}j};22.0,4.0)$

(20) $\equiv \frac{1}{\pi }\frac{4.0}{{(T_{\mathrm{F}}^{\mathrm{E}j}-22.0)}^2+{4.0}^2},$(20)

(21) $T_{\mathrm{F}}^i\sim \mathit{N}(T_{\mathrm{F}}^{\mathrm{E}j},{0.3}^2),$(21)

where $\mathit{f}(\mathit{T};x_0,\mathit{\gamma })$ is the Cauchy distribution varying $\mathit{T}$ with median $x_0$ and scale parameter $\mathit{\gamma }$. Compared with the normal distribution, the heavier tails of the Cauchy distribution are more robust to outliers, thereby providing a more flexible and tolerant approach to modeling errors in data with potential anomalies. The range limits for $T_{\mathrm{F}}^{\mathrm{E}j}$ and $T_{\mathrm{F}}^i$ are set to

(22) $16.0\mathrm{K}<T_{\mathrm{F}}^{\mathrm{E}j}<35.0\mathrm{K},$(22)

(23) $16.0\mathrm{K}<T_{\mathrm{F}}^i<35.0\mathrm{K}.$(23)

The probability distributions for $T_{\mathrm{F}}^{\mathrm{E}j}$ and $T_{\mathrm{F}}^i$ are broad, but $T_{\mathrm{F}}^i$ is closely related to $T_{\mathrm{F}}^{\mathrm{E}j}$ in Equation (21)(21) $T_{\mathrm{F}}^i\sim \mathit{N}(T_{\mathrm{F}}^{\mathrm{E}j},{0.3}^2),$(21) . Under these conditions, the model Hamiltonian is as follows: if $\mathit{T}\ge T_{\mathrm{F}}^i$ (the paraelectric state),

(24) $Y^i(\mathit{n})\sim \mathit{N}\left(\hat{a}+{\hat{b}}_1\times \left(\mathit{X}(\mathit{n})-T_{\mathrm{F}}^i\right),{{\mathit{\sigma }}_{\mathrm{obs}}}^2\right),$(24)

${\sigma }_{\mathrm{obs}}\sim \left|\mathit{f}({\sigma }_{\mathrm{obs}};0,2)\right|,$

while if $\mathit{T}<T_{\mathrm{F}}^i$ (the ferroelectric state),

(25) $Y^i(\mathit{n})\sim \mathit{N}\left(\hat{a}+{\hat{b}}_2\times \left(\mathit{X}(\mathit{n})-T_{\mathrm{F}}^i\right),{{\mathit{\sigma }}_{\mathrm{obs}}}^2\right),$(25)

${\sigma }_{\mathrm{obs}}\sim \left|\mathit{f}({\sigma }_{\mathrm{obs}};0,2)\right|,$

where $\mathit{X}(\mathit{n})$ is the temperature at the $\mathit{n}$-th $\mathit{T}$ point, i.e. $\mathit{X}(1)=40.0$ K, $\cdots$, $\mathit{X}(159)=14.1$ K. Unlike Equation (13)(13) $\mathit{\delta }=\left(1-\frac{\theta }{{360}^{\circ }}\right)\times \mathit{\lambda },$(13) , Equations (24)(24) $Y^i(\mathit{n})\sim \mathit{N}\left(\hat{a}+{\hat{b}}_1\times \left(\mathit{X}(\mathit{n})-T_{\mathrm{F}}^i\right),{{\mathit{\sigma }}_{\mathrm{obs}}}^2\right),$(24) and (25(25) $Y^i(\mathit{n})\sim \mathit{N}\left(\hat{a}+{\hat{b}}_2\times \left(\mathit{X}(\mathit{n})-T_{\mathrm{F}}^i\right),{{\mathit{\sigma }}_{\mathrm{obs}}}^2\right),$(25) ) use linear statistics to address the errors in $Y^i(\mathit{n})$ because we estimate the intersection of two straight lines on the linear graph.

Below 40.0 K, the displacement of the sample position is negligible because the stress apparatus underwent sufficient thermal shrinkage. Therefore, as shown in Figure 10(b), the measurement error corresponding to ${\sigma }_{\mathrm{obs}}$ could be small. The parameters ${\hat{b}}_1$ and ${\hat{b}}_2$ in Equations (24)(24) $Y^i(\mathit{n})\sim \mathit{N}\left(\hat{a}+{\hat{b}}_1\times \left(\mathit{X}(\mathit{n})-T_{\mathrm{F}}^i\right),{{\mathit{\sigma }}_{\mathrm{obs}}}^2\right),$(24) and (25(25) $Y^i(\mathit{n})\sim \mathit{N}\left(\hat{a}+{\hat{b}}_2\times \left(\mathit{X}(\mathit{n})-T_{\mathrm{F}}^i\right),{{\mathit{\sigma }}_{\mathrm{obs}}}^2\right),$(25) ) were constrained to ${\hat{b}}_1>0$ (a positive slope) and ${\hat{b}}_2<0$ (a negative slope), respectively, but the parameter $\hat{a}$ was not constrained. The MCMC simulations were performed with different random initial values for 20,000 iterations each over eight chains, with the first 1,000 being the warm-up period. As detailed in the Supplementary materials, the $n_{\mathrm{eff}}$ derived from the MCMC simulations exceeded $\sim$40% of the total sample size, which consisted of 152,000 iterations. For convergence diagnostics, the Gelman – Rubin statistics were computed for seven parameters, including the log probability of the unconstrained parameters associated with the estimation of $T_{\mathrm{F}}^{\mathrm{E}j}$ and $T_{\mathrm{F}}^i$, and $\hat{R}$ of all variables and pixels remained below 1.05, indicating satisfactory convergence (see Supplementary materials).

To estimate $T_{\mathrm{c}}$, the structural phase transition temperature at each pixel, $T_{\mathrm{c}}^i$ ($\mathit{i}=1\sim 42,280$), and in the E$\mathit{j}$ cluster, $T_{\mathrm{c}}^{\mathrm{E}j}$ ($\mathit{j}=1\sim 4$), were determined similar to $T_{\mathrm{F}}$. Figure 10(c) shows that the $\mathit{\delta }(\mathit{T})$ peak at $\sim$120 K has a large width and significant jump and bend, making it difficult to use two straight lines to identify the intersection. Despite the computational burden, we could detect a trend change from increasing to decreasing by comparing the values of $\mathit{\delta }$ at the current $\mathit{T}$ and a previous $\mathit{T}$. This approach could be called ‘Bayesian temperature-series change point detection’ as it is derived from the time-series analysis [33]. In this method, the order of the $\mathit{T}$-series data and the uniformity of the $\mathit{T}$ intervals may affect the estimation results; however, the data were entered in their measured order without any correction. In Figure 10(c), several types of slopes can be observed around the peak top. In this study, by setting a wide distribution of measurement error, we detected temperatures that entered a strongly decreasing trend (see Supplementary materials). The $\mathit{\delta }(\mathit{T})$ curve between 149.7 and 100.1 K for each pixel was smoothed by averaging every eight $\mathit{T}$ points, resulting in a series of 74 points and an interval of approximately 0.68 K. Although the $\mathit{T}$ interval of the $T_{\mathrm{c}}$ estimate is larger than that of the $T_{\mathrm{F}}$ estimate, this discrepancy is not a problem because the $\mathit{T}$ intervals are approximately the same when normalized by the respective phase transition temperatures based on the scaling laws [1]. Furthermore, the smoothed $\mathit{\delta }(\mathit{T})$ curve between 149.7 and 100.1 K at the $\mathit{i}$-th pixel was standardized and denoted as $Y^i(\mathit{n})$. To explore the parameter space more independently, we defined the structural phase transition temperatures as follows:

(26) $T_{\mathrm{c}}^i=T_{\mathrm{c}}^{\mathrm{E}j}+T_{\mathrm{c}}^{\mathrm{o}\mathrm{f}\mathrm{f}i},$(26)

where $T_{\mathrm{c}}^{\mathrm{o}\mathrm{f}\mathrm{f}i}$ is the deviation between $T_{\mathrm{c}}^i$ and $T_{\mathrm{c}}^{\mathrm{E}j}$. Post optimization, the prior probability for the MCMC simulations was defined as follows:

(27) $107.0\mathrm{K}<T_{\mathrm{c}}^{\mathrm{E}j}<130.0\mathrm{K},$(27)

$T_{\mathrm{c}}^{\mathrm{E}j}=\mathit{N}(117.0,{1.0}^2),$

(28) $-0.5\mathrm{K}<T_{\mathrm{c}}^{\mathrm{o}\mathrm{f}\mathrm{f}i}<+0.5\mathrm{K},$(28)

$T_{\mathrm{c}}^{\mathrm{o}\mathrm{f}\mathrm{f}i}=\mathit{N}(0,5^2).$

We set a wider range for $T_{\mathrm{c}}^{\mathrm{E}j}$ but a narrower range for $T_{\mathrm{c}}^{\mathrm{o}\mathrm{f}\mathrm{f}i}$. Preliminary runs indicated that exploring ${\hat{T}}_{\mathrm{c}}^{\mathrm{off}\mathit{i}}$ does not significantly affect the estimation of $T_{\mathrm{c}}^i$, as ${\hat{T}}_{\mathrm{c}}^{\mathrm{off}\mathit{i}}$ is uniformly distributed over the substrate (see Supplementary materials). To minimize the risk of converging to local minima, we systematically reduced the range limit while closely monitoring the posterior probability distribution. The final exploration range was then determined using Equation (28)(28) $-0.5\mathrm{K}<T_{\mathrm{c}}^{\mathrm{o}\mathrm{f}\mathrm{f}i}<+0.5\mathrm{K},$(28) , considering the balance between the computational time and accuracy. Under these conditions, the model Hamiltonian is as follows: if $\mathit{T}\ge T_{\mathrm{c}}^i$ (the cubic phase),

(29) $Y^i(\mathit{n})=\mathit{N}\left(Y^i\left(\mathit{n}-1\right)+{\hat{c}}_1,{{\mathit{\sigma }}_{\mathrm{obs}}}^2\right),$(29)

${\sigma }_{\mathrm{obs}}=\left|\mathit{f}({\sigma }_{\mathrm{obs}};0,5)\right|,$

while if $\mathit{T}<T_{\mathrm{c}}^i$ (the tetragonal phase),

(30) $Y^i(\mathit{n})=\mathit{N}\left(Y^i\left(\mathit{n}-1\right)+{\hat{c}}_2,{{\mathit{\sigma }}_{\mathrm{obs}}}^2\right),$(30)

${\sigma }_{\mathrm{obs}}=\left|\mathit{f}({\sigma }_{\mathrm{obs}};0,5)\right|,$

where $Y^i(\mathit{n})$ and $Y^i(\mathit{n}-1)$ are the standardized $\mathit{\delta }(\mathit{T})$ curves at the $\mathit{i}$-th pixel for the $\mathit{n}$-th and ($\mathit{n}-1$)-th $\mathit{T}$ points, respectively. The measurement error ${\sigma }_{\mathrm{obs}}$ was set to a distribution that was sufficiently wide to be regarded as an essentially uniform distribution. The parameters ${\hat{c}}_1$ and ${\hat{c}}_2$ in Equations (29)(29) $Y^i(\mathit{n})=\mathit{N}\left(Y^i\left(\mathit{n}-1\right)+{\hat{c}}_1,{{\mathit{\sigma }}_{\mathrm{obs}}}^2\right),$(29) and (30(30) $Y^i(\mathit{n})=\mathit{N}\left(Y^i\left(\mathit{n}-1\right)+{\hat{c}}_2,{{\mathit{\sigma }}_{\mathrm{obs}}}^2\right),$(30) ) are assumed as ${\hat{c}}_1>0$ (an increasing trend) and ${\hat{c}}_2<0$ (a decreasing trend), respectively. The MCMC simulations were performed with different random initial values for 20,000 iterations each over four chains, with the first 1,000 as the warm-up period. As detailed in the Supplementary materials, the $n_{\mathrm{eff}}$ derived from the MCMC simulations exceeded $\sim$20% of the total sample size, which consisted of 76,000 iterations. The convergence diagnostics used the Gelman – Rubin statistics for 81 parameters, including the log probability of the unconstrained parameters associated with the estimation of $T_{\mathrm{c}}^{\mathrm{E}j}$ and $T_{\mathrm{c}}^{\mathrm{o}\mathrm{f}\mathrm{f}i}$, and $\hat{R}$ of all variables and pixels were below 1.05, indicating satisfactory convergence (see Supplementary materials).

4.4. Analysis for each pixel

Figure 11 shows the results of the exploration values for $T_{\mathrm{c}}$, expressed as ${\hat{T}}_{\mathrm{c}}^{\mathrm{E}\mathit{j}}$ and ${\hat{T}}_{\mathrm{c}}^{\mathit{i}}$, to visually represent the posterior distributions for the MCMC simulations. Figure 11(a) shows a box and whisker plot of ${\hat{T}}_{\mathrm{c}}^{\mathrm{E}\mathit{j}}$ ($\mathit{j}=1\sim 4$). Table 3 presents a five-number summary and the mean required for constructing this plot. The exploration of the parameter space for $T_{\mathrm{c}}^{\mathrm{E}j}$ is set wide in Equation (27)(27) $107.0\mathrm{K}<T_{\mathrm{c}}^{\mathrm{E}j}<130.0\mathrm{K},$(27) ; therefore, the minimum and maximum values are widely scattered, but their ranges are sufficiently narrower than their limits. No significant differences in ${\hat{T}}_{\mathrm{c}}^{\mathrm{E}\mathit{i}}$ were observed between clusters. Figure 11(b) shows a probability density histogram for ${\hat{T}}_{\mathrm{c}}^{\mathit{i}}$, summed over the exploration values at all pixels in each cluster. Although this plot is unsuitable for understanding the characteristics of ${\hat{T}}_{\mathrm{c}}^{\mathit{i}}$ on a pixel-by-pixel basis, the distribution of ${\hat{T}}_{\mathrm{c}}^{\mathit{i}}$ showed no significant differences between clusters. We calculated the Pearson correlation coefficient ($\mathit{r}$) to examine the relationship between ${\hat{T}}_{\mathrm{c}}^{\mathrm{E}\mathit{j}}$ (the higher level) and ${\hat{T}}_{\mathrm{c}}^{\mathit{i}}$ (the lower level). Because the coefficient $\mathit{r}$ measures the relative relationship between two groups, changes in the scale or position of the variables do not affect the strength or direction of the correlation. In general, a linear relationship is considered strong when $|\mathit{r}|>0.7$, moderate when $0.3<|\mathit{r}|\le 0.7$, and weak or indicative of almost no correlation when $|\mathit{r}|\le 0.3$ [61]. For the correlation between ${\hat{T}}_{\mathrm{c}}^{\mathrm{E}\mathit{j}}$ and ${\hat{T}}_{\mathrm{c}}^{\mathit{i}}$, the coefficient $\mathit{r}$ exceeds 0.94 for all pixels, because the exploration of ${\hat{T}}_{\mathrm{c}}^{\mathrm{off}\mathit{i}}$ did not affect the estimation of ${\hat{T}}_{\mathrm{c}}^{\mathit{i}}$. This indicates a uniform distribution of ${\hat{T}}_{\mathrm{c}}^{\mathit{i}}$ over the substrate (see Supplementary materials). Figure 11(c) shows the distribution of the mean of ${\hat{T}}_{\mathrm{c}}^{\mathit{i}}$ for all pixels, which appears almost uniform. Although the standard deviation is often used when discussing the degree of distribution, it does not accurately represent the characteristics of the distribution if the dataset has an asymmetric distribution or numerous outliers. Furthermore, to compare the distribution of $T_{\mathrm{c}}$ and $T_{\mathrm{F}}$, we must use indicators that can be analyzed at different scales. Therefore, a quantitative assessment of the variation in the mean of ${\hat{T}}_{\mathrm{c}}^{\mathit{i}}$ is facilitated by Table 4, which presents the calculations of median absolute deviation (MAD) and coefficient of variation (CV) [65]. In addition, the bootstrap 95% confidence interval is also shown after 100,000 resamplings. The MAD metric, which indicates the median deviation from the median, provides a robust assessment of data scatter that is less sensitive to outliers. The CV metric, which is the standard deviation divided by the mean, reflects the relative variation. These metrics are useful for comparing datasets of different scales. From Table 4, the MAD shows that the variation is smaller for E1 and E2 than for E3 and E4, but the difference is small. Furthermore, the CV shows that the standard deviation is estimated to be approximately 0.1$\sim$0.2% of the mean; thus, the overall variation in the mean of ${\hat{T}}_{\mathrm{c}}^{\mathit{i}}$ is almost quantitatively uniform.

Figure 11. (a) A box and whisker plot of exploration values of structural phase transition temperatures in clusters, expressed as ${\hat{T}}_{\mathrm{c}}^{\mathrm{E}1}$, ${\hat{T}}_{\mathrm{c}}^{\mathrm{E}2}$, ${\hat{T}}_{\mathrm{c}}^{\mathrm{E}3}$, and ${\hat{T}}_{\mathrm{c}}^{\mathrm{E}4}$. (b) Probability density histogram of the sum of the exploration values of ${\hat{T}}_{\mathrm{c}}^{\mathit{i}}$ for each pixel when structural phase transition temperatures at a pixel are expressed as $T_{\mathrm{c}}^i$. All histograms overlap. (c) Distribution plot of the mean of ${\hat{T}}_{\mathrm{c}}^{\mathit{i}}$.

Table 3. Five-number summary and mean of exploration values of the structural phase transition temperatures for respective clusters expressed as ${\hat{T}}_{\mathrm{c}}^{\mathrm{E}1}$, ${\hat{T}}_{\mathrm{c}}^{\mathrm{E}2}$, ${\hat{T}}_{\mathrm{c}}^{\mathrm{E}3}$, and ${\hat{T}}_{\mathrm{c}}^{\mathrm{E}4}$.

Cluster	Mean (K)	Min (K)	1st (K)	Median (K)	3rd (K)	Max (K)
E1	116.94	111.25	116.21	116.92	117.65	122.81
E2	116.94	111.38	116.20	116.91	117.65	122.66
E3	116.89	111.30	116.13	116.88	117.63	123.17
E4	116.89	111.67	116.13	116.87	117.64	122.99

Table 4. Median absolute deviation (MAD) and coefficient of variation (CV) in the distribution of the mean of ${\hat{T}}_{\mathrm{c}}^{\mathit{i}}$ and ${\hat{T}}_{\mathrm{F}}^{\mathit{i}}$ for respective clusters. The upper panel of each column displays the values of MAD and CV, whereas the lower panel provides the bootstrap 95% confidence intervals for them.

	Mean of ${\hat{T}}_{\mathrm{c}}^{\mathit{i}}$		Mean of ${\hat{T}}_{\mathrm{F}}^{\mathit{i}}$
Cluster	MAD	CV	MAD	CV
Total	0.110	0.00196	0.561	0.0655
	$[0.108,0.111]$	$[0.00193,0.00198]$	$[0.552,0.570]$	$[0.0647,0.0663]$
E1	0.075	0.00122	1.181	0.0845
	$[0.072,0.077]$	$[0.00119,0.00125]$	$[1.125,1.241]$	$[0.0832,0.0859]$
E2	0.090	0.00139	1.013	0.0780
	$[0.087,0.093]$	$[0.00135,0.00143]$	$[0.961,1.055]$	$[0.0764,0.0796]$
E3	0.117	0.00206	0.470	0.0423
	$[0.114,0.120]$	$[0.00202,0.00209]$	$[0.458,0.482]$	$[0.0413,0.0432]$
E4	0.121	0.00217	0.406	0.0357
	$[0.119,0.124]$	$[0.00212,0.00221]$	$[0.397,0.416]$	$[0.0351,0.0364]$

Figure 12 shows the results of the exploration values for $T_{\mathrm{F}}$, represented as ${\hat{T}}_{\mathrm{F}}^{\mathrm{E}\mathit{j}}$ and ${\hat{T}}_{\mathrm{F}}^{\mathit{i}}$. Figure 12(a) shows a box and whisker plot of ${\hat{T}}_{\mathrm{F}}^{\mathrm{E}\mathit{j}}$ ($\mathit{j}=1\sim 4$). Table 5 provides a five-number summary and the mean required for constructing this plot. For E1 and E2, ${\hat{T}}_{\mathrm{F}}^{\mathrm{E}1}$ and ${\hat{T}}_{\mathrm{F}}^{\mathrm{E}2}$ are scattered on the high-$\mathit{T}$ side; thus, the mean is more than 0.5 K above the median and the IQR is more than twice as large as for E3 and E4. Figure 12(b) shows a probability density histogram for ${\hat{T}}_{\mathrm{F}}^{\mathit{i}}$, summed over the exploration values at all pixels in each cluster. In E1 and E2, the histograms show 3–4 peaks, indicating a mixture of several groups. Conversely, for E3 and E4, the histograms are concentrated on the lower-$\mathit{T}$ side. For the correlation between ${\hat{T}}_{\mathrm{F}}^{\mathrm{E}\mathit{j}}$ and ${\hat{T}}_{\mathrm{F}}^{\mathit{i}}$ in each cluster, the coefficients $\mathit{r}$ are widely scattered within $0.1<\mathit{r}<0.6$, although the means of $\mathit{r}$ in the respective clusters are 0.25–0.31 (see Supplementary materials). This indicates that the lower-level estimation is successfully performed for each pixel, while being restricted by the higher-level estimation of the entire cluster. Figure 12(c) shows the distribution of the mean of ${\hat{T}}_{\mathrm{F}}^{\mathit{i}}$. The MAD and CV, and their bootstrap 95% confidence intervals after 100,000 resamplings were calculated to evaluate the distributions of the mean of ${\hat{T}}_{\mathrm{F}}^{\mathit{i}}$ and are shown in Table 4. The MAD results show that the distributions of the mean of ${\hat{T}}_{\mathrm{F}}^{\mathit{i}}$ from the median are more than twice as wide for E1 and E2 than for E3 and E4. The CV results show that the standard deviation of the mean of ${\hat{T}}_{\mathrm{F}}^{\mathit{i}}$ is approximately 8% of the mean for E1 and E2, and approximately 4% of the mean for E3 and E4. For both metrics, E1 and E2 are a pair and E3 and E4 are another pair. Table 6 shows Cohen’s $\mathit{d}$ with Hedges’ correction of distributions of ${\hat{T}}_{\mathrm{F}}^{\mathit{i}}$ and the bootstrap 95% confidence interval after 10,000 resamplings of $|\mathit{d}|$ for each pair of clusters to assess significant differences. Consequently, significant differences exist between the clusters, except for the pair E1 and E2 and pair E3 and E4. This is consistent with the ${\delta }^{\mathrm{A}v}(\mathit{T})$ curves shown in Figure 7(a).

Figure 12. (a) A box and whisker plot of exploration values of ferroelectric phase transition temperatures in clusters, expressed as ${\hat{T}}_{\mathrm{F}}^{\mathrm{E}1}$, ${\hat{T}}_{\mathrm{F}}^{\mathrm{E}2}$, ${\hat{T}}_{\mathrm{F}}^{\mathrm{E}3}$, and ${\hat{T}}_{\mathrm{F}}^{\mathrm{E}4}$. (b) Probability density histogram of the sum of the exploration values of ${\hat{T}}_{\mathrm{F}}^{\mathit{i}}$ for each pixel when structural phase transition temperatures at a pixel are expressed as $T_{\mathrm{F}}^i$. The histograms for E1 and E2 overlap, and those for E3 and E4 overlap. (c) Distribution plot of the mean of ${\hat{T}}_{\mathrm{F}}^{\mathit{i}}$.

Table 5. Five-number summary and mean of exploration values of the ferroelectric phase transition temperatures for respective clusters expressed as ${\hat{T}}_{\mathrm{F}}^{\mathrm{E}1}$, ${\hat{T}}_{\mathrm{F}}^{\mathrm{E}2}$, ${\hat{T}}_{\mathrm{F}}^{\mathrm{E}3}$, and ${\hat{T}}_{\mathrm{F}}^{\mathrm{E}4}$.

Cluster	Mean (K)	Min (K)	1st (K)	Median (K)	3rd (K)	Max (K)
E1	23.36	18.99	21.58	22.75	24.96	30.84
E2	23.44	19.00	21.87	22.90	24.75	31.31
E3	21.56	18.42	20.77	21.43	22.12	31.02
E4	21.56	18.31	20.88	21.49	22.11	30.14

Table 6. Cohen’s $\mathit{d}$ with Hedges’ correction of the distributions of the exploration values ${\hat{T}}_{\mathrm{F}}^{\mathit{i}}$ in two combinations within E1, E2, E3, and E4. Because the sign of the $\mathit{d}$ value is insignificant, the upper panel of each column displays the absolute value of $\mathit{d}$, whereas the lower panel details the bootstrap 95% confidence intervals for $|\mathit{d}|$.

	E1	E2	E3
E2	0.0412
	$[0.0402,0.0422]$
E3	1.4137	1.5656
	$[1.4127,1.4147]$	$[1.5645,1.5666]$
E4	1.5051	1.6810	0.0025
	$[1.5041,1.5061]$	$[1.6800,1.6820]$	$[0.0020,0.0031]$

From Ref [18]. the larger $T_{\mathrm{F}}$ regions indicate the stress and/or strain concentration. On the basis of linear relationship between stress and $T_{\mathrm{c}}$ expressed as Equation (1)(1) ${\sigma }_{\mathrm{Cv}}=30\times \left(T_{\mathrm{c}}-105\right)\left(\mathrm{MPa}\right).$(1) , the stress distribution was also considered. By uniformly subtracting 105 K from the means of ${\hat{T}}_{\mathrm{c}}^{\mathit{i}}$ in clusters, the standard deviations remain unchanged, while the mean and median are changed. Therefore, the MADs remain constant, but the CVs change, and the relative relationship between the CVs could change. However, if the means of ${\hat{T}}_{\mathrm{c}}^{\mathit{i}}$ in clusters are multiplied by the same factor, the relationship between the MADs does not change as the same coefficients are multiplied for them. Comparing the MADs, the relationship of the means of ${\hat{T}}_{\mathrm{c}}^{\mathit{i}}$ between clusters is opposite to that of the means of ${\hat{T}}_{\mathrm{F}}^{\mathit{i}}$. In Ref [16], we have already proposed evaluating the stress and/or strain distribution using $\mathit{\delta }$ caused by the photoelastic effect. The distributions of $\mathit{\delta }(\mathit{T})$ and $\mathit{\psi }(\mathit{T})$ at 130.9 K in Figure 1 and those at 14.1 K in Figure 1 of Ref [16]. are similar to the distribution of the mean of ${\hat{T}}_{\mathrm{F}}^{\mathit{i}}$ in Figure 12(c). To evaluate the similarity quantitatively, Figure 13 shows the heat maps of the coefficient $\mathit{r}$ between $\mathit{\delta }(\mathit{T})$, $\mathit{\psi }(\mathit{T})$, means of ${\hat{T}}_{\mathrm{c}}^{\mathit{i}}$, and means of ${\hat{T}}_{\mathrm{F}}^{\mathit{i}}$. Concerning the correlations between the means of ${\hat{T}}_{\mathrm{c}}^{\mathit{i}}$ and ${\hat{T}}_{\mathrm{F}}^{\mathit{i}}$, there are moderate correlations for E1 but weak or no correlations for E2–E4. Furthermore, for the relationship between $\mathit{\delta }(\mathit{T})$ and the means of ${\hat{T}}_{\mathrm{F}}^{\mathit{i}}$, there are moderate correlations for E1 and E2, but almost no correlations for E3 and E4 at 130.9 K, whereas there are strong correlations for E1 and E2, and moderate correlations for E3 and E4 at 14.1 K. Contrary to our expectation, we could not find any pair of variables showing a strong correlation for E3 and E4, despite almost uniform stress and/or strain. Consequently, the stress and/or strain distribution could be better estimated from $\mathit{\delta }(14.1\mathrm{K})$ or the mean of ${\hat{T}}_{\mathrm{F}}^{\mathit{i}}$ than from the mean of ${\hat{T}}_{\mathrm{c}}^{\mathit{i}}$.

Figure 13. Heat maps for each cluster of Pearson correlation coefficients ($\mathit{r}$) between $\mathit{\delta }$, $\mathit{\psi }$, mean of ${\hat{T}}_{\mathrm{c}}^{\mathit{i}}$ and mean of ${\hat{T}}_{\mathrm{F}}^{\mathit{i}}$ at (a) 130.9 K and (b) 14.1 K.

Finally, we speculate why the correlation is stronger for E1 and E2 than for E3 and E4. If the Pockels effect causes an increase in the $\mathit{\delta }(\mathit{T})$ curve below $T_{\mathrm{F}}$, the coefficient $\mathit{r}$ should be stronger for E3 and E4 than for E1 and E2 because of the less disordered crystal orientation. Meanwhile, considering the photoelastic effect, we can assume a strong linear relationship between the local stress and/or strain and $T_{\mathrm{F}}$. Although it is not easy to distinguish between these effects experimentally, we hypothesize that flexoferroelectricity is induced by local strain, which is an interesting concept [66], but it borders on speculative territory in this study. At zero electric field, the ensemble average of the normal spontaneous polarization within a pixel ($\ge$ 10 $\times$ 10 $\mathit{\mu }$m${ }^2$) is very small, whereas that of the flexoferroelectric polarization may be large because the strain has already biased it in a certain direction. If this scenario is true, flexoferroelectric polarization would be generated locally for E1 and E2 in addition to the normal ferroelectric polarization distributed throughout the substrate. This hypothesis can be tested by applying electric fields at various $\mathit{\sigma }$. For example, the flexoferroelectricity is bound to the lattice distortions; therefore, different electric field dependencies may emerge within the substrate. In the future, to confirm the nonlinear relationship between ${\hat{T}}_{\mathrm{c}}^{\mathit{i}}$, ${\hat{T}}_{\mathrm{F}}^{\mathit{i}}$, and stress and/or strain, we will also apply machine learning using nonlinear models such as decision trees and random forests to unravel the complex dynamics governing the flexoferroelectricity and normal ferroelectricity.

5. Conclusion

In this study, the stress-induced ferroelectric domains elongated along $⟨111{⟩}_{\mathrm{c}}$ in SrTiO${ }_3$ were investigated by combining machine learning and statistics to extract deep insights from large $\mathit{T}$-series birefringence imaging data. After applying the advanced preprocessing to characterize the spatial information of the $\mathit{T}$-series data, we first performed the multivariate $\mathit{K}$-shape clustering method focused on the shape similarity of the $\mathit{T}$ dependencies. This result showed that a large spontaneous polarization appears for E1 and E2 owing to stress and/or strain concentration. Furthermore, hierarchical Bayesian $\mathit{T}$-series estimation was used to ensure consistent analysis of 42,280-pixel $\mathit{T}$-series data to estimate $T_{\mathrm{c}}^i$ and $T_{\mathrm{F}}^i$ for each pixel as constrained by the overall $T_{\mathrm{c}}^{\mathrm{E}j}$ and $T_{\mathrm{F}}^{\mathrm{E}j}$ in the cluster. For the means of ${\hat{T}}_{\mathrm{c}}^{\mathit{i}}$, the uniform distribution over the substrate is strongly supported. Furthermore, the means of ${\hat{T}}_{\mathrm{F}}^{\mathit{i}}$ for E1 and E2 are widely scattered to the high-$\mathit{T}$ side. The higher region coincides well with the stress and/or stress concentration regions. Consequently, the hierarchical Bayesian $\mathit{T}$-series estimation strongly supports the $\mathit{K}$-shape clustering results. Based on the correlation coefficients for E1 and E2, a strong correlation exists between $\mathit{\delta }(14.1\mathrm{K})$ and the means of ${\hat{T}}_{\mathrm{F}}^{\mathit{i}}$, but the correlations become moderate for E3 and E4. Because of the weak or no correlation between the means of ${\hat{T}}_{\mathrm{c}}^{\mathit{i}}$ and ${\hat{T}}_{\mathrm{F}}^{\mathit{i}}$, it is difficult to deduce the distribution of stress and/or strain from these means. Our multistep $\mathit{T}$-series analysis method has proven considerably useful for studying critical phenomena near phase transition temperatures. Further development of $\mathit{T}$-series analysis is essential for the widespread use of machine learning in condensed matter physics.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/27660400.2024.2342234

Additional information

Funding

This research was partially supported by a Grant-in-Aid for Scientific Research (KAKENHI) Grant Numbers [JP21K04897] and [JP23K03283] from the Japan Society for the Promotion of Science.