††thanks: These authors equally contributed to this work.
Authors to whom correspondence should be addressed:
<EMAIL_ADDRESS><EMAIL_ADDRESS>These authors equally contributed to this work.
Authors to whom correspondence should be addressed:
<EMAIL_ADDRESS><EMAIL_ADDRESS>
# Stoichiometric growth of SrTiO3 films
via Bayesian optimization with adaptive prior mean

Yuki K. Wakabayashi  NTT Basic Research Laboratories, NTT Corporation, Atsugi,
Kanagawa 243-0198, Japan  Takuma Otsuka  NTT Communication Science
Laboratories, NTT Corporation, Seika-cho, Kyoto 619-0237, Japan  Yoshiharu
Krockenberger  NTT Basic Research Laboratories, NTT Corporation, Atsugi,
Kanagawa 243-0198, Japan  Hiroshi Sawada  NTT Communication Science
Laboratories, NTT Corporation, Seika-cho, Kyoto 619-0237, Japan  Yoshitaka
Taniyasu  Hideki Yamamoto  NTT Basic Research Laboratories, NTT Corporation,
Atsugi, Kanagawa 243-0198, Japan

###### Abstract

Perovskite insulator SrTiO3 is expected to be applied to the next generation
of electronic and photonic devices as high-$k$ capacitors and photocatalysts.
However, reproducible growth of highly insulating stoichiometric SrTiO3 films
remains challenging due to the difficulty of the precise stoichiometry control
in perovskite oxide films. Here, to grow stoichiometric SrTiO3 thin films by
fine-tuning multiple growth conditions, we developed a new Bayesian
optimization (BO)-based machine learning method that encourages the
exploration of the search space by varying the prior mean to get out of
suboptimal growth condition parameters. Using simulated data, we demonstrate
the efficacy of the new BO method, which reproducibly reaches the global best
conditions. With the BO method implemented in machine-learning-assisted
molecular beam epitaxy (ML-MBE), highly insulating stoichiometric SrTiO3 film
with no absorption in the band gap was developed in only 44 MBE growth runs.
The proposed algorithm provides an efficient experimental design platform that
is not as dependent on the experience of individual researchers and will
accelerate not only oxide electronics but also various material syntheses.

##  I Introduction

Perovskite insulator SrTiO3 (STO) (cubic structure with the lattice constant
of 3.905 Å), having a band gap of 3.2 eV, is one of the most promising
materials for oxide electronics.Shi  _et al._ (2020); Pai  _et al._ (2018);
Fujimoto and Kingery (1985); Sakuma  _et al._ (1998) It is expected to be
applied to high-$k$ capacitorsSakuma  _et al._ (1998); Galt  _et al._ (1993)
and photocatalystMavroides, Kafalas, and Kolesar (1976); Konta  _et al._
(2004); Ng  _et al._ (2010) owing to its high dielectric constant of 100–200
at room temperature,Sakuma  _et al._ (1998); Lee  _et al._ (2008) chemical
stability,Shi  _et al._ (2020) almost 100% quantum efficiency of
photocatalytic water splitting under ultraviolet light (UV),Mavroides,
Kafalas, and Kolesar (1976); Takata  _et al._ (2020) and compatibility with
other perovskite oxides.Pai  _et al._ (2018); Ohtomo and Hwang (2004); Bowen
_et al._ (2008); Wakabayashi  _et al._ (2021a); Takada  _et al._ (2021);
Wakabayashi  _et al._ (2021b); Takada  _et al._ (2022) In addition, when it is
doped by cation substitution, adding oxygen vacancies or cation vacancies,
many interesting physical properties or phenomena emerge, such as
superconducting states,Koonce  _et al._ (1967); Bert  _et al._ (2011); Li  _et
al._ (2011); Ahadi  _et al._ (2019) ferroelectricity,Jang  _et al._ (2010)
high mobility carriers,Son  _et al._ (2010); Matsubara  _et al._ (2016) and
blue light emission.Kan  _et al._ (2005) However, mid-gap states originating
from off-stoichiometry defects, such as oxygen and cation vacancies, are known
to cause leakage current in STO capacitors,Popescu  _et al._ (2014); Baek  _et
al._ (2020) and also cause mid-gap absorption that may decrease the
photocatalytic activity of STO.Chen  _et al._ (2018); Kumar  _et al._ (2020)
Therefore, to utilize the potential of STO as a high-$k$ capacitor or
photocatalyst, it is essential to grow stoichiometric STO epitaxial films
without mid-gap states. Since the growth of stoichiometric STO entails fine-
tuning of multiple growth conditions, including the supplied flux ratio of Ti
and Sr, the growth temperature, and oxidation strength in the case of
molecular beam epitaxy (MBE), only a few papers have reported highly
insulating stoichiometric STO films having the same lattice constants as bulk
STO and no absorption in the band gap.Lee  _et al._ (2016a); Flores  _et al._
(2017)

The conventional trial-and-error approach to optimizing the growth conditions
is time-consuming and costly, and the reproducibility of optimization depends
on the individual researcher. In contrast, data-driven decision-making
approaches have attained high-throughput in experiments where machine learning
models, such as Bayesian optimization (BO) and artificial neural networks, are
incrementally updated by newly measured data.Mueller, Kusne, and Ramprasad
(2015); Lookman, Alexander, and Rajan (2016); Burnaex and Panov (2015);
Agrawal and Choudhary (2016); Ueno  _et al._ (2018); Ren  _et al._ (2018);
Wakabayashi  _et al._ (2018); Li  _et al._ (2018); Hou  _et al._ (2019); Xue
_et al._ (2016); G.Baird, Liu, and Sparks (2022) BO is a sample-efficient
approach for global optimization,Snoek, Larochelle, and Adams (2012) which has
proven itself useful for streamlining the optimization of the thin film growth
conditions.Wakabayashi  _et al._ (2019a); Takiguchi  _et al._ (2020); Shimizu
_et al._ (2020); Wakabayashi  _et al._ (2022a) However, a technical challenge
for growth optimization has remained. Namely, the search procedure needs to
find a suitable parameter reliably since experiments are costly in terms of
time, labor, and expense. To this end, exploration needs to be encouraged,
especially when the suitable parameter region lies in a complex shape in the
search space. Such a complex-shaped growth parameter space may force the
search method to become stuck at suboptimal growth condition parameters.

In this study, to obtain highly insulating stoichiometric STO films, we
develop a new BO method that encourages the exploration by adapting the
hyperparameter of the prediction model to get out of suboptimal parameters. We
demonstrate the efficacy of our adaptation first by using simulated data and
then through implementation for real materials growth: machine-learning-
assisted molecular beam epitaxy (ML-MBE) of STO films. Results of MBE growth
and crystallographic analyses of grown samples are accumulated to produce the
next growth conditions with BO (Fig. 1). As a result, we developed highly
insulating stoichiometric STO films with lattice constants identical to that
of bulk STO. Visible-to-UV light spectroscopy shows no optical absorption in
the band gap, and the films were achieved in only $44$ MBE growth runs. The
reproducible highly insulating stoichiometric STO films will contribute to the
development of the next generation of electronic and photonic devices.

##  II Methods

###  II.1 Bayesian optimization with adaptive prior mean

This section outlines how BO tackles the optimization problem and its
adaptation of the prior mean function. Detailed formulations are presented in
Appendix A. BO is a method for optimizing a black box function in the form of
$y=f({\bf x})$, where function $f:\mathbb{R}^{D}\rightarrow\mathbb{R}$ is
unknown and expensive to evaluate given $D$-dimensional input ${\bf
x}\in\mathpzc{X}\subset\mathbb{R}^{D}$ of a specified search space
$\mathpzc{X}$. In materials growth optimization, ${\bf x}$ and $y$ represent
the growth parameters and physical properties used to evaluate grown
materials, respectively. Examples of physical properties include electrical
resistance and X-ray diffraction intensity. BO searches the parameter space by
repeating the following steps: (i) Construct a prediction model based on the
Gaussian process (GP)Rasmussen and Williams (2006) given the past $n$
observations $\mathpzc{D}_{n}=\left\\{({\bf x}_{i},y_{i})\right\\}_{i=1}^{n}$,
(ii) evaluate an acquisition function to find a promising ${\bf x}^{\prime}$
that is likely to give a good function output, and (iii) evaluate the new
point ${\bf x}_{n+1}={\bf x}^{\prime}$ and acquire its function value
$y_{n+1}$ to update the prediction model using
$\mathpzc{D}_{n+1}=\mathpzc{D}_{n}\cup\left\\{(y_{n+1},{\bf
x}_{n+1})\right\\}$.

The GP uses a prior mean function $\eta({\bf x})$ for predicting the value of
$f({\bf x}^{\prime})$ at unseen ${\bf x}^{\prime}$. Nevertheless, this
function is typically fixed as $\eta({\bf x})=0$ or a certain constant
Shahriari  _et al._ (2016); Frazier (2018); Snoek, Larochelle, and Adams
(2012) throughout the BO iterations because of the black box nature of $f$.
Our method still employs a constant function in the form of $\eta({\bf
x})=m_{0}$, but it adaptively calculates a hyperparameter $m_{0}$ using past
observations $\mathpzc{D}_{n}$. Figure 2 visualizes an example of difference
in prediction of outcome at unseen parameters caused by different choices of
$m_{0}$ (see also Appendix A.2). This example uses the two-dimensional Ackley
function as $f({\bf x})$ [Fig. 2 (a)], where the search progress has focused
on the left-top region. Since this function has four isolated peaks with
different heights, BO needs to explore the search space instead of persisting
in one of the peaks. In the example in Fig. 2, the $m_{0}=0$ case keeps on
searching the skirt of left-top sub-optimal peak [Fig. 2 (d)]. In contrast,
the $m_{0}\approx 0.2$ case jumps into the unexplored left-bottom area [Fig. 2
(e)]. While greater $m_{0}$ encourages the exploration, fixing $m_{0}$ at a
large value may not always be efficient: a large $m_{0}$ tends to produce an
optimistic prediction in an unexplored parameter space. This can trigger
unnecessary explorations leading to a plateau of the sub-optimal function
value. Thus, the choice of $m_{0}$ needs to take account of the balance
between the exploration and exploitation in the parameter search.

To mitigate this dilemma, we developed three methods called adaptive leveling
(AL), empirical Bayes (EB) and empirical Bayes uniform (EBu). The AL method
draws $m_{0}$ from the uniform distribution between the minimum and maximum of
the past observations to obtain a balance between the exploitation and
exploration by the stochastic choice. We also developed the EB and EBu methods
to justify the modification of hyperparameter of prior using the observations.
Their further descriptions are deferred to Appendix A.3.

###  II.2 ML-MBE growth and sample characterizations

Epitaxial STO films with a thickness of 60 nm were grown on STO (001)
substrates in a custom-designed MBE system with multiple e-beam evaporators
[Fig. 1 (a)]. Detailed information about the MBE system is described
elsewhere.Naito and Sato (1995); Yamamoto, Krockenberger, and Naito (2013);
Wakabayashi  _et al._ (2022b, 2019b, 2021c) We precisely controlled the Sr and
Ti elemental fluxes by monitoring the flux rates with an electron-impact-
emission-spectroscopy sensor and feeding the results back to the power
supplies for the e-beam evaporators. The oxidation during growth was carried
out with ozone (O3) gas ($\sim 15\%$ O3 $+85\%$ O2) introduced through an
alumina nozzle pointed at the substrate. For the stoichiometric STO growth, it
is important to fine tune the growth conditions (the ratio of the Ti flux to
the Sr flux, growth temperature, and local ozone pressure at the growth
surface).Kumar  _et al._ (2020); Naito, Yamamoto, and Sato (1998); Jalan  _et
al._ (2009); Ohnishi  _et al._ (2005) To systematically change the Ti flux
ratio to the Sr flux, we changed the Ti flux while keeping the Sr flux at 0.98
Å/s. The growth temperature was controlled by the heater shown in Fig. 1 (a).
We can adjust the local ozone pressure at the growth surface by changing the
O3-nozzle-to-substrate distance [Fig. 1 (a)] while keeping the flow rate of O3
gas at $\sim 0.25$ sccm.

We executed the BO algorithm in a three-dimensional space. The search windows
for the Ti flux rate, growth temperature, and O3-nozzle-to-substrate distance
were 0.20–0.33 Å/s, 600–900∘C, and 10–80 mm, respectively. We searched equally
spaced grid points for each parameter. The number of points of the respective
quantities was 100. Since the three-dimensional parameter space consisted of
1,000,000 ($100^{3}$) points, performing a trial for the entire space in a
point-by-point manner is unrealistic, as only several runs can be carried out
per day with a typical MBE system. In order to evaluate the stoichiometry of
the films, we measured $\theta$–$2\theta$ scanned x-ray diffraction (XRD)
[Fig. 1 (b)] since the increase in the lattice constant is a good indicator of
the magnitude of the off-stoichiometry of STO by changes in the cation and/or
oxygen concentration.Naito, Yamamoto, and Sato (1998); Jalan  _et al._ (2009);
Ohnishi  _et al._ (2008); Lebeau  _et al._ (2009) Therefore, we adopted the
difference in the $c$-axis lattice constant of the film and the substrate
($\Delta c$) as the evaluation value. Thermal conductivity might be considered
as a useful metric for evaluating the crystalline quality of SrTiO3 films with
high-quality films of sufficient thickness reproducing the thermal
conductivity observed in bulk single crystals.Oh  _et al._ (2011); Brooks  _et
al._ (2015) However, to make the thermal conductivity of a SrTiO3 film
distinguishable from that of the SrTiO3 substrate, a film thickness of several
hundred nm is required,Oh  _et al._ (2011); Brooks  _et al._ (2015) and
therefore thermal conductivity is not suitable for the evaluation of the
SrTiO3 films with a thickness of 60 nm used in this study. If films are thick
enough and the thermal conductivity measurements are reliable, optimization
using BO methods with thermal conductivity as the evaluation value should also
be possible. When XRD diffractions from the STO phase were indiscernible
and/or diffractions from SrO, TiO2, or Srn+1TinO3n+1 ($n$: integer;
$n$$\neq$$\infty$) Ruddlesden-Popper seriesLee  _et al._ (2013) precipitates
(impurity phases) appeared, we defined the evaluation value of those samples
to be the worst experimental $\Delta c$ value by that time. This imputation of
the missing data generated when the designated phase is not formed enabled a
direct search of the wide three-dimensional parameter space.Wakabayashi  _et
al._ (2022a)

Here, a black box function $\Delta c=f({\bf x})$ is the target function
specific to our STO films, and ${\bf x}$ represents the growth parameters (Ti
flux rate, growth temperature, and O3-nozzle-to-substrate distance). We used
data $\mathpzc{D}_{n}={({\bf x}_{i},y_{i})}_{i=1}^{n}$ [Fig. 1 (c)] obtained
from past $n$ MBE growths and XRD measurements [Fig. 1 (b)] of STO films to
construct a model to predict the value of $f({\bf x})$ at an unseen ${\bf x}$.
To this end, we used the GP to estimate the mean $m$ and variance $s^{2}$ at
an arbitrary parameter value ${\bf x}$ (see Appendix A.1 for details).
Specifically, the GP predicts the value of $f({\bf x})$ as a Gaussian-
distributed variable $\mathpzc{N}(m({\bf x}),s^{2}({\bf x}))$, where $m$ and
$s^{2}$ depend on ${\bf x}$ and $\mathpzc{D}_{n}$. In short, $m({\bf x})$ and
$s^{2}({\bf x})$ represent the expected value and uncertainty of $\Delta c$ at
${\bf x}$. To consider the inherent noise in the $\Delta c$ of STO films grown
under nominally the same conditions, the variance of the observation noise
$\sigma_{\varepsilon}^{2}$ of the GP model was set to $0.01$. In our
implementation, we used the Matérn $\frac{5}{2}$ kernel since it is good at
fitting functions with steep gradients.Snoek, Larochelle, and Adams (2012) We
iterated the routine after the initial MBE growth with five random initial
growth parameters and XRD measurements. First, the GP was updated using the
data set at the time [Fig. 1 (c)]. Subsequently, to assign the value of the
growth parameter in the next run, we calculated the expected improvement (EI)
[Fig. 1 (d)].Močkus, Tiesis, and Žilinskas (1978)

##  III Results and Discussion

###  III.1 Experiments with simulated data

This section investigates the optimization performance for simulated functions
by five methods: the baseline with $m_{0}=0$, a simple adaptation of $m_{0}$
by taking the average of observed data $m_{0}=\frac{1}{n}\sum_{i=1}^{n}y_{i}$
referred to as DA, which stands for data averaging, and the methods that
adjust $m_{0}$: AL, EB, and EBu. We use two functions—the Ackley
functionAckley (1987); Molga and Smutnicki (2005) and the Rosenbrock
functionRosenbrock (1960); Molga and Smutnicki (2005). The boundary of search
space $\mathpzc{X}$ was set at $-0.5$ and $2$ for each element. These
functions allow us to set an arbitrary number of dimensions $D$. In this
study, we used $D=2,4,$ and 6. Generally speaking, larger dimensionality $D$
makes black box optimization more challenging. The Ackley function with $D=2$
is illustrated in Figure 2 (a) and the Rosenbrock function with $D=2$ is
displayed in Figure 3, respectively. For each configuration, all methods
(baseline, DA, AL, EB, and EBu) were iterated until 100 observations were
obtained. These optimization processes were repeated five times with five
randomly chosen initial observations. Each evaluation contained noise with
$\sigma_{\varepsilon}^{2}=0.001$. These functions are described in Appendix B.

Figure 4 shows the optimization results for the Ackley and Rosenbrock
functions with $D=2,4$, and 6. Each curve indicates the best observation value
$\displaystyle\left(\max_{1\leq i\leq n}y_{i}\right)$ averaged the over five
runs as a function of the number of observations $n$. The shaded area
indicates the best and worst observations among the five runs. A curve that
rises with fewer observations indicates a better search algorithm, one that
requires fewer resources before reaching a high value of $y$. A vertically
narrow shaded area means that the method performs robustly against randomness
in the search process, such as the initial choice of parameters and random
seeds.

For the Ackley function with $D=2$ [Fig. 4 (a)], all methods reached the top
peak on average. In particular, the baseline, DA and AL found the optimal
value after 50 observations in all five trials of the optimization process.
Since EB and EBu were stuck at the second-best peak occasionally, their
average performance was slightly inferior to that of the baseline, DA and AL.
In larger dimensionalities $D=4$ and $6$ for the Ackley function [Figs. 4 (b)
and 4 (c)], the baseline method struggled to improve. This is because the
baseline tends to be trapped at one of the suboptimal peaks. In contrast, DA,
AL, EB, and EBu with adaptive $m_{0}$ clearly outperformed the baseline. This
supports the efficacy of the variable prior mean for BO. Nevertheless, none of
the methods attained the maximum value of 0.5. This explains the general
difficulty of black box optimization in high-dimensional search space,
especially when the objective function has multiple and disconnected peaks.
The performance of DA was worse in $D=6$ [Fig. 4 (c)]. Since DA method tends
to stabilize the value of $m_{0}$ with more observations, the algorithm failed
to explore novel regions.

Figures 4 (d)–4 (f) show the optimization results for the Rosenbrock function
with $D=2,4$, and 6. Unlike the case with the Ackley function, all methods
performed comparably well under all conditions. These results mean that the
optimization of the Rosenbrock function is easy due to its concave surface
[Fig. 3]. Since the configuration of search space $\mathpzc{X}$ was restricted
to a limited region between -0.5 and 2, a high evaluation value $f({\bf
x})\geq 0.3$ was present for a large portion of the search space. This allowed
all methods to work reasonably well in our experiment. With that said, DA
method showed a slower improvements when $D=4$ [Fig. 4 (e)] and $D=6$ [Fig. 4
(f)]. This is because an occasional observation of low $y_{i}$ decreases
$m_{0}$ of DA, which caused conservative and pessimistic prediction that led
to slower improvements.

Among AL, EB and EBu that change $m_{0}$ at each step, the performance was
similar in most configurations. With that said, some trials of EB and EBu were
inferior to that of AL in, for example, the optimization result of the Ackley
functions with $D=2$ and $6$ and the Rosenbrock function with $D=4$. Owing to
the stability and reliability of the performance, we adopted the AL method for
the ML-MBE growth of STO films.

###  III.2 Application to ML-MBE of STO film

To obtain stoichiometric STO films with no absorption in the band gap, we grew
STO films by the AL method implemented in ML-MBE. Figure 5 shows how the BO
algorithm predicts $\Delta c$ values with unseen parameter configurations and
acquires new data points. The process starts with five random initial growth
parameters and gains experimental $\Delta c$ values for the updated GP model
with 10 [Fig. 5 (a)], 27 [Fig. 5 (b)] and 44 [Fig. 5(c)] observations. Two-
dimensional plots of the predicted $\Delta c$, $s({\bf x})$ and EI values at
the O3-nozzle-to-substrate distance, at which the highest EI value was
obtained, are shown in the lower panels [Figs. 5 (d)–5 (l)]. Within the first
ten samples from the start, the STO phase had not formed at two growth
conditions [green spheres in Fig. 5 (a)]. Thus, we defined the $\Delta c$
value of these samples as the worst experimental one at that time. This
imputation of experimental failure enabled a direct search of the wide three-
dimensional parameter space.Wakabayashi  _et al._ (2022a) According to the GPR
prediction from the ten samples, the highest EI was obtained at the Ti flux
$=0.26$ Å/s, growth temperature $=777~{}^{\circ}$C, and O3-nozzle-to-substrate
distance $=39.5$ mm [Fig. 5 (a)]. This predicted growth condition yielded the
$\Delta c$ of $0.036$ Å, smaller than the minimum value of $0.045$ Å at that
time [Fig. 5 (b)]. The more surrounding data points there are, the smaller $s$
becomes, and the less there are, the larger it becomes. Therefore, the region
with relatively small $s$ of the predicted $\Delta c$ became larger as the
number of experimental samples increased from 10 to 44 [Figs. 5 (g)– 5 (i)],
meaning that the prediction accuracy had increased, resulting in the lower EI
values [Figs. 5 (j)–5 (l)]. Through this optimization process, in which the
exploration is encouraged by varying the prior mean, the lowest $\Delta c$
value decreased and reached an ideal value of 0 in only 44 MBE growth runs
(Fig. 6). The ideal $\Delta c$ of 0 was achieved at the Ti flux $=0.32$ Å/s,
growth temperature $=852~{}^{\circ}$C, and O3-nozzle-to-substrate distance
$=13.5$ mm [Fig. 5 (c)]. The achievement of the target material with the
desired properties in such a small number of optimizations demonstrates the
efficacy of the AL method for high-throughput materials growth.

The optimized Ti flux ($0.32$ Å/s) corresponds to the Ti/Sr ratio of 1.026,
which is near the stoichiometric value of 1. This may be due to the low vapor
pressure of the most stable strontium oxide SrO and titanium oxide
TiO2,Guguschev  _et al._ (2014) for which the sticking coefficients of both Sr
and Ti become almost 1—even if SrO and TiO2 are concomitantly formed during
the growth of SrTiO3, they do not desorb from the growth surface and are
eventually transformed to SrTiO3. Generally, in complex oxide films, if the
sticking coefficients of each constituent cation were known as functions of
growth temperature and oxidation strength, the optimum supplied flux ratio
could be predicted, at least in principle. However, one cannot optimize the
whole growth conditions by predicting or just using reasoning from a crystal
growth perspective. Instead, the combined impact of the growth temperature,
the oxidation strength, and the supplied flux ratio on the oxide film growth
can only be obtained empirically through experiments. This is because actual
crystal growth dynamics are complicated and cannot be comprehended even when
thermodynamic phase diagrams are available. The proposed method enables high-
quality and efficient materials growth, independent of the researcher’s
knowledge and experience, even in cases where such prior knowledge is limited.

###  III.3 Crystallographic and optical properties of stoichiometric STO films

We experimentally characterized the physical properties of the stoichiometric
STO film with $\Delta c$ of 0. For comparison, we also examined the physical
properties of the off-stoichiometric STO film with $\Delta c$ of 0.045 Å grown
under one of the first random growth conditions (Ti flux $=0.29$ Å/s, growth
temperature $=796~{}^{\circ}$C, and O3-nozzle-to-substrate distance $=52$ mm).
The sheet resistance was measured by a standard two-point method with Ag
electrodes deposited on the STO surface. The stoichiometric STO is highly
insulating, exceeding the measurable range (sheet resistance $>50$ M$\Omega$),
while the off-stoichiometric STO film shows a relatively small sheet
resistance of $1.4$ k$\Omega$. This result indicates that precise
stoichiometry adjustment is necessary to obtain highly resistive STO. The
crystallinity of the STO films was examined by XRD, atomic force microscopy
(AFM), and scanning transmission electron microscopy (STEM). Figure 7 (a)
shows XRD $\theta$–$2\theta$ scans of the stoichiometric STO film around the
(002) STO Bragg peak. The XRD $\theta$–$2\theta$ scans of the STO film grown
on (001) (LaAlO3)0.3(SrAl0.5Ta0.5O3)0.7 (LSAT) under the same growth
conditions are also shown. It has been reported that non-stoichiometric STO
films have larger lattice constants than stoichiometric ones. Kumar  _et al._
(2020); Naito, Yamamoto, and Sato (1998); Jalan  _et al._ (2009); Ohnishi  _et
al._ (2005, 2008); Brooks  _et al._ (2009) The XRD pattern for the
stoichiometric STO on STO shows a good overlap between film and substrate
peaks without XRD fringes. The lack of fringes in the XRD data—which would be
observed for finite repetition of the unit cell—is a typical feature of
stoichiometric STO films Kumar  _et al._ (2020); Ohnishi  _et al._ (2008)
since the films merge with the substrate and become indistinguishable. In
contrast, the fringes are clearly observed for the films heteroepitaxially
grown on the LSAT substrate, which allows for thickness estimation of the SRO
films. The film thickness estimated from the periods of the Laue fringes (62
nm) agrees very well with that calculated by the Sr flux rate (60 nm), whose
sticking coefficient is 1. Figures 7 (b) and 7 (c) show the AFM images of the
stoichiometric STO film. The root-mean-square roughness is 0.25 nm, indicating
that the stoichiometric STO film has smooth surfaces.

Figure 8 shows high-angle annular dark-field (HAADF)- and annular bright-field
(ABF)-STEM images of the stoichiometric and off-stoichiometric STO films taken
with a JEOL JEM-ARM 200F microscope. Since the intensity in the HAADF-STEM
image is proportional to $\sim\mathrm{Z}^{n}$ ($n\sim$ 1.7–2.0, and Z is the
atomic number),Pennycook and Jesson (1991) the brighter spheres and darker
ones in Figs. 8 (a) and 8 (c)] are assigned to Sr- ($\mathrm{Z}=38$) and Ti-
($\mathrm{Z}=22$) occupied columns, respectively. The ABF-STEM images [8 (b)
and 8 (d)] represent atomic arrangement of oxygen since oxygen is emphasized
in annular bright-field ABF-STEM images.Okunishi  _et al._ (2009) The film and
the substrate are nearly indistinguishable in the HAADF-STEM image for the
stoichiometric film [Fig. 8 (a)], indicating the ideal cationic arrangement at
the interface. In contrast, the threading dislocations perpendicular to the
film surface are observed in the ABF-STEM image for the off-stoichiometric
film [Fig. 8 (d)], which are merely observed in the HAADF-STEM image as well
[Fig. 8 (c)]. Such threading dislocations have been reported in Sr-rich STO
films and are thought to be Ruddlesden–Popper planar faults.Brooks  _et al._
(2009) In addition, the magnified ABF-STEM image [inset in Fig. 8 (d)] reveals
strong contrast due to local atomic dechanneling.Lee  _et al._ (2016b); Li
_et al._ (2021) Since oxygen is emphasized in ABF-STEM unlike HAADF-STEM
images [Fig. 8 (c)], the contrasts in the ABF-STEM image should come from the
oxygen vacancies. The oxygen vacancies in the off-stoichiometric STO film are
consistent with the growth conditions with an oxidation strength lower than
that for the stoichiometric STO film (O3-nozzle-to-substrate distances are
$13.5$ and $52$ mm for the stoichiometric STO and off-stoichiometric STO
films, respectively).

To determine the Ti/Sr composition ratios, we carried out energy dispersive
x-ray spectroscopy (EDS) measurements taken also with a JEOL JEM-ARM 200F
microscope. Figure 9 shows the EDS spectra for the stoichiometric and off-
stoichiometric STO films. Only Sr, Ti, and O peaks are observed in both films,
indicating no observable impurities in the films. The Ti/Sr composition ratio
was estimated by the Ti $K\alpha$/Sr $L\alpha$ integrated intensity ratios
normalized by that of the STO substrate, i.e., the Ti/Sr ratio of the STO
substrate is assumed to be 1. The estimated Ti/Sr ratios in the stoichiometric
and off-stoichiometric films are 1.00 and 0.94, respectively. Note that
typical accuracy of the EDS for the Sr and Ti $L\alpha$ integrated intensity
is $\pm$0.01–$\pm$0.04.Gao  _et al._ (2018); Alonso  _et al._ (2019)

Figure 10 shows the optical absorption of the stoichiometric and off-
stoichiometric STO films at room temperature. The sudden increase of the
absorption at 3.2 eV originates from the O 2$p$-to-Ti 3$d$ charge transfer
transition in STO films and substrates, indicating the bandgap of STO.Konta
_et al._ (2004) The absorption spectrum of the stoichiometric STO is identical
to that of the STO substrate, indicating that the stoichiometric STO could be
an ideal mother material for photocatalysis applications. In contrast, non-
stoichiometric STO shows the Drude (free-electron-carrier) absorption (photon
energy $<1$ eV) and absorption from deep-impurity states ($1$ eV $<$ photon
energy $<3.2$ eV).Kumar  _et al._ (2020) The absorptions at 2.4 and 2.9 eV may
originate from the excitation of electrons trapped by oxygen vacancies since
it is widely observed in reduced STO.Pennycook and Jesson (1991)

##  IV Conclusion

We demonstrated the stoichiometric growth of STO films via Bayesian
optimization with an adaptive hyperparameter of a prior mean function. To
obtain highly insulating stoichiometric STO films, we developed a new BO
method that encourages the exploration by adjusting the prior mean to get out
of suboptimal parameters. Using simulated data, we found the efficacy of all
the methods that vary the prior mean value, reproducibly reaching the global
best conditions. Among these methods, we employed the AL method for ML-MBE. In
only 44 MBE growth runs, our approach attained highly insulating
stoichiometric STO films having no absorption in the band gap, which will
contribute to the next generation of electronic and photonic devices. The
proposed algorithm provides an efficient experimental design platform that is
not as dependent on the experience and skills of individual researchers. It
will enhance the efficiency of not only oxide electronics but also various
material syntheses and autonomous syntheses.Burger  _et al._ (2020); Ziatdinov
_et al._ (2022); Stach  _et al._ (2021)

## Authors’ Contributions

Y. K. Wakabayashi: Conceptualization (equal); Methodology (equal); Software
(supporting); Validation (equal); Investigation (lead); Supervision (equal);
Writing – Original Draft (equal); Writing – Review & editing (equal). T.
Otsuka: Conceptualization (equal); Methodology (equal); Software (lead);
Validation (equal); Investigation (supporting); Supervision (equal); Writing –
Original Draft (equal); Writing – Review & editing (equal). Y. Krockenberger:
Investigation (supporting); Writing – Review & editing (supporting). H.
Sawada, Y. Taniyasu, H. Yamamoto: Writing – Review & editing (supporting).

## Data Availability

The data and code that support the findings of this study are openly available
in GitHub at https://github.com/nttcslab/adaptive-leveling-BO.

Figure 1:  Flow of ML-MBE growth using BO. (a) Schematic illustration of our
multisource oxide MBE setup. EIES: Electron Impact Emission Spectroscopy. (b)
X-ray diffraction (XRD) $\theta$–$2\theta$ scan of the STO film with a $\Delta
c$ of 0.045 Å, as an example. STO sub. means the peak from the substrate. (c)
Growth conditions for four samples, as an example. (d) Two-dimensional plots
of EI values at the O3-nozzle-to-substrate distance of 65.5 mm obtained from
the collected data for 27 samples, as an example.  Figure 2:  (a) Ackley
function on two-dimensional parameter space [$x$1,$x$2]. Among four peaks, the
left-bottom one at $[0,0]^{\top}$ gives the largest value. (b), (c) Predicted
mean $m({\bf x})$ and (d), (e) acquisition function $a_{\mathrm{EI}}({\bf x})$
with 20 observations depicted with $\times$ marks. In (d) and (e), the $\star$
mark represents the position of maximum of $a_{\mathrm{EI}}({\bf x})$. (b),
(d) Results with $m_{0}=0$. (c), (e) Results with $m_{0}\approx 0.2$.

Figure 3:  Rosenbrock function with $D=2$ employed for simulated optimization
experiment.  Figure 4:  (a)–(c) Optimization results for the Ackley function
with (a) $D=2$, (b) $D=4$, and (c) $D=6$. (d)–(f) Optimization results for the
Rosenbrock function with (d) $D=2$, (e) $D=4$, and (f) $D=6$. Solid lines
indicate the best evaluation value found so far as a function of number of
observations averaged over five trials. Shaded areas means the best and worst
performance among the five trials.  Figure 5:  (a)–(c) Experimental $\Delta c$
values in the three-dimensional growth parameter space for 10 (a), 27 (b), and
44 (c) samples. The green spheres indicate the NaN points at which the STO
phase was not obtained. The red spheres indicate the most promising conditions
with the highest EI values, which should be examined in the next growth run.
The red planes indicate the cutting plane of the O3-nozzle-to-substrate
distance, at which the highest EI value was obtained. (d)–(l): Two-dimensional
plots of predicted $\Delta c$ values (d)–(f), $s$ values (g)–(i), and EI
values (j)–(l) at O3-nozzle-to-substrate distance of 39.5 mm [(d),(g),(j)],
65.5 mm [(e),(h),(k), and 10 mm [(f),(i),(l)], which were obtained from the
collected data for 10 [(d),(g),(j)], 27 [(e),(h),(k)], and 44 [(f),(i),(l)]
observations, respectively. The O3-nozzle-to-substrate distance was that at
which the highest EI value was obtained.  Figure 6:  Actual $\Delta c$ values
and lowest experimental $\Delta c$ plotted as a function of growth number.
Figure 7:  (a) XRD $\theta$–$2\theta$ scans of the stoichiometric STO film
grown on (001) STO and LSAT substrates. Inset shows the XRD $\theta$–$2\theta$
scan of the off-stoichiometric STO film. (b) AFM image of the stoichiometric
STO film on (001) STO. (c) Magnified image of (b).  Figure 8:  (a), (c) HAADF-
STEM and (b), (d) ABF-STEM images of the stoichiometric [(a), (b)] and off-
stoichiometric [(c), (d)] STO films along the [100] direction. Dashed lines
indicate the interfaces between the grown STO layers and the substrates.
Insets show magnified images at the center of the films.  Figure 9:  EDS
spectra for the stoichiometric and off-stoichiometric STO films. The inset
shows a magnified view at the Ti $K\alpha$. The spectra are normalized at the
Sr $L\alpha$ peak intensities for easy comparison.  Figure 10:  Optical
absorptions of the STO substrate and the stoichiometric and non-stoichiometric
STO films on STO at room temperature.

##  Appendix A Bayesian optimization with adaptive prior mean function

Bayesian optimization (BO) is a method for optimizing a black box function
$f({\bf x})$:

$\displaystyle\max_{{\bf x}\in\mathpzc{X}}f({\bf x}),$ (1)

where we can evaluate the function value given a $D$-dimensional parameter
${\bf x}\in\mathpzc{X}\subset\mathbb{R}^{D}$ as:

$\displaystyle y$ $\displaystyle=f({\bf x})+\varepsilon,$ (2)

but the underlying function $f({\bf x})$ is unknown. In materials growth
optimization, ${\bf x}$ and $y$ represent growth parameters and physical
properties used to evaluate grown materials, respectively. Examples of
physical properties include electrical resistance and X-ray diffraction
intensity. Here, we assume observation y contains an additive Gaussian noise
$\varepsilon\sim\mathpzc{N}(0,\sigma_{\varepsilon}^{2})$. While we will
present the maximization problem, note that we can equivalently cope with the
minimization problem by negating the observation. Furthermore, we assume the
search space $\mathpzc{X}$ is bounded so that each element in ${\bf x}$ may be
normalized between 0 and 1, as described in existing work.Wakabayashi  _et
al._ (2022a) In addition, we may combine the floor padding trickWakabayashi
_et al._ (2022a) to handle experimental failures where the observation of Eq.
(2) is unavailable. BO optimizes the black box function by constructing a
prediction model to find an unexplored parameter with a good chance to improve
its function value. Then, the predicted model is iteratively updated after
obtaining the actual observation corresponding to the parameter. This process
is typically iterated until a certain function value is attained or the number
of observations reaches a certain budget for the experiment.

###  A.1 GP prediction and acquisition function

The Gaussian process is adopted for the prediction model. This model predicts
the outcome of unseen parameters as a normal distribution.

$\displaystyle p(y^{\prime}|{\bf x}^{\prime},\mathpzc{D}_{n})$
$\displaystyle=\mathpzc{N}\left(m({\bf x}^{\prime}),s^{2}({\bf
x}^{\prime})\right),$ (3)

where the past $n$ observations are denoted as $\mathpzc{D}_{n}=\left\\{({\bf
x}_{1},y_{1}),\dots,({\bf x}_{n},y_{n})\right\\}$. The predictive mean $m({\bf
x}^{\prime})$ and variance $s^{2}({\bf x}^{\prime})$ are calculated using the
kernel $k(\cdot,\cdot)$ as follows:

$\displaystyle m({\bf x}^{\prime})$ $\displaystyle=\eta({\bf x}^{\prime})+{\bf
k}_{n}^{\top}{\bf K}_{n}^{-1}\left({\bf y}_{n}-\boldsymbol{\eta}_{n}\right),$
(4) $\displaystyle s^{2}({\bf x}^{\prime})$ $\displaystyle=k({\bf
x}^{\prime},{\bf x}^{\prime})-{\bf k}_{n}^{\top}{\bf K}_{n}^{-1}{\bf k}_{n},$
(5)

where the $i$th element of vector ${\bf k}_{n}\in\mathbb{R}^{n}$ is given by
$k({\bf x}_{i},{\bf x}^{\prime})$ and the element at the $i$th row and $j$th
column of Gram matrix ${\bf K}_{n}\in\mathbb{R}^{n\times n}$ is $k({\bf
x}_{i},{\bf x}_{j})$. We used the Matérn $\frac{5}{2}$ kernel for our GP.
Vector ${\bf y}_{n}$ represents the past observations; ${\bf
y}_{n}=\left[y_{1},\dots,y_{n}\right]^{\top}$ with $\cdot^{\top}$ means the
transpose operator. Function $\eta({\bf x})$ is the prior mean function that
shifts the predictive mean by subtracting
$\boldsymbol{\eta}_{n}=\left[\eta({\bf x}_{1}),\dots,\eta({\bf
x}_{n})\right]^{\top}$ from observation ${\bf y}_{n}$ and adding $\eta({\bf
x}^{\prime})$ back to predictive mean $m({\bf x}^{\prime})$ Intuitively, when
parameter ${\bf x}^{\prime}$ is distant from any of observed data ${\bf
x}_{1},\dots,{\bf x}_{n}$, the predictive mean is dominated by the prior mean
function, i.e., $m({\bf x}^{\prime})\approx\eta({\bf x}^{\prime})$ because
elements of ${\bf k}_{n}$ approach zero with distance-based kernels, such as
the Matérn or radial basis function kernels.

Given the predictions, the acquisition function is evaluated to decide which
parameter to employ in the next trial. We use the expected improvement (EI)
criterion

$\displaystyle a_{\mathrm{EI}}$
$\displaystyle=\mathbb{E}_{p\left(y^{\prime}|{\bf
x}^{\prime},\mathpzc{D}_{n}\right)}\left[\left(y^{\prime}-\bar{y}_{n}\right)\mathbb{I}_{y^{\prime}\geq\bar{y}_{n}}\right],$
(6)

where $\displaystyle\bar{y}_{n}=\max_{1\leq i\leq n}y_{i}$ is the largest
observation and $\mathbb{I}_{y^{\prime}\geq\bar{y}_{n}}$ is the indicator
function that equals 1 when $y^{\prime}\geq\bar{y}_{n}$ and 0 otherwise. This
acquisition function evaluates the expectation of improvement over the best
observation $\bar{y}_{n}$ at unseen parameter ${\bf x}^{\prime}$.

###  A.2 Impact of prior mean function on optimization

The choice of $\eta({\bf x})$ will affect the acquisition function as well as
the parameter search efficiency. The typical choice of the prior mean function
is $\eta({\bf x})=0$, particularly because we lack the knowledge of underlying
black box function $f({\bf x})$. We replace the prior mean with a constant
function in the form of $\eta({\bf x})=m_{0}$ since a more flexible functional
form of $\eta$ to potentially approximate $f$ is inaccessible.

Figure 2 shows the difference in the prediction mean and the acquisition for
the two-dimensional Ackley function with $m_{0}=0$ and $m_{0}\approx 0.2$. See
Appendix B and Eq. (13) for the description of the Ackley function. As shown
in Fig. 2 (a), the Ackley function $f({\bf x})$ has four peaks with the left-
bottom one at ${\bf x}=[0,0]^{\top}$ being the highest. In this example, the
left-top suboptimal peak has been intensively searched [Figs. 2 (b) and 2
(c)]. Comparison with $m_{0}=0$ [Fig. 2 (b)] and $m_{0}\approx 0.2$ [Fig. 2
(c)] shows that the predictive mean $m({\bf x}^{\prime})$ with $m_{0}\approx
0.2$ is larger than that with $m_{0}=0$ in the half bottom region. This
difference in the predicted mean results in the difference in the acquisition
function [Figs. 2 (d) and 2 (e)]. The use of $m_{0}>0$ increases the
predictive mean of unexplored region and leads to finding other peaks.

While a greater $m_{0}$ worked favorably in the example in Fig. 2, large
$m_{0}$ may not always be efficient: a large $m_{0}$ tends to produce an
optimistic prediction in unexplored regions. This can trigger unnecessary
explorations leading to a plateau at the suboptimal function value. Thus, the
choice of hyperparameter $m_{0}$ needs to take account of the balance between
the exploration and exploitation in the parameter search.

###  A.3 Methods for adapting prior mean function hyperparameter

The discussion in Section A.2 motivates us to adapt $m_{0}$. A simple way is
to average observed data $m_{0}=\frac{1}{n}\sum_{i=1}^{n}y_{i}$. We call this
approach DA (data averaging). Despite its simplicity, DA works well in
practice [Fig. 4], but tends to saturate $m_{0}$ to a certain value as more
observations are accumulated. This can compensate the exploration of parameter
search. In the following, we present three methods with varying $m_{0}$ for
efficient BO.

The first method is adaptive leveling (AL) that randomly sets $m_{0}$ between
the best and worst observations every time a new observation is acquired:

$\displaystyle m_{0}\sim\mathpzc{U}\left(\min_{1\leq i\leq n}y_{i},\max_{1\leq
i\leq n}y_{i}\right),$ (7)

where $\mathpzc{U}(l,h)$ denotes the uniform distribution on interval $[l,h]$.
We expect that $m_{0}$ is not always too optimistic or pessimistic by the
stochastic choice.

While the design of AL is intuitive, this may be unnatural in light of
Bayesian inference: $m_{0}$, the hyperparameter of the prior mean, is
specified after observing data. The empirical Bayes methodCasella (1985);
Carlin and Louis (2000) can justify the choice of $m_{0}$ based on
observations ${\bf y}_{n}$ by maximizing the marginal likelihood. The marginal
log likelihood, or the evidence, of the Gaussian process for data
$\mathpzc{D}_{n}$ is given as follows:Rasmussen and Williams (2006)

$\displaystyle\ln p({\bf y}_{n}|{\bf x}_{1},\dots,{\bf x}_{n})$
$\displaystyle=-\frac{1}{2}\left({\bf y}_{n}-{\bf
1}_{n}m_{0}\right)^{\top}{\bf K}_{n}^{-1}\left({\bf y}_{n}-{\bf
1}_{n}m_{0}\right)+\mathrm{const}.$ (8)

Here, ${\bf 1}_{n}$ is an $n$-dimensional vector with all elements being 1 and
the constant term in Eq. (8) includes quantities irrelevant to $m_{0}$. The
marginal likelihood is maximized when

$\displaystyle m_{0}$ $\displaystyle=\frac{{\bf 1}_{n}^{\top}{\bf
K}_{n}^{-1}{\bf y}_{n}}{{\bf 1}_{n}^{\top}{\bf K}_{n}^{-1}{\bf 1}_{n}}.$ (9)

This quantity can be interpreted as a weighted average of elements in ${\bf
y}_{n}$, where the weight is given by ${\bf K}_{n}^{-1}{\bf 1}_{n}$. We refer
to this choice of $m_{0}$ as EB that stands for empirical Bayes. We may
optionally add some randomness to $m_{0}$ by considering the condition when
the marginal log likelihood improves over the baseline configuration
$m_{0}=0$. The empirical Bayes uniform (EBu) draws $m_{0}$ as follows:

$\displaystyle m_{0}$
$\displaystyle\sim\left\\{\begin{array}[]{lr}\mathpzc{U}\left(0,2\frac{{\bf
1}_{n}^{\top}{\bf K}_{n}^{-1}{\bf y}_{n}}{{\bf 1}_{n}^{\top}{\bf
K}_{n}^{-1}{\bf 1}_{n}}\right)&\mathrm{if}\;\;\frac{{\bf 1}_{n}^{\top}{\bf
K}_{n}^{-1}{\bf y}_{n}}{{\bf 1}_{n}^{\top}{\bf K}_{n}^{-1}{\bf 1}_{n}}\geq
0,\\\ \mathpzc{U}\left(2\frac{{\bf 1}_{n}^{\top}{\bf K}_{n}^{-1}{\bf
y}_{n}}{{\bf 1}_{n}^{\top}{\bf K}_{n}^{-1}{\bf
1}_{n}},0\right)&\mathrm{otherwise}.\end{array}\right.$ (12)

EBu has a mechanism similar to AL for switching between the exploration and
exploitation, as well as a guarantee for a better marginal likelihood compared
with the baseline $m_{0}=0$. In contrast, AL may use a value of $m_{0}$ with
less marginal likelihood; its range is restricted to the range of past
observations. This can circumvent an extreme choice of $m_{0}$ that may cause
unstable search.

Indeed, EB and EBu may choose $m_{0}\in\left[\min_{1\leq i\leq
n}y_{i},\max_{1\leq i\leq n}y_{i}\right]$. This is because some elements in
${\bf K}_{n}^{-1}{\bf 1}_{n}$ can be negative, and thus their weighted average
can produce the extrapolation of observations in ${\bf y}_{n}$. Since EBu uses
an extended interval between 0 and $\displaystyle 2\frac{{\bf
1}_{n}^{\top}{\bf K}_{n}^{-1}{\bf y}_{n}}{{\bf 1}_{n}^{\top}{\bf
K}_{n}^{-1}{\bf 1}_{n}}$, $m_{0}$ of EBu has more chance to be outside of the
past observations.

##  Appendix B Objective functions used in simulated data experiments

The Ackley function is defined as follows:

$\displaystyle f_{\mathrm{A}}({\bf x})$ $\displaystyle\propto
C_{1}\left(\exp\left(-\sqrt{\frac{\sum_{d=1}^{D}x_{d}^{2}}{25D}}\right)+\exp\left(\cos\left(4\frac{\sum_{d=1}^{D}x_{d}}{D}\right)\right)\right)+C_{2},$
(13)

where $x_{d}$ denotes the $d$th element of ${\bf x}$ and $C_{1}$ and $C_{2}$
are set such that $f_{\mathrm{A}}({\bf x})$ ranges from -0.5 to 0.5 on
$\mathpzc{X}$. This function has $2^{D}$ peaks in $\mathpzc{X}$ and takes its
maximum value at ${\bf x}={\bf 0}$.

The Rosenbrock function is defined as follows:

$\displaystyle f_{\mathrm{R}}$ $\displaystyle\propto-
C_{3}\left(\sum_{d=1}^{D}(x_{d}-1)^{2}+5\sum_{d=1}^{D-1}\left(x_{d}-x_{d+1}^{2}\right)^{2}\right)+C_{4},$
(14)

where $C_{3}$ and $C_{4}$ are configured such that $f_{\mathrm{R}}({\bf x})$
ranges from -0.5 to 0.5 on $\mathpzc{X}$. This function has a valley-shaped
surface with its unique peak located at ${\bf x}={\bf 1}$.

##  Appendix C Optimization of deviation in lattice constant

In the STO film growth experiment, our objective was to minimize $\Delta c$,
the absolute difference of lattice constants. While $\Delta c\geq 0$ by
definition, the prediction of GP assumes $\Delta c$ may take a negative value
by fitting the normal distribution. We may consider two types of fix for this
mismatch. First one is to handle the logarithm of difference: $y=-\log\Delta
c$. Note that we negated the difference since BO is presented as a
maximization of $y$. While this naturally maps the nonnegative value to a real
value, a large difference $\Delta c$ is distorted and less emphasized. The
second approach is to truncate the probability of $\Delta c<0$, that is,
$p(y>0|{\bf x},\mathpzc{D}_{n})=0$, where $y=-\Delta c$. We can modify the EI
criterion to derive an analytic form of

$\displaystyle a_{\mathrm{EI}}({\bf x})$
$\displaystyle=\mathbb{E}_{\tilde{p}\left(y^{\prime}|{\bf
x}^{\prime},\mathpzc{D}_{n}\right)}\left[\left(y^{\prime}-\bar{y}_{n}\right)\mathbb{I}_{y^{\prime}\geq\bar{y}_{n}}\right],$
(15)

where the modified predictive probability is

$\displaystyle\tilde{p}(y|{\bf x},\mathpzc{D}_{n})$
$\displaystyle\propto\left\\{\begin{array}[]{cc}p(y|{\bf
x},\mathpzc{D}_{n})&\mathrm{if}\;\;y\leq 0,\\\
0&\mathrm{otherwise}.\end{array}\right.$ (18)

We carried out an experiment for minimizing the difference between the
function value and its target. We compared four methods: three of them
minimizes the absolute difference between the target and function value, while
the other uses a special acquisition function that evaluates the closeness of
the prediction to the target value. The first method is the most
straightforward one; it simply handles the absolute difference between the
target and function value: $y=-|t-f({\bf x})|$, where $t$ denotes the target
value. We call this method absolute here. The second method (log absolute)
maximizes $y=-\log|t-f({\bf x})|$ to let $y$ take a negative value. The third
method (truncation) modifies the acquisition function as in Eq. (15) using
$y=-|t-f({\bf x})|$. The last method Urenholt and Jensen (2019) fits the GP to
the function directly using $y=f({\bf x})$ to form the predictive mean $m({\bf
x}^{\prime})$ and variance $s^{2}({\bf x}^{\prime})$ as in Eq. (3). Letting
the squared distance denote $d({\bf x}^{\prime})=\left(t-f({\bf
x}^{\prime})\right)^{2}$, we model the distance using the non-central chi-
squared distribution $\chi_{\mathrm{NC}}^{2}$ as

$\displaystyle\left.\frac{d}{s^{2}({\bf x}^{\prime})}\right|{\bf
x}^{\prime},\mathpzc{D}_{n}$
$\displaystyle\sim\chi_{\mathrm{NC}}^{2}\left(\lambda({\bf
x}^{\prime})\right),$ (19)

where the degree of freedom of $\chi_{\mathrm{NC}}^{2}$ is 1 and its non-
centrality is $\lambda({\bf x}^{\prime})=\left(\frac{t-m({\bf
x}^{\prime})}{s({\bf x}^{\prime})}\right)^{2}$. The acquisition function
evaluates how $d$ can be close to zero by, for example, assessing the distance
at predetermined percentile of the distribution in Eq. (19). We refer to the
last method as chi2.

Our experiments used the Ackley and Rosenbrock functions with $D=2$ and set
the target as $t=0$. Namely, we minimized the absolute value of these
functions. Figure 11 (a) and (b) show the value $|f({\bf x})-t|$ of Ackley
function and Rosenbrock function, respectively. Figure 11 (c) and (d) show the
optimization results of each function. For both functions, the optimization
process with each method was repeated five times until 100 observations
acquired. Solid lines indicate the smallest difference between the function
value and target as a function of the number of observations averaged over the
five runs. Shaded areas mean the best and worst performance among the five
runs. The results in Figure 11 (c) and (d) were given by the AL method for GP
prediction. In the case of the Ackley function [Fig. 11 (c)], all approaches
eventually reached $f({\bf x})=0$ after 40 observations. In contrast, the
truncation approach failed once out of five trials [Fig. 11 (d)]. The rest of
the methods showed stable performance for finding $f({\bf x})=0$. Figure 11
(e) and (f) investigate the effect of the choice of $m_{0}$ (AL vs. baseline)
on the chi2 approach. While AL made little difference for the Rosenbrock
function [Fig. 11 (f)], AL accelerated the optimization in the case of the
Ackley function [Fig. 11 (e)]. The observation in an unexplored region was
encouraged when $m_{0}\approx 0$ was chosen. This behavior enhanced the chance
of finding $f({\bf x})=0$ on the wavy surface of the absolute Ackley function
[Fig. 11 (a)].

Our results showed that all the approaches except the truncation one gave
comparable performance for both functions. While chi2 reduced the difference
between the function and target values slightly quicker than the other
approaches, we adopted the simplest absolute approach for the STO film growth
experiment in expectation of its robust behavior and easier interpretation of
progress in the wild environment. Since the MBE growth experiment requires
extensive time, we leave the use of chi2 for optimizing the growth parameters
as future work.

Figure 11:  (a) Absolute value of Ackley function with $D=2$. The original
function is displayed in Fig. 2 (a). (b) Absolute value of Rosenbrock function
with $D=2$. The original function is displayed in Fig. 3. (c) Optimization
results of minimizing the difference between the function value and $0$, i.e.,
the absolute value of the Ackley function. (d) Optimization results of
minimizing the absolute value of the Rosenbrock function. (c), (d) Solid lines
indicate the best evaluation value found so far as a function of the number of
observations averaged over five trials. Shaded areas mean the best and worst
performance for the five runs. (e) Optimization results of absolute Ackley
function using $\chi_{\mathrm{NC}}^{2}$ with $m_{0}=0$ (base) and AL. (f)
Optimization results of absolute Rosenbrock function using
$\chi_{\mathrm{NC}}^{2}$ with $m_{0}=0$ (base) and AL.

## References

  * Shi  _et al._ (2020) X. L. Shi, H. Wua, Q. Liu, W. Zhousa, S. Lu, Z. Shao, M. Dargusch, and Z. G. Chen, “SrTiO3-based thermoelectrics: Progress and challenges,” Nano Energy 78, 105195 (2020).
  * Pai  _et al._ (2018) Y. Y. Pai, A. T. Yyler, P. Irvin, and J. Levy, “Physics of SrTiO3-based heterostructures and nanostructures: A review,” Rep. Prog. Phys. 81, 036503 (2018).
  * Fujimoto and Kingery (1985) M. Fujimoto and W. D. Kingery, “Microstructures of SrTiO3 internal boundary layer capacitors during and after processing and resultant electrical properties,” J. Am. Ceram. Soc. 68, 169 (1985).
  * Sakuma  _et al._ (1998) T. Sakuma, S. Yamamichi, S. Matsubara, H. Yamaguchi, and Y. Miyasaka, “Barrier layers for realization of high capacitance density in SrTiO3 thin-film capacitor on silicon,” Appl. Phys. Lett. 57, 2431 (1998).
  * Galt  _et al._ (1993) D. Galt, J. C. Price, J. A. Beall, and R. H. Ono, “Characterizatoin of a tunable thin film microwave YBa2Cu3O7-x/SrTiO3 coplanar capacitor,” Appl. Phys. Lett. 63, 3078 (1993).
  * Mavroides, Kafalas, and Kolesar (1976) J. G. Mavroides, J. A. Kafalas, and D. F. Kolesar, “Photoelectrolysis of water in cells with SrTiO3 anodes,” Appl. Phys. Lett. 28, 241 (1976).
  * Konta  _et al._ (2004) R. Konta, T. Ishii, H. Kato, and A. Kudo, “Photocatalytic activities of noble metal ion doped SrTiO3 under visible light irradiation,” J. Phys. Chem. B 108, 8992 (2004).
  * Ng  _et al._ (2010) B. J. Ng, S. Xu, X. Zhang, H. Y. Yang, and D. D. Sun, “Hybridized nanowires and cubes : A novel architecture of a heterojunctioned TiO2/SrTiO3 thin film for efficient water splitting,” Adv. Funct. Mater. 20, 4287 (2010).
  * Lee  _et al._ (2008) S. W. Lee, O. S. Kwon, J. H. Han, and C. S. Hwang, “Enhanced electrical properties of thin films grown by atomic layer deposition at high temperature for dynamic random access memory applications,” Appl. Phys. Lett. 92, 222903 (2008).
  * Takata  _et al._ (2020) T. Takata, J. Jiang, Y. Sakata, M. Nakabayashi, N. Shibata, V. Nandal, K. Seki, T. Hisatomi, and K. Domen, “Photocatalytic water splitting with a quantum efficiency of almost unity,” Nature 581, 411 (2020).
  * Ohtomo and Hwang (2004) A. Ohtomo and H. Y. Hwang, “A high-mobility electron gas at the LaAlO3/SrTiO3 heterointerface,” Nature 427, 423 (2004).
  * Bowen  _et al._ (2008) M. Bowen, M. Bibes, A. Barthelemy, J. P. Contour, A. Anane, Y. Lemaitre, and A. Fert, “Nearly total spin polarization in La2/3Sr1/3MnO3 from tunneling experiments,” Appl. Phys. Lett. 82, 233 (2008).
  * Wakabayashi  _et al._ (2021a) Y. K. Wakabayashi, S. K. Takada, Y. Krockenberger, Y. Taniyasu, and H. Yamamoto, “Wide-range epitaxial strain control of electrical and magnetic properties in high-quality SrRuO3 films,” ACS Appl. Electron. Mater. 3, 2712 (2021a).
  * Takada  _et al._ (2021) S. K. Takada, Y. K. Wakabayashi, Y. Krockenberger, S. Ohya, M. Tanaka, Y. Taniyasu, and H. Yamamoto, “Thickness-dependent quantum transport of Weyl fermions in ultra-high-quality SrRuO3 films,” Appl. Phys. Lett. 118, 092408 (2021).
  * Wakabayashi  _et al._ (2021b) Y. K. Wakabayashi, M. Kobayashi, Y. Takeda, K. Takiguchi, H. Irie, S. Fujimori, T. Takeda, R. Okano, Y. Krockenberger, Y. Taniyasu, and H. Yamamoto, “Single-domain perpendicular magnetization induced by the coherent O $2p$-Ru $4d$ hybridized state in an ultra-high-quality SrRuO3 film,” Phys. Rev. Mater. 5, 124403, 1 (2021b).
  * Takada  _et al._ (2022) S. K. Takada, Y. K. Wakabayashi, Y. Krockenberger, T. Nomura, Y. Kohama, H. Irie, K. Takiguchi, S. Ohya, M. Tanaka, Y. Taniyasu, and H. Yamamoto, “High-mobility two-dimensional carriers from surface Fermi arcs in magnetic Weyl semimetal films,” npj Quantum Mater. 7, 102 (2022).
  * Koonce  _et al._ (1967) C. S. Koonce, M. L. Cohen, J. F. Schooley, W. R. Hosler, and E. R. Pfeiffer, “Superconducting transition temperatures of semiconducting SrTiO3,” Phys. Rev. 163, 380 (1967).
  * Bert  _et al._ (2011) J. A. Bert, B. Kalisky, C. Bell, M. Kim, Y. Hikita, H. Y. Hwang, and K. A. Moler, “Direct imaging of the coexistence of ferromagnetism and superconductivity at the LaAlO3/SrTiO3 interface,” Nat. Phys. 7, 767 (2011).
  * Li  _et al._ (2011) L. Li, C. Richter, J. Mannhart, and R. C. Ashoori, “Coexistence of magnetic order and two-dimensional superconductivity at LaAlO3/SrTiO3 interfaces,” Nat. Phys. 7, 762 (2011).
  * Ahadi  _et al._ (2019) K. Ahadi, L. Galletti, Y. Li, S. S. Rezaie, W. Wu, and S. Stemmer, “Enhancing superconductivity in SrTiO3 films with strain,” Sci. Adv. 5, eeaw0120 (2019).
  * Jang  _et al._ (2010) H. W. Jang, A. Kumar, S. Denev, M. D. Biegalski, P. Maksymovych, C. W. Bark, C. T. Nelson, C. M. Folkman, S. H. Baek, N. Balke, C. M. Brooks, D. A. Tenne, D. G. Schlom, L. Q. Chen, X. Q. Pan, S. V. Kalinin, V. Gopalan, and C. B. Eom, “Ferroelectricity in strain-free SrTiO3 thin films,” Phys. Rev. Lett. 104, 197601 (2010).
  * Son  _et al._ (2010) J. Son, P. Moetakef, B. Jalan, O. Bierwagen, N. J. Wright, R. E. Herbert, and S. Stemmer, “Epitaxial SrTiO3 films with electron mobilities exceeding 30,000 cm2v-1s-1,” Nat. Mater. 9, 482 (2010).
  * Matsubara  _et al._ (2016) Y. Matsubara, K. S. Takahashi, M. S. Bahramy, Y. Kozuka, D. Maryenko, J. Falson, A. Tsukazaki, Y. Tokura, and M. Kawasaki, “Observation of the quantum hall effect in $\delta$-doped SrTiO3,” Nat. Commun. 7, 11631 (2016).
  * Kan  _et al._ (2005) D. Kan, T. Terashima, R. Kanda, A. Masuno, K. Tanaka, S. Chu, H. Kan, A. Ishizumi, Y. Kanemitsu, Y. Shimakawa, and M. Takano, “Blue-light emission at room temperature from Ar+-irradiated SrTiO3,” Nat. Mater. 4, 816 (2005).
  * Popescu  _et al._ (2014) D. Popescu, B. Popescu, G. Jegert, S. Schmelzer, U. Boettger, and P. Lugli, “Feasibility study of SrRuO3/SrTiO3/SrRuO3 thin film capacitors in DRAM applications,” IEEE Trans. Electron Devices 61, 2130 (2014).
  * Baek  _et al._ (2020) J. Baek, L. Thai, S. Yeon, and H. Seo, “Aluminum doping for optimization of ultrathin and high-$k$ dielectric layer based on SrTiO3,” J. Mater. Sci. Technol. 42, 28 (2020).
  * Chen  _et al._ (2018) H. Chen, F. Zhang, W. Zhang, Y. Du, and G. Li, “Negative impact of surface Ti3+ defects on the photocatalytic hydrogen evolution activity of SrTiO3,” Appl. Phys. Lett. 112, 013901 (2018).
  * Kumar  _et al._ (2020) M. Kumar, P. Basera, S. Saini, and S. Bhattacharya, “Role of defects in photocatalytic water splitting: monodoped vs codoped SrTiO3,” J. Phys. Chem. C 124, 10272 (2020).
  * Lee  _et al._ (2016a) H. N. Lee, S. S. A. Seo, W. S. Choi, and C. M. Rouleau, “Growth control of oxygen stoichiometry in homoepitaxial SrTiO3 films by pulsed laser epitaxy in high vacuum,” Sci. Rep. 6, 19941 (2016a).
  * Flores  _et al._ (2017) A. M. H. Flores, J. Chen, L. M. T. Martínez, A. Ito, and E. Moctezuma, “Laser assisted chemical vapor deposition of nanostructured NaTaO3 and SrTiO3 thin films for efficient photocatalytic hydrogen evolution,” Fuel 197, 174 (2017).
  * Mueller, Kusne, and Ramprasad (2015) T. Mueller, A. G. Kusne, and R. Ramprasad, “Machine learning in materials science,” in  _Reviews in Computational Chemistry_ , Vol. 29 (Wiley, Hoboken, 2015).
  * Lookman, Alexander, and Rajan (2016) T. Lookman, F. J. Alexander, and K. Rajan,  _Information Science for Materials Discovery and Design_ (Springer, Cham, 2016).
  * Burnaex and Panov (2015) E. Burnaex and M. Panov,  _Statistical Learning and Data Sciences_ (Springer, Cham, 2015).
  * Agrawal and Choudhary (2016) A. Agrawal and A. N. Choudhary, “Materials informatics and big data: Realization of the “fourth paradigm” of science in materials science,” APL Mater. 4, 053208 (2016).
  * Ueno  _et al._ (2018) T. Ueno, H. Hino, A. Hashimoto, Y. Takeichi, M. Sawada, and K. Ono, “Adaptive design of an X-ray magnetic circular dichroism spectroscopy experiment with Gaussian process modelling,” npj Comput. Mater. 4, 4 (2018).
  * Ren  _et al._ (2018) F. Ren, L. Ward, T. Williams, K. Laws, C. Wolverton, J. Hattrick-Simpers, and A. Mehta, “Accelerated discovery of metallic glasses through iteration of machine learning and high-throughput experiments,” Sci. Adv. 4, eaaq1556 (2018).
  * Wakabayashi  _et al._ (2018) Y. K. Wakabayashi, T. Otsuka, Y. Taniyasu, H. Yamamoto, and H. Sawada, “Improved adaptive sampling method utilizing Gaussian process regression for prediction of spectral peak structures,” Appl. Phys. Express 11, 112401 (2018).
  * Li  _et al._ (2018) X. Li, Z. Hou, S. Gao, Y. Zeng, J. Ao, Z. Zhou, B. Da, W. Liu, Y. Sun, and Y. Zhang, “Efficient optimization of the performance of Mn2+-doped kesterite solar cell: Machine learning aided synthesis of high efficient Cu2(Mn,Zn)Sn(S,Se)4 solar cells,” Sol. RRL 2, 1800198 (2018).
  * Hou  _et al._ (2019) Z. Hou, Y. Takagiwa, Y. Shinohara, Y. Xu, and K. Tsuda, “Machine-learning-assisted development and theoretical consideration for the Al2Fe3Si3 thermoelectric material,” ACS Appl. Mater. Interfaces 11, 11545 (2019).
  * Xue  _et al._ (2016) D. Xue, P. V. Balachandran, R. Yuan, T. Hu, X. Qian, E. R. Dougherty, and T. Lookman, “Accelerated search for BaTiO3-based piezoelectrics with vertical morphotropic phase boundary using Bayesian learning,” Proc. Natl. Acad. Sci. U. S. A. 113, 11301 (2016).
  * G.Baird, Liu, and Sparks (2022) S. G.Baird, M. Liu, and T. D. Sparks, “High-dimensional Bayesian optimization of 23 hyperparameters over 100 iterations for an attention-based network to predict materials property: A case study on CrabNet using Axplatform and SAASBO,” Comput. Mater. Sci. 211, 111505 (2022).
  * Snoek, Larochelle, and Adams (2012) J. Snoek, H. Larochelle, and R. P. Adams, “Practical Bayesian optimization of machine learning algorithms,” Advances in Neural Information Processing Systems 25 (2012).
  * Wakabayashi  _et al._ (2019a) Y. K. Wakabayashi, T. Otsuka, Y. Krockenberger, H. Sawada, Y. Taniyasu, and H. Yamamoto, “Machine-learning-assisted thin-film growth: Bayesian optimization in molecular beam epitaxy of SrRuO3 thin films,” APL Mater. 7, 101114 (2019a).
  * Takiguchi  _et al._ (2020) K. Takiguchi, Y. K. Wakabayashi, H. Irie, Y. Krockenberger, T. Otsuka, H. Sawada, S. A. Nikolaev, H. Das, M. Tanaka, Y. Taniyasu, and H. Yamamoto, “Quantum transport evidence of Weyl fermions in an epitaxial ferromagnetic oxide,” Nat. Commun. 11, 4969 (2020).
  * Shimizu  _et al._ (2020) R. Shimizu, S. Kobayashi, Y. Watanabe, Y. Ando, and T. Hitosugi, “Autonomous materials synthesis by machine learning and robotics,” APL Mater. 8, 111110 (2020).
  * Wakabayashi  _et al._ (2022a) Y. K. Wakabayashi, T. Otsuka, Y. Krockenberger, H. Sawada, Y. Taniyasu, and H. Yamamoto, “Bayesian optimization with experimental failure for high-throughput materials growth,” npj Comput. Mater. 8, 180 (2022a).
  * Rasmussen and Williams (2006) C. E. Rasmussen and C. K. I. Williams,  _Gaussian Processes for Machine Learning_ (MIT Press, 2006).
  * Shahriari  _et al._ (2016) B. Shahriari, K. Swersky, Z. Wang, R. P. Adams, and N. de Freitas, “Taking the human out of the loop: A review of Bayesian optimization,” Proc. of the IEEE 104, 148 (2016).
  * Frazier (2018) P. I. Frazier, “A tutorial on Bayesian optimization,” arXiv preprint , 1807.02811 (2018).
  * Naito and Sato (1995) M. Naito and M. H. Sato, “Stoichiometry control of atomic beam fluxes by precipitated impurity phase detection in growth of (Pr,Ce)2CuO4 and (La,Sr)2CuO4 films,” Appl. Phys. Lett. 67, 2557 (1995).
  * Yamamoto, Krockenberger, and Naito (2013) H. Yamamoto, Y. Krockenberger, and M. Naito, “Multi-source MBE with high-precision rate control system as a synthesis method sui generis for multi-cation metal oxides,” J. Cryst. Growth 378, 184 (2013).
  * Wakabayashi  _et al._ (2022b) Y. K. Wakabayashi, Y. Krockenberger, T. Otsuka, H. Sawada, Y. Taniyasu, and H. Yamamoto, “Intrinsic physics in magnetic Weyl semimetal SrRuO3 films addressed by machine-learning-assisted molecular beam epitaxy,” Jap. J. Appl. Phys. (2022b).
  * Wakabayashi  _et al._ (2019b) Y. K. Wakabayashi, Y. Krockenberger, N. Tsujimoto, T. Boykin, S. Tsuneyuki, Y. Taniyasu, and H. Yamamoto, “Ferromagnetism above 1000 K in a highly cation-ordered double-perovskite insulator Sr3OsO6,” Nat. Commun. 10, 535 (2019b).
  * Wakabayashi  _et al._ (2021c) Y. K. Wakabayashi, S. K. Takada, Y. Krockenberger, K. Takiguchi, S. Ohya, M. Tanaka, Y. Taniyasu, and H. Yamamoto, “Structural and transport properties of highly ru-deficient SrRu0.7O3 thin films prepared by molecular beam epitaxy: Comparison with stoichiometric SrRuO3,” AIP Adv. 11, 035226 (2021c).
  * Naito, Yamamoto, and Sato (1998) M. Naito, H. Yamamoto, and H. Sato, “Reflection high-energy electron diffraction and atomic force microscopy studies on homoepitaxial growth of SrTiO (001),” Physica C 305, 233 (1998).
  * Jalan  _et al._ (2009) B. Jalan, R. Engel-Herbert, N. Wright, and S. Stemmer, “Growth of high-quality SrTiO3 films using a hybrid molecular beam epitaxy approach,” J. Vac. Sci. Technol. 27, 461 (2009).
  * Ohnishi  _et al._ (2005) T. Ohnishi, M. Lippmaa, T. Yamamoto, S. Meguro, and H. Koinuma, “Improved stoichiometry and misfit control in perovskite thin film formation at a critical fluence by pulsed laser deposition,” Appl. Phys. Lett. 87, 241919 (2005).
  * Ohnishi  _et al._ (2008) T. Ohnishi, K. Shibuya, T. Yamamoto, and M. Lippmaa, “Defects and transport in complex oxide thin films,” J. Appl. Phys. 103, 103703 (2008).
  * Lebeau  _et al._ (2009) J. M. Lebeau, R. E. Herbert, B. Jalan, J. Cagnon, P. Moetakef, S. Stemmer, and G. B. Stephenson, “Stoichiometry optimization of homoepitaxial oxide thin films using X-ray diffraction,” Appl. Phys. Lett. 95, 142905 (2009).
  * Oh  _et al._ (2011) D. W. Oh, J. Ravichandran, C. W. Liang, W. Siemons, B. Jalan, C. M. Brooks, M. Huijben, D. G. Schlom, S. Stemmer, L. W. Martin, A. Majumdar, R. Ramesh, and D. G. Cahill, “Thermal conductivity as a metric for the crystalline quality of SrTiO3 epitaxial layers,” Appl. Phys. Lett. 98, 221904 (2011).
  * Brooks  _et al._ (2015) C. M. Brooks, R. B. Wilson, A. Schäfer, J. A. Mundy, M. E. Holtz, D. A. Muller, J. Schubert, D. G. Cahill, and D. G. Schlom, “Tuning thermal conductivity in homoepitaxial SrTiO3 films via defects,” Appl. Phys. Lett. 107, 051902 (2015).
  * Lee  _et al._ (2013) C. H. Lee, N. J. Podraza, Y. Zhu, R. F. Berger, S. Shen, M. Sestak, R. W. Collins, L. F. Kourkoutis, J. A. Mundy, H. Wang, Q. Mao, X. Xi, L. J. Brillson, J. B. Neaton, D. A. Muller, and D. G. Schlom, “Effect of reduced dimensionality on the optical band gap of SrTiO3,” Appl. Phys. Lett. 102, 122901 (2013).
  * Močkus, Tiesis, and Žilinskas (1978) J. Močkus, V. Tiesis, and A. Žilinskas, “The application of Bayesian methods for seeking the extremum,” Towards Glob. Optim. 2, 117 (1978).
  * Ackley (1987) D. H. Ackley,  _A Connectionist Machine for Genetic Hillclimbing_ (Kluwer Academic Publishers, 1987).
  * Molga and Smutnicki (2005) M. Molga and C. Smutnicki, “Test functions for optimization needs,” (2005), https://marksmannet.com/RobertMarks/Classes/ENGR5358/Papers/functions.pdf.
  * Rosenbrock (1960) H. H. Rosenbrock, “An automatic method for finding the greatest or least value of a function,” The Comput. J. 3, 175 (1960).
  * Guguschev  _et al._ (2014) C. Guguschev, D. Klimm, F. Langhans, Z. Galazka, D. Kok, U. Juda, and R. Uecker, “Top-seeded solution growth of SrTiO3 crystals and phase diagram studies in the SrO-TiO2 system,” CrystEngComm 16, 1735–1740 (2014).
  * Brooks  _et al._ (2009) C. M. Brooks, L. F. Kourkoutis, T. Heeg, J. Schubert, D. A. Muller, and D. G. Schlom, “Growth of homoepitaxial SrTiO3 thin films by molecular-beam epitaxy,” Appl. Phys. Lett. 94, 162905 (2009).
  * Pennycook and Jesson (1991) S. J. Pennycook and D. E. Jesson, “High-resolution Z-contrast imaging of crystals,” Ultramicroscopy 37, 14 (1991).
  * Okunishi  _et al._ (2009) E. Okunishi, I. Ishikawa, H. Sawada, F. Hosokawa, M. Hori, and Y. Kondo, “Visualization of light elements at ultrahigh resolution by stem annular bright field microscopy,” Microscopy and Microanalysis 15, 164–165 (2009).
  * Lee  _et al._ (2016b) S. A. Lee, H. Jeong, S. Woo, J. Y. Hwang, S. Y. Choi, S. D. Kim, M. Choi, S. Roh, H. Yu, J. Hwang, S. W. Kim, and W. S. Choi, “Phase transitions via selective elemental vacancy engineering in complex oxide thin films,” Sci. Rep. 6, 23649 (2016b).
  * Li  _et al._ (2021) T. Li, S. Deng, H. Liu, S. Sun, H. Li, S. Hu, S. Liu, X. Xing, and J. Chen, “Strong room-temperature ferroelectricity in strained SrTiO3 homoepitaxial film,” Adv. Mater. 33, 2008316 (2021).
  * Gao  _et al._ (2018) P. Gao, R. Ishikawa, B. Feng, A. Kumamoto, N. Shibata, and Y. Ikuhara, “Atomic-scale structure relaxation, chemistry and charge distribution of dislocation cores in SrTiO3,” Ultramicroscopy 184, 217 (2018).
  * Alonso  _et al._ (2019) L. M. Alonso, C. Mochales, L. Nascimento, and W. D. Müller, “Electrochemical deposition of Sr and Sr/Mg-co-substituted hydroxyapatite on Ti-40Nb alloy,” Mater. Lett. 248, 65 (2019).
  * Burger  _et al._ (2020) B. Burger, P. M. Maffettone, V. V. Gusev, C. M. Aitchison, Y. Bai, X. Wang, X. Li, B. M. Alston, R. C. B. Li, N. Rankin, B. Harris, R. S. Sprick, and A. I. Cooper, “A mobile robotic chemist,” Nature 583, 237 (2020).
  * Ziatdinov  _et al._ (2022) M. A. Ziatdinov, Y. Liu, A. N. Morozovska, E. A. Eliseev, X. Zhang, I. Takeuchi, and S. V. Kalinin, “Hypothesis learning in automated experiment : Application to combinatorial materials libraries,” Adv. Mater. 34, 2201345 (2022).
  * Stach  _et al._ (2021) E. Stach, B. DeCost, A. G. Kusne, J. H. Simpers, K. A. Brown, K. G. Reyes, J. Schrier, S. Billinge, T. Buonassisi, I. Foster, C. P. Gomes, J. M. Gregoire, A. Mehta, J. Montoya, E. Olivetti, C. Park, E. Rotenberg, S. K. Saikin, S. Smullin, V. Stanev, and B. Maruyama, “Autonomous experimentation systems for materials development : A community perspective,” Mater. 4, 2702 (2021).
  * Casella (1985) G. Casella, “An introduction to empirical Bayes data analysis,” The Am. Stat. 39, 83 (1985).
  * Carlin and Louis (2000) B. P. Carlin and T. A. Louis, “Empirical Bayes: Past, present and future,” J. of the Am. Stat. Assoc. 95, 1286 (2000).
  * Urenholt and Jensen (2019) A. K. Urenholt and B. S. Jensen, “Efficient Bayesian optimization for target vector estimation,” in  _Proc. of the 22nd International Conference on Artificial Intelligence and Statistics_ (2019) p. 2661.