these phenomena manifest in doped materials, where local atomic disorder can significantly affect competition between the order parameters, as exemplified, e.g., by ferroelectric relaxors or cation-substituted manganites, even though in certain cases, the phase coexistence can be induced by an external stimulus such as magnetic or electric fields, or uniaxial or hydrostatic pressure or strain. 28,29,75,76 These materials offer a challenge in terms of developing models capable of describing and predicting their behavior. For homogeneous materials, the natural description on the mesoscale level can be explored within Landau theory, where the properties of interest are associated with the appropriate order parameter and the energy of the system is represented as its functional. For ferroelectrics, ferroelastics, and multiferroics, the relevant description can be given by Ginzburg-Landau theory with a real order parameter, [30][31][32] whereas complex order parameter fields can be used for description of superconductors. [33][34][35][36] The mesoscopic order parameter-based approach naturally allows for classification of phase transitions, identification of topological defects, etc. 37 However, the applicability of these mesoscopic models for atomically disordered systems remains a topic of continuous controversy and has been for several decades. [38][39][40] Alternatively, the physics of these materials can be explored using lattice models, with the materials structure or functionality represented via a collection of local spins interacting via local interactions, 41,42 potentially mediated by long-range depolarization fields. These models can be explored numerically, or in special cases analytically, to yield corresponding phase diagrams, susceptibilities, response functions, and other macroscopic functionalities. 41,42 On a local level, these models can give rise both to ordered and disordered ground states, offering considerably deeper insight into the physics of the material. 12 Yet, the key issue both for the mesoscale and lattice-type models is comparison with an experiment, as a way to explore the applicability of the model, extract relevant materials parameters, and enable forward prediction and thus materials optimization. For spatially uniform ground states, the combination of diffraction techniques with macroscopic measurements offers comprehensive insight into materials functionality, since the mesoscopic order parameter fields provide a sufficient description of the system. This, however, is not the case for disordered systems or systems with competing and spatially inhomogeneous ground states, since macroscopically averaged scattering data necessarily lose the relevant information. 43,44 Until recently, the analysis of such systems was possible only via complex diffraction techniques, where the analysis of the correlated disorder allowed reconstruction of possible microscopic models. [44][45][46] However, the associated inverse problems generally lead to multiple possible solutions, necessitating additional structural and compositional information. Recently, the emergence of high-resolution aberration-corrected electron microscopy has provided an opportunity to visualize these systems at the atomic level and detect minute distortions of local atomic structures that can serve as a proxy for local order parameter fields. Namely, site variations of symmetry can be extracted from the atomic coordination neighborhood to identify changes in local symmetry class, i.e., different phases, 47,48 or quantify the strength of symmetry-lowering distortions which provide measures for several order parameter fields. For instance, affine-type transforms of the local atomic lattice (i.e., expansions, compressions, and shears) serve as proxies for strain fields. [49][50][51][52][53] Additional examples include scalar lattice measurements, such as oxygen octahedral tilt 54 or tetragonality, 47,49,51,54,55 and vector quantities, such as inversionsymmetry breaking from electrical polarization order parameter fields in displacive ferroelectrics. 49,[55][56][57] While these techniques have been primarily used for qualitative studies, the information can be used to infer quantitative physics of materials via mesoscopic 58,59 or discrete models. 60,61 However, this approach further opens a question as to the role of non-observed degrees of freedom or unknown physics of interactions between the observed units, which can act as confounders or latent variables for observed behaviors. As an example, observation of the larger than expected Pt-Pt distance by Sohlberg et al. 62 was interpreted as a presence of a capping OH group in a Pt trimer, providing insight into its catalytic activity. In this case, the discrepancy between theoretical prediction and quantitative experimental measurements was used to infer the presence of (an unobservable) structural element. A more complex case is the transition from the electronic carrier screening to oxygen vacancy screening at ferroelectric interfaces. 63 Another example is the observation of a small region of atomic scale distortion (similar to nanoislands) in a cleaved manganite crystal in an otherwise featureless scanning tunneling microscopy image, which was interpreted as a real-space image of a trapped polaron. 64 However, these and many other examples rely on theoretical models that can fully explain the observed phenomena, and often require lengthy experiment and hypothesis testing cycles proposing and ruling out possible explanations for observed behavior. Often, this approach requires development of novel theoretical models or introduces difficult to access approximations. In some cases, the direct information on the latent and confounding factors (beyond their existence) is unavailable, leading to the considerable uncertainty in interpretation of experimental results, as exemplified by multiple reports of ferroelectricity in ultrathin oxides that can be attributed both to ferroelectricity 65,66 and electrochemical effects, 67 or the interplay between the two. [68][69][70] Generally, the discrepancy between the theory and observations (Bayesian surprise 71 ) is used to guide scientific research via proposing the existence of new entities or new physical mechanisms.

Here, we explore whether predictability and uncertainty of physical behaviors can be used to gain insight into associated physical mechanisms in the presence of potential confounding, bias, and latent factors, via machine learning methods. We argue that experimentally observed correlations active over the large volumes of experimental data can be used to make a general prediction, whereas the regions where these correlative relationships are violated suggest the presence of new physics, either in the form of new laws or presence of latent variables as illustrated in Fig. 1. For example, for an ideal ferroelectric material, the local polarization can be expected to be strongly affected by the polarization in neighboring unit cells, but only weakly sensitive to the chemical composition. At the same time, in the vicinity of the phase boundaries (in composition space) or domain boundaries (in real space), the effects of the compositional fluctuations will be more pronounced, resulting in larger uncertainties in prediction of relevant properties. At the same time, the functional form of these relationships is generally unknown. Thus, an ability to correlate the descriptors to the observations with a universal function approximator that includes uncertainty quantification will enable determination of when such descriptors are adequate (and, therefore, that our inclinations are supported) and when they appear inadequate to the task of functionality prediction. The latter can then lead directly to a search for new physical descriptors and further our understanding of these systems. This approach complements the classical physical paradigm of comparison to the known functional laws, or numerically derived models. Here, we use the Gaussian processes as universal interpolators that allow uncertainty quantification in predictions via the posterior predictive density and explore the relationship between the predictability and a priori physical mechanisms. As a model system, we choose paradigmatic multiferroic material BiFeO 3 (BFO) with Sm substitutional doping of Bi grown as a continuous composition spread library. Pure BFO is a rhombohedral ferroelectric (R3c) with the polarization oriented along the h111i axis of a prototypical cubic unit cell. Upon Sm substitution of $14%, the material transforms into an orthorhombic non-ferroelectric phase. 72,73 The two are separated by the morphotropic phase boundary between symmetryincompatible phases. 74,75 The material structure and behavior in the vicinity of the morphotropic phase boundary (MPB) are highly complex, with many functional properties showing pronounced maxima and microstructure exhibiting the presence of the complex multiscale domain patterns, nanoscale regions, etc. 76,77 While macroscopically many of these phases identify as monoclinic or triclinic, in the proximity to the MPB, the corresponding atomistic structures remain unknown. Notably, the structure and functionality of materials at an MPB remain one of the most challenging questions in the physics of ferroic systems due to the multiple interactions between phases, domain structures, and disorder.

A continuous composition spread of Bi 1-x Sm x FeO 3 in the range of 0 x 0.2 was fabricated on a SrTiO 3 (001) substrate by pulsed laser deposition. A uniform gradient was produced by alternate ablation of BiFeO 3 and SmFeO 3 targets with a moving shadow mask. A chemical composition range, from 0 to 20% Sm, was determined by wavelength dispersive spectroscopy (WDS) measurements. This composition range has been structurally characterized by x-ray diffraction and piezoresponse force microscopy (PFM), 78,79 and it displays a transition from highly ordered striped ferroelectric domains, to a ferroelectric mosaic, to a non-ferroelectric piezoresponse with increasing Sm%. 80 Cross-sectional atomic resolution STEM imaging was performed for nominal x ¼ 0%, 7%, and 20% Sm concentrations by focused ion beam (FIB) lift out from a single composition-spread film. Data were collected along the [001] pseudocubic zone axis at 200 kV on a Nion UltraSTEM with a high-angle annular dark-field (HAADF) detector, providing visualization of atomic columns according to their massthickness. To preserve local spatial relationships against raster-scan artifacts, the HAADF datasets were reconstructed from an orthogonal pair of HAADF images. 81 An example of HAADF dataset for the 0% Sm BiFeO 3 endmember is shown in Figs. 2(a) and 2(b). We quantify the local spatial relationships based on atom positions, determined by fit as 2D Gaussians, within a local neighborhood of the 5-cation perovskite unit cell as outlined in Ref. 82. Pertinent to the subsequent analysis, we define descriptors for each B-site (Fe) centered unit cell which contain structural {V, h, a, b}, chemical {I 1 , I 5 }, and polarization {P} informational content, as defined in Table I. They are, respectively, the unit cell volume (V), internal angle (h), in-plane lattice vector (a), outof-plane lattice vector (b), mean HAADF intensity (I 1 ), A-to B-site intensity difference (I 5 ), and the vector of Fe-cation displacement from the A-site defined centro-symmetric positions (P), which we hereafter refer to as the "polar displacement vector." Diagrams of the V, h, a, b, and P descriptors, derived from atom positions, are shown for an example unit cell in Fig. 2(b). Native units for these descriptors are pixels (datasets have a uniform scan resolution), with a real space conversion of 0.156 A ˚/pixel. I 1 , and I 5 derive from the HAADF atomiccolumn intensities measured as the Gaussian weighted 9-pixel intensity centered at the atom fit position. These two selected HAADF intensity descriptors correspond to the average (I 1 ) and A-to B-site difference (I 5 ) of cation mass-thickness and are indicators of sample thickness and cation composition variation, respectively.

To explore the predictability of specific physical behaviors, we adopt the Gaussian process (GP) method. [83][84][85][86] GP regression learns a function f over all source-target pairs D ¼ {(x 1 , y 1 ), …(x N , y N )}, with each pair related by y ¼ f(x) þ e, where e is Gaussian observation noise, by performing Bayesian inference in a function space, assuming (a) Gaussian process (GP) regression can, like a physical law with known functional form, be used to establish a relationship between two sets of parameters: in this case, between the atomic coordinate within a lattice and the local polarization. In the classical case, the relationship between the two is established based on known physical model, e.g., polarization profile across the domain wall. In GP case, the relationship is established based on multiple observations of the system providing interpolation function.

Note that an important aspect of GP approach is that prediction also yields uncertainty, allowing us to establish the significance of prediction and providing an independent channel of information on the process. Deviation from this behavior in certain locations can indicate the presence of hidden variables. For example, the presence of a dopant atom at a ferroelectric wall [(b) and (c)] will result in larger uncertainty in prediction of polarization at the certain lattice sites from that of the neighbors.

that function f has a prior distribution f $ GPð0; K f x; x 0 ð ÞÞ, where K f is a covariance function (kernel). 83 Naturally, GP is used as a powerful interpolator, where the covariance matrix of the GP posterior distribution serves as an estimator of the uncertainty in the interpolation. It should be noted that GP is a "lazy" machine learning algorithm that is evaluated at runtime and does not require specific "training." The optimization of the GP entails maximizing the log marginal likelihood, by means of tuning the hyperparameters h of the kernel function K, i.e.,

In the above equation, we note that the first term is essentially a penalty term for failure to correctly predict function values, whereas the second term penalizes the model complexity. Here, we used radial basis function

, where the variance r and the length scale l are kernel hyperparameters, which are learned from the data by maximizing the log-marginal likelihood. For large datasets (>10 3 ), computing the log marginal likelihood becomes nearly intractable, and instead, a small set of m function points is used as support or inducing variables, which are inferred along with kernel hyperparameters. 87 Here, we adapted the inducing points-based sparse implementation of GP from the Pyro probabilistic programing library, 88 which allows training a GP model with modern GPU accelerators, reducing the analysis time for datasets with >10 4 points from hours to minutes. Once the model is trained, we can calculate the mean prediction and the associated variance for each point (see the accompanying notebook at https://git.io/Jv1Bn for the full code of GP regression). Here, the high variance reflects higher uncertainty in GP prediction. To gain insight into the predictability of the physical parameters in these materials, we first explore the polarization behavior in the system. Here, we assume that the material possesses an electrical polarization field with a proxy representation by the polar displacement vector descriptor P, defined as the off-centering of the B-site cation. The following treatment is on the two components of P, noted as P x (i,j) and P y (i,j), where i and j are the lattice positions of the unit cell. We explore whether the polarization at the individual lattice site can be predicted from the polarization in the surrounding area, i.e., regressions

between an 8-component vector of polarization values in the 8 adjacent lattice sites and polarization in the given site. We further explore the predictability of the polarization component from the other polarization component, i.e., regression

The spatial degeneracy of the ferroelectric polarization frequently results in formation of polydomain configurations to minimize stray field energies and is a product of growth kinetics. Such a polydomain structure is observed for the 0% Sm ferroelectric BiFeO 3 endmember shown in Figs. 2(c)-2(e). This image contains three domains indicated by directional arrow labels in Fig. 2(c). The sharp boundaries, i.e., domain walls, are classified according to their polarization rotation across the wall. The domain walls present in the image are 109 and 180 , as labeled in Fig. 2(c), identified by their 2D projected rotations (where the z-axis rotation must be inferred) and conforming to known preferential (100) and (110) planes, respectively. In both instances, the domain walls deviate from their charge-neutral planes to ones with divergent polarization orientations, labeled with "þ" symbols representing their bound electrical charges, classified as charged domain walls.

Not surprisingly, the predicted polarization, denoted as Pyy , from the P y polarization neighborhood shows values and distribution very close to the original values [Fig. 2(f)], with the visible effect of GP being reduction and smoothing of the noise, i.e., essentially a convolution effect. However, non-trivial behavior emerges in the uncertainty maps [Fig. 2(g)], where significant uncertainty of prediction is observed in the vicinity of domain walls and in a number of distributed locations throughout the material. Notably, the prediction uncertainty is zero inside the substrate, i.e., local laws are well-defined.

The predicted polarization Pxy from the P y polarization neighborhood as compared to the true polarization distribution [Fig. 2(g)] shows considerably more interesting behavior. Note that in most locations of the image, the predicted polarization is close to the true polarization, as can be expected since GP is local interpolator. The prediction is reliable in the vicinity of the 109 domain wall. However, in the vicinity of the charged domain wall, the polarization shows strong deviations from the true polarization distribution, with clearly visible regions of opposite polarization direction. Careful examination of all polarization components shows that there may be fluctuations in the true polarization field; however, these are relatively minor. Therefore, the predicted "anomalous" polarization in the vicinity of the charged domain wall may be an indicator of confounding physical factors.

Interestingly, the uncertainty map in Fig. 2(i) shows overall higher uncertainty in prediction compared to GP ny!y , with the larger uncertainty values in the vicinity of domain walls and substrate interface. However, the prediction uncertainty has spatial distribution generally different from the difference between predicted and observed polarization.

Similar analyses for the other two compositions corresponding to the intermediate phase (x ¼ 0.07) and the orthorhombic phase (x ¼ 0.2) are shown in Fig. 3. Here, the intermediate phase [Figs. 3(a) and 3(b)] shows the presence of a ferroelectric phase with non-zero polarization values in the vicinity of the interface and a disordered phase in the film. The nature of this disordered phase is unclear and can correspond to a nanodomain system averaged in the z-axis direction of the electron beam, or a true disordered intermediate phase. Similar to the rhombohedral x ¼ 0 sample, the Pyy image shows features similar to the true image P y [Fig. 3(a)], albeit with considerably lower noise level and a slight blur. The corresponding uncertainty map shows regions with high uncertainty at the boundaries between the ferroelectric and disordered phase. Interestingly, these features seem to be new, and while they partially overlap with the domain walls, there is not a one-to-one correspondence. Again, the uncertainty for prediction within the substrate is uniform and small. The true P x image, predicted Pxy , and corresponding uncertainty are shown in Fig. 3(b). Notably, the uncertainty in GP is similar for Pyy and Pxy .

Finally, similar analysis is performed for the x ¼ 0.2 orthorhombic phase [Figs. 3(c) and 3(d)]. In this case, the polarization TABLE I. Unit cell descriptors. Calculations correspond to A-and B-site 2D Gaussian atom fits per Fe-centered HAADF unit cell according to site labels in Fig. 2(b).

In-plane lattice vector A

Out-of-plane lattice vector

Internal angle Angle between a and b V Internal volume Volume of convex hull A 1;xy , A 2;xy , A 3;xy , and A 4;xy

Cation intensity asymmetry A-site vs. B-site

xy distribution shows only weak contrast within the image and demonstrates small inclusions with period doubling (discussed later). For this composition, both Pyy and Pxy are very close to the true distributions P y and P x , and corresponding uncertainty maps show increase in uncertainty in the vicinity of the modulated phase regions.

To summarize, for all compositions explored, the predicted and uncertainty maps demonstrate additional spatial features. Some of these are concentrated at the regions where polarization changes rapidly, including domain walls and the interface. Yet, others are concentrated in the "cloud-like" regions within the material and do not bear obvious resemblance with the individual atomic descriptor maps. Hence, we conclude that the prediction and uncertainty maps contain non-trivial information or at least highlight the discrepancies where locally observed polarization differs from that expected given the average polarization behavior in the sample, pointing to new physical phenomena or the presence of latent variables.

We further proceed to explore the predictability and prediction uncertainty between the structural descriptors and polarization. For this, we perform the regressions

where S ij is the parameter vector of local structural descriptors.

In principle, such analysis can be performed using all the permutations of the available structural descriptors and their gradients (e.g., bag of features 89 approach), etc., to obtain comprehensive insight into functional relationships between the descriptors. Given the computational constraints and being informed by the general physics of ferroic materials, here we explore five groups of descriptors, namely, V; h; a; b; I 1 ; I 5 f g , a; b f g, V; h f g, V; h; a; b f g , and I 1 ; I 5 f g. This choice of groups is made given that I 1 and I 5 are parameters describing local chemical composition, V is related to the local phase state and composition (via Vegard terms), h is the measure of the unit cell deformation from cubic/tetragonal into modulate phase, and a and b define tetragonality of the unit cell. These variables are not independent; however, they can be considered as semiquantitative descriptors for composition, ferroelectricity, and proximity to the orthorhombic phase, respectively. Notably, all these descriptors are invariant with respect to the C 2 rotations, whereas polarization components change sign. Hence, the regression Eq. ( 3) explores whether polarization can be predicted from this set of structural and chemical descriptors and provide the associated uncertainties.

The analysis for rhombohedral BiFeO 3 is shown in Fig. 4. Here, the top row visualizes the true polarization distribution P y [Fig. 3(a)] and the predictions Pys for the different groups of structural descriptors [Figs. 3(b)-3(f)]. Upon examining the characteristic features in the images, we note that none of the groups allows for reliable reconstruction of polarization field, as can be expected from the physics of the problem. At the same time, the groups containing V; h f g as descriptors allow us to predict the changes in polarization in the charged domain and separate the ferroelectric phase and the substrate. The prediction based on a; b f g is completely non-discriminative, and only the anomalies at the domain walls are visible. Finally, chemical descriptors I 1 ; I 5 f g yield the maps resembling chemical disorder in the system, but do not contain any elements of the true polarization distribution or phase contrast between ferroelectric and substrate.

In comparison, the uncertainty maps shown in the bottom row of Fig. 4 show considerably richer behavior. Here, the prediction uncertainties for descriptor groups containing a; b f g show clear anomaly at the 109 wall, which is invisible for V; h f g and I 1 ; I 5 f g descriptors. At the same time, chemical descriptors I 1 ; I 5 f ggive rise to an uncertainty maximum at the ferroelectric-substrate interface. While this feature can be noticed in three other maps, it is much weaker. The charged domain wall is visible only in the maps containing V; h f g. Finally, I 1 ; I 5 f gseem to provide regions of higher uncertainty uncorrelated to the true potential distributions and potentially related to the presence of chemical disorder.

A similar analysis for the intermediate composition, x ¼ 0.07, is shown in Fig. 5. Similar to the ferroelectric phase, in general, the predicted polarization field is very dissimilar from the true polarization, with only the maps containing V; However, very interesting behavior is associated with the V; h f g group of variables that clearly shows the emergence of the modulated phase throughout the top two-thirds of the image. This phase is essentially unnoticeable in the original polarization image. What is particularly noteworthy is that this predicted ordering differs from the original one in the dataset.

The Fourier transform of each predicted P y image can be seen in Figs. 7(a)-7(h), where each FFT was calculated with a Hanning window to account for edge effects. To account for differences in signal strength, each FFT is normalized to the area under the cross section of the central spot. We can see that the maximum FFT intensity is observed in the original image [Fig. 7(a)] with some superlattice spots appearing just above the background signal. Similar diffraction spots can be seen in many of the prediction FFTs including those predicted from P y8 (the 8 neighbor P y values as used for the Pyy GP prediction), P x8 (the 8 neighbor P x values), V; h; a; b; I 1 ; I 5 f g , V; h; a; b f g , and V; h f g. In order to compare the intensities of the diffracted spots, line profiles of each FFT were taken such that the diffraction spots are included, as seen in Figs. 7(a) and 7(i). The normalized intensities of the diffracted spots for each descriptor group are shown in Fig. 7(j). The maximum intensity of the diffracted spots occurs when using the full set of descriptors V; h; a; b; I 1 ; I 5 f g , with dominant contributions coming from V; h; a; b f gand V; h f g. Here, further analysis established that h is the preponderant origin of the modulation. In comparison, I 1 ; I 5 f g and a; b f g correspond to almost zero predicted modulation.

While the diffracted spots appear in the FFT, the corresponding ordered region in the real space ground truth (P y ) image is difficult to identify, though they can clearly be seen in the real space images of the predictions from V; h; a; b; I 1 ; I 5 f g , V; h; a; b f g , and V; h f g. Completely unexpectedly, the FFT for I 1 ; I 5 f gshows a modulation of a different type, corresponding to ( 12 , 1 2 ); however, they are very weak. The detailed analysis of the data suggests that these might be originating from the bottom right corner of the sample, as can be confirmed by inverse FFT analysis. Note that the prediction of the period doubling phase by GP is a highly non-trivial result. We pose that a possible explanation for this behavior is that the compositional and structural descriptors are that of the orthorhombic phase, but its emergence is suppressed due to the confinement and mechanical depolarization effects.

To explore whether behaviors and uncertainties derived using Gaussian process regression can be derived via linear analysis, we explored the linear correlation structure between the predictors, predictions, and the ground truth images using canonical correlation analysis (CCA). 90 In this analysis, two sets of multi-dimensional variables (X, Y) are connected by finding linear combinations of variables that maximize their correlation. New variables are defined by the linear combinations with weights w a and w b such that U ¼ w a 0 X and V ¼ w b 0 Y. The weights (w a , w b ) are determined by maximizing Eq. ( 4), subject to unit variance, w a 0 Xw a ¼ w b 0 Yw b ¼ 1, to account for scaling

The Pearson correlation coefficient, which quantifies the magnitude of association between two variables, 91 is then calculated between the learned pair of canonical variables. Specifically, the correlation coefficients are calculated given the (X, Y) pairs: ({descriptor group}, P y ), ({predictions from descriptor group}, P y ), and ({descriptor group}, {predictions from descriptor group}), where ({}, P y ) is each descriptor group or prediction from a descriptor group. In the last case, each descriptor group is only compared to the prediction derived from that group, i.e., no cross-comparisons are made between descriptor groups and predictions. , and (j) the normalized FFT intensity of the diffraction spot, outlined in red in (i), as a function of prediction group. Note that (i) predicted FFT intensity in group 1 is two times higher than in original image, and group 5 predicts a different ordering than original dataset, in h11i rather than h01i direction.

The CCA analysis of the pairwise combinations between the ground truth, descriptors, and GP predictions is shown in Fig. 8 for the chosen groups of descriptors and for three studied compositions. As a general trend, the descriptors show relatively high correlations with the prediction, which is maximal for the pure BiFeO 3 and decreases for higher Sm concentrations, across all descriptor groups. However, this correlation is well below unity and can be as small as $0.5 for high Sm content BFO. It is also important to note that the linear correlation explored here provides average information across image. Hence, the behaviors observed at the localized structural elements (domain walls, interfaces, etc.) must be interpreted independently, as above.

The correlation between the prediction and the truth is generally above the correlation between the descriptors and predictors, clearly illustrating that non-linear GP regression is capable of better prediction than linear estimates. This improvement is particularly pronounced for the high Sm composition, where corresponding correlation coefficients are higher by as much as a factor of 2.

The analysis of the quality of prediction (i.e., correlation between prediction and truth) for different descriptor groups across the compositions reveals several universal observations. In all cases, it appears that the P y8 group predicts P y extremely well. This behavior is expected, given the strong correlations between the polarization for adjacent unit cells. The correlation between P x8 and P y is fairly strong in BFO and becomes considerably weaker for the high Sm composition. Finally, extremely interesting behaviors are observed for the structural and chemical descriptors. Here, for the strongly ferroelectric material, the best prediction is achieved with the structural descriptors including V, h, a, and b, whereas inclusion of chemical descriptors I 1 and I 5 considerably lowers the prediction. The opposite behavior is observed for the intermediate and high Sm concentrations, where the quality of prediction improves considerably once the chemical effects are included.

The behavior of GP prediction uncertainty across the compositions and descriptor groups is shown in Fig. 9. Note that uncertainty is available only for the Gaussian (or, more generally, Bayesian) methods. As expected, for polarization descriptors P y8 ! P y , the uncertainty is lowest, i.e., prediction is highly reliable. At the same time, highly non-trivial behavior is observed for the predictions based on the different polarization component P x8 and structural and chemical descriptors. Here, the uncertainty in prediction is highest for BFO. This makes sense, since in strongly ferroelectric BFO polarization is less sensitive to the chemical and structural factors because the corresponding potential well is deepest and depolarization effects are strongest. The orthorhombic phase shows the intermediate uncertainty, presumably because polarization is close to zero. Interestingly, the lowest uncertainty in prediction is observed for intermediate phase. We ascribe this behavior to the fact that it is a mixture of ferroelectric and non-ferroelectric phases which have different lattice parameters, and ferroelectric phase is polarized in one direction. Hence, polarization is most sensitive to the local compositional fluctuations, allowing for good predictability. Interestingly, this lability is also the origin of the high functional responses of morphotropic materials. At the same time, this information is absent in a linear correlation analysis.

Similar to correlation analysis, including additional parameters slightly reduces the uncertainty, e.g., prediction from V; h; a; b f g reduces uncertainty compared to a; b f g and V; h f g. Inclusion of I 1 and I 5 considerably reduces uncertainty for intermediate, somewhat less high Sm material, but not pure one. Again, this behavior comports well with the expectations because compositional effects can be especially pronounced for regions with phase competition.

To summarize, here we explored whether the predictability and quantified uncertainty of physical behaviors can be used to determine manifest and latent physical mechanisms. Classically, this approach is implemented via, e.g., comparison of theory to experiment, to establish the presence of latent physical variables or unknown mechanisms, and subsequently refine their functional form and associated model parameters. Here, we implement this paradigm via analysis of predictability and uncertainty of Gaussian process regression between multiple spatially distributed degrees of freedoms. In other words, we explored whether local deviations of the characteristic behaviors established globally can be used to identify the latent local physical behaviors.

This approach is implemented for multiferroic Sm-doped BiFeO 3 across the broad composition range as a model. Using the atomic resolution STEM data, we directly visualize the atomic structure and structural, physical, and chemical order parameter fields for the material. GP regression is used to establish the predictability of local polarization field from the different groups of parameters, including the adjacent polarization values and several combinations of physical and chemical descriptors including lattice parameters, column intensities, etc. We observe that certain elements of materials domain structure, including charged and uncharged domain walls and interface with substrate, are best predicted with certain specific combinations of descriptors, and this predictability and associated uncertainties are consistent across the composition series. For example, the 180 domain walls are best predicted from unit cell sizes, whereas charged walls are best predicted from the molar volume and angles. Interestingly, the GP regression is shown to considerably exceed predictions based on linear correlative models. The associated physical mechanisms generally agree with the expected behaviors for ferroelectric materials across the morphotropic phase boundary with the non-ferroelectric phase. Remarkably, we show that uncertainty in prediction provides a new channel of information, presumably allowing to identify the regions with latent mechanism, in this case compositional fluctuations. Finally, and surprisingly, the GP regression predicts the experimentally unobserved ordering in the Sm-doped BFO phase.

Overall, we believe that the proposed approach can be broadly used in exploration of potential novel physical mechanisms in both experimental and modeling datasets. For instance, such tools may be extremely useful for interpretation of molecular dynamics data of disordered systems, or for analysis of nanoscale segregation in alloys by inspection of data from atom probe tomography. Future work can focus on more advanced functional interpolators, such as neural network-induced Gaussian process models, or Bayesian neural nets, as well as more complex descriptors that take into consideration longer range effects (e.g., graph-based methods). Similarly, one may apply this methodology to determine mechanisms in light-induced dynamics with local probes-for instance, exploration of photoconductive effects and ionic transport in perovskite photovoltaics, or for isolating factors involved in ferroelectric fatigue from nanoscale imaging and chemical spectroscopy.