Notice: This manuscript has been authored by UT-Battelle, LLC, under Contract No.

DE-
AC0500OR22725 with the U.S. Department of Energy. The United States Government retains and
the publisher, by accepting the article for publication, acknowledges that the United States
Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or
reproduce the published form of this manuscript, or allow others to do so, for the United States
Government purposes. The Department of Energy will provide public access to these results of
federally sponsored research in accordance with the DOE Public Access Plan
(http://energy.gov/downloads/doe-public-access-plan).

1
 Hypothesis learning in automated experiment: application to combinatorial
materials libraries
Maxim Ziatdinov,1,2,a Yongtao Liu,2 Anna N. Morozovska,3 Eugene A. Eliseev,4
Xiaohang Zhang,5 Ichiro Takeuchi,5 and Sergei V. Kalinin2,b
1
Computational Sciences and Engineering Division, Oak Ridge National Laboratory, Oak
Ridge, TN 37831
2
Center for Nanophase Materials Sciences, Oak Ridge National Laboratory, Oak Ridge, TN 37831
3
Institute of Physics, National Academy of Sciences of Ukraine, 46, pr. Nauky, 03028 Kyiv,
Ukraine
4
Institute for Problems of Materials Science, National Academy of Sciences of Ukraine,
Krjijanovskogo 3, 03142 Kyiv, Ukraine
5
Department of Materials Science and Engineering, University of Maryland, College Park, MD
20742

Machine learning is rapidly becoming an integral part of experimental physical discovery via
automated and high-throughput synthesis, and active experiments in scattering and electron/probe
microscopy. This, in turn, necessitates the development of active learning methods capable of
exploring relevant parameter spaces with the smallest number of steps. Here we introduce an active
learning approach based on co-navigation of the hypothesis and experimental spaces. This is
realized by combining the structured Gaussian Processes containing probabilistic models of the
possible system’s behaviors (hypotheses) with reinforcement learning policy refinement
(discovery). This approach closely resembles classical human-driven physical discovery, when
several alternative hypotheses realized via models with adjustable parameters are tested during an
experiment. We demonstrate this approach for exploring concentration-induced phase transitions
in combinatorial libraries of Sm-doped BiFeO3 using Piezoresponse Force Microscopy, but it is
straightforward to extend it to higher-dimensional parameter spaces and more complex physical
problems once the experimental workflow and hypothesis-generation are available.

a
ziatdinovma@ornl.gov
b
sergei2@ornl.gov

2
 Recent advances in accelerated and combinatorial synthesis have allowed rapid fabrication
of immense number of materials’ compositions in the form of combinatorial spread libraries in
pulsed laser deposition,1-7 high-throughput microfluidic8-12 and pipetting robot sets.13-19 In parallel,
conventional sample fabrication has been drastically accelerated via laboratory robotization and
autonomous and combined human-automation workflows. 20-24 Complementary to the large-
throughput synthesis are rapid characterization methods that allow establishing structure-property
relations in the specific systems and in certain cases target functionalities of interest.
However, even with this outstanding progress, simple grid-based exploration of the
immense compositional spaces is impractical. For example, for a quarternary alloy with the
properties of interest peaking within 1% from the morphotropic phase boundary or quantum
critical point, a standard grid search would require sampling 10 6 compositions for discovery.
Correspondingly, numerous strategies for active learning in such settings were proposed, using
statistical learning methods to navigate compositional space. These were based on direct Bayesian
optimization (BO) of the functionality of interest 25-27 or relied on the presence of cheap proxy
signals to select candidates for expensive target measurements. In several cases, theoretical
calculations were used to augment the BO-based strategies. 28-30
However, fundamentally, this discovery process has been preponderantly based on the
Gaussian process (GP) regression,31, 32 a non-parametric method that learns the behavior of
function of interest over relatively low-dimensional parameter space. The kernel function defining
the degree of correlation within this parameter space is typically defined within a certain functional
form and the parameters of this function are learned from the available data at each
optimization/exploration step. However, the simple GP-based methods do not incorporate any
physical model of behavior within the system. Consequently, they can be prone to trivial solutions,
often require a fairly large number of steps to reconstruct the complicated data distributions, and
are limited to low-dimensional parameter spaces.
An alternative approach for the active learning is based on the reinforcement learning (RL)
methods.33-36 The RL field has experienced rapid growth over last several years and comprises
multiple subfields. However, common to simple RL algorithms is extreme data requirements,
stemming from the necessity to build either a low-dimensional model of the process, or discover
low-dimensional representations of the state and policies of the system. Parenthetically, it is
important to note that many classical RL problems become trivial if the generative physical model
of the system is known explicitly.
Recently, we have introduced the structured GP (sGP) approach combining the flexibility
of GP models with the expressive power of the physical priors, and demonstrated significant
acceleration of the active learning/optimization based on the sGP.37 Here, we demonstrate that this
approach can be further extended towards physics discovery via active learning of competing
hypotheses, effectively combining the sGP and reinforcement learning. We refer to this approach
as hypothesis learning. This approach is illustrated for exploring concentration-induced phase
transition in Sm-doped BiFeO3 using a Piezoresponse force microscopy (PFM) but can be expected
to be applicable universally.

3
Figure 1. Hypothesis learning. The workflow of the scientific discovery process based on co-
navigation of the experimental and hypothesis space.

The general concept of the proposed hypothesis learning approach is illustrated in Figure
1. Here, we assume a classical experimental set-up in the context of sequential materials synthesis
or sequential characterization of compositional- or parameter spread systems. The examples of the
former are classical wet chemical solid-state synthesis/ceramics processing, or film growth via
techniques such as pulsed laser deposition, atomic layer deposition, or chemical deposition.
Common to all these techniques is that the material is fabricated for a fixed set of parameters such
as component ratios, stoichiometries, and deposition conditions. Jointly, these form a (relatively)
low-dimensional parameter space spanning chemical and physical control variables. Once the film
is grown (or material is synthesized), its properties are explored and based on the results, a new
exploration point (e.g., different composition or growth condition) is chosen. Notably, a similar
approach can be applied for composition spread libraries (where the composition changes across
the sample plane), thickness or growth temperature gradients, 4, 38 or samples under heat, electric,
or magnetic field gradients, if the local measurements are performed sequentially. These in turn
can be based on the modalities of scanning probe microscopies (piezoresponse, conductive,
microwave, acoustic), electron microscopy, micro-Raman imaging, or on the classical
semiconductor characterization performed via microfabricated electrode arrays.
Traditionally, these synthesis-measurement cycles are performed using human intuition
and past experience for selecting initial parameter settings, and often proceed via dense grid
sampling of the compositional or parameter spaces. This can include characterization of the 1D
compositional library, or synthesis of the epitaxial solid solution films. However, such a strategy
is not efficient, since functionalities of interest are often concentrated in specific parts of the
parameter space, necessitating the development of active learning methods for experimental
planning, namely selection of the next experimental point based on the results of previous
experiments.
In the last several years, the emergence of the automated synthesis and control systems
have stimulated the implementation of the Gaussian Process-based active learning (GP-AL) and
Bayesian Optimization (GP-BO) methods. By definition, GP seeks to reconstruct a scalar function

4
𝑓(𝑥) over N-dimensional parameter space x ∈ RN from the sparse measurements (yi, xi), i = 1, …
n, where the measurements are the noisy values of the function. Crucial for this discussion is that
GP also provides uncertainty estimates for the predicted function values. Hence, it can be directly
used as a basis for the AL or BO, which both represent a sequential measurement protocol where
the subsequent point for measuring is selected based on the previous ones. In the GP-AL, we
operate in the pure exploratory mode seeking to minimize the uncertainty by selecting at each step
a point associated with the largest uncertainty. In the GP-BO, the predicted function values and
uncertainty are combined in the acquisition function and used to guide the optimization of the
(unknown) underlying function. The principles and applications of the GP-BO are described in
details in Ref[39].
The fundamental limitation of the GP-based AL and BO for physical sciences is that the
state of the system is learned in the form of the kernel function describing the correlations between
functional behavior. The latter can be chosen from a given functional form, most typically radial
basis function or Matern.40 In special cases, more complex kernel functions such as asymmetric
Gaussian or spectral kernels are chosen. However, while this approach allows building a non-
parametric model of system behavior, the kernel function itself is not directly linked to the physics
of the system. Curiously, in most cases the mean function of GP is taken to be zero. In this manner,
variability of system behavior is effectively reduced to correlation structure.
Note that the alternative approach for the analysis of physical behaviors can be found in
the field of classical Bayesian inference,41, 42 where a probabilistic model is constructed using the
known functional form(s) reflecting a physical model, and some prior knowledge on model
parameters.43-45 Given the experimental observations, the posteriors of the parameter values are
inferred, representing an increase of knowledge based on observations.
Here, we introduce a combined approach, referred to as hypothesis learning, integrating
active discovery and Bayesian Inference, and further integrate this approach with a structured GP.
Our method is summarized in Algorithm 1 and is based on the idea that during active learning a
correct model of the system’s behavior decreases the overall Bayesian uncertainty about the system
under study. The method consists of the (very) short warm-up phase and exploration phase and
assumes the existence of several measurements at randomly (or uniformly) initialized coordinates
of the parameter space (we call them ‘seed’ measurements). The role of the warm-up phase is to
produce a ‘momentum’ for the exploration phase. This is necessary since in the active learning
setup we aim at discovering the overall data distribution with a minimal number of steps which is
usually much smaller than the number of exploration steps in standard RL problems. We also must
supply to the algorithm a list of probabilistic models that, based on our understanding of the
system, can describe the system’s behavior. These models can come as standalone parametric
models with appropriate priors over their parameters or they can be wrapped into sGPs. The
practical differences between the two approaches will be explained further in the text. It is assumed
that there is one ‘correct’ model in the list, but we do not know which one. We also note that if the
‘correct’ model is absent, the flexibility of the sGP approach allows for convergence comparable
to vanilla GP.

5
 Algorithm 1 (Hypothesis Learning)
Inputs: Initial measurements D, List of unmeasured points 𝑋∗, List of models M (standalone or
wrapped into sGP) and corresponding sample-averaged rewards Ra (initialized at zeros),
Exploration steps Nsteps, Warm-up steps Nwarmup, Reward function R,  in the -greedy policy.
for i=1, …, Nwarmup do
for each model in M do
Run BI to obtain posterior samples for model parameters given D
Compute posterior predictive uncertainty, 𝕍[𝑓∗ ], over 𝑋∗
( ∗)
Store the total uncertainty, 𝕍 = ∑ 𝕍[𝑓∗ ], of the predicted function values
end
Reward the model that produced the lowest uncertainty in prediction
Use the rewarded model to derive the next measurement point, xnext = argmax(𝕍[𝑓∗ ])
Perform measurement in xnext; update D and 𝑋∗
end
Average model rewards over the warm-up steps, update Ra
for i=1, …, Nsteps do
Use -greedy policy to sample model from M based on Ra[i] or at random
Run BI to obtain posterior samples for parameters of selected model given D
Compute posterior predictive uncertainty, 𝕍[𝑓∗ ], over 𝑋∗
Reward/penalize according to R (𝕍 , 𝕍 ), update Ra
Derive the next measurement point, xnext = argmax(𝕍[𝑓∗ ])
Perform measurement in xnext; update D and 𝑋∗
end

Once we have our ‘seed’ measurements and a list of models, we proceed to the warm-up phase
where we conduct Bayesian inference (BI) for each model. In the BI, the prior, 𝑃(𝜃), represents
our knowledge about the system before the measurement. The measurement produces the data, D,
based on which the posterior distribution, 𝑃(𝜃|𝐷), is computed via Bayes formula as
( | ) ( )
𝑃(𝜃|𝐷) = (1)
( )
where 𝑃(𝐷|𝜃) represents the likelihood that this data can be generated by the model with
parameters 𝜃 and the denominator 𝑃(𝐷) = ∫ 𝑃(𝐷|𝜃)𝑃(𝜃) 𝑑𝜃 defines the space of possible
outcomes. This step yields “degree of trust” in the model via the derived posterior distributions.
In practice, the BI is performed using sampling algorithms based on the Markov Chain
Monte Carlo (MCMC) techniques.46-48 Here, specifically, we are using the No-U-Turn Sampler
(NUTS) algorithm as implemented in the NumPyro probabilistic programming library. 49 The
posterior predictive mean and variance for a new point, 𝑥∗ , given the measured data, 𝐷, are then
computed as
𝑃(𝑥∗ |𝐷) = ∫ 𝑃(𝑥∗ |𝜃)𝑃(𝜃|𝐷)𝑑𝜃 ≈ ∑ 𝑃(𝑥∗ |𝜃 , 𝐷) = 𝑓∗ , (2a)
𝕍[𝑓∗ ] = ∑ (𝑓∗ − 𝑓∗ ) , (2b)
where 𝜃 ~ 𝑃(𝜃|𝐷) are samples drawn from the posterior. We then select a model with the lowest
total (or median) uncertainty over all the unmeasured coordinates, 𝑋∗, and reward it according to
a pre-defined reward function. The uncertainty map produced by the rewarded model is then used
to select the coordinate(s) of the next measurement point the same way as in the standard AL setup.

6
This “warm-up” procedure can be repeated multiple times (Nwarmup in Algorithm 1) to ensure that
the model selection was not affected by bad initialization of the seed points. We note that
conducting BI for each model is computationally costly and time-consuming (the cost and time
quickly rise as new points are added) and in practice, depending on the number of models, we limit
the number of warm-up steps to ~1-5.
After the warm-up phase is completed, we use a standard -greedy policy to sample only
one model at each exploration step. The total uncertainty computed from the posterior of a sampled
model is compared to the total uncertainty recorded in the previous step resulting in either positive
(the uncertainty decreased) or negative (the uncertainty increased) reward. The next measurement
point is then derived from the uncertainty map produced by the sampled model.

Figure 2. (a) Synthetic data describing a 1D system with a discontinuous “phase transition”. The
full data (‘noisy observations’) is shown but is never seen by our models. (b-d) The results of
active learning after 15 steps (3 warm-up + 12 exploration) with standalone probabilistic models:
model 1 (b), model 2 (e), model 3 (d), with sample-averaged rewards denoted as Ra.

As a practical example chosen here, we are interested in the active learning of phase
diagram that has a transition between different phases. The phase transition manifests in a
discontinuity of a measurable system’s property, such as heat capacity. An example using synthetic
data is shown in Figure 2(a) where x-axis is a tunable parameter (such as temperature) and y-axis
is a system’s property that is measured experimentally. However, we usually do not know where
a phase transition occurs precisely, nor are we aware of the exact behavior of the property of
interest in different phases. We note that using a standard GP-AL is not an optimal choice in these

7
cases as simple GP struggles around the discontinuity point. Instead, we are going to utilize
structured probabilistic models of expected system’s behavior.
We hypothesize that the system of interest can be described by a piecewise function as

𝑔 (𝑥), 𝑥 < 𝑥
𝜇(𝑥) = (3)
𝑔 (𝑥), 𝑥 ≥ 𝑥

where 𝑔 (𝑥) are the functions describing the system’s behavior before and after the transition
point 𝑥𝑐 . We consider three different models based on Eq (3), namely i) a power-law behavior
before transition point and linear behavior after transition point (model 1), ii) a linear behavior
before and after transition point (model 2), iii) a power law behavior before and after transition
point (model 3). We write them down as three standalone probabilistic models and put Bayesian
priors over their parameters. Specifically, the parameter 𝑥𝑐 is modeled using a uniform prior while
the rest of the parameters (slopes and power-law coefficients) are described by normal priors (for
more implementation details, see the accompanying source code).
We apply the hypothesis learning (HL) algorithm to the system with a single-phase
transition shown in Figure 2(a). The algorithm is initialized with four seed datapoints and our
measurements at the coordinates suggested by the algorithm uncover new data points from the
existing ground truth. We defined a simple reward function,

+1, 𝕍 <𝕍
𝑅 𝕍 ,𝕍 = (4)
−1, 𝕍 ≥𝕍

where 𝕍 is model’s total uncertainty for the prediction over the unmeasured part of the parameter
space at exploration step i. Here, the  parameter in the -greedy policy was set to 0.4 in all
experiments, reflecting the need for exploration over small (tens of data points) number of
attempts. However, similar to classical RL methods, smaller  values and annealing of  can be
trivially implemented. Furthermore, the exploration policy can be based on the Thompson
sampling of the posteriors or use policy embeddings. Here, given the small number of theoretical
models, we did not explore these strategies.
Figure 2(b-d) shows the result of the Bayesian fit over the entire parameter space for the
three models after running HL for 15 steps. The model that received the highest reward (i.e., is
favored by our algorithm) clearly provided the best fit. Hence, we were able both to learn a correct
data distribution with a small number of sparse measurements while also identifying a correct
model that describes the system’s behavior.
The evolution of the total uncertainty during the exploration for a single warm-up step and
for three warm-up steps are shown in Figure 3 (a) and 3 (b), respectively, for the larger number of
exploration steps. Again, in both cases, the algorithm favored the third model where the system’s
behavior before and after the transition point is described by the power-law functions. To confirm
the reproducibility and robustness of our algorithm, we ran it for 30 different initializations. The
histograms with averaged rewards received by three models are shown in Figure 3 (c) and 3 (d)
for a single warm-up step and for three warm-up steps, respectively. The algorithm showed a good
performance in both cases although the delineation of the correct model (model 3) was better when
we ran the warm-up procedure three times.

8
Figure 3. Hypothesis-driven active learning with standalone structured probabilistic models. (a,
b) Evolution of the total uncertainty during the active learning with standalone probabilistic models
with one ‘warm-up’ step (a) and with three ‘warm-up’ steps (b). See Algorithm 1 for the definition
of the ‘warm-up’ steps.’ At each exploration step, the model is selected according to ε-greedy
policy with ε=0.4. For the warm-up steps, only the model that produced the lowest uncertainty is
shown. (c, d) Statistics of the sample-averaged rewards for 30 different random initializations of
the four seed measurements for one ‘warm-up’ step (c) and for three ‘warm-up’ steps (d).

One disadvantage of the current approach is that it fails when the list of probabilistic
models does not contain the ‘correct’ model, as one can quickly conclude from Fig. 2(bc, ). At the
same time, most theoretical models describe the experimental reality only to a certain
approximation, and the complete/correct model may simply be not available.
To address this issue, we utilize a structured Gaussian process (sGP) method where a fully
Bayesian GP is augmented by a structured probabilistic model of expected system’s behavior. 37
Recall that in standard GP, we place a multivariate normal (MVN) prior over a target function 𝑓
such that 𝑓~𝑀𝑉𝑁(𝑚(𝑥), 𝐾(𝑥, 𝑥 )) and 𝑦 = 𝑓(𝑥) + 𝜀, where 𝜀 is a normally distributed noise.
The prior mean function 𝑚 is usually set to 0 and the GP is completely defined by its kernel
function 𝐾. In the sGP, we replace the GP’s constant prior mean function with one of the
probabilistic models defined earlier. The parameters of the GP’s new probabilistic mean function
𝑚 and those of the GP’s kernel 𝐾 (here chosen to be Matern40) are inferred simultaneously via
NUTS. Hence, the posterior predictive uncertainty associated with the parameters of sGP’s mean
function is automatically taken into account in our active learning setup and provides additional

9
control via the choice of the model’s priors. The sGP’s probabilistic mean function is expected to
capture general or even partial trends in the data, but it does not have to be precise. Note that we
replace standalone models M in the Algorithm 1 with sGPs containing the same models without
changing any other part of the algorithm.

Figure 4. Hypothesis-driven active learning with structured Gaussian processes (sGPs). (a-c) The
results of active learning after 15 steps (3 warm-up + 12 exploration) with sGPs augmented by
three different structured probabilistic models (same settings as in Fig. 2b-d), with sample-
averaged rewards denoted as Ra. Results of the active learning with vanilla GP (i.e., with prior
mean function set to 0) is shown in (d) for comparison.

Figures 4 and 5 show the sGP results for the same set of conditions as used for the
standalone probabilistic models in Fig. 2 and 3. In the sGP, the flexibility of the GP kernel can
partially correct for the wrong model choice. This offers an advantage for the reconstruction from
sparse data but also introduces a potential hurdle for learning the correct model in our setup, since
even with a wrong model the total uncertainty will likely keep decreasing. Therefore, a bad
initialization may result in being “locked-up” with an incorrect model (Fig. 5a).
Fortunately, this can be easily addressed by performing several warm-up steps (Fig. 5b)
where we compare the absolute values of the total (or median) uncertainties between different sGP
models. Indeed, as one can see from the histograms of sample-averaged rewards for 30 different
initializations in Fig. 5c and 5d, the increased number of warmup steps allows for a more robust
identification of the correct model. The Bayesian fits over the entire parameter space for the three
models after running the HL for 15 steps demonstrate that even with ‘wrong’ models, the sGP can
still recover the data distribution reasonably well (Fig. 4a, b). At the same time, the best
10
reconstruction was produced by sGP with the third model as its prior mean function (Fig. 4c) that
received the largest reward (0.82), in agreement with the results of the standalone probabilistic
models and the known ground truth. For comparison, we also showed the results of the active
learning with a vanilla GP (i.e., GP with a prior mean function set to 0), which completely failed
around the discontinuity point (Fig. 4d). In the Supplementary Note I, we describe the application
of sGP to several more systems with discontinuous behavior, including a scenario where none of
the proposed models correctly describes a behavior of the full system.

Figure 5. Hypothesis-driven active learning with structured Gaussian processes (sGPs). (a, b)
Evolution of the total uncertainty during the active learning with sGPs augmented by three
different structured probabilistic models (the models are the same as in Fig. 3) with one ‘warm-
up’ step (a) and with three ‘warm-up’ steps (b). See Algorithm 1 for the definition of the ‘warm-
up steps.’ At each exploration step, the model is selected according to ε-greedy policy with ε=0.4.
For the warm-up steps, only the model that produced the lowest uncertainty is shown. (c, d)
Statistics of the sample-averaged rewards for 30 different random initializations of the four seed
measurements for one ‘warm-up’ step (c) and for three ‘warm-up’ steps (d)

Having confirmed the performance of our approach on synthetic data, we move to the real-
time PFM experiments on the combinatorial library of Sm-doped BiFeO 3, as shown in Figure 6.
This material system has been extensively studied in the context of the evolution of ferroelectric
properties, with the pure BiFeO3 being rhombohedral ferroelectric materials and 20% Sm-doped
BiFeO3 being the orthorhombic non-ferroelectric.50 The intermediate phases at the boundary
between crystallographically-incompatible phases are known to exhibit complex nanoscale

11
ordering patterns and can be associated with enhanced electromechanical responses. 51 However,
these behaviors are strongly dependent on the mechanical and electrical boundary conditions, 52-54
resulting in strong variability of behaviors between bulk ceramics and thin films. Generally, the
mechanisms describing the emergence of these properties are poorly understood and are of
continuous interest for physical community.55-61 The compositional spread library hence encodes
the full constant-temperature cross-section of the phase diagram in a given concentration range.
Previously, we have extensively explored the property of this combinatorial library using high-
resolution electron microscopy.45, 62
Here, we have explored the mesoscale electromechanical responses in this material system
using PFM experiments.63-66 This technique is based on the detection of the local
electromechanical response induced by application of small amplitude (~1 Vac) periodic bias (100-
300 kHz) to the scanning probe. The detection of the response as a function of slowly (~Hz)
varying large amplitude (3-10 V) waveforms allows measuring the local hysteresis loops. 67-71 The
shape of the hysteresis loop is closely linked to tip-induced domain dynamics, 72, 73 and thus
represents local polarization switching mechanisms. 74, 75 Here, we use the switching spectroscopy
PFM approach, where hysteresis loops are measured over the rectangular grid of points, offering
better statistics and information on the spatial variability of switching behavior. 76, 77 Important for
subsequent discussion is that in linear approximation, the PFM signal is independent on tip-surface
contact radius and hence surface topography, making PFM signal quantitative. 78-81 Shown in
Figure 6a are band excitation PFM results showing the sample morphology and domain structures
from three representative locations with varying Sm concentration, along with the picture of the
Sm concentration gradient sample in Figure 6b. Figure 6c shows hysteresis loops from six
locations corresponding to various Sm concentration, indicating the dependence of switching
behavior on Sm concentration.

Figure 6. Piezoresponse force microscopy (PFM) results of Sm doped BFO sample. (a)
Topography and corresponding band excitation PFM amplitude images from three representative

12
locations: 0% Sm side, center, and 20% Sm side, which correspond to the three optical images in
(b), respectively. (b) Picture of the Sm-BFO sample, which has a gradually increased Sm ratio
from the left side to the right side. The optical images under a Cypher microscopy optical camera
show different color due to Sm ratio variation, this color gradient can also be observed by naked
eyes; however, note that the color shown by the microscope camera deviates due to optical path
effects and illumination conditions. (c). Hysteresis loops from six selected locations obtained in
Band excitation piezoresponse spectroscopy measurements; note that 56 locations representing
gradual Sm concentration change were uniformly marked along the sample.

To determine the possible models of the hysteresis loop behavior across the compositional
phase transition, we develop a simple Ginburg-Landau based model. Here, we assume that
piezoresponse (PR) loops are determined by the ferroelectric polarization reversal via piezoelectric
coupling as:

𝑃𝑅 (V) ≅ 𝑃𝑅 + 𝑑 𝑃 (V), (5a)

𝛤 + 𝛼(𝑇, 𝑥)𝑃 + 𝛽𝑃 + 𝛾𝑃 + 𝛿 𝑙𝑛 = 𝐸 (𝑡), (5b)

where 𝑑 is an effective piezocoefficient related to a pure BiFeO3. Applied voltage is 𝑉 =

𝑉 (𝑡) + 𝑉 𝑠𝑖𝑛(𝜔𝑡), where 𝑉 (𝑡) is a slow-varying amplitude of the voltage applied to the tip.
( )
The acting field 𝐸 (𝑡) = − , where ℎ is the BiFeO3 thickness, and we assume that an ultra-
thin gap of width  separates the BiFeO3 surface from the tip. Constant 𝜀 is a relative background
permittivity of BFO, 𝜀 10, and 𝜀 ~(1 − 10) is a relative permittivity in the gap.82

The coefficient 𝛼(𝑇, 𝑥) linearly depends on the Sm concentration: 83, 84

𝛼(𝑇, 𝑥) = 𝛼 𝑇 − 𝑇 (𝑥) , 𝑇 (𝑥) = 𝑇 1− , (6)

where 𝑇 is the temperature, 𝑇 is a Curie temperature of a pure BiFeO3, and 𝑥 is the relative content
of rare-earth Sm atoms. The critical concentration 𝑥 is a fitting parameter. Note that Sm doping
also affects the parameters 𝑃𝑅 , 𝛽, 𝛾 and 𝛿, and the effect can be nonlinear and complex. Only the
shift 𝑇 (𝑥) contains a linear component.

Several cases of possible materials behaviors are analyzed in the Supplementary Note II.
Here, we considered the two possible models that relate the measured loop area S to the dopant
concentration x. The first model is defined as
𝑆 1− + 𝐶, 𝑥 ≤ 𝑥 ,
𝑆= (7)
𝐶, 𝑥 > 𝑥
where xc is the transition point. This model does not have a discontinuous jump at the transition
point. The second model is defined as

13
 𝑆= 𝑆 1− +𝐶 , 𝑥 ≤𝑥 , (8)
𝐶, 𝑥>𝑥

and is characterized by a discontinuity at xc. Note that this analysis is thermodynamic in nature
and, in particular assumes structurally homogeneous phases. As such, it does not account for
potential formation of morphotropic phases with enhanced electromechanical responses.
We note that, unlike in the examples on synthetic data, we do not know the true function
in the actual experiment. We, therefore, cannot exclude the possibility that none of these models
is correct. We place a uniform prior over the transition point location and log-normal priors on all
the remaining parameters of the models as well as GP kernel hyperparameters and do not change
them during the experiment.

Figure 7. Application of Algorithm 1 with two different sGP models to real experiment. (a)
Evolution of the uncertainty during the active learning with sGPs. (b, c) sGP predictions sampled
from the posterior distribution at step 11 (b) and 12 (c), before performing a new measurement.

The experiment starts by measuring two extreme points on the concentration axis. These
seed points are used to initialize Algorithm 1 with two sGP models. We used five warmup states
and seven exploration steps. The median uncertainty as a function of the active learning steps is
shown in Fig. 7a. At the end of the exploration, the first model received a cumulative sample-
averaged reward of 0.13. The second model received a reward of 0.56. This suggests that the
second model is better at describing the current system. However, the second model experienced
significant difficulties in convergence at each exploration step, as judged by the effective number
of MCMC samples and the Gelman-Rubin convergence statistics (see Supplementary Note III).
The convergence problems indicate that, even though the second model received a larger reward,
it may not be the correct model for describing the system’s behavior (see also a discussion in the
Supplementary Note I for synthetic data). To get a better understanding, we examined the sGP
predictions sampled from the posterior distribution at the two last steps of the exploration. Shown
in Fig. 7b and 7c are the sampled sGP predictions together with the experimental values measured
in the points suggested by the algorithm. We can see there’s a discontinuity-like behavior at
x~17%, which explains why the algorithm favored the second model, preceded by a prominent
peak-like feature. The latter was not a part of any prior model, which explains the problems with
convergence. The inconsistency between the prior model and observations resulted in larger
uncertainty around the transition region, allowing us to quickly discover a potentially new
behavior. This discovery, in turn, will require the development of an improved model of the

14
system’s behavior to account for the peak before the transition, i.e., the behavior associated with
the morphotropic behavior.

Figure 8. The “automated scientist” concept. Here, the initial hypothesis space is formed as a
random symbolic sequence up to a given length. The hard physical constraints (dimensionality,
conservation laws, symmetry) are used to select the physical hypothesis and provide hard priors
(e.g., non-negativity). The prior “soft” knowledge is used to ascribe relevant probability and form
parameter priors. These models form inputs to the (active) hypothesis learning algorithm. Here,
the part in red dashed rectangular corresponds to Figure 1.
We further consider the opportunity to combine the hypothesis learning algorithm with the
automated hypothesis generation, as illustrated in Figure 8. Over the last decade, several notable
frameworks have been proposed to discover of the laws and functional expressions from the
experimental data. They include Eureka originally developed for the dynamic data 85 and more
recent frameworks such as PySindy,86 PySR,87 SISSO88, and many others. The overview of
symbolic regression methods in materials science is given in Ref [ 89]. These methods operate with
the static data, but generally allow the generation and simplification of the symbolic functional
form consistent with the observations. Recent work illustrates how this approach can be
generalized to incorporate physical laws.90 Other developments in this area includes frameworks
such as Feynman AI.91
We argue that the same approach can be extended for experimental hypothesis learning,
interfacing hypothesis generation with experiment. Clearly, the hypothesis spaces formed by all
possible symbolic expressions are intractable for experimental verification, as are all physically
possible models. Hence, an additional stage will be down the selection of the hypotheses using the
prior knowledge derived from domain expertise or literature mining. 92-95 In the context of HL, this
includes assigning probabilities to the possible theories and identifying the priors on the relevant
parameters. The formed list of the hypotheses and corresponding priors can be used as an input for
the HL, substituting the model list in Fig. 1. It is important to note that this approach dynamically
allows further matching the degree of compression of the hypothesis space and expected
experimental budget. For example, when the experiments are high throughput, and fast, hypothesis
space can be expanded, whereas for expensive experiments, the effort in down selecting
hypotheses is justified. Similarly, hypothesis selection can be based on different selection
strategies, varying degrees of risk, level of assumption, i.e., within its own optimization strategies.

15
Overall, this approach combines hypotheses generation and testing and allows constraining the
experimental budget.
Finally, it is important to note that the hypotheses generation is not necessarily limited to
symbolic models and can also include, for example, the lattice Hamiltonians and force field
models. However, in this case, the inverse problems are considerably more complex and are
addressed only for special cases.96, 97 Similarly, this will necessitate transition from discrete models
with Bernoulli distribution to more complex representations such as continuous embeddings, etc.
To summarize, we introduced an active learning approach based on co-navigation of
hypothesis and experimental spaces via a combination of fully Bayesian structured Gaussian
Processes with reinforcement learning policy refinement. This approach closely resembles
classical human-driven physical discovery, when multiple alternative hypotheses realized via
models with adjustable parameters are tested during the experiment. We demonstrated this method
for exploring concentration-induced phase transitions in the combinatorial libraries of Sm-doped
BiFeO3 using PFM measurements.
Looking forward, we believe that the proposed approach makes a strong case for the
synergistic combination between scanning probe and other imaging methods with combinatorial
material science. While traditionally combinatorial studies have been limited by the need to
quantify materials structure and functionality across the libraries, t recent advances in the focused
X-Ray methods and SPM offer a solution. At the same time, the proposed hypothesis-driven
Gaussian Processes framework developed here further allows incorporation and selection between
physical models. Note that while the approach is implemented here for the 1D case, it is
straightforward to extend it to higher-dimensional parameter spaces and more complex physical
problems.

Acknowledgements: This effort (M.Z.) was performed and supported at Oak Ridge National
Laboratory’s Center for Nanophase Materials Sciences (CNMS), a U.S. Department of Energy,
Office of Science User Facility. Y.L. and S.V.K. was supported as part of the center for 3D
Ferroelectric Microelectronics (3DFeM), an Energy Frontier Research Center funded by the U.S.
Department of Energy (DOE), Office of Science, Basic Energy Sciences under Award Number
DE-SC0021118. A.N.M. was supported by the National Academy of Sciences of Ukraine (the
Target Program of Basic Research of the National Academy of Sciences of Ukraine "Prospective
basic research and innovative development of nanomaterials and nanotechnologies for 2020 -
2024", Project № 1/20-Н, state registration number: 0120U102306) and received funding from the
European Union’s Horizon 2020 research and innovation programme under the Marie
Skłodowska-Curie grant agreement No 778070. The work at the University of Maryland was
supported in part by the National Institute of Standards and Technology Cooperative Agreement
70NANB17H301. The authors gratefully acknowledge Dr. Bobby Sumpter for careful reading and
editing of the manuscript.

16
Experimental details

A Sm-doped BiFeO3 (Sm-BFO) thin film, grown on SrTiO3 substrate, with Sm concentration
gradually increases from 0% to 20% along the sample was used as a model system. The detailed
growth conditions and characterization are described in previous publications. 45, 62

Band excitation piezoresponse spectroscopy (BEPS) measurements were performed on an Oxford

Instrument Asylum Research MFP3D atomic force microscopy system. A National Instruments
DAQ card and a LabView framework are equipped for band excitation measurement. All
measurements were carried out using Budget Sensor Multi75E-G Cr/Pt coated AFM probes (∼3
N/m).

56 locations were marked on the side of sample, which correspond to the Sm concentration
spanning from 0% to 20%. Band excitation piezoresponse spectroscopy measurements were
performed on the locations determined by the hypothesis learning algorithm.

Code availability
Code is available without restrictions at https://github.com/ziatdinovmax/hypoAL

17
 References
1. Curtarolo, S.; Hart, G. L. W.; Nardelli, M. B.; Mingo, N.; Sanvito, S.; Levy, O., The high-
throughput highway to computational materials design. Nat. Mater. 2013, 12 (3), 191-201.
2. Takeuchi, I.; Famodu, O. O.; Read, J. C.; Aronova, M. A.; Chang, K. S.; Craciunescu, C.; Lofland,
S. E.; Wuttig, M.; Wellstood, F. C.; Knauss, L.; Orozco, A., Identification of novel compositions of
ferromagnetic shape-memory alloys using composition spreads. Nature Materials 2003, 2 (3), 180-184.
3. Ohkubo, I.; Christen, H. M.; Kalinin, S. V.; Jr., G. E. J.; Rouleau, C. M.; Lowndes, D. H., High-
throughput growth temperature optimization of ferroelectric SrxBa1−xNb2O6 epitaxial thin films using a
temperature gradient method. Applied Physics Letters 2004, 84 (8), 1350-1352.
4. Christen, H. M.; Rouleau, C. M.; Ohkubo, I.; Zhai, H. Y.; Lee, H. N.; Sathyamurthy, S.; Lowndes,
D. H., An improved continuous compositional-spread technique based on pulsed-laser deposition and
applicable to large substrate areas. Review of Scientific Instruments 2003, 74 (9), 4058-4062.
5. Wang, J.; Yoo, Y.; Gao, C.; Takeuchi, I.; Sun, X.; Chang, H.; Xiang, X.-D.; Schultz, P. G.,
Identification of a Blue Photoluminescent Composite Material from a Combinatorial Library. Science
1998, 279 (5357), 1712-1714.
6. Nikolaev, P.; Hooper, D.; Perea-López, N.; Terrones, M.; Maruyama, B., Discovery of Wall-
Selective Carbon Nanotube Growth Conditions via Automated Experimentation. ACS Nano 2014, 8 (10),
10214-10222.
7. Nikolaev, P.; Hooper, D.; Webber, F.; Rao, R.; Decker, K.; Krein, M.; Poleski, J.; Barto, R.;
Maruyama, B., Autonomy in materials research: a case study in carbon nanotube growth. npj
Computational Materials 2016, 2 (1), 16031.
8. Steiner, S.; Wolf, J.; Glatzel, S.; Andreou, A.; Granda, J. M.; Keenan, G.; Hinkley, T.; Aragon-
Camarasa, G.; Kitson, P. J.; Angelone, D.; Cronin, L., Organic synthesis in a modular robotic system
driven by a chemical programming language. Science 2019, 363 (6423), eaav2211.
9. Angelone, D.; Hammer, A. J. S.; Rohrbach, S.; Krambeck, S.; Granda, J. M.; Wolf, J.; Zalesskiy,
S.; Chisholm, G.; Cronin, L., Convergence of multiple synthetic paradigms in a universally programmable
chemical synthesis machine. Nature Chemistry 2021, 13 (1), 63-69.
10. Abdel-Latif, K.; Epps, R. W.; Bateni, F.; Han, S.; Reyes, K. G.; Abolhasani, M., Self-Driven
Multistep Quantum Dot Synthesis Enabled by Autonomous Robotic Experimentation in Flow. Advanced
Intelligent Systems 2021, 3 (2), 2000245.
11. Epps, R. W.; Bowen, M. S.; Volk, A. A.; Abdel-Latif, K.; Han, S.; Reyes, K. G.; Amassian, A.;
Abolhasani, M., Artificial Chemist: An Autonomous Quantum Dot Synthesis Bot. Advanced Materials
2020, 32 (30), 2001626.
12. Epps, R. W.; Felton, K. C.; Coley, C. W.; Abolhasani, M., Automated microfluidic platform for
systematic studies of colloidal perovskite nanocrystals: towards continuous nano-manufacturing. Lab on
a Chip 2017, 17 (23), 4040-4047.
13. Higgins, K.; Valleti, S. M.; Ziatdinov, M.; Kalinin, S. V.; Ahmadi, M., Chemical Robotics Enabled
Exploration of Stability in Multicomponent Lead Halide Perovskites via Machine Learning. ACS Energy
Lett. 2020, 5 (11), 3426-3436.
14. Higgins, K.; Valleti, S. M.; Ziatdinov, M.; Kalinin, S. V.; Ahmadi, M., Chemical Robotics Enabled
Exploration of Stability in Multicomponent Lead Halide Perovskites via Machine Learning. ACS Energy
Letters 2020, 3426-3436.
15. Heimbrook, A.; Higgins, K.; Kalinin, S. V.; Ahmadi, M., Exploring the physics of cesium lead
halide perovskite quantum dots via Bayesian inference of the photoluminescence spectra in automated
experiment Nanophotonics 2021.
16. Gu, E.; Tang, X.; Langner, S.; Duchstein, P.; Zhao, Y.; Levchuk, I.; Kalancha, V.; Stubhan, T.;
Hauch, J.; Egelhaaf, H. J.; Zahn, D.; Osvet, A.; Brabec, C. J., Robot-Based High-Throughput Screening of
Antisolvents for Lead Halide Perovskites. Joule 2020, 4 (8), 1806-1822.

18
17. Choubisa, H.; Askerka, M.; Ryczko, K.; Voznyy, O.; Mills, K.; Tamblyn, I.; Sargent, E. H., Crystal
Site Feature Embedding Enables Exploration of Large Chemical Spaces. Matter 2020, 3 (2), 433-448.
18. Zhao, Y.; Zhang, J.; Xu, Z.; Sun, S.; Langner, S.; Hartono, N. T. P.; Heumueller, T.; Hou, Y.; Elia,
J.; Li, N.; Matt, G. J.; Du, X.; Meng, W.; Osvet, A.; Zhang, K.; Stubhan, T.; Feng, Y.; Hauch, J.; Sargent,
E. H.; Buonassisi, T.; Brabec, C. J., Discovery of temperature-induced stability reversal in perovskites
using high-throughput robotic learning. Nature Communications 2021, 12 (1), 2191.
19. Li, Z.; Najeeb, M. A.; Alves, L.; Sherman, A. Z.; Shekar, V.; Cruz Parrilla, P.; Pendleton, I. M.;
Wang, W.; Nega, P. W.; Zeller, M.; Schrier, J.; Norquist, A. J.; Chan, E. M., Robot-Accelerated
Perovskite Investigation and Discovery. Chemistry of Materials 2020, 32 (13), 5650-5663.
20. MacLeod, B. P.; Parlane, F. G. L.; Morrissey, T. D.; Hase, F.; Roch, L. M.; Dettelbach, K. E.;
Moreira, R.; Yunker, L. P. E.; Rooney, M. B.; Deeth, J. R.; Lai, V.; Ng, G. J.; Situ, H.; Zhang, R. H.; Elliott,
M. S.; Haley, T. H.; Dvorak, D. J.; Aspuru-Guzik, A.; Hein, J. E.; Berlinguette, C. P., Self-driving laboratory
for accelerated discovery of thin-film materials. Science Advances 2020, 6 (20), 8.
21. MacLeod, B. P.; Parlane, F. G. L.; Morrissey, T. D.; Häse, F.; Roch, L. M.; Dettelbach, K. E.;
Moreira, R.; Yunker, L. P. E.; Rooney, M. B.; Deeth, J. R.; Lai, V.; Ng, G. J.; Situ, H.; Zhang, R. H.; Elliott,
M. S.; Haley, T. H.; Dvorak, D. J.; Aspuru-Guzik, A.; Hein, J. E.; Berlinguette, C. P., Self-driving laboratory
for accelerated discovery of thin-film materials. Science Advances 2020, 6 (20), eaaz8867.
22. Coley, C. W.; Thomas, D. A.; Lummiss, J. A. M.; Jaworski, J. N.; Breen, C. P.; Schultz, V.; Hart,
T.; Fishman, J. S.; Rogers, L.; Gao, H.; Hicklin, R. W.; Plehiers, P. P.; Byington, J.; Piotti, J. S.; Green,
W. H.; Hart, A. J.; Jamison, T. F.; Jensen, K. F., A robotic platform for flow synthesis of organic
compounds informed by AI planning. Science 2019, 365 (6453), eaax1566.
23. Burger, B.; Maffettone, P. M.; Gusev, V. V.; Aitchison, C. M.; Bai, Y.; Wang, X.; Li, X.; Alston,
B. M.; Li, B.; Clowes, R.; Rankin, N.; Harris, B.; Sprick, R. S.; Cooper, A. I., A mobile robotic chemist.
Nature 2020, 583 (7815), 237-241.
24. Langner, S.; Häse, F.; Perea, J. D.; Stubhan, T.; Hauch, J.; Roch, L. M.; Heumueller, T.; Aspuru-
Guzik, A.; Brabec, C. J., Beyond Ternary OPV: High-Throughput Experimentation and Self-Driving
Laboratories Optimize Multicomponent Systems. Advanced Materials 2020, 32 (14), 1907801.
25. Nikolaev, P.; Hooper, D.; Perea-Lopez, N.; Terrones, M.; Maruyama, B., Discovery of Wall-
Selective Carbon Nanotube Growth Conditions via Automated Experimentation. Acs Nano 2014, 8 (10),
10214-10222.
26. Nikolaev, P.; Hooper, D.; Webber, F.; Rao, R.; Decker, K.; Krein, M.; Poleski, J.; Barto, R.;
Maruyama, B., Autonomy in materials research: a case study in carbon nanotube growth. npj Comput.
Mater. 2016, 2, 6.
27. Kusne, A. G.; Yu, H.; Wu, C.; Zhang, H.; Hattrick-Simpers, J.; DeCost, B.; Sarker, S.; Oses, C.;
Toher, C.; Curtarolo, S.; Davydov, A. V.; Agarwal, R.; Bendersky, L. A.; Li, M.; Mehta, A.; Takeuchi, I.,
On-the-fly closed-loop materials discovery via Bayesian active learning. Nat. Commun. 2020, 11 (1),
5966.
28. Stach, E.; DeCost, B.; Kusne, A. G.; Hattrick-Simpers, J.; Brown, K. A.; Reyes, K. G.; Schrier, J.;
Billinge, S.; Buonassisi, T.; Foster, I.; Gomes, C. P.; Gregoire, J. M.; Mehta, A.; Montoya, J.; Olivetti, E.;
Park, C.; Rotenberg, E.; Saikin, S. K.; Smullin, S.; Stanev, V.; Maruyama, B., Autonomous
experimentation systems for materials development: A community perspective. Matter 2021, 4 (9),
2702-2726.
29. Pollice, R.; Gomes, G. D.; Aldeghi, M.; Hickman, R. J.; Krenn, M.; Lavigne, C.; Lindner-
D'Addario, M.; Nigam, A.; Ser, C. T.; Yao, Z. P.; Aspuru-Guzik, A., Data-Driven Strategies for Accelerated
Materials Design. Accounts Chem. Res. 2021, 54 (4), 849-860.
30. Sun, S. J.; Tiihonen, A.; Oviedo, F.; Liu, Z.; Thapa, J.; Zhao, Y. C.; Hartono, N. T. P.; Goyal, A.;
Heumueller, T.; Batali, C.; Encinas, A.; Yoo, J. J.; Li, R. P.; Ren, Z.; Peters, I. M.; Brabec, C. J.; Bawendi,

19
M. G.; Stevanovic, V.; Fisher, J.; Buonassisi, T., A data fusion approach to optimize compositional
stability of halide perovskites. Matter 2021, 4 (4), 1305-1322.
31. Dai, C. Y.; Bruss, I. R.; Glotzer, S. C., Efficient Phase Diagram Sampling by Active Learning. J. Phys.
Chem. B 2020, 124 (7), 1275-1284.
32. Noack, M. M.; Doerk, G. S.; Li, R. P.; Fukuto, M.; Yager, K. G., Advances in Kriging-Based
Autonomous X-Ray Scattering Experiments. Scientific Reports 2020, 10 (1), 17.
33. Silver, D.; Huang, A.; Maddison, C. J.; Guez, A.; Sifre, L.; van den Driessche, G.; Schrittwieser,
J.; Antonoglou, I.; Panneershelvam, V.; Lanctot, M.; Dieleman, S.; Grewe, D.; Nham, J.; Kalchbrenner,
N.; Sutskever, I.; Lillicrap, T.; Leach, M.; Kavukcuoglu, K.; Graepel, T.; Hassabis, D., Mastering the game
of Go with deep neural networks and tree search. Nature 2016, 529 (7587), 484-+.
34. Schmidhuber, J., Deep learning in neural networks: An overview. Neural Netw. 2015, 61, 85-117.
35. Mnih, V.; Kavukcuoglu, K.; Silver, D.; Rusu, A. A.; Veness, J.; Bellemare, M. G.; Graves, A.;
Riedmiller, M.; Fidjeland, A. K.; Ostrovski, G.; Petersen, S.; Beattie, C.; Sadik, A.; Antonoglou, I.; King,
H.; Kumaran, D.; Wierstra, D.; Legg, S.; Hassabis, D., Human-level control through deep reinforcement
learning. Nature 2015, 518 (7540), 529-533.
36. Taylor, M. E.; Stone, P., Transfer Learning for Reinforcement Learning Domains: A Survey. J.
Mach. Learn. Res. 2009, 10, 1633-1685.
37. Ziatdinov, M. A.; Ghosh, A.; Kalinin, S. V., Physics makes the difference: Bayesian optimization
and active learning via augmented Gaussian process. arXiv:2108.10280 2021.
38. Ohkubo, I.; Christen, H. M.; Kalinin, S. V.; Jellison, G. E.; Rouleau, C. M.; Lowndes, D. H., High-
throughput growth temperature optimization of ferroelectric SrxBa1-xNb2O6 epitaxial thin films using a
temperature gradient method. Appl. Phys. Lett. 2004, 84 (8), 1350-1352.
39. Shahriari, B.; Swersky, K.; Wang, Z.; Adams, R. P.; Freitas, N. d., Taking the Human Out of the
Loop: A Review of Bayesian Optimization. Proceedings of the IEEE 2016, 104 (1), 148-175.
40. Rasmussen, C. E.; Williams, C. K. I., Gaussian Processes for Machine Learning (Adaptive
Computation and Machine Learning). The MIT Press: 2005.
41. Martin, O., Bayesian Analysis with Python: Introduction to statistical modeling and probabilistic
programming using PyMC3 and ArviZ, 2nd Edition. Packt Publishing: 2018.
42. Lambert, B., A Student’s Guide to Bayesian Statistics. SAGE Publications Ltd; 1 edition: 2018.
43. Somnath, S.; Law, K. J. H.; Morozovska, A. N.; Maksymovych, P.; Kim, Y.; Lu, X.; Alexe, M.;
Archibald, R.; Kalinin, S. V.; Jesse, S.; Vasudevan, R. K., Ultrafast current imaging by Bayesian inversion.
Nat. Commun. 2018, 9.
44. Vasudevan, R. K.; Kelley, K. P.; Eliseev, E.; Jesse, S.; Funakubo, H.; Morozovska, A.; Kalinin, S.
V., Bayesian inference in band excitation scanning probe microscopy for optimal dynamic model
selection in imaging. J. Appl. Phys. 2020, 128 (5), 10.
45. Nelson, C. T.; Vasudevan, R. K.; Zhang, X. H.; Ziatdinov, M.; Eliseev, E. A.; Takeuchi, I.;
Morozovska, A. N.; Kalinin, S. V., Exploring physics of ferroelectric domain walls via Bayesian analysis of
atomically resolved STEM data. Nat. Commun. 2020, 11 (1), 12.
46. Andrieu, C.; de Freitas, N.; Doucet, A.; Jordan, M. I., An Introduction to MCMC for Machine
Learning. Machine Learning 2003, 50 (1), 5-43.
47. Brooks, S.; Gelman, A.; Jones, G.; Meng, X.-L., Handbook of markov chain monte carlo. CRC
press: 2011.
48. McDannald, A.; Frontzek, M.; Savici, A. T.; Doucet, M.; Rodriguez, E. E.; Meuse, K.; Opsahl-
Ong, J.; Samarov, D.; Takeuchi, I.; Kusne, A. G., On-the-fly Autonomous Control of Neutron Diffraction
via Physics-Informed Bayesian Active Learning. arXiv preprint arXiv:2108.08918 2021.
49. Phan, D.; Pradhan, N.; Jankowiak, M., Composable Effects for Flexible and Accelerated
Probabilistic Programming in NumPyro. arxiv preprint: ArXiv:1912.11554 2019.

20
50. Kan, D.; Palova, L.; Anbusathaiah, V.; Cheng, C. J.; Fujino, S.; Nagarajan, V.; Rabe, K. M.;
Takeuchi, I., Universal Behavior and Electric-Field-Induced Structural Transition in Rare-Earth-
Substituted BiFeO(3). Adv. Funct. Mater. 2010, 20 (7), 1108-1115.
51. Troyanchuk, I. O.; Karpinsky, D. V.; Bushinsky, M. V.; Mantytskaya, O. S.; Tereshko, N. V.; Shut,
V. N., Phase Transitions, Magnetic and Piezoelectric Properties of Rare-Earth-Substituted BiFeO3
Ceramics. Journal of the American Ceramic Society 2011, 94 (12), 4502-4506.
52. Zeches, R. J.; Rossell, M. D.; Zhang, J. X.; Hatt, A. J.; He, Q.; Yang, C. H.; Kumar, A.; Wang, C.
H.; Melville, A.; Adamo, C.; Sheng, G.; Chu, Y. H.; Ihlefeld, J. F.; Erni, R.; Ederer, C.; Gopalan, V.;
Chen, L. Q.; Schlom, D. G.; Spaldin, N. A.; Martin, L. W.; Ramesh, R., A Strain-Driven Morphotropic
Phase Boundary in BiFeO(3). Science 2009, 326 (5955), 977-980.
53. Borisevich, A. Y.; Eliseev, E. A.; Morozovska, A. N.; Cheng, C. J.; Lin, J. Y.; Chu, Y. H.; Kan, D.;
Takeuchi, I.; Nagarajan, V.; Kalinin, S. V., Atomic-scale evolution of modulated phases at the
ferroelectric-antiferroelectric morphotropic phase boundary controlled by flexoelectric interaction. Nat.
Commun. 2012, 3.
54. Grinberg, I.; Suchomel, M. R.; Davies, P. K.; Rappe, A. M., Predicting morphotropic phase
boundary locations and transition temperatures in Pb- and Bi-based perovskite solid solutions from
crystal chemical data and first-principles calculations. J. Appl. Phys 2005, 98 (9).
55. Woodward, D. I.; Knudsen, J.; Reaney, I. M., Review of crystal and domain structures in the
PbZrxTi1-xO3 solid solution. Phys. Rev. B 2005, 72 (10).
56. Damjanovic, D., Ferroelectric, dielectric and piezoelectric properties of ferroelectric thin films
and ceramics. Rep. Prog. Phys. 1998, 61 (9), 1267-1324.
57. Dagotto, E., Complexity in strongly correlated electronic systems. Science 2005, 309 (5732), 257-
262.
58. Dagotto, E.; Hotta, T.; Moreo, A., Colossal magnetoresistant materials: The key role of phase
separation. Phys. Rep.-Rev. Sec. Phys. Lett. 2001, 344 (1-3), 1-153.
59. Imada, M.; Fujimori, A.; Tokura, Y., Metal-insulator transitions. Reviews of Modern Physics 1998,
70 (4), 1039-1263.
60. Kimura, T.; Goto, T.; Shintani, H.; Ishizaka, K.; Arima, T.; Tokura, Y., Magnetic control of
ferroelectric polarization. Nature 2003, 426 (6962), 55-58.
61. Tokura, Y., Critical features of colossal magnetoresistive manganites. Rep. Prog. Phys. 2006, 69
(3), 797-851.
62. Ziatdinov, M.; Nelson, C. T.; Zhang, X. H.; Vasudevan, R. K.; Eliseev, E.; Morozovska, A. N.;
Takeuchi, I.; Kalinin, S. V., Causal analysis of competing atomistic mechanisms in ferroelectric materials
from high-resolution scanning transmission electron microscopy data. npj Comput. Mater. 2020, 6 (1), 9.
63. Kolosov, O.; Gruverman, A.; Hatano, J.; Takahashi, K.; Tokumoto, H., NANOSCALE
VISUALIZATION AND CONTROL OF FERROELECTRIC DOMAINS BY ATOMIC-FORCE MICROSCOPY. Phys.
Rev. Lett. 1995, 74 (21), 4309-4312.
64. Eng, L. M.; Bammerlin, M.; Loppacher, C.; Guggisberg, M.; Bennewitz, R.; Luthi, R.; Meyer, E.;
Huser, T.; Heinzelmann, H.; Guntherodt, H. J., Ferroelectric domain characterisation and manipulation: a
challenge for scanning probe microscopy. Ferroelectrics 1999, 222 (1-4), 153-162.
65. Gruverman, A.; Kholkin, A., Nanoscale ferroelectrics: processing, characterization and future
trends. Rep. Prog. Phys. 2006, 69 (8), 2443-2474.
66. Balke, N.; Bdikin, I.; Kalinin, S. V.; Kholkin, A. L., Electromechanical Imaging and Spectroscopy of
Ferroelectric and Piezoelectric Materials: State of the Art and Prospects for the Future. Journal of the
American Ceramic Society 2009, 92 (8), 1629-1647.
67. Eng, L. M.; Grafstrom, S.; Loppacher, C.; Schlaphof, F.; Trogisch, S.; Roelofs, A.; Waser, R.,
Three-dimensional electric field probing of ferroelectrics on the nanometer scale using scanning force

21
microscopy. In Advances in Solid State Physics 41, Kramer, B., Ed. Springer-Verlag Berlin: Berlin, 2001;
Vol. 41, pp 287-298.
68. Gruverman, A.; Alexe, M.; Meier, D., Piezoresponse force microscopy and nanoferroic
phenomena. Nature communications 2019, 10 (1), 1-9.
69. Suzuki, K.; Hosokura, T.; Okamoto, T.; Steffes, J.; Murayama, K.; Tanaka, N.; Huey, B. D., Effect
of surface charges on the polarization of BaTiO3 thin films investigated by UHV-SPM. Journal of the
American Ceramic Society 2018, 101 (10), 4677-4688.
70. Denning, D.; Guyonnet, J.; Rodriguez, B. J., Applications of piezoresponse force microscopy in
materials research: from inorganic ferroelectrics to biopiezoelectrics and beyond. International
Materials Reviews 2016, 61 (1), 46-70.
71. Bark, C.; Sharma, P.; Wang, Y.; Baek, S. H.; Lee, S.; Ryu, S.; Folkman, C.; Paudel, T. R.; Kumar,
A.; Kalinin, S. V., Switchable induced polarization in LaAlO3/SrTiO3 heterostructures. Nano Letters 2012,
12 (4), 1765-1771.
72. Bdikin, I. K.; Shvartsman, V. V.; Kholkin, A. L., Nanoscale domains and local piezoelectric
hysteresis in Pb(Zn1/3Nb2/3)O-3-4.5%PbTIO3 single crystals. Appl. Phys. Lett. 2003, 83 (20), 4232-4234.
73. Bdikin, I. K.; Kholkin, A. L.; Morozovska, A. N.; Svechnikov, S. V.; Kim, S. H.; Kalinin, S. V.,
Domain dynamics in piezoresponse force spectroscopy: Quantitative deconvolution and hysteresis loop
fine structure. Appl. Phys. Lett. 2008, 92 (18).
74. Morozovska, A. N.; Eliseev, E. A.; Kalinin, S. V., Domain nucleation and hysteresis loop shape in
piezoresponse force spectroscopy. Appl. Phys. Lett. 2006, 89 (19).
75. Morozovska, A. N.; Kalinin, S. V.; Eliseev, E. A.; Svechnikov, S. V., Local polarization switching in
piezoresponse force microscopy. Ferroelectrics 2007, 354, 198-207.
76. Jesse, S.; Lee, H. N.; Kalinin, S. V., Quantitative mapping of switching behavior in piezoresponse
force microscopy. Rev. Sci. Instrum. 2006, 77 (7), 073702.
77. Jesse, S.; Baddorf, A. P.; Kalinin, S. V., Switching spectroscopy piezoresponse force microscopy
of ferroelectric materials. Appl. Phys. Lett. 2006, 88 (6), 062908.
78. Kalinin, S. V.; Karapetian, E.; Kachanov, M., Nanoelectromechanics of piezoresponse force
microscopy. Phys. Rev. B 2004, 70 (18), 184101.
79. Karapetian, E.; Kachanov, M.; Kalinin, S. V., Stiffness relations for piezoelectric indentation of
flat and non-flat punches of arbitrary planform: Applications to probing nanoelectromechanical
properties of materials. Journal of the Mechanics and Physics of Solids 2009, 57 (4), 673-688.
80. Kalinin, S. V.; Eliseev, E. A.; Morozovska, A. N., Materials contrast in piezoresponse force
microscopy. Applied Physics Letters 2006, 88 (23).
81. Eliseev, E. A.; Kalinin, S. V.; Jesse, S.; Bravina, S. L.; Morozovska, A. N., Electromechanical
detection in scanning probe microscopy: Tip models and materials contrast. J. Appl. Phys 2007, 102 (1).
82. Woo, C. H.; Zheng, Y., Depolarization in modeling nano-scale ferroelectrics using the Landau free
energy functional. Applied Physics a-Materials Science & Processing 2008, 91 (1), 59-63.
83. Morozovska, A. N.; Eliseev, E. A.; Chen, D. Y.; Shvetz, V.; Nelson, C. T.; Kalinin, S. V., Melting of
spatially modulated phases at domain wall/surface junctions in antiferrodistortive multiferroics. Phys.
Rev. B 2020, 102 (7), 14.
84. Morozovska, A. N.; Karpinsky, D. V.; Alikin, D. O.; Abramov, A.; Eliseev, E. A.; Glinchuk, M. D.;
Yaremkevich, A. D.; Fesenko, O. M.; Tsebrienko, T. V.; Pakalniskis, A.; Kareiva, A.; Silibin, M. V.; Sidski,
V. V.; Kalinin, S. V.; Kholkin, A. L., A combined theoretical and experimental study of the phase
coexistence and morphotropic boundaries in ferroelectric-antiferroelectric-antiferrodistortive
multiferroics. Acta Mater. 2021, 213, 14.
85. Schmidt, M.; Lipson, H., Distilling Free-Form Natural Laws from Experimental Data. Science 2009,
324 (5923), 81-85.

22
86. Kaptanoglu, A. A.; de Silva, B. M.; Fasel, U.; Kaheman, K.; Callaham, J. L.; Delahunt, C. B.;
Champion, K.; Loiseau, J.-C.; Kutz, J. N.; Brunton, S. L., PySINDy: A comprehensive Python package for
robust sparse system identification. arXiv preprint arXiv:2111.08481 2021.
87. Cranmer, M.; Sanchez-Gonzalez, A.; Battaglia, P.; Xu, R.; Cranmer, K.; Spergel, D.; Ho, S. J. a. p.
a., Discovering symbolic models from deep learning with inductive biases. 2020.
88. Ouyang, R.; Curtarolo, S.; Ahmetcik, E.; Scheffler, M.; Ghiringhelli, L. M., SISSO: A compressed-
sensing method for identifying the best low-dimensional descriptor in an immensity of offered
candidates. Physical Review Materials 2018, 2 (8), 083802.
89. Wang, Y.; Wagner, N.; Rondinelli, J. M. J. M. C., Symbolic regression in materials science. 2019, 9
(3), 793-805.
90. Ashok, D.; Scott, J.; Wetzel, S. J.; Panju, M.; Ganesh, V., Logic Guided Genetic Algorithms.
Proceedings of the AAAI Conference on Artificial Intelligence 2021, 35 (18), 15753-15754.
91. Udrescu, S.-M.; Tegmark, M., AI Feynman: A physics-inspired method for symbolic regression.
Science Advances 6 (16), eaay2631.
92. Kim, E.; Huang, K.; Saunders, A.; McCallum, A.; Ceder, G.; Olivetti, E., Materials Synthesis
Insights from Scientific Literature via Text Extraction and Machine Learning. Chemistry of Materials
2017, 29 (21), 9436-9444.
93. Vasudevan, R. K.; Ziatdinov, M.; Chen, C.; Kalinin, S. V., Analysis of citation networks as a new
tool for scientific research. MRS Bull. 2016, 41 (12), 1009-1015.
94. Young, S. R.; Maksov, A.; Ziatdinov, M.; Cao, Y.; Burch, M.; Balachandran, J.; Li, L. L.;
Somnath, S.; Patton, R. M.; Kalinin, S. V.; Vasudevan, R. K., Data mining for better material synthesis:
The case of pulsed laser deposition of complex oxides. J. Appl. Phys. 2018, 123 (11), 11.
95. Tshitoyan, V.; Dagdelen, J.; Weston, L.; Dunn, A.; Rong, Z. Q.; Kononova, O.; Persson, K. A.;
Ceder, G.; Jain, A., Unsupervised word embeddings capture latent knowledge from materials science
literature. Nature 2019, 571 (7763), 95-+.
96. Vlcek, L.; Ziatdinov, M.; Maksov, A.; Tselev, A.; Baddorf, A. P.; Kalinin, S. V.; Vasudevan, R. K.,
Learning from Imperfections: Predicting Structure and Thermodynamics from Atomic Imaging of
Fluctuations. ACS Nano 2019, 13 (1), 718-727.
97. Vlcek, L.; Maksov, A.; Pan, M. H.; Vasudevan, R. K.; Kalinin, S. V., Knowledge Extraction from
Atomically Resolved Images. Acs Nano 2017, 11 (10), 10313-10320.
98. Vasudevan, R. K.; Balke, N.; Maksymovych, P.; Jesse, S.; Kalinin, S. V., Ferroelectric or non-
ferroelectric: Why so many materials exhibit "ferroelectricity" on the nanoscale. Applied Physics Reviews
2017, 4 (2).
99. Ducharme, S.; Fridkin, V. M.; Bune, A. V.; Palto, S. P.; Blinov, L. M.; Petukhova, N. N.; Yudin, S.
G., Intrinsic ferroelectric coercive field. Phys. Rev. Lett. 2000, 84 (1), 175-178.

23
 Supplementary Materials

I. Application of hypoAL for active learning of different discontinuous functions.

Here we illustrate application of hypoAL (Algorithm 1) to two different discontinuous
functions. In the first case, there is one correct model in the input list M. In the second case, none
of the models in the input list is correct.
I.1. In the Supplementary Figure 1, we show data generated by the discontinuous
function of the form:
𝐴 + 𝐵 ⁄(𝑥 − 𝑥 ) , 𝑥 < 𝑥
𝑓(𝑥) = (S1)
𝐶 ⁄(𝑥 − 𝑥 ) , 𝑥 ≥ 𝑥

The noisy observations generated by this function together with the ‘true function’ are
shown in Supplementary Figure 1.

Supplementary Figure 1. The noisy observations of discontinuous function in Eq. S1

Same as in the main text example, the full data (“noisy observations”) in Supplementary
Figure 1 is never seen by the algorithm. Instead, we start with four ‘seed’ measurement to
initialize the algorithm. The algorithm is provided with the two possible models of system’s
behavior, wrapped into structured Gaussian processes. The first model is a Lorentzian peak
function of the form 𝐴 + 𝐵⁄ (𝑥 − 𝑥 ) + 𝐶 , with a Uniform(0, 1) prior on 𝑥 and Normal (0,
1) priors on A, B, and C parameters. The second model is the one from Eq. S1, with the same
priors on its parameters.
The hypoAL results after 20 steps (3 warm-up + 17 exploration) with the same settings as
used in the main text are shown in Supplementary figure 1, together with the reward values for
each model.
Supplementary Figure 2. Results of the hypotheses-driven active learning after 20 steps (3
warm-up + 17 exploration) with sGPs augmented by two different probabilistic models, with
sample-averaged rewards denoted as Ra.

I.2. In Supplementary Figure 3, we show noisy observations of a discontinuous function

which, unlike the one in the main text, has two transition points. We are going to use the same
HypoAL settigns with the same models, wrapped into structured Gaussian process, as in the
main text. This means that none of the models is correct (as they all assume only a single
transition).

Supplementary Figure 3. The noisy observations of discontinuous function with two “phase
transitions”.

The hypoAL results are shown in Supplementary Figure 4. The first two models that
assume a linear behavior after the transition point received high positive rewards whereas the third
model that assumes a power law behavior before and after the transition point received a negative
reward. Although none of the models is correct, the preference for the first two models can be
easily understood since it is easier to approximate the behavior after the first transition point with
a linear behavior than a power law one. In both cases, the flexibility of GP kernel partially
compensates for the incorrect model choice. We note that the models with high rewards described
correctly one of the transitions and provided a hint for the existence of the second transition.
Supplementary Figure 4. Results of the hypotheses-driven active learning after 25 steps (3 warm-
up + 12 exploration) with sGPs augmented by the same probabilistic models as in the main text,
with sample-averaged rewards denoted as Ra.

In the example above, the correct piece of information in all the models was the presence
of a transition point. The obvious question is: how the hypoAL algorithm distinguishes between
a model that is only partially correct and a model that is completely wrong? To test this, we run
hypoAL using two models: the third model in the example above and a wrong model of the form
ax2 + b with Normal(0, 1) priors over its parameters.

Supplementary Figure 5. Results of the hypotheses-driven active learning after 25 steps (3 warm-
up + 12 exploration) with sGPs augmented by a partially correct model (a) and by a wrong model
(b), with sample-averaged rewards denoted as Ra.

The results shown in Supplementary Figure 5 clearly demonstrate that when presented with
a wrong model (Fig. S5b) and a partially correct model (Fig. S5a), the algorithm chooses the
partially correct model. At the same time, a model that is only partially correct typically results in
a slow MCMC convergence at each step. In this case, the MCMC diagnostics may show a low
effective number of samples (n_eff) and high values (>1.05) of Gelman-Rubin criteria (r_hat), as
demonstrated below for the final AL step that performed Bayesian inference on the first model
parameters using NUTS with 5000 samples:
Such diagnostics can serve as an indicator that the model needs to be updated.

II. Phenomenological models of the local piezoresponse formation in Bi 1-xSmxFeO3

One of the most widely used model for the formation of vertical PR loops is the semi-
empirical equation of Landau type, which reflects the fact that PR loops are determined by the
ferroelectric polarization reversal via piezoelectric coupling, rather than by electrostrictive
coupling as it follows from the classical Landau theory. The mathematical formulation of the above
statement is:
𝑃𝑅 (V) ≅ 𝑃𝑅 + 𝑑 𝑃 (V), (1a)
𝛤 + 𝛼(𝑇, 𝑥)𝑃 + 𝛽𝑃 + 𝛾𝑃 + 𝛿 𝑙𝑛 = 𝐸 (𝑡), (1b)
where 𝑑 is an effective piezocoefficient related to a pure BiFeO3. Applied voltage is 𝑉 =
𝑉 (𝑡) + 𝑉 𝑠𝑖𝑛(𝜔𝑡), where 𝑉 (𝑡) is a slow-varying amplitude of the voltage applied to the tip.
( )
The acting field 𝐸 (𝑡) = − , where ℎ is the BiFeO3 thickness, and we assume that an ultra-
thin gap of width  separates the BiFeO3 surface from the tip. Constant 𝜀 is a relative background
permittivity of BFO, 𝜀 10, and 𝜀 ~(1 − 10) is a relative permittivity in the gap.82
The coefficient 𝛼(𝑇, 𝑥) linearly depends on the Sm concentration: 83, 84
𝛼(𝑇, 𝑥) = 𝛼 𝑇 − 𝑇 (𝑥) , 𝑇 (𝑥) = 𝑇 1 − , (2)
where 𝑇 is the temperature, 𝑇 is a Curie temperature of a pure BiFeO3, and 𝑥 is the relative content
of rare-earth Sm atoms. The critical concentration 𝑥 is a fitting parameter. Note that Sm doping
affect also on the parameters 𝑃𝑅 , 𝛽, 𝛾 and 𝛿, but the effect can be nonlinear and complex. Only
the shift 𝑇 (𝑥) contains a linear component.

Case I. The coefficients 𝛽 > 0 and 𝛾 = 𝛿 = 0 for the second order ferroelectrics of the
displacement type For the case the width of the PR loop is proportional to 𝐸 = − , and the
√

loop height is proportional to the spontaneous polarization 𝑃 = − . The PR loop area 𝑆, which
has a quasi-rectangular shape in a quasi-static case, is roughly proportional to their product, but
has sense only when 𝛼 < 0. So:
⎧𝑆 − 1− , 𝑇≤𝑇 1− ,
𝑆(𝑇, 𝑥) = (3a)
⎨
⎩ 0, 𝑇 ≥ 𝑇 1− .
Here 𝑆 and 𝑥 are fitting parameters, and 𝑇 is tabulated.
 For the purposes of the implementation in the code, we can introduce 𝑥 as the
concentration for which phase transition happens at room temperature, 𝑇 = 𝑇 1 − . With
this re-designation Eq.(3a) can be rewritten as:
𝑆 1− , 𝑥≤𝑥 ,
𝑆(𝑇 , 𝑥) = (3b)
0, 𝑥 > 𝑥 .
There is no jump at 𝑥 = 𝑥 and zero afterwards. To account for the possible electrostatic effects, 98
we additionally add the constant offset that we assume to be independent of concentration.

Case II. The coefficients 𝛽 < 0, 𝛾 > 0 and 𝛿 = 0 for the first order ferroelectrics of displacement
type. For the case, the width of the PR loop is proportional to the spontaneous polarization 𝑃 =
, the loop height is proportional to the coercive field 𝐸 = 2𝛽 +
⁄
9𝛽 − 20𝛼𝛾 .99 The PR loop area 𝑆, which has a quasi-rectangular shape
in a quasi-static case, is roughly proportional to their product, and has sense only in the polar phase:
⎧𝑆 2𝛽 + 9𝛽 − 20𝛼𝛾 , 𝑇≤𝑇 1− + ,
𝑆= .
⎨
⎩ 0, 𝑇≥𝑇 1− + ,
(3b)
where 𝛼(𝑇, 𝑥) is a linear function of x according to Eq.(2). Equation (3b) is too cumbersome for
the machine learning. Allowing for 𝑃 = at the transition temperature, 𝑇 = 𝑇 1− +

, it can be simplified further as:

⎧
⎪𝑆 1− − , 𝑇≤𝑇 1− +𝜃,
𝑆(𝑇, 𝑥) ≈ (4a)
⎨
⎪ 0, 𝑇≥𝑇 1− + 𝜃,
⎩
Here 𝑆 and 𝑥 are fitting parameters, 𝜃 = and 𝑇 are tabulated. Note the expression (4a) is
valid near the transition point only.
Introducing the concentration 𝑥 for which phase transition happens at room temperature, 𝑇 =
𝑇 1− , Eq.(4a) can be rewritten as:

𝑆(𝑇 , 𝑥) = 𝑆 1− +𝑆 , 𝑥 ≤𝑥 , (4b)
0, 𝑥 > 𝑥 .
There is a jump of the height 𝑆 at 𝑥 = 𝑥 and zero afterwards. Similar to case I, we add a constant
concentration-independent offset to account for the possible electrostatic effects

Case III. For the order-disorder type ferroelectrics we can put 𝛼 < 0, 𝛽 = 𝛾 = 0 and 𝛿 > 0.
Important than within the model 𝛼 is x-dependent but should be temperature independent.
Moreover, we can assume that the unknown parameter 𝑃 can be x-dependent. Eq.(1b) can be
rewritten as
𝛼(𝑥)𝑃 + 𝑘 𝑇𝑎𝑟𝑐𝑡𝑎𝑛ℎ − 𝑃 𝐸 = 0. (4a)
Using a series expansion 𝑎𝑟𝑐𝑡𝑎𝑛ℎ(𝑦) ≈ 𝑦 + 𝑦 + 𝑦 in Eq.(4a) we can reduce it either to the
Case I or to the Case II in dependence on the cutting term in the expansion. In the simplest case:
𝛼 + + =𝐸 , (4b)

From Eq.(4b), the coercive field 𝐸 = − 𝛼+ , and the loop height is proportional
√

to the spontaneous polarization 𝑃 = 𝑃 − 𝛼+ . The PR loop area 𝑆 is roughly

proportional to their product, but has sense only when 𝛼 + < 0. So:

𝑆 𝛼(𝑥) + , 𝛼(𝑥) + < 0,

𝑆(𝑇, 𝑥) = (5a)
0, 𝛼 + ≥ 0.
Here 𝑆 , 𝑃 and 𝛼 are fitting parameters, by the x-dependence can be assumed:
𝛼(𝑥) = 𝛼 1 − . (5b)
Note that this functional form is equivalent to case I, and hence these forms cannot be separated
in this approximation.

Table I. Material parameters for bulk BFO at room temperature 293K

coefficient BiFeO3
(collected and recalculated from Ref.[a])
Symmetry rhombohedral
b 9
7 -2
𝛼 (10 C ·mJ) 7.53, 𝑇 =1100 K, 𝛼 = 8.37 105C-2·mJ/K
8 -4 5 𝛽 11= 48, 𝛽 12= 8
𝛽 (10 C ·m J)
8 -6 9 Poorly known, maybe zero
𝛾 (10 C ·m J)

a) J. X. Zhang, Y. L. Li, Y. Wang, Z. K. Liu, L. Q. Chen, Y. H. Chu, F. Zavaliche, and R. Ramesh.

J. Appl. Phys. 101, 114105 (2007).

III. The MCMC diagnostics for hypoAL steps in Figure 7.

The effective number of samples (n_eff) and the Gelman-Rubin convergence criteria (r_hat) for
the two hypoAL steps shown in Figure 7a and 7b are shown below. In each case the Bayesian
inference was performed using NUTS algorithm with 5000 warmup steps (not to be confused with
the warm-up phase in the hypoAL) and 5000 samples.
Similar to the example in Supplementary Note I, the problems with MCMC convergence may
indicate that a model is only partially correct (for the completely wrong model, the GP kernel will
take over resulting in a good convergence but a trivial solution).