Exploring physics of ferroelectric domain walls via Bayesian analysis of atomically resolved STEM data

The physics of ferroelectric domain walls is explored using the Bayesian inference analysis of atomically resolved STEM data. We demonstrate that domain wall profile shapes are ultimately sensitive to the nature of the order parameter in the material, including the functional form of Ginzburg-Landau-Devonshire expansion, and numerical value of the corresponding parameters. The preexisting materials knowledge naturally folds in the Bayesian framework in the form of prior distributions, with the different order parameters forming competing (or hierarchical) models. Here, we explore the physics of the ferroelectric domain walls in BiFeO3 using this method, and derive the posterior estimates of relevant parameters. More generally, this inference approach both allows learning materials physics from experimental data with associated uncertainty quantification, and establishing guidelines for instrumental development answering questions on what resolution and information limits are necessary for reliable observation of specific physical mechanisms of interest.

U nique phenomena emerging at ferroelectric domain walls [1][2][3][4][5][6][7] have attracted the attention of researchers for many decades. Since the early days of ferroelectricity, it was recognized that the minimization of electrostatic energy of depolarization fields necessitates formation of ferroelectric domains, with the domain walls containing excess free energy compared to the bulk phase [8][9][10][11][12] . For several decades, the attention of the condensed matter physics and materials sciences communities alike was focused preponderantly on the domain properties and dynamics, whereas the walls were essentially treated as 2D objects. Correspondingly, even though remarkably advanced Ginzburg-Landau theory based theoretical models of wall structures were available as early as the late 1950s 13,14 , these theories were experimentally unverifiable beyond macroscopic thermodynamic descriptors due to the lack of the high-resolution imaging tools capable of probing wall structures on the nanometer and atomic levels.
This situation has changed drastically in the last fifteen years, in which the improvements in characterization tools made such studies possible and the interest of physics community shifted to the intrinsic physics and applications of the ferroelectric domain walls. After the discovery of enhanced local conductivity at ferroelectric domain walls 1 , tunable electronic properties have been demonstrated 4,15 and given rise to continuous research efforts towards domain wall electronics [2][3][4]7,15,16 . In conjunction with these experimental advances, a number of groups have theoretically explored the physics of the ferroic domain walls using the mesoscopic [17][18][19] and DFT models [20][21][22][23] , and have demonstrated that suppression of the primary order parameter at the wall core can give rise to additional magnetic or polar functionalities 5,[24][25][26][27] . The internal wall structure and hence conductivity are further strongly affected by the presence of flexoelectric interactions 28,29 , and can thus be used to establish the strength of the latter 30 . In addition to purely physical functionalities, the domain walls were also shown to interact with the chemical subsystem in materials 31 , giving rise to phenomena ranging from ferroelectric aging to vacancy segregation. While many of these theoretical advances suggest potential emergent physics at domain walls, experimental verification often remains a challenge even for atomic-scale, realspace imaging tools like (Scanning) Transmission Electron Microscopy (STEM).
Indeed, the equilibrium domain walls in proper ferroelectrics are usually very narrow, of the order of the atomic lattice parameter. Thus, STEM characterization of their internal structure has been enabled in the last decade by the introduction and proliferation of atomic-resolution spherical aberration corrected microscopes, allowing direct observations of the ferroelectric domain wall (and other interfaces) on the atomic level [32][33][34][35][36][37][38][39][40][41][42] . Quantitative information on the wall structure has been compared to Ginzburg-Landau-Devonshire (GLD) based models to yield materials parameters 40,[43][44][45] , and to validate DFT calculations 21 . However, despite the advance in imaging capabilities, the amount of material-specific information remains highly limited. Indeed, from early works of Ivanchik and Zhirnov it has been known that the order parameter profile across the domain wall is determined by the type of ferroelectric (second or first order), polarization behavior at the wall (Bloch, Ising, or Néel), and presence of the secondary order parameters or flexoelectric interactions 46,47 . Yet, while this information can theoretically be extracted from the wall profile, the experimental manifestation in the atomic structure can be subtle against instrumental noise or artifacts, is discretized at the level of atom positions, and is viewed in projection, precluding information from the 3rd spatial dimension (e.g., observation of Bloch character). This raises the statistical question of certainty in comparing multiple models to STEM profiles, or conversely, an estimation of the level of spatial resolution/information limit of the imaging system required to distinguish separate models.
Bayesian inference provides the probability of a model/ hypothesis given a set of experimental observations, a powerful statistical technique rising in popularity with corresponding computational and statistical tools (we refer interested readers to more comprehensive texts such as refs. [48][49][50] ). Here we develop this Bayesian inference framework for the analysis of material physics from structural imaging data. We demonstrate its application to the exploration of domain wall physics in the prototypical ferroelectric BiFeO 3 deriving the posterior probabilities for three GLD domain wall models from atomic resolution STEM experimental observations. Predicted domain wall shapes are dependent on parameters of the underpinning Landau theory. Under the assumptions of the validity of the theory, this then suggests that the domain wall profiles observed allow inversion to yield the parameters of the underlying Landau model, and thus, infer the order of the phase transition. We incorporate the effects of imaging noise and lattice discretization, and demonstrate how these can affect inferred materials properties. The required resolution limits to explore progressively fine details of domain wall physics are established, providing an answer to questions such as "how good of a microscope is necessary to address specific aspects of physics in a given materials system".

Results
As the first step in this analysis, we discuss the relationship between the physics of the ferroelectric material and the order parameter profile across the domain wall. For multiferroic materials with the general form of the long-range order parameters, the wall profiles can be found using the classical LGD approach. Two vector long-range order parameters, polarization components P i and oxygen octahedral tilts Φ i , were used for the description of the antiferrodistortive (AFD), ferroelectric (FE), and antiferroelectric (AFE) long-range orders in the rare-earth doped R x Bi 1−x FeO 3 , where R = Sm, La, Pr, Eu, etc. 27,[51][52][53] . For completeness, we also add the antiferroelectric (AFE) long-range orders to the description. The bulk part of LGD thermodynamic potential consists of several contributions, which are listed in Supplementary Methods of Supplementary Information. For further explanations, we list only the compact form of the FE and AFE contributions as: where the FE and AFE order parameters, , are introduced, P a i and P b i are the polarization components of two equivalent sublattices "a" and "b" [54][55][56] . As usual for proper and incipient ferroelectrics, the coefficients a k are temperature dependent and obey the linear law, a k (T) = α T [T-T C ], where T C is the Curie temperature, and T is the absolute temperature; negative a k (T) supports FE or AFE state. The sign and value of γ ab ij determines the AFE and FE phases coexistence. For a general case of domain structured or spatially modulated system, one should solve the coupled Euler-Lagrange (EL) equations of states, which are expressed via the variational derivatives of the functional (1), δG External and depolarization fields, E ext i and E d i , which contribute to the electric field, can be found from the electrostatic equation for electric displacement D, div D = 0, with boundary conditions at the surfaces, interfaces and/or electrodes. Elastic fields, which are, in fact, the secondary order parameters, satisfy equation of state and mechanical equilibrium equations, whereas the strains and/or stresses should be defined at the system boundaries.
As a rule, the biquadratic coupling to the polar subsystem is small, and depolarization field is absent for uncharged walls. In this case, it is possible to reduce the non-local free energies to the decoupled system of equations for the order parameters: where g P ijkl ¼ g aa ijkl þ g ab ijkl and g A ijkl ¼ g aa ijkl À g ab ijkl Since "aa" constants are for next-nearest neighbors sublattices, while "ab" constants are for the nearest neighbors sublattices, one can assume that, g aa ijkl ( g ab ijkl , and so g P ijkl % g ab ijkl and g A ijkl % Àg ab ijkl . Here, we derive analytical solutions for Eq. (1) for several specific cases as described below.
Model 1. For the ferroelectrics with the second order phase transition the order parameter P across the uncharged domain wall can be found in the one component and one-dimensional approximations. Namely: where P 2 S ¼ Àa 1 =a 11 is the spontaneous polarization, x 3 -x 0 is the distance from center of the domain wall plane, and L c ¼ is the correlation length 57 . Eq. (3a) is valid at a 1 < 0, a 11 > 0, a 111 = 0 [see Fig. 1a, d].
Model 3. For the ferroelectrics with the second order phase transition in the presence of possible polarization rotation, Pprofile can be found for the specific case, a 12 = 6a 11 46 , and the solution is the superposition of two tanh-profiles: the correlation length is L c ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi , and R 0 = x b -x a is an arbitrary constant. Exact expressions (4) describe rotational Ising-Bloch-type uncharged domain wall with Landau free energy density g LGD ¼ a 1 The ranges of the possible values for the 4 parameters in Eqs. (3)(4) are P S = (0.1−1) C/m 2 , L c = (0.5-5) nm, η = (0-100), and R 0 /L c = (0-10). Experimentally, an additional variable is the wall position, x 0 . Note that within the LGD approach, the "true" independent parameters are the coefficients in the GLD expansion, a 1 , a 11 , a 111 , and g 44 . However, for the analysis of the experimental data we treat the phenomenological wall parameters as independent.
Note that similar analysis can be performed for more complex physical mechanisms, albeit in these cases the numerical solutions are generally required. Further note that Model 1, Eq. (3a) is a special case of Model 2, Eq. (3b) for η = 0, as well as of Model 3, Eqs. (4) at R 0 = 0, as it can be expected from the physics of the problem. At the same time, Models 2 and 3 are alternative models, corresponding to the dissimilar physics of the material.
While the solutions (Eqs. 3-4) are limited for specific numerical values of free energy expansion, finite element analysis (FEM) confirms that a direct variational method with slightly more complex trial functions can be used to describe the rotation domain wall at a 12 ≠ 6a 11 : where P a,b,c , x a,b,c , a, b, c, and μ are variational parameters. These analyses set the context for the problem of physics determination from experimental data. Namely, the question we seek to answer is in which cases the more complex behaviors defined by Model 2, Eq. (3b), can be reliably differentiated from the simplest behavior of Model 1, Eq. (3a), thus differentiating materials with first and second order phase transitions. Similarly, how well can Models 1, 2, and 3, or the more complex models in Eq. (5), be separated given the noise level of the measurement system and the discretization in measurements induced by the underlying lattice. How can prior knowledge of the material system be used to narrow down the answers? And finally, given the possible range of materials properties, can we establish the requirements on microscope resolution/information limit required for these to be determinable?
The answers for these questions can be naturally explored in the context of Bayesian inference. Bayes formula relates the prior and posterior probabilities as where D represents the experimental data, p(D|θ i ) is the likelihood of the data given the model, i, with model parameters, θ i , and p(θ i ) is the prior, i.e., probability of the model. Finally, p(D) is the denominator that defines the total possible outcomes.
Despite the long history of Bayesian theory, its adoption by the basic science fields has been slow. Part of the reason is that only a few special distribution classes allow analytical solutions of Eq. (6), whereas many realistic model classes require numerical methods and extremely cumbersome integrations. Secondly, the classical argument against a Bayesian approach is the need for the prior distributions, and dependence of the answers on the priors. Here we note that while problemic in medicine, sociology, or economics, in physical areas domain knowledge is often sufficiently developed to provide meaningful priors, as explored below. The main point to note here is that in Bayesian inference we are aiming to compute the posterior distributions of the parameters of some model, conditioned on data, which is usually done with an appropriate sampling method such as the Metropolis algorithm. This allows one to numerically estimate the posterior in (Eq. 6) for all model parameters. Note that the variance in the data can itself be modeled by a parameter, as noted below.
Here, we formulate a Bayesian regression models based on Models 1, 2, and 3. We note that STEM data does not provide an absolute calibration for polarization magnitude and the wall position is a priori unknown. Correspondingly, we chose weakly informative priors for these parameters. Similarly, for a second order ferroelectric described by Eq. (3a), the correlation length, L c , is the sole parameter defining wall structure, and hence the corresponding prior can also be weak. For a first order ferroelectric described by Eq. (3b), the correlation length, L c , and η are parameters to be inferred. Note that model (3a) is the special case of (3b) for η = 0, and hence the target of Bayesian analysis is to establish whether η is practically equivalent to zero (and hence the material is second order), or nonzero, and hence the material is first order. Finally, model Eq. (4a) has a different functional form than Eq. (3b), and hence determination of the parameters L c , R 0 and separation of models 2 and 3 is the task for Bayesian inference.
The behavior of the posterior distributions was extensively explored using synthetic data made via a-priori known models (see Python notebook Supplementary Data 1) as a function of parameter values, noise level, and digitization step. It was established that for very thin domain walls with the thickness comparable with the lattice spacing, the lattice discretization becomes the most limiting factor in the analysis. Correspondingly, the inference allows only to distinguish the main features of the observed physical behaviors. Hence here as priors, we choose weakly informative priors for all associated parameters, as shown in Table 1.
Domain walls in the rhombohedral proper-ferroelectric archetype BiFeO 3 were taken as the model system for this analysis. A pseudocubic (pc) perovskite, BiFeO 3 adopts spontaneous polarization spatially degenerate along the eight <111> pc Polarization  1 2 Fig. 1 LGD ferroelectric domain wall models. a-c Free energy vs. order parameter-polarization components P i /P S and (d-f) distribution of these P i /P S across the domain wall for Models 1-3. Note that we fixed the coefficients a 1 , a 11 , and a 111 , and recalculated the spontaneous polarization P S and parameter η from them for Model 2. For Model 3, component P 1 is not observable by STEM. Table 1 Priors used in sampling.
x 0~U niform(0,L) x 0~U niform(0,L) P S~U niform(0.5 P est ,2 P est ) P S~U niform(0.5 P est ,2 P est ) P S~U niform(0.5 P est ,2 P est ) L C~U niform(0,5) L C~U niform(0,5) L C~U niform(0,5) var~HalfNormal (1) var~HalfNormal (1) R 0~H alfNormal (10) Ranges are shown in parenthesis, L is size of the image in unit cells, P est is estimated saturated polarization. ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-19907-2 directions with three possible rotation angles between adjacent domains: 71°, 109°, and 180°. A BiFeO 3 thin film was grown on a SrTiO 3 substrate by pulsed laser deposition, the electrostatic energy from the insulating substrate interface promoting polydomain formation. A cross section of the film was imaged along the <100> pc direction by high-angle annular dark field (HAADF) STEM, producing a mass-thickness sensitive picture of the atomic columns (Fig. 2a). This does not provide a direct measure of electric polarization, however BiFeO 3 exhibits a strong lattice coupling in the form of non-centrosymmetric Bi-site displacements, which can be used as a proxy, and we hereafter refer to them with the same P notation. This cation polar displacement was calculated for the 4-atom nearest neighbors after a 2orthogonal image scan-artifact reconstruction 58 and Gaussian fitting. A colorized vector map of this polar displacement is shown in Fig. 2b, clearly illustrating the polydomain structure of the film. In this case, typical equilibrium 109°and 180°domain walls are found forming on the [100] pc and ½ 101 pc planes, respectively. Subregions indicated by the white boxes are used as the input experimental data for Bayesian network analysis. A magnified view of the 109°domain wall from the region of interest (ROI) in Fig. 2b is shown in Fig. 2c Fig. 2d. A lattice coordinate space is utilized for subsequent analysis, defining the normal distance from the domain wall in lattice parameter units. The polarization components are defined in reference to the domain wall plane, the perpendicular (P[100] pc ) and parallel (P[001] pc ) axes for the 109°domain wall. A latticecoordinate spatial plot of these two components (Fig. 2e, f) and view of their profile data (Fig. 2g, h) illustrate the polarization transition occurring parallel to the plane. This component is used as the input for the Bayesian network analysis.
The model parameters for the 109°domain wall are shown in Fig. 3. The 90% highest posterior density interval is illustrated in the vicinity of the domain wall for all three models, overlaying mean values of the experimental data for each unit lattice distance, Fig. 3a. A least squares fit curve is also shown along with corresponding fit parameters. However, the power of the Bayesian analysis over a minimum value optimization is the probability information it provides. Further insight into the wall behavior can be inferred from the selected pairwise joint probability distributions as shown in Fig. 3b-d. Here, for most parameters the pairwise distributions are clearly marginalizable, i.e., the joint probability distribution can be well approximated as a product of the marginal distributions for individual parameters. This behavior implies that the variables are statistically independent, and the knowledge of one variable does not improve understanding of the other one. This behavior is clearly shown for wall position and polarization, P s -x 0 , pair. Similarly, P s -L c and x 0 -η are close to marginalizable.
At the same time, the posterior distribution for the L c and η shows very different behavior. The joint distribution is not marginalizable and shows strong parameter dependence. This implies that the knowledge of one of the parameters can significantly affect the amount of information we learn about the other one. As an example, the fit with the narrow parameter for L c in the (0.9-1.1) interval is compared to the elements in Fig. 3b This analysis clearly illustrates the natural way in which the Bayesian inference allows to explore the available data given the past knowledge of the physics of the system. We have further performed the analysis for the third model Eq. (4a) and associated parameters and distribution functions are provided in the Supplementary Data 1.
To compare the models, we apply the widely applicable information criteria (WAIC) 59 . Most model selection criteria are based on two terms: one term that describes how well the model fits to the observed data, and a second term that penalizes models with greater degrees of freedom. As such, the WAIC is equal to (LL-p WAIC ), where LL ¼   where LL is the log likelihood, which is calculated by drawing S samples from the posteriors of the parameters of a model, and the second term represents the effective number of parameters of the model, which is calculated by summing the variance of the loglikelihood. This is done for each model to compute the WAIC, and then the weights are calculated using a pseudo-Bayesian model averaging approach with AIC-type weighting 60 . The WAIC scores for the models are shown in Table 2. Based on the Bayesian analysis Model 1, the second-order ferroelectric, is the most likely. However, the relative probabilities of all three models are very close and demonstrate that given the observations for weakly informative priors the models cannot be distinguished. We argue (and confirm later) that this indistinguishability of the models in this case is due to the small domain width (compared to lattice spacing), which precludes the subtle differences in domain wall structure between models from being reliably identified.
We repeat this analysis for the case of the 180°domain wall ROI in Fig. 2a, the experimental data shown in Fig. 4. The 180°d omain wall corresponds to a reversal of the polarization along the same axis, also appearing as a 180°transition in [010] pc projection, albeit the P[010] pc component is again not being observed. The polar displacement is again broken into perpendicular and parallel components. For the ½ 101 pc domain wall, this corresponds to ½ 101 pc and [101] pc , respectively, and the latter is used as the input to the Bayesian analysis. There is an observable asymmetry with P not centered around zero, a known artifact of small off-axis mistilts which manifests in a similar polar cation asymmetry 61 . Since we are concerned with the polarization delta of the domain wall, this component is treated with a fixed offset term, P 0 , added to the fitting function.
The model parameters for the 180°domain wall are shown in Fig. 5. The 90% highest posterior density interval for the three models and least squares fits (Fig. 1a) are, like the 109°case, very similar, with Models 2 and 3 adopting parameters that approach Model 1, the second-order ferroelectric Eq. 3a. WAIC scores for the three models again also suggest Model 1 is most likely ( Table 2). Selected pairwise joint probability distributions from Model 2 are shown in Fig. 5b-g with results similar to the 109°c ase, albeit with an additional P 0 parameter. P 0 -P S (Fig. 5b), x 0 -P S (Fig. 5c), and η-x 0 (Fig. 5f) are clearly marginalizable, whereas L c -P S (Fig. 5d), and P 0 -x 0 (Fig. 5e) are somewhat close, and η-L c (Fig. 5g) is not. Although parameter correlations may not be linear, and clearly aren't in cases such as η-L c , Pearson correlation coefficients are helpful to highlight this distinction, with P 0 -P S (0.11), x 0 -P S (−0.04), η-x 0 (0.14), L c -P S (0.20), P 0 -x 0 (−0.27), and η-L c (−0.80).
The analysis above illustrates the determination of the physically relevant parameters of the material given the experimental observations as STEM atomic coordinates, and past knowledge in the form of the Bayesian priors on relevant materials parameters. As expected for ferroelectric materials with the extremely narrow domain walls, ultimately this consideration becomes the limiting factor in these studies. In other words, while domain wall shape is a measure of the physics of the order parameter in the system, practically the STEM observations are limited by the discreteness of the lattice and the noise in the system.
To explore the effect of the noise and sampling on the differentiability of specific physical behaviors from the observational data, we perform a range of numerical experiments on synthetic data. Here, the synthetic data set is generated using Model 2, Eq. (3b), for a first order ferroelectric with a specific set of ground truth parameters, chosen here to be (P s , L c , η) = (0.    Fig. 6 is the simulated domain wall profile for N = 40 and noise level σ = 0.001, 0.03, and 0.1. Weight values from comparison of WAIC scores from Models 1 and 2 as a function of noise level are shown in Fig. 6b. Here, 0 corresponds to identification of the correct Model (#2) of the generated data. The correct physical model can be determined from data only for noise levels below~0.02. For higher noise values the weight centers round 50% (chance), and the model cannot be established from experiment alone. However, as shown in Supplementary  Fig. 1, if the model and its parameters are partially known, they can be used as physics-based prior to determine materials properties more precisely.

Shown in
The corresponding inferred parameters are shown in Fig. 6c-e. Here, in all cases the inferred noise levelσ is very close to the ground truth level, σ. The point estimates of the domain wall parameters b P s ,L c ,η coincide with the ground truth values for small noises, and start to deviate for large noise level (abovẽ 0.01-0.02). For parameters b P s andx 0 (reconstructed wall position), the uncertainty grows~linearly with the noise, as can be expected from the functional form of the model. For parameterŝ L c andη the dependence is more complicated, and particularly for L c the inferred value is centered at the ground truth value and is weakly affected by noise.
Finally, we explore the combined effect of the sampling and the noise level on the separability of the physical models given experimental data and only weak physics-based priors. In this approach, we aim to answer the fundamental questions as to what level of microscope resolution/information limit is required to reliably determine the generative physics model from the data. Note that the issue of the information limit in the measured atomic positions and its relationship to microscope parameters is not explored here and we refer to other publications where these studies are performed [62][63][64] . Here, we seek to explore only the question as to which extent the correct (here, a priori known) physical model can be determined from experimental data given the sampling induced by lattice, and to which extent uncertainty in atomic positions (here, noise) affect this inference. However, we do not address the origin and mitigation strategies for the uncertainties for atomic positions, To explore this issue, we introduce the ground truth model corresponding to Model 2 with (P s , L c , η) = (0.5, 2.0, 2.0). We further create the calculated profiles at different samplings and noise levels. Bayesian inference is used to calculate the WAIC for Models 1 and 2. Where the WAIC comparison weight is close to 0, the models can be reliably separated, i.e., the physics of material can be established from the data. Where the value is 0.5 or larger the correct model cannot be determined from the data. The modeling results are shown in Fig. 7, where for convenience the units for pixel spacing are chosen to be comparable to domain wall width, i.e., 30/(N L c ). Figure 7a clearly illustrates that models can be distinguished as long as pixel size is comparable to the wall width, with the threshold value~1.9. At the same time, for the noise levels above 0.02 the model cannot be established. This threshold seems to be only weakly dependent on the pixel size.
We further explore the effect of prior physical knowledge on physics extraction. Figure 7b, c illustrates the effect of transition from weak to strong physical priors, where distributions of possible parameter values are much better defined. During the inference, such strong priors tend to produce uniform posteriors or posteriors sharply concentrated on the boundary or interval, as opposed to the Gaussian-like posteriors for weak priors. Note that the effect of strong prior can be roughly compared to the three-fold reduction in the noise level (compare Fig. 7a, c). Notably, the absolute knowledge of priors ( Fig. 7c) does not considerably improve analysis compared to strong priors (Fig. 7b). However, incorrect priors, Fig. 7d, has strongly deleterious effect on analysis, effectively precluding model inference for all but extremely high-quality data. Overall, the additional physical knowledge can provide significant improvement in the analysis. However, incorrect knowledge provides a much stronger effect, calling for care with analysis.

Discussion
To summarize, the physics of ferroelectric domain walls in BFO is explored via Bayesian inference analysis of atomically resolved STEM data. This approach allows for determination of materials parameters in the form of a relevant posterior distribution, based on prior materials knowledge in the form of prior distributions and available experimental data. The Bayesian inference can further be extended to analyze the likelihood of alternative models for materials physics. Here, we show that for noninformative or weakly informative priors (equivalent to classical point estimates in least-square fits) we can establish the model parameters and their posterior distributions, as well as attempt to distinguish the models. For the specific case of 109°and 180°d omain walls in BiFeO 3 , the combination of sampling and noise preclude a reliable differentiation of physical mechanisms.
However, incorporation of materials knowledge in the form of prior distributions of parameters can significantly narrow posterior distributions and allow for differentiation of generative mechanisms.
As a future perspective, we note that rather than using analytical models, Bayesian inference can be based on numerical solvers or even atomistic methods. However, this approach will considerably increase computational complexity, and may require approximate models based on Gaussian Processing 48,50 or deep learning-based interpolators. Similarly, this approach can be applied to other physical models, including those with more complex order parameters, in the presence of electronic or ionic screening. This includes non-equilibrium cases, e.g., strain during preparation or heating profile, as long as the numerical schemes for forward modeling are available.
More generally, Bayesian methods allow for a natural framework to distinguish possible physical mechanisms in observations, and allows us to very clearly ascertain to what extent we learn more from knowledge of the microscope and materials. This analysis clearly illustrates that additional physics or knowledge leads to the increase of the value of physical experiment. However, incorrect knowledge can strongly obviate any potential information gain. Taken over the multiple domains, it provides clear and quantifiable stimulus towards development of high-resolution microscopies and high information content probes such as four-dimensional (4D) STEM, and allows exploring cost-benefit considerations in the microscope development. , and deposition pressure of~100 mTorr from Oxygen flow after a base pressure of~2 × 10 −8 Torr.

Methods
The STEM sample was fabricated by a Focused Ion Beam cross-sectional liftout and subsequent Argon ion milling in a Fischione NanoMill with a final energy of 0.5 keV. STEM was performed on a NION UltraSTEM operating at 200 keV, the data utilized here corresponding to the High-Angle Annular Dark Field (HAADF) detector. Cation polar displacements were calculated for the 4-atom nearest neighbors after a 2-orthogonal image scan-artifact reconstruction 58 and Gaussian atom-fitting.

Data availability
The supporting data is publicly available in a Zenodo repository with identifier https:// doi.org/10.5281/zenodo.4122471.

Code availability
Analysis code for this work is available in a Python notebook as Supplementary Data 1.
Received: 14 April 2020; Accepted: 28 October 2020;  The prior distribution of parameters are labeled, black text indicates comparison values, blue represents additional knowledge that allows more informed prior, and red represents incorrect knowledge. a Weakly informed priors on polarization and correlation length. b Strongly informed priors. c Exact polarization and weakly informed correlation length priors. d Incorrect polarization and weakly informed correlation length priors.