# Bayesian Optimization for Stable Properties Amid Processing Fluctuations in
Sputter Deposition

Ankit Shrivastava<EMAIL_ADDRESS>Sandia National Laboratories  Matias
Kalaswad<EMAIL_ADDRESS>Sandia National Laboratories  Joyce O. Custer
<EMAIL_ADDRESS>Sandia National Laboratories  David P. Adams
<EMAIL_ADDRESS>Sandia National Laboratories  Habib N. Najm
<EMAIL_ADDRESS>Sandia National Laboratories

###### Abstract

We introduce a Bayesian optimization approach to guide the sputter deposition
of molybdenum thin films, aiming to achieve desired residual stress and sheet
resistance while minimizing susceptibility to stochastic fluctuations during
deposition. Thin films are pivotal in numerous technologies, including
semiconductors and optical devices, where their properties are critical.
Sputter deposition parameters, such as deposition power, vacuum chamber
pressure, and working distance, influence physical properties like residual
stress and resistance. Excessive stress and high resistance can impair device
performance, necessitating the selection of optimal process parameters.
Furthermore, these parameters should ensure the consistency and reliability of
thin film properties, assisting in the reproducibility of the devices.
However, exploring the multidimensional design space for process optimization
is expensive. Bayesian optimization is ideal for optimizing inputs/parameters
of general black-box functions without reliance on gradient information. We
utilize Bayesian optimization to optimize deposition power and pressure using
a custom-built objective function incorporating observed stress and resistance
data. Additionally, we integrate prior knowledge of stress variation with
pressure into the objective function to prioritize films least affected by
stochastic variations. Our findings demonstrate that Bayesian optimization
effectively explores the design space and identifies optimal parameter
combinations meeting desired stress and resistance specifications.

##  I Introduction

Recent advances in computational, machine learning and statistical methods
have significantly impacted the capabilities for optimally designing
metallurgical thin films. Thin films are used in various microelectronics,
optics, and mechanical applications ranging from integrated circuits to
reflective or protective coatings. Among metallurgical thin films, refractory
metals are beneficial due to their remarkable resistance to heat and wear.
Their high melt temperatures, low diffusivity, and good electrical
conductivity often make them ideal for applications involving elevated
temperatures and electronics. Thermally stable contacts, Ohmic contacts, and
interconnects are prevalent in the semiconductor industry, and all capitalize
on such properties of refractory metal thin films and coatings Ashtiani
(2023); Messica, Meirav, and Shtrikman (1995); Erofeev, Fedin, and Fedina
(2018). Notably, refractory thin films can be fabricated using various
methods, such as chemical vapor deposition, sputter deposition, and vacuum
thermal evaporation. The requisite process steps, conditions, and tolerances
play a crucial role in achieving target properties, such as low resistance
during the fabrication of thin films for specific applications.

Determining the best process conditions to fabricate thin films with desired
properties requires optimizing the process-property models. Often, information
about such models is limited, necessitating exploring the process conditions
space. Exploration based on expert intuition and grid search Yasrebi  _et al._
(2014) has been used for optimizing process conditions for thin films.
However, these methods may require many experiments to be carried out,
depending on the dimension of the process and property parameters (design
space), making them time-consuming and costly. Due to experimental budget
constraints, an efficient design of experiments (DoE) to search for the
optimum is imperative. In addition to producing thin films with target
properties, the optimal properties achieved must remain robust to unavoidable
process fluctuations. This robustness is essential to ensure consistently good
process target properties, sometimes at the expense of some sub-optimality. A
highly optimized but non-robust process output is not a desirable setup for a
manufacturing context.

Over the years, various DoE methods Yücel, Güler, and Yücel (2014); Dreer,
Krismer, and Wilhartitz (1999); Ramana  _et al._ (2012) have been proposed to
guide toward optimum process conditions sequentially. These DoE methods
involve iteratively exploring sampled processing conditions and learning a
process-property surrogate model (response surface) to determine the next set
of process conditions for experiments. DoE methods typically sample process
conditions using Latin hypercube sampling, full factorial Fisher (1936);
Madhava (1956), or central compositeBox and Wilson (1951) to learn a response
surface. For robust design, methods such as the Taguchi methodTaguchi (1985)
have been proposed to handle experimental noise. The extension of DoE methods
for constraint optimizations has also been suggested in work such as that by
Jones et al. Jones, Schonlau, and Welch (1998). However, these methods learn a
predefined parametric form, such as linear and quadratic surfaces, which
limits their applicability to complex systems. Additionally, these methods
tend to be exploitative and often converge to a local optimum, requiring
extensive parameter space exploration to reach the global optimum. Extensive
exploration increases the number of trials exponentially, making classical DoE
methods costly for resource-intensive experiments such as sputtering
deposition.

Data-driven Bayesian optimization (BayesOpt) algorithms have been found to
model complex systems and explore the process conditions more
efficientlyPackwood (2017); Garnett (2023). Unlike classical DoE methods,
BayesOpt employs a probability distribution over functions using nonparametric
models such as the Gaussian process. The flexibility of this nonparametric
model enables learning complex nonlinear relationships. Furthermore, the
probability distribution quantifies uncertainties, which help differentiate
the regions of parameter space as explored and unexplored, helping to balance
exploitation and exploration. Balancing exploitation and exploration helps
prevent locally optimal solutions and suboptimal searches, leading to faster
convergence to a global optimum. Due to its probabilistic structure, BayesOpt
can also handle experimental noise with fewer trials than the Taguchi method.
Further, various BayesOpt algorithmsDaulton  _et al._ (2022); Sanders  _et
al._ (2019); Le and Branke (2020); Nguyen  _et al._ (2021) have been proposed
to obtain robust process output efficiently.

In materials manufacturing, Bayesian optimization has guided experiments and
optimized process conditions for desired properties. It has been demonstrated
that BayesOpt converges faster than traditional DoE methods in multiple
material science applications Liang  _et al._ (2021); Lookman  _et al._
(2019). In particular, BayesOpt has optimized thin film growth targeting
desired properties. Studies have successfully optimized crystal structures
Ohkubo  _et al._ (2021), optical contrast Kusne  _et al._ (2020), residual
resistivity ratio Wakabayashi  _et al._ (2019), and contact resistance
Miyagawa  _et al._ (2021). In recent years, various strategies based on
BayesOpt have been suggested to address challenging optimization problems in
material science, such as multi-objective problems Gopakumar  _et al._ (2018);
Khatamsaz  _et al._ (2022), and mixed-variable problems with robust design
Zhang, Apley, and Chen (2020).

In this work, we present a Bayesian optimization construction to identify
process conditions that satisfy multiple property constraints for thin films
while simultaneously maximizing optimized process robustness against minor
variations in sputter deposition conditions. To our knowledge, existing work
for optimizing thin film deposition has not addressed the challenge of
concurrently identifying process conditions that aim to satisfy both
optimality and robustness.

The proposed approach is tested on molybdenum (Mo) thin films, a well-studied
refractory metal, with results readily validated against existing literature.
Coatings with Mo have practical uses, such as in Mo/Si multilayers for
extreme-UV optics Barbee, Mrowka, and Hettrick (1985); Mirkarimi (1999) and
back contacts in CIGS solar cells Orgassa, Schock, and Werner (2003); Li  _et
al._ (2016). In the latter application, and often generally, the Mo layer
should have low electrical resistance and residual stress to ensure good
adhesion to other layers. We used these criteria to build an objective
function for the Bayesian optimization of Mo films presented here.

In our BayesOpt guided experiment, we fabricated Mo thin films using sputter
deposition, a well-established technique compatible with refractory metal
microelectronics. We deposited a range of films with various sputter powers
and argon (Ar) pressures – two process parameters that typically significantly
influence film properties. We then used data from measurements of these
fundamental properties (e.g., resistance, stress) to update the process-
property model and calculate the next set of deposition conditions to explore.
We repeated this process until a stable set of conditions was found that
produced Mo thin films with the desired properties.

We observed that BayesOpt can guide the sputter deposition experiment
efficiently. With a few additional experiments and prior knowledge, the
algorithm explored feasible power and pressure values, ultimately discovering
the optimal configuration that satisfied the desired criteria for Mo films. It
was also observed that along with being robust to experimental noise during
deposition, our BayesOpt construction was also robust to unquantified
environmental uncertainties due to deposition chamber cleaning.

The remainder of this article is organized as follows. In section II, we
discuss the design objectives for Mo thin film. In section III, we elucidate
the Bayesian optimization methodology employed to achieve the design
specifications. Section IV delves into the details of our experimental setup.
Section V presents the obtained results, and in section VI, we draw our
conclusions from the study.

##  II Design specifications for Thin Films

We aim to discover power and pressure conditions for sputter deposition of Mo
thin films such that they meet the following design specifications (criteria):

  1. 1.

The first criterion requires that the films exhibit residual stresses in the
range of ($-300$ MPa, $300$ MPa). For Mo thin films, this range is considered
acceptable for low residual stress. A low residual stress criterion is crucial
for ensuring the integrity of the resulting thin film against various failure
mechanisms. Like almost all thin film applications, films must adhere to a
substrate or underlying surface. Films with large residual stresses are more
susceptible to failure mechanisms such as buckling, cracking, and delamination
– phenomena detrimental to the device or coating’s overall functionAbadias
_et al._ (2018).

  2. 2.

The second criterion specifies that the films exhibit resistance below $3$
$\Omega/\text{sq}$, which falls within an acceptable range for low sheet
resistance. In most microelectronic applications, it is also crucial that
metal layers have low resistance to minimize losses and resistive heating.

  3. 3.

The third criterion requires that the films be as dense as possible in a
qualitative sense. Density plays a crucial role in critical applications, such
as microelectronics, where low-density metal layers lead to higher resistance
(pores/voids restrict electrical paths). We enforce this criterion by favoring
design configurations where the derivative of stress with Ar pressure is
positive. This criterion is supported by the behavior of magnetron-sputtered
metal films. These films demonstrate a compressive-to-tensile stress
transition with increasing process pressure. Films with stress near this
transition (characterized by a large, positive slope) are generally denser
than those deposited at higher pressures, assuming other factors are held
similar. At Ar pressures greater than this transition point — where the
derivative of stress becomes negative — films are typically underdense and
porousAbadias  _et al._ (2018).

  4. 4.

Finally, the fourth objective requires that the films’ stress be least
susceptible to pressure fluctuations. As already discussed, sputtered metal
films can undergo a large change in stress over a narrow range of process
conditions. Film stress can be particularly sensitive to sputter pressure. We
target a robust set of process conditions where a slightly different pressure
causes a minimal change in the film stress. We achieve this objective by
minimizing the derivative of the stress with Ar pressure. This objective helps
address the manufacturing challenge of finding deposition conditions where the
properties of a sputter-deposited metal are robust to process fluctuations.

##  III Optimization for Design Specifications

We use BayesOpt to iteratively search for the optimal power and pressure pairs
within the bounded domain of feasible sputter deposition configurations,
denoted by $\mathcal{X}\subset\mathbb{R}_{\geq 0}^{2}$ such that they satisfy
the design specifications as discussed in Section II.

In the following, we introduce our notation. During each iteration of
experiments, thin films are deposited using a given process condition (design
point) $\bm{x}=(x_{\text{po}},x_{\text{pr}})\in\mathcal{X}$, where the
subscript “po" and “pr" denotes the power and pressure dimension respectively.
Subsequently, their stress $S(\bm{x})$ and sheet resistance $R(\bm{x})$ are
measured. The set of $n$ power and pressure pairs at which thin film
measurements have been collected is denoted as
$\mathcal{D}=[\bm{x}_{1},\bm{x}_{2},\dots,\bm{x}_{n}]$, where
$\mathcal{D}\subset\mathcal{X}$. The set of the outputs of any function
$y(\bm{x})$ at the observed pairs are denoted as
$y(\mathcal{D})=[y(\bm{x}_{1}),y(\bm{x}_{2})\dots y(\bm{x}_{n})]^{T}$. The set
of collected property measurements with corresponding process conditions
$\\{\mathcal{D},S(\mathcal{D}),R(\mathcal{D})\\}$ is referred to as the
observed dataset.

The primary challenge lies in the fact that the stress $S(\bm{x})$ and sheet
resistance $R(\bm{x})$ functions are unknown (black box). The only information
available is the set of measurements of stress, $S(\mathcal{D})$ and sheet
resistance $R(\mathcal{D})$ at observed design points $\mathcal{D}$. Using the
available information, we first obtain an objective function (Section III.1),
which we then optimize using the Bayesian optimization algorithm (Section
III.2) to satisfy the design specifications. The method is summarized in
algorithm LABEL:alg:bayesian_optimization.

###  III.1 Obtaining a Unified Design Objective

The choice of the objective function is made such that it inherently accounts
for all design specifications. First, we start by selecting four distinct
smooth functions, $f_{1},f_{2},f_{3},f_{4}:\mathcal{X}\rightarrow\mathbb{R}$,
that satisfy the four design specifications individually. Afterward, we
combine the separate functions into a single unified smooth function, $f$,
that serves as an objective for our optimization. Selecting smooth functions
helps effectively fit the objective function by the surrogate model in the
Bayesian algorithm. The four functions are defined in the following.

With $T_{s}$ defined as a threshold stress magnitude, we employ the first
objective $f_{1}(\bm{x})$ to enforce the requirement that the stress satisfies
$-T_{s}\leq S(\bm{x})\leq T_{s}$. Thus, $f_{1}(\bm{x})$ is defined as

$\displaystyle f_{1}(\bm{x})$
$\displaystyle=\tau(T_{s}+S(\bm{x}))+\tau(T_{s}-S(\bm{x}))$ (1)

where

$\displaystyle\tau(z)$ $\displaystyle=\frac{1}{1+\exp^{-z}}\quad
z\in\mathbb{R}$

is the standard logistic function, which maps its input to a range between $0$
and $1$ with a smooth step-like transition. This objective maps the stress
values to $(0,1)$ such that $f_{1}(\bm{x})\simeq 1$ for
$S(\bm{x})\in[-T_{s},T_{s}]$ and $0$ elsewhere.

The second objective requires the films to have $R(\bm{x})<T_{r}$, where
$T_{r}$ is a specified threshold sheet resistance value. This objective is
enforced using $f_{2}(\bm{x})$, a modified form of the logistic function,
given as

$\displaystyle f_{2}(\bm{x})$
$\displaystyle=\frac{1}{1+\exp^{-m(T_{r}-R(\bm{x}))}},$ (2)

where the scalar $m$ controls the transition slope in range $(0,1)$. This
function maps $R(\bm{x})$ values below $T_{r}$ to $1$, and $0$ elsewhere.

The third objective, requiring the thin films to be dense, is enforced by
ensuring that the derivative of stress with pressure is positive,
$\frac{\partial S(\bm{x})}{\partial x_{\text{pr}}}>0$, as discussed in Section
II. Since the function of stress $S(\bm{x})$ is unknown, to obtain the
derivative, we first fit a Gaussian process (GP) mean interpolant $g(\bm{x})$
over stress observations and then compute its derivative with respect to
pressure. Selecting the GP mean as an interpolant eliminates the need to make
assumptions about stress’s functional form. The details about the GP
construction are provided in Section III.2.1. The GP mean function is given by

$\displaystyle S(\bm{x})\approx g(\bm{x})$
$\displaystyle=k_{\phi}(\bm{x},\mathcal{D})(k_{\phi}(\mathcal{D},\mathcal{D})+\sigma^{2}_{\text{noise}}I_{n})^{-1}S(\mathcal{D}),$
(3)

where $I_{n}=\text{diag}(1,\ldots,1)\in\mathbb{R}^{n\times n}$, we presume a
zero GP prior mean function, additive Gaussian independent identically
distributed data noise with zero mean and $\sigma^{2}_{\text{noise}}$
variance, and $k_{\phi}(\cdot,\cdot)$ is a squared exponential radial basis
function (RBF) kernel with hyperparameters $\phi$. The kernel structure
follows the formulation defined in Eq. 11, Section III.2.2, and its
hyperparameters are tuned similarly. It is important to note that $g(\bm{x})$
differs from the surrogate model utilized in the context of Bayesian
optimization as described in Section III.2. To ensure clarity, we will use the
term ’interpolant’ to represent $g(\bm{x})$ for the rest of the article. With
this surrogate for $S(\bm{x})$ in hand, the 3rd objective is met using a
hyperbolic tangent function, and $f_{3}(\bm{x})$, defined as

$\displaystyle f_{3}(\bm{x})$
$\displaystyle=\frac{\tanh(g_{\text{pr}}(\bm{x}))+1}{2};\quad
g_{\text{pr}}(\bm{x}):=\frac{\partial g(\bm{x})}{\partial x_{\text{pr}}}.$ (4)

where, $f_{3}(\bm{x})$ takes the derivative of stress with pressure as input
and maps to $1$ for the nonnegative input and $0$ elsewhere.

The final objective requires minimal film stress dependence on pressure
fluctuation. To achieve this, we maximize the function $f_{4}(\bm{x})$,
defined as

$\displaystyle f_{4}(\bm{x})$
$\displaystyle=\frac{-g_{\text{pr}}(\bm{x})}{B}+1,$ (5)

where $B$ is a constant scale factor. This function output increases linearly
with decreasing $g_{\text{pr}}(\bm{x})$. We enforce the $f_{4}(\bm{x})$ range
to be roughly between $0$ and $1$, aligning it in scale with the other three
functions, with a suitable choice of $B$.

Finally, we define the composite function $f(\bm{x})$ as

$\displaystyle f(\bm{x})$
$\displaystyle=f_{1}(\bm{x})f_{2}(\bm{x})f_{3}(\bm{x})f_{4}(\bm{x}).$ (6)

The overall output is influenced by all functions, effectively taking care of
all four criteria. The composite function $f(\bm{x})$ is the objective
function for our optimization problem.

We have illustrated the behavior of these functions in Fig. 1, where we offer
an example of stress and sheet resistance in a one-dimensional pressure-
variation context for a fixed power condition. The stress and sheet resistance
dependence on pressure is shown in Figs. 1(a) and 1(b). The corresponding
profiles of $f_{1},f_{2},f_{3},f_{4}$ are shown in Figs. 1(c), 1(d), 1(e), and
1(f), respectively. The final objective function $f$ is shown in Fig. 1(g).
Notably, for this example, the maximum of $f$ lies at a pressure value of $2$
mTorr for the given power setting, indicating the deposition conditions that
satisfy the design criteria.

(a) stress, $S(\bm{x})$ vs pressure

(b) resistance, $R(\bm{x})$ vs pressure

(c) stress criterion, $f_{1}(\bm{x})$ vs pressure

(d) resistance criterion, $f_{2}(\bm{x})$ vs pressure

(e) postive derivative criterion, $f_{3}(\bm{x})$ vs pressure

(f) minimum derivative criterion, $f_{4}(\bm{x})$ vs pressure

(g) unified objective function, $f(\bm{x})$ vs pressure

Fig. 1:  Figures 1(a) and 1(b) illustrate the stress and resistance variations
with pressure for a fixed power, respectively. For the given stress and
resistance profiles, the behavior of the four design criteria functions
$f_{1}$, $f_{2}$, $f_{3}$, and $f_{4}$ are shown in Figs. 1(c), 1(d), 1(e),
1(f), respectively. The pressure values at which these functions reach their
maximum satisfy the corresponding design criteria. The pressure that maximizes
the unified objective function $f$, shown in Fig. 1(g), satisfies all design
criteria.

Note that the four design criteria can be conceptualized as elements of a
constrained optimization problem, with the first three criteria formulated as
constraints and the fourth used as the objective function. One way to solve
this problem involves using a Lagrangian approach, where constraints are
included with penaltiesGramacy  _et al._ (2016) within the objective function,
transforming the problem into an unconstrained optimization task. This
adjustment removes the strict criteria of the original constraints, opting
instead to penalize deviations from them, effectively relaxing the constraint
criteria. Several extensions of BayesOpt algorithms have also been proposed to
solve constrained optimization problems where the methods leverage variants of
acquisition functions such as expected improvementGardner  _et al._ (2014);
Feliot, Bect, and Vazquez (2017) and predictive entropy search with
constraintsHern  _et al._ (2016). However, these variants of acquisition
functions for constrained problems can be costly for resource-intensive
experiments, as they can result in the rejection of exploitation-based
proposals. Given the limited experimental budget, we opted to formulate the
problem in a manner inspired by filter-based optimization
approachesPourmohamad and Lee (2020); Sainvitu and Toint (2007). The proposed
approach prioritizes the search toward the feasible space, reducing the number
of trials required to achieve the optimum.

###  III.2 Bayesian Optimization Algorithm

The BayesOpt algorithm aims to find a design that optimizes the objective
function $f(\cdot)$ by iteratively evaluating $f(\cdot)$ at existing design
points and using these evaluations, $f(\mathcal{D})$ to propose new design
points. The process involves two tasks: (1) building a probabilistic surrogate
model from the existing observations; see Section III.2.1, and (2) using an
acquisition function that guides experiments where to sample next in search
space; see Section III.2.3. The iteration process continues until a
termination criterion is achieved; see Section III.3.

####  III.2.1 Building Surrogate Model of Objective Function

We use a Gaussian process Rasmussen and Williams (2005) (GP) as the surrogate
model, which, when fitted to $(\mathcal{D},f(\mathcal{D}))$, provides a
probabilistic prediction of $f(\bm{x})$ at all points $\bm{x}\in\mathcal{X}$.
This information is used in the acquisition function to sample the next point.
GPs are nonparametric models, which makes them highly flexible to fit complex
mappings.

A GP is a probability distribution over random functions that are jointly
Gaussian at any finite set of points. It is defined by its mean function
$\mu(\bm{x})$ and covariance function $k(\bm{x},\bm{x}^{\prime})$. In our
present setting, we define the covariance function as a squared exponential
kernel $k_{\theta}(\bm{x},\bm{x}^{\prime})$, specified later in this section
so that we write the GP as

$\displaystyle\hat{f}(\bm{x})$
$\displaystyle\sim\mathcal{GP}(\mu(\bm{x}),k_{\theta}(\bm{x},\bm{x}^{\prime}))\quad\forall\bm{x}\in\mathcal{X}.$
(7)

Here, $\hat{f}$ is a surrogate function of $f$, sampled from the Gaussian
process $\mathcal{GP}$. Using observed data independent coordinates
$\mathcal{D}=[\bm{x}_{1},\ldots,\bm{x}_{n}]$, we define a prior distribution

$\mathcal{GP}(\mu(\mathcal{D}),k_{\theta}(\mathcal{D},\mathcal{D})),$ (8)

where $\mu(\mathcal{D})=[\mu(\bm{x}_{1}),\mu(\bm{x}_{2})\dots\mu(\bm{x}_{n})]$
is the prior mean vector and
$k_{\theta}(\mathcal{D},\mathcal{D})\in\mathbb{R}^{n\times n}$ is the prior
covariance matrix. The estimate of the function $f$ at any unobserved design
point $\bm{x}^{\ast}$ is obtained using the following posterior distribution

$\displaystyle
P(\hat{f}(\bm{x}^{\ast})|f(\mathcal{D}),\theta)=\mathcal{N}(\mu_{\text{pos}}(\bm{x}^{\ast}),\sigma_{\text{pos}}^{2}(\bm{x}^{\ast})),$
(9)

where, with $k_{\theta}(\bm{x}^{\ast},\mathcal{D})\in\mathbb{R}^{1\times n}$,
$k_{\theta}(\mathcal{D},\bm{x}^{\ast})\in\mathbb{R}^{n\times 1}$, presuming
additive Gaussian independent identically distributed data noise with zero
mean and $\sigma^{2}_{\text{noise}}$ variance, and a prior GP mean function
$\mu(\bm{x})\equiv\bm{0}$, we have

$\displaystyle\begin{split}&\mu_{\text{pos}}(\bm{x}^{\ast})=k_{\theta}(\bm{x}^{\ast},\mathcal{D})(k_{\theta}(\mathcal{D},\mathcal{D})+\sigma^{2}_{\mathrm{noise}}I_{n})^{-1}f(\mathcal{D})\\\
&\sigma_{\text{pos}}^{2}(\bm{x}^{\ast})=k_{\theta}(\bm{x}^{\ast},\bm{x}^{\ast})-\\\
&\quad
k_{\theta}(\bm{x}^{\ast},\mathcal{D})(k_{\theta}(\mathcal{D},\mathcal{D})+\sigma^{2}_{\mathrm{noise}}I_{n})^{-1}k_{\theta}(\mathcal{D},\bm{x}^{\ast})\end{split}$
(10)

The $\mu_{\text{pos}}(\bm{x}^{\ast})$ and
$\sigma_{\text{pos}}^{2}(\bm{x}^{\ast})$ are the posterior mean and variance
of the function $f$ at $\bm{x}^{\ast}$ conditioned on the observed data,
respectively. We employ an anisotropic squared exponential RBF kernel, where,
for any two points $(\bm{x},\bm{x}^{\prime})$,

$\displaystyle k_{\theta}(\bm{x},\bm{x}^{\prime})$
$\displaystyle=A\cdot\exp\left(-(\bm{x}-\bm{x}^{\prime})^{T}M^{-1}(\bm{x}-\bm{x}^{\prime})\right)$
(11)

where $M=\text{diag}(l_{\text{po}},l_{\text{pr}})\in\mathbb{R}^{2}$, and
$(l_{\text{po}},l_{\text{pr}})$ are the correlation lengths in the power and
pressure directions, respectively. The values of hyperparamters
$\theta=\\{A,l_{\text{po}},l_{\text{pr}}\\}$ are selected using hyperparamter
tuning, discussed in Section III.2.2. Finally, the posterior mean and variance
at all points in search space $\mathcal{X}$ is used to estimate the
acquisition function; see Section III.2.3.

####  III.2.2 Selecting Hyperparameters of Gaussian Process

The hyperparameters $\theta$ used in the Gaussian process for building the
surrogate model during Bayesian optimization and the hyperparameters $\phi$
for the interpolant $g(\bm{x})$ are selected by maximizing the log-likelihood
of the observed data $f(D)$, as shown below with
$\lambda:=(\theta,\sigma_{\text{noise}})$

$\displaystyle\hat{\theta}$
$\displaystyle=\operatorname*{arg\max}_{\theta}\log P(f(\mathcal{D})|\lambda)$
(12) $\displaystyle\log P(f(\mathcal{D})|\lambda)$
$\displaystyle=-\frac{1}{2}f(\mathcal{D})^{T}(k_{\theta}(\mathcal{D},\mathcal{D})+\sigma_{\text{noise}}^{2}I_{n})^{-1}f(\mathcal{D})$
$\displaystyle\quad-\frac{1}{2}\log|k_{\theta}(\mathcal{D},\mathcal{D})+\sigma_{\text{noise}}^{2}I_{n}|$
$\displaystyle\quad-\frac{n}{2}\log{2\pi}.$ (13)

The optimization problem, Eq. 12, is solved using gradient descent on Eq. 13.
Since the unified objective function $f(\cdot)$ is deterministic, we choose a
predetermined small value of $\sigma_{\text{noise}}=0.01$ for regularization
purposes. As for the GP regression context for the stress surrogate in Eq. (3)
above, we used $\sigma_{\text{noise}}=0.058$ based on the observed data noise
from the experiments conducted prior to BayesOpt.

####  III.2.3 Acquisition Function

The acquisition function provides the means for BayesOpt to balance
exploration and exploitation. We use the upper confidence bound (UCB) as the
acquisition function to sample new power and pressure values where the sputter
deposition is conducted. We chose UCB given its simplicity and broad utility.
The UCB acquisition function characterizes the input space $\mathcal{X}$ into
regions for exploration with high uncertainty and regions for exploitation
that likely contain the objective function’s optimum. The UCB function is the
sum of the posterior mean and scaled standard deviation obtained from the
surrogate model, given by

$\displaystyle\text{ucb}(\bm{x})$
$\displaystyle=\mu_{\text{pos}}(\bm{x})+\beta\cdot\sigma_{\text{pos}}(\bm{x}),$
(14)

where the hyperparameter $\beta$ controls the balance between exploration and
exploitation. A larger value of $\beta$ will drive the BayesOpt algorithm to
enhance its priority for exploration over exploitation. The solution of the
UCB optimization problem is the next design point in the search space
$\mathcal{X}$, given by

$\displaystyle\bar{\bm{x}}$
$\displaystyle=\operatorname*{arg\max}_{\bm{x}\in\mathcal{X}}\text{ucb}(\bm{x}).$
(15)

The optimal solution of the UCB function is obtained using gradient-based
optimizer L-BFGS-B with several random starts to estimate good hyperparameter
values. The optimum value $\bar{\bm{x}}$ is the next point where subsequent
sputter deposition is conducted.

###  III.3 Algorithm termination strategy

Our approach for algorithm termination hinges on achieving the convergence of
the optimal solution $\bm{x}^{*}$. During BayesOpt guided search, the optimal
solution at each iteration is identified as the configuration within the
current dataset $\mathcal{D}$ that maximizes the objective function $f(x)$ as
shown below.

$\displaystyle\bm{x}^{*}=\operatorname*{arg\max}_{\bm{x}\in\mathcal{D}}f(\bm{x})$
(16)

The optimal solution is determined by simply identifying the configuration for
which the function attains its maximum value in the set $f(\mathcal{D})$. We
monitor the optimal solution at each iteration until its value remains
stagnated for a maximum predefined count. Subsequently, we assess all four
objective criteria, objective function, and uncertainty and conclude the
termination of the experiment.

##  IV Experimental Setup

For this work, the stress and sheet resistance threshold values were chosen as
$T_{s}=300$ MPa and $T_{r}=3$ $\Omega/\text{sq}$. The BayesOpt-guided
experiment was started with the initial measurements of stress and sheet
resistance of the thin film, deposited using $27$ different configurations of
power and pressure settings, as shown in Fig. 2.

At each iteration, a pair of power and pressure settings for the sputter
deposition of Mo films was recommended by the BayesOpt algorithm.
Subsequently, the deposited films’ stress and sheet resistance were measured
and fed back into the BayesOpt algorithm for the next iteration. Section IV.1
describes the settings under which the Bayesian optimization algorithm was
executed. Section IV.2 describes the sputter deposition of Mo films, and
sections IV.4 and IV.5 describe the measurements of stress and sheet
resistance, respectively.

(a)

(b)

Fig. 2:  The plots depict the variations of residual stress, shown in Fig.
2(a), and variations of sheet resistance, shown in Fig. 2(b) with pressure and
powers in the dataset before Bayesian optimization. The observations from the
plots assisted in constructing the objective functions.

###  IV.1 Bayesian optimization and objective function setup

The BayesOpt algorithm was initialized using stress and sheet resistance
measurements for the $27$ initial configurations. The BayesOpt search space,
$\mathcal{X}$, was defined between the pressure values of $2$ mTorr to $23$
mTorr and the power values ranging from $50$ W to $750$ W. The objective
function constants were set as $m=5$ in Eq. 2 and $B=1900$ in Eq. 5. The
Gaussian process hyperparameters $\phi$, for the interpolant $g(\bm{x})$ and
$\theta$ for the surrogate model, were tuned as explained in Section III.2.2.
Finally, $\beta=1$ was used for the UCB function.

The code for the BayesOpt algorithm was implemented in Python 3.9.2 and
executed on a Ubuntu 20.04.6 LTS 64-bit machine with 32 Gb RAM and 11th Gen
Intel Core i7-1185G7 @ 3.00GHz × 8 processor. During each iteration, the
algorithm required roughly $10-12$ minutes to finish the computation and
provide the next sputter deposition configuration.

###  IV.2 Thin Film Deposition

(a) Process Chamber

(b) Sputtering Process

Fig. 3:  Vacuum chamber geometry, as shown in Fig. 3(a), is used for physical
vapor deposition employing a top-down sputter geometry with true planetary
sample stage motion (i.e., orbit + spin). Ultra-high-purity argon gas is
consistently introduced into the chamber and regulated by a mass flow
controller (MFC). Fig. 3(b) shows the cross-section drawing of key atomistic
processes involved in sputtering, transport, and film growth.

Molybdenum films were sputter deposited in an Innotec $240$ SD system, as
illustrated in Fig. 3. Before deposition, the chamber was cryo-pumped to a
base pressure of $\leq 5\times 10^{-7}$ Torr. A planar target $20.32$ cm ($8$
inches) in diameter and made of high-purity Mo $(99.95\%)$ was used to deposit
the material in a sputter-down configuration. The substrate rotation stage
consists of six platens moving in planetary motion, with a working distance of
$\approx$9.25 cm between the substrates and target. All films were deposited
onto $5.08$ cm ($2$ inches) diameter Si ($100$) substrates with $400$ nm of
thermal oxide. These substrates were chosen because they possess a single
biaxial modulus, significantly simplifying internal stress value calculations
using Stoney’s equation. Three nominally identical Si wafers were loaded onto
three platens for redundancy for each deposition. Ultra-high purity argon
flowed into the chamber at $50$ sccm, and its working pressure was maintained
via a downstream pressure control (throttled gate valve). Ar pressure was held
constant and monitored to ensure any changes in pressure were $<1\%$ for each
deposition. The same deposition time was used for all sputter pressures so
that the samples spent similar time in the flux of the sputter gun.

###  IV.3 Substrate and Film Thickness

Substrate thickness was measured at five locations on each wafer using a drop
gauge. Each deposition included an additional Si piece of $\approx~{}1$ sq. cm
for film thickness measurements. A portion of this ride-along sample was
masked with Kapton tape so that its removal after deposition would create a
step edge amenable to thickness measurements. To maintain congruity with other
Si substrates, witness pieces were of similar thickness, placed in the center
of a deposition platen at the same working distance, and exposed to the flux
for the same time as other samples. A DEKTAK XT surface profilometer measured
the film height at five step-edge locations, from which an average film
thickness was obtained.

###  IV.4 Stress Measurement of Films

Before the Mo deposition, an FLX2320-S wafer scanning instrument, which
measures the deflection from an incident laser beam, was used to determine an
initial radius of curvature for each wafer. After deposition, each coated
wafer was re-measured with the instrument to find the post-deposition wafer
radius of curvature. Using these values – together with the film thickness,
the wafer thickness measured previously, and the biaxial modulus of Si (100) –
the average in-plane stress of the film was calculated according to Stoney’s
equation Stoney (1909).

###  IV.5 Sheet Resistance Measurement of Films

The Mo film’s room temperature sheet resistance measurements were made using a
four-point probe and a constant current source and meter. Samples were probed
with a single tip assembly having four $100~{}\mu$m radius tips arranged with
1.0 mm spacing. Sheet resistance was measured six times for each of the Mo-
coated wafers at three locations relative to the wafer flat (top, bottom, and
center) and in forward bias and reverse bias at each position. Insulating
substrates (i.e., with thermal oxide) helped isolate the measurement to be
only of the film.

##  V Results

The BayesOpt-guided sputter deposition experiment continued for $10$
iterations until the change in optimal configuration $\bm{x}^{*}=(2\text{
mTorr},620\text{ W})$ no longer exceeded the predefined thresholds in both
power ($\epsilon_{\text{pr}}=0.5$) and pressure
($\varepsilon_{\text{po}}=10$).

During each iteration of the experiments, multiple thin films were deposited
using the configuration proposed by BayesOpt, denoted as $\bar{\bm{x}}$, and
their properties were measured. Subsequently, these measured properties were
used to update the objective function $f(x)$ and estimate the optimal
configuration, $\bm{x}^{*}$. For additional information about the experiments
and property measurements, please refer to Appendix A.

The trends of the BayesOpt algorithm proposals are discussed in Section V.1,
and the rationale behind the conclusion about the optimal solution is detailed
in Section V.2. For the convenience of the discussion, we denote the subregion
of the search space $\mathcal{X}$, with pressure $<3$ mTorr as
$\mathcal{X}_{1}$ and that with pressure $\geq 3$ mTorr as $\mathcal{X}_{2}$.

###  V.1 Bayesian Optimization iteration trends

(a) Pressure v/s iterations

(b) Power v/s iterations

Fig. 4:  The figures illustrate trends observed during the BayesOpt guided
search, including the proposed sputter deposition configuration for the next
iteration and the optimal configuration in that iteration. Figs. 4(a) and 4(b)
show pressure and power trends over iterations, respectively. See Table II in
appendix A for the specific values.

The trends in the configurations proposed by BayesOpt, denoted as
$\bar{\bm{x}}$, for the next sputter deposition experiment at each iteration
are shown in Fig. 4. It was observed that for the first six iterations,
BayesOpt proposed configurations with a fixed pressure value of $2$ mTorr and
varying power (region $\mathcal{X}_{1}$). Then, in later iterations, it
proposed configurations in the complementary region, $\mathcal{X}_{2}$. To
offer insight into the rationale behind the suggested trends, the contour
plots of the posterior mean $\mu_{\text{pos}}(\bm{x})$, standard deviation
$\sigma_{\text{pos}}(\bm{x})$ and $\text{ucb}(\bm{x})$ at the first, seventh,
and tenth iterations are shown in Figs. 5, 6 and 7 respectively.

During the first six iterations, as shown _e.g._ for the first iteration in
Fig. 5, contour plots show that the posterior mean, $\mu_{\text{pos}}(\bm{x})$
in $\mathcal{X}_{1}$ dominated the standard deviation term,
$\beta\cdot\sigma_{\text{pos}}(\bm{x})$ in the entire search space in terms of
influencing $\text{ucb}(\bm{x})$. The high posterior mean in $\mathcal{X}_{1}$
was due to high values in the first three design criteria,
$f_{1},f_{2},f_{3}$, suggesting a likely location for the objective function
maxima. Accordingly, the algorithm exploited the region $\mathcal{X}_{1}$ to
find the optimum. Later, during iterations $7$, $8$ and $9$, as shown _e.g._
for the seventh iteration in Fig. 6, contour plots show the standard
deviation, $\beta\cdot\sigma_{\text{pos}}(\bm{x})$ in the region
$\mathcal{X}_{2}$ dominating the posterior mean in terms of influencing
$\text{ucb}(\bm{x})$. Hence, the algorithm explored the region
$\mathcal{X}_{2}$ due to the lack of observations causing high local
uncertainty. In the $10$th iteration, as shown in Fig. 7, the algorithm
revisited exploitation within the $\mathcal{X}_{1}$ region, only to shift back
to exploration in the $\mathcal{X}_{2}$ region during the $11$th iteration. In
summary, the algorithm goes back and forth between exploration and
exploitation to efficiently search for the optimal solution by balancing the
need to gather information about the objective function and leveraging that
information to refine the search.

(a) Exploitation term, $\mu_{\text{pos}}(\bm{x})$

(b) Exploration term, $\beta\cdot\sigma_{\text{pos}}(\bm{x})$

(c) Acqusition function, $\text{ucb}(\bm{x})$

Fig. 5: The figure shows the contour plots of exploitation term,
$\mu_{\text{pos}}$ in Fig. 5(a), exploration term
$\beta\cdot\sigma_{\text{pos}}$ in Fig. 5(b), and UCB in Fig. 5(c) for the
first iteration with $\beta=1$. At this iteration, BayesOpt is exploiting the
region $\mathcal{X}_{1}$.

(a) Exploitation term, $\mu_{\text{pos}}(\bm{x})$

(b) Exploration term, $\beta\cdot\sigma_{\text{pos}}(\bm{x})$

(c) Acqusition function, $\text{ucb}(\bm{x})$

Fig. 6: The figure shows the contour plots of exploitation term,
$\mu_{\text{pos}}$ in Fig. 6(a), exploration term
$\beta\cdot\sigma_{\text{pos}}$ in Fig. 6(b), and UCB in Fig. 6(c) for the
seventh iteration with $\beta=1$. At this iteration, BayesOpt is exploring the
region $\mathcal{X}_{2}$.

(a) Exploitation term, $\mu_{\text{pos}}(\bm{x})$

(b) Exploration term, $\beta\cdot\sigma_{\text{pos}}(\bm{x})$

(c) Acqusition function, $\text{ucb}(\bm{x})$

Fig. 7: The figure shows the contour plots of exploitation term,
$\mu_{\text{pos}}$ in Fig. 7(a), exploration term
$\beta\cdot\sigma_{\text{pos}}$ in Fig. 7(b), and UCB in Fig. 7(c) for the
tenth iteration with $\beta=1$. At this iteration, BayesOpt is exploring the
region $\mathcal{X}_{2}$.

During the experiments, we evaluated the optimal solution $\bm{x}^{*}$ at each
iteration, as shown in Fig. 4. It was observed that the optimal solution
changed during the initial iterations until it remained consistent at
$(2\text{ mTorr},620\text{ W})$ for the last four iterations. The Bayesian
optimization algorithm was terminated based on the convergence of the optimal
solution $\bm{x}^{*}$.

Note that the final/converged optimal values of stress and resistance,
$(2\text{ mTorr},620\text{ W})$ were in fact proposed at iteration $5$, but
the “optimal“ solution identified by our iterative method remained consistent
only after iteration $8$. This delay occurred because our algorithm also
attempts to learn accurate values of the objective function $f$ at the
observed configurations $\mathcal{D}$. The evaluation of the objective
function depends not only on the values of stress and resistance but also on
the derivative of stress with pressure, refer Eq. 6. Since we cannot directly
obtain the derivative values from the experiments, we rely on the interpolant,
$g(\bm{x})$ for calculations. As more experiments are conducted, $g(\bm{x})$
converges to a more accurate representation of the stress function $S(\bm{x})$
in the region $\mathcal{X}_{1}$, as does its derivative. This leads to the
convergence of our estimated objective function values $f(\mathcal{D})$. Thus,
the delayed convergence of the optimal solution highlights the algorithm’s
process of learning and adapting to the characteristics of the objective
function, ultimately leading to more accurate estimates over time.

###  V.2 Validating the optimal solution

Based on the observations in Fig. 4 and termination criterion discussed in
Section III.3, the optimal configuration was concluded as $(2\text{
mTorr},620\text{ W})$. For verification that the optimal configuration for the
sputter deposition satisfied the design objectives, we analyzed the sheet
resistance values, the interpolant $g(\bm{x})$ for stress, and the derivative
$g_{\text{pr}}(\bm{x})$ at the observed configurations $\mathcal{D}$. These
results are illustrated in Fig. 8. The optimal configuration was first
observed in iteration 5, with resulting stresses $(-180.9,-218.0,-215.5)$ MPa
and sheet resistances $(0.68,0.68,0.77)$ $\Omega/\text{sq}$, refer appendix A
for details. These values satisfy stress requirements in $[-300,300]$ MPa, and
sheet resistance $<3$ $\Omega/\text{sq}$. The optimum is indicated in the
field plots in Fig. 8(d), putting it in the context of the overall variation
of the fields of interest. The derivative of stress with pressure is positive
at the optimum, and the overall objective function is high.

(a) residual stress interpolant, $g(\bm{x})$ in MPa

(b) derivative of stress with pressure, $g_{\text{pr}}(\bm{x})$

(c) sheet resistance, $R(\bm{x})$ in $\Omega/\text{sq}$

(d) objective function, $f(\bm{x})$

Fig. 8: The figure shows the contour plots of stress in Fig. 8(a), the
derivative of stress with pressure in Fig. 8(b), the sheet resistance in Fig.
8(c), and the objective function $f$ in Fig. 8(d) after $10$ iterations. The
optimal conditions are achieved at $(2\text{ mTorr},620\text{ W})$, ensuring
that the observed stress falls within the range of $[-300,300]$ MPa, the sheet
resistance is below $3$ $\Omega/\text{sq}$, and the derivative is both
positive and at its minimum, fulfilling all specified criteria.

For further clarity, we plot the variation of $g_{\mathrm{pr}}(\bm{x})$ and
$f(\bm{x})$ with power, at the optimal pressure of 2 mTorr, in Fig. 9. Here,
we see, indeed, the optimal power at $620$ W. Note that the derivative is, in
fact, minimal around $200$ W; however, clearly, the region preceding $300$ W
fails to meet the other three criteria, such that full objective function
$f(\bm{x})$ is relatively low in that region. Evidently, the configuration
with the smallest derivative value that fulfills all criteria is found at
$620$ W, exhibiting notably low uncertainty given the locally high data
density. We note that the other cluster with high $f(\bm{x})$, although not as
high as the optimum, is the region around $(2\text{ mTorr},320\text{ W})$. It
is possible that this region, with further exploration, might be competitive
with the optimum. However, since we stopped at ten iterations, given
experimental budget constraints, we have the optimal solution at
$\bm{x}^{\ast}=($2 mTorr, 620 W$)$.

(a)

(b)

Fig. 9:  The plot in Fig. 9(a) illustrates the plots of the derivative of
stress with pressure, $g_{\text{pr}}(\bm{x})$ against power. Similarly, Fig.
9(b) illustrates the GP fit of the objective function $f(\bm{x})$ based on
observations at $\mathcal{D}$. Fig. 9(b) presents the variations of the
posterior mean, $\mu_{\text{pos}}(\bm{x})$, (as a dashed line) and the
standard deviation, $\sigma_{\text{pos}}(\bm{x})$, (as a shaded area) with
power. with the posterior mean represented and the standard deviation The peak
value of the objective function is identified at the point
$(2~{}\text{mTorr},620~{}\text{W})$. Both plots are shown for the pressure $2$
mTorr, after $10$ iterations.

The optimal configuration (2 mTorr, 620 W) met the stress and resistance
criteria, with stress values of $(-180.9,-218.0,-215.5)$ MPa and sheet
resistances of $(0.68,0.68,0.77)$ $\Omega/\text{sq}$. Further, to ensure that
this configuration is least affected by minor pressure fluctuations,
additional experiments were conducted in the region $\mathcal{X}_{1}$ to
analyze the sensitivity of film stress to pressure. We conducted experiments
around the optimal power setting of $620$ W, exploring values at $590$ W and
$650$ W, each differing by $30$ W from the optimal setting. Additionally, we
varied the pressure incrementally by $0.5$ mTorr, starting at the optimal
pressure of $2$ mTorr. The observations of mean stress for this experimental
setup are shown in Table I. More detailed information on these additional
experiments is provided in Appendix B. For sensitivity analysis, we computed
the derivative of stress with pressure corresponding to each power setting, as
presented in the table. Results revealed that stress exhibited the least
sensitivity to pressure at $(2$ mTorr, 620 W$)$, with a derivative value of
$339.23$, which increased away from the optimal conditions in either power or
pressure while satisfying all the other criteria, further supporting the
conclusions drawn from the BayesOpt algorithm.

Table I: Table displaying the mean of residual stresses in thin films deposited using the utilized configuration $\bm{x}$ during supplementary experiments. The sensitivity of stress to pressure was analyzed by computing the derivative of stress with respect to pressure for each power setting. The derivative is smallest at $(2\text{mTorr},620\text{W})$, further supporting the optimum observed in Fig. 9.  utilized configuration |  mean stress (MPa) |  derivative of stress with pressure
---|---|---
$\bm{x}$ |  $\bar{S}$ |  $\frac{\Delta\bar{S}}{\Delta x_{\text{pr}}}$
$(2,590)$ |  $-85.58$ |  $975.5$
$(2.5,590)$ |  $402.17$ |
$(2,620)$ |  $-44.00$ |  $399.23$
$(2.5,620)$ |  $155.62$ |  $1030.16$
$(3,620)$ |  $670.70$ |
$(2,650)$ |  $-256.62$ |  $1115.6$
$(2.5,650)$ |  $301.06$ |

##  VI Conclusion

We used Bayesian optimization to identify optimal operating conditions for the
sputter deposition of Mo thin films on Si, targeting desired robust film
properties. Specifically, the objective was to find a combination of power and
pressure such that film stress is in the range of $(-300,300)$ MPa, the sheet
resistance is below $3$ mTorr, and the film is dense. We used a positive
derivative of stress vs. pressure requirement to identify sputter
configurations that lead to dense films. Although density was not directly
measured in this work, we reasonably conclude that the positive stress
gradient requirement promoted denser films. This is supported by our previous
studies of Mo films using the same apparatus described herein, which show by
Rutherford Backscattering Spectrometry that lower process pressures (near the
positive stress gradient) produce denser filmsKalaswad  _et al._ (2023).
Further, to promote stability in stress against pressure fluctuation, we
targeted conditions with a low gradient of stress vs. pressure. Four distinct
smooth functions were employed to achieve these objectives, using a unified
objective function obtained by multiplying the individual objective functions.

The Bayesian optimization algorithm performs an adaptive search over the
parameter space for optimal solutions, iteratively proposing the next sputter
configuration for thin film fabrication and subsequent stress and sheet
resistance measurement. In each iteration, stress values for thin films
obtained from known configurations were utilized to construct a Gaussian
process mean interpolant, whose derivative can be evaluated analytically. The
objective function values at the known configurations were derived from sheet
resistance, stress, and the derivative of stress with respect to pressure. A
Gaussian process was fitted to these objective function values, providing
predictive posterior mean and variance. These, in turn, were used to compute
the upper confidence bound function, guiding the algorithm to propose the next
sputter deposition configuration. The Bayesian optimization guided search
continued for $10$ iterations until the calculated optimal solution at each
step remained consistent. It was observed that the sputter configuration that
maximized the objective function was ($2$ mTorr, $620$ W).

The value of stress, sheet resistance, and derivative of stress with pressure
were verified at the optimal solution. It was observed that all the criteria
were satisfied. The stress values at the optimal configuration were
$(-180.9,-218.0,-215.5)$ MPa, falling within the desired range of $(-300,300)$
MPa, satisfying the first criteria. Similarly, the sheet resistance
measurements were $(0.68,0.68,0.77)\ \Omega/\text{sq}$ remaining below $3\
\Omega/\text{sq}$ which satisfied the second criteria. Finally, the derivative
of stress with pressure at the optimal configuration was positive and
sufficiently low, satisfying the last two criteria. Moreover, further
experiments conducted for sensitivity analysis near the BayesOpt proposed
optimum configuration provided supporting evidence that this configuration was
indeed at a local optimum.

We found that Bayesian optimization efficiently explored the configuration
space, identifying regions likely to yield optimal results. This work also
showed that a unified objective function can simultaneously achieve multiple
desired properties of Mo thin films while ensuring their robustness against
minor variations in sputter deposition conditions. While we did not explore
different weighting of the various objectives, _e.g._ where weight exponents
$(w_{1},\ldots,w_{4})$ are used, resulting in
$f(\bm{x}):=\prod_{i}f_{i}^{w_{i}}(\bm{x})$, this can easily be done, allowing
evaluation of selective weighting strategies. Moreover, comparing the
convergence and efficiency of constrained optimization using a filter-based
objective function with other Bayesian optimization approaches that
incorporate constraints would be useful. The present approach can be extended
to more complex problems involving high dimensional multimodal datasets
involving _e.g._ X-ray diffraction and scanning electron microscope image
data.

###### Acknowledgements.

The authors acknowledge funding under the BeyondFingerprinting Sandia Grand
Challenge Laboratory Directed Research and Development (GC LDRD) program.
Sandia National Laboratories is a multi-mission laboratory managed and
operated by National Technology & Engineering Solutions of Sandia, LLC
(NTESS), a wholly owned subsidiary of Honeywell International Inc., for the
U.S. Department of Energy’s National Nuclear Security Administration
(DOE/NNSA) under contract DE-NA0003525. This written work is authored by an
employee of NTESS. The employee, not NTESS, owns the right, title and interest
in and to the written work and is responsible for its contents. Any subjective
views or opinions that might be expressed in the written work do not
necessarily represent the views of the U.S. Government. The publisher
acknowledges that the U.S. Government retains a non-exclusive, paid-up,
irrevocable, world-wide license to publish or reproduce the published form of
this written work or allow others to do so, for U.S. Government purposes. The
DOE will provide public access to results of federally sponsored research in
accordance with the DOE Public Access Plan.

## Author contributions

The design of the objective function, implementation of the Bayesian
optimization algorithm, and performance analysis were conducted by AS. MK and
JC conducted the sputter deposition of Mo thin films and measured stress and
sheet resistance values. HNN and DPA proposed the research idea and supervised
the work. All the authors verified the data and results and contributed to
writing and reviewing the manuscript.

## Competing interests

The authors declare no competing interests.

## Materials & Correspondence

Requests and correspondence should be addressed to Ankit Shrivastava.

## Data and Code Availability

The Bayesian optimization code for this work is available at
https://github.com/ashriva16/bayesian-optimization-sputter-deposition.

##  Appendix A Data from Bayesopt guided experiments

This section provides comprehensive information about the BayesOpt guided
experiments. The values of the proposed configurations $\bar{\bm{x}}$, and
optimal configurations $\bm{x}^{*}$ during Bayesian optimization, discussed
Fig. 4, are listed in Table II. The proposed configurations were rounded off
for the experiments, as feasible power values were restricted to multiples of
$10$ W, and feasible pressure values were constrained to increments of $0.5$
mTorr. The corresponding feasible configurations used during the experiments
are also listed in the table.

At each iteration, several films were deposited using the respective utilized
configuration, and their residual stress and sheet resistance were measured,
as detailed in Tables III and IV, respectively. Note that more than three
depositions were performed for iterations $3$ and $4$. During the series of
experiments, chamber cleaning became necessary due to the accumulation of
previously sputtered materials inside the chamber. Over time, this buildup can
peel off from chamber components and contaminate samples or cause arcing
during a deposition. The cleaning process can also introduce moisture and
other vapors into the chamber, which is known to affect the film properties of
depositions immediately following the cleaningMattox (1998). This effect
likely explains a larger variation in stress values of iterations $3$ and $4$,
which include films deposited before and after chamber cleaning.

Table II:  Tabulated data showing trends in configurations observed during the BayesOpt guided search at each iteration. The feasible configurations utilized for the experiments are shown in the fourth column.  Iteration |  Optimal |  Proposed |  Utilized
---|---|---|---
$\mathbf{i}$ |  $\bm{x}^{*}$ |  $\bar{\bm{x}}$ |  $\bm{x}$
1 |  $(2.0,750)$ |  $(2.0,648.28)$ |  $(2.0,650)$
2 |  $(2.0,750)$ |  $(2.0,433.32)$ |  $(2.0,430)$
3 |  $(2.0,750)$ |  $(2.0,558.49)$ |  $(2.0,560)$
4 |  $(2.0,650)$ |  $(2.0,667.06)$ |  $(2.0,670)$
5 |  $(2.0,650)$ |  $(2.0,615.34)$ |  $(2.0,620)$
6 |  $(2.0,620)$ |  $(2.0,322.11)$ |  $(2.0,320)$
7 |  $(2.0,320)$ |  $(4.35,323.98)$ |  $(4.5,320)$
8 |  $(2.0,620)$ |  $(18.46,302.29)$ |  $(18.5,300)$
9 |  $(2.0,620)$ |  $(6.08,644.67)$ |  $(6.0,650)$
10 |  $(2.0,620)$ |  $(2.0,236.79)$ |  $(2.0,240)$
11 |  $(2.0,620)$ |  $(12.24,305.94)$ |  Stopped
Table III: Residual stress measurements, along with their computed mean and standard deviation at each iteration of the BayesOpt guided experiments.  $\mathbf{i}$ | Utilized |  Measurements in MPa |  mean |  stdev
---|---|---|---|---
1 | $(2.0,650)$ |  $-123.5,-134.8,-191$ |  $-149.8$ |  $36.2$
2 | $(2.0,430)$ |  $330.0,353.5,349.8$ |  $344.4$ |  $12.6$
3 | $(2.0,560)$ |  $97.9,40.1,105.5,$ $28.3,-15.3,3.2$ |  $43.3$ |  $49.2$
4 | $(2.0,670)$ |  $111.2,-239.6,-171.1,$ $-354,-411.1,-73.9,$ $-276.6,-309.3$ |  $-243.4$ |  $117.9$
5 | $(2.0,620)$ |  $-180.9,-218.0,-215.5$ |  $-204.8$ |  $20.7$
6 | $(2.0,320)$ |  $148.1,218.6,164.0$ |  $176.9$ |  $37.0$
7 | $(4.5,320)$ |  $1648.4,1660.3,1626.8$ |  $1645.2$ |  $17.0$
8 | $(18.5,300)$ |  $39.7,38.4,36.4$ |  $38.2$ |  $1.7$
9 | $(6.0,650)$ |  $1534.3,1549.20,1528.03$ |  $1537.2$ |  $10.9$
10 | $(2.0,240)$ |  $766.4,756.8,722.1$ |  $748.4$ |  $23.3$
Table IV: Sheet resistance measurements, along with their computed mean and standard deviation at each iteration of the BayesOpt guided experiments.  $\mathbf{i}$ | Utilized |  Measurements in $\Omega/\text{sq}$ |  mean |  stdev
---|---|---|---|---
1 | $(2.0,650)$ |  $0.87,0.86,0.85$ |  $0.86$ |  $0.01$
2 | $(2.0,430)$ |  $0.74,0.74,0.73$ |  $0.74$ |  $0.01$
3 | $(2.0,560)$ |  $0.64,0.64,0.63,$ $0.84,0.78,0.81$ |  $0.72$ |  $0.10$
4 | $(2.0,670)$ |  $0.72,0.66,0.70,$ $0.65,0.65,0.66,$ $0.74,0.72$ |  $0.69$ |  $0.04$
5 | $(2.0,620)$ |  $0.68,0.68,0.77$ |  $0.71$ |  $0.05$
6 | $(2.0,320)$ |  $1.44,1.44,1.42$ |  $1.43$ |  $0.01$
7 | $(4.5,320)$ |  $1.68,1.68,1.7$ |  $1.69$ |  $0.01$
8 | $(18.5,300)$ |  $8.04,7.9,8.14$ |  $8.03$ |  $0.12$
9 | $(6.0,650)$ |  $1.68,1.68,1.70$ |  $1.69$ |  $0.01$
10 | $(2.0,240)$ |  $0.81,0.82,0.85$ |  $0.83$ |  $0.02$

##  Appendix B Data from experiments for optimum verification

This appendix provides comprehensive information about the additional
experiments conducted for verification of the optimum solution obtained by
Bayesian optimization. The experiments were conducted with configurations near
the optimal configuration to observe the sensitivity of residual stress
dependence on pressure to process variations. For each configuration, three
films were deposited, and the residual stress and sheet resistance were
measured as shown in Tables V and VI, respectively.

Note that small differences in the stress measurements were observed at
$(2.0~{}\text{mTorr},~{}620~{}\text{W})$ as compared to the corresponding
measurements taken during the BayesOpt-guided experiments. Given the
significant time gap between experiments, several factors may have influenced
the discharge voltage of the sputter target. These include small target
impedance changes when the target is installed and re-installed, the number of
kWh the target has been used, and small pressure differences that are outside
what the capacitance manometer on the system can measure. Nevertheless, thin-
film stress and resistivity exhibit similar trends with power and pressure
despite minor deviations in absolute values. Additionally, all design criteria
continue to be satisfied at the optimum configuration proposed by BayesOpt.

Table V:  Measurements of residual stress, along with their computed mean and standard deviation during the additional experiments.  Utilized |  Measurements in MPa |  mean |  stdev
---|---|---|---
$(2.0,590)$ |  $-64.50,-84.73,-107.52$ |  $-85.58$ |  $21.52$
$(2.5,590)$ |  $382.83,409.49,414.19$ |  $402.17$ |  $16.91$
$(2.0,620)$ |  $-52.34,-35.65,-19.57$ |  $-43.99$ |  $16.39$
$(2.5,620)$ |  $140.26,175.44,151.19$ |  $155.63$ |  $18.01$
$(3.0,620)$ |  $679.36,691.20,641.25$ |  $670.60$ |  $26.10$
$(2.0,650)$ |  $-303.79,-254.90,-211.19$ |  $-256.63$ |  $46.32$
$(2.5,650)$ |  $292.57,314.78,295.83$ |  $301.06$ |  $11.99$
Table VI:  Measurements of sheet resistance, along with their computed mean and standard deviation during the additional experiments.  Utilized |  Measurements in $\Omega/\text{sq}$ |  mean |  stdev
---|---|---|---
$(2.0,590)$ |  $0.63,0.61,0.62$ |  $0.62$ |  $0.01$
$(2.5,590)$ |  $0.74,0.71,0.72$ |  $0.72$ |  $0.01$
$(2.0,620)$ |  $0.78,0.77,0.76$ |  $0.77$ |  $0.01$
$(2.5,620)$ |  $0.62,0.61,0.64$ |  $0.62$ |  $0.01$
$(3.0,620)$ |  $0.66,0.65,0.67$ |  $0.66$ |  $0.01$
$(2.0,650)$ |  $0.58,0.57,0.57$ |  $0.57$ |  $0.01$
$(2.5,650)$ |  $0.61,0.59,0.60$ |  $0.60$ |  $0.01$

## References

## References

  * Ashtiani (2023) K. Ashtiani, Semiconductor Digest. 5, 40 (2023).
  * Messica, Meirav, and Shtrikman (1995) A. Messica, U. Meirav, and H. Shtrikman, Thin Solid Films 257, 54–57 (1995).
  * Erofeev, Fedin, and Fedina (2018) E. V. Erofeev, I. V. Fedin, and V. V. Fedina, J. Phys.: Conf. Ser. 1115, 032068 (2018).
  * Yasrebi  _et al._ (2014) N. Yasrebi, B. Bagheri, P. Yazdanfar, B. Rashidian, and P. Sasanpour, Acta Metall. Sin. (Engl. Lett.) 27, 324–330 (2014).
  * Yücel, Güler, and Yücel (2014) E. Yücel, N. Güler, and Y. Yücel, J. Alloy. Compd. 589, 207–212 (2014).
  * Dreer, Krismer, and Wilhartitz (1999) S. Dreer, R. Krismer, and P. Wilhartitz, Surf. Coat. Technol. 114, 29–38 (1999).
  * Ramana  _et al._ (2012) G. V. Ramana, P. Saravanan, S. Kamat, and Y. Aparna, Appl. Surf. Sci. 261, 110–117 (2012).
  * Fisher (1936) R. A. Fisher, Brit. Med. J. 1, 554 (1936).
  * Madhava (1956) K. Madhava, Rev. Inst. Int. Stat. , 64–72 (1956).
  * Box and Wilson (1951) G. Box and K. Wilson, J. Roy. Stat. Soc. B. Met. 13, 1–45 (1951).
  * Taguchi (1985) G. Taguchi, Commun. Stat.-Theor. M. 14, 2785–2801 (1985).
  * Jones, Schonlau, and Welch (1998) D. R. Jones, M. Schonlau, and W. J. Welch, J. Global. Optim. 13, 455–492 (1998).
  * Packwood (2017) D. Packwood,  _Bayesian optimization for materials science_ (Springer, 2017).
  * Garnett (2023) R. Garnett,  _Bayesian optimization_ (Cambridge University Press, 2023).
  * Daulton  _et al._ (2022) S. Daulton, S. Cakmak, M. Balandat, M. A. Osborne, E. Zhou, and E. Bakshy, in  _International Conference on Machine Learning_ (PMLR, 2022) pp. 4831–4866.
  * Sanders  _et al._ (2019) N. D. Sanders, R. M. Everson, J. E. Fieldsend, and A. A. Rahat, arXiv preprint arXiv:1904.11416 (2019).
  * Le and Branke (2020) H. P. Le and J. Branke, in  _International Conference on Uncertainty Quantification & Optimisation_ (Springer, 2020) pp. 245–259.
  * Nguyen  _et al._ (2021) Q. P. Nguyen, Z. Dai, B. K. H. Low, and P. Jaillet, in  _International Conference on Machine Learning_ (PMLR, 2021) pp. 8063–8072.
  * Liang  _et al._ (2021) Q. Liang, A. E. Gongora, Z. Ren, A. Tiihonen, Z. Liu, S. Sun, J. R. Deneault, D. Bash, F. Mekki-Berrada, S. A. Khan, _et al._ , Npj Comput. Mater. 7, 188 (2021).
  * Lookman  _et al._ (2019) T. Lookman, P. V. Balachandran, D. Xue, and R. Yuan, Npj Comput. Mater. 5, 21 (2019).
  * Ohkubo  _et al._ (2021) I. Ohkubo, Z. Hou, J. Lee, T. Aizawa, M. Lippmaa, T. Chikyow, K. Tsuda, and T. Mori, Mater. Today Phys. 16, 100296 (2021).
  * Kusne  _et al._ (2020) A. G. Kusne, H. Yu, C. Wu, H. Zhang, J. Hattrick-Simpers, B. DeCost, S. Sarker, C. Oses, C. Toher, S. Curtarolo, _et al._ , Nat. Commun. 11, 5966 (2020).
  * Wakabayashi  _et al._ (2019) Y. K. Wakabayashi, T. Otsuka, Y. Krockenberger, H. Sawada, Y. Taniyasu, and H. Yamamoto, APL. Mater. 7 (2019).
  * Miyagawa  _et al._ (2021) S. Miyagawa, K. Gotoh, K. Kutsukake, Y. Kurokawa, and N. Usami, Sol. Energ. Mat. Sol. C. 230, 111251 (2021).
  * Gopakumar  _et al._ (2018) A. M. Gopakumar, P. V. Balachandran, D. Xue, J. E. Gubernatis, and T. Lookman, Sci. Rep. 8, 3738 (2018).
  * Khatamsaz  _et al._ (2022) D. Khatamsaz, B. Vela, P. Singh, D. D. Johnson, D. Allaire, and R. Arróyave, Acta Mater. 236, 118133 (2022).
  * Zhang, Apley, and Chen (2020) Y. Zhang, D. W. Apley, and W. Chen, Sci. Rep. 10, 4924 (2020).
  * Barbee, Mrowka, and Hettrick (1985) T. W. Barbee, S. Mrowka, and M. C. Hettrick, Appl. Optics 24, 883–886 (1985).
  * Mirkarimi (1999) P. B. Mirkarimi, Opt. Eng. 38, 1246–1259 (1999).
  * Orgassa, Schock, and Werner (2003) K. Orgassa, H. W. Schock, and J. Werner, Thin Solid Films 431, 387–391 (2003).
  * Li  _et al._ (2016) W. Li, X. Yan, A. G. Aberle, and S. Venkataraj, Procedia Engineer. 139, 1–6 (2016).
  * Abadias  _et al._ (2018) G. Abadias, E. Chason, J. Keckes, M. Sebastiani, G. B. Thompson, E. Barthel, G. L. Doll, C. E. Murray, C. H. Stoessel, and L. Martinu, J Vac. Sci. Technol. A 36 (2018).
  * Gramacy  _et al._ (2016) R. B. Gramacy, G. A. Gray, S. Le Digabel, H. K. Lee, P. Ranjan, G. Wells, and S. M. Wild, Technometrics 58, 1–11 (2016).
  * Gardner  _et al._ (2014) J. R. Gardner, M. J. Kusner, Z. E. Xu, K. Q. Weinberger, and J. P. Cunningham, ICML , 937–945 (2014).
  * Feliot, Bect, and Vazquez (2017) P. Feliot, J. Bect, and E. Vazquez, Journal of Global Optimization 67, 97–133 (2017).
  * Hern  _et al._ (2016) J. M. Hern, M. A. Gelbart, R. P. Adams, M. W. Hoffman, Z. Ghahramani, _et al._ , J. Mach. Learn. Res. 17, 1–53 (2016).
  * Pourmohamad and Lee (2020) T. Pourmohamad and H. K. Lee, Technometrics 62, 303–312 (2020).
  * Sainvitu and Toint (2007) C. Sainvitu and P. L. Toint, Optim. Methods Softw. 22, 835–848 (2007).
  * Rasmussen and Williams (2005) C. E. Rasmussen and C. K. I. Williams,  _Gaussian Processes for Machine Learning_ (The MIT Press, 2005).
  * Stoney (1909) G. G. Stoney, Proc. R. Soc. Lond. 82, 172–175 (1909).
  * Kalaswad  _et al._ (2023) M. Kalaswad, J. O. Custer, S. Addamane, R. M. Khan, L. Jauregui, T. F. Babuska, A. Henriksen, F. W. DelRio, R. Dingreville, B. L. Boyce, and D. P. Adams, Integr. Mater. Manuf. Innov. 12, 118–129 (2023).
  * Mattox (1998) D. M. Mattox, “Preparation and cleaning of vacuum surfaces,” in  _Handbook of Vacuum Science and Technology_ (Elsevier, 1998) pp. 553–606.