[table]capposition=top

# On the Effects of Heterogeneous Errors on Multi-fidelity Bayesian
Optimization

Zahra Zanjani Foumani  Department of Mechanical and Aerospace Engineering,
University of California, Irvine  Amin Yousefpour  Department of Mechanical
and Aerospace Engineering, University of California, Irvine  Mehdi Shishehbor
Department of Mechanical and Aerospace Engineering, University of California,
Irvine  Ramin Bostanabad Corresponding Author<EMAIL_ADDRESS>Department of
Mechanical and Aerospace Engineering, University of California, Irvine

Abstract
Bayesian optimization (BO) is a sequential optimization strategy that is
increasingly employed in a wide range of areas including materials design. In
real world applications, acquiring high-fidelity (HF) data through physical
experiments or HF simulations is the major cost component of BO. To alleviate
this bottleneck, multi-fidelity (MF) methods are used to forgo the sole
reliance on the expensive HF data and reduce the sampling costs by querying
inexpensive low-fidelity (LF) sources whose data are correlated with HF
samples. However, existing multi-fidelity BO (MFBO) methods operate under the
following two assumptions that rarely hold in practical applications: $(1)$ LF
sources provide data that are well correlated with the HF data on a global
scale, and $(2)$ a single random process can model the noise in the fused
data. These assumptions dramatically reduce the performance of MFBO when LF
sources are only locally correlated with the HF source or when the noise
variance varies across the data sources. In this paper, we dispense with these
incorrect assumptions by proposing an MF emulation method that $(1)$ learns a
noise model for each data source, and $(2)$ enables MFBO to leverage highly
biased LF sources which are only locally correlated with the HF source. We
illustrate the performance of our method through analytical examples and
engineering problems on materials design.

Keywords: Bayesian optimization; multi-fidelity modeling; emulation; Gaussian
process; interval score.

##  1 Introduction

Bayesian optimization (BO) is a sequential and sample-efficient global
optimization technique that is increasingly used in the optimization of
expensive-to-evaluate (and typically black-box) functions [1]. BO has two main
ingredients: an emulator which is typically a Gaussian process (GP) and an
acquisition function (AF) [2]. The first step in BO is to train an emulator on
some initial data. Then, an auxiliary optimization is solved to determine the
new sample that should be added to the training data. The objective function
of this auxiliary optimization is the AF whose evaluation relies on the
emulator. Given the new sample, the training data is updated and the entire
emulation-sampling process is repeated until the convergence conditions are
met [3].

Although BO is a highly efficient technique, the total cost of optimization
can be substantial if it solely relies on the accurate but expensive high-
fidelity (HF) data source. To mitigate this issue, multi-fidelity (MF)
techniques are widely adopted [4, 5, 6] where one uses multiple data sources
of varying levels of accuracy and cost in BO. The fundamental principle behind
MF techniques is to exploit the correlation between low-fidelity (LF) and HF
data to decrease the overall sampling costs [7, 8].

Over the past two decades many multi-fidelity BO (MFBO) strategies have been
proposed which primarily differ in terms of their emulator and AF. Most
existing strategies rely on the Co-Kriging method [9], Kennedy and O’Hagan’s
bi-fidelity approach [10], and the BoTorch package [11]. These MFBO methods
have some major drawbacks such as inability to simultaneously leverage
multiple LF sources, sensitivity to the sampling costs (where highly
inexpensive LF sources are heavily sampled which, in turn, causes numerical
and convergence issues), and presuming simple bias forms (e.g., an additive
function [10]) for the LF sources.

Some of these limitations are recently addressed in [12] where the authors
propose to (1) use latent map Gaussian processes (LMGPs) for emulation, and
(2) quantify the information value of LF and HF samples differently. Their AF
is cost-aware in that it considers the sampling cost in quantifying the value
of HF and LF data points. Henceforth, we refer to this method as MFBO.

(a) MFBO

(b) Fidelity Manifold

(c) MFBOUQ

(d) Convergence History

Figure 1: Effect of heterogeneous noise and model fidelity on MFBO: HF data
are noisy ($\sigma_{noise}=1$) and expensive while the LF data are
deterministic and cheap. In this example, LF1 is more correlated with the HF
source for $x<5$ while LF2 correlates better with HF data for $x>5$. The
sampling cost of the HF and two LF sources are $10/1/1$, respectively. (a)
demonstrates the sampling history of MFBO which excludes LF2 from the sampling
process based on the fidelity manifold of LMGP shown in (b) which indicates
that LF2 is globally more correlated with the HF data compared to LF1. The
emulator in (a) learns a single noise process for LF2 and HF data and the
manifold in (b) is learnt via the initial data. (c) MFBOUQ is the approach we
propose in this paper effectively explores the space (as it samples the HF
source more in $x>5$) and leverages LF1 in $x<5$. As for LF2, MFBOUQ mostly
samples from $x>5$ since this region includes the optimum of LF2 and is more
correlated with the HF. (d) MFBOUQ outperforms MFBO in finding the optimum of
HF source ($y^{*}$) for various initial conditions (the large noise variance
of HF data cause both approaches to have some errors upon convergence).
Initial data are not shown (a) and (c).

While MFBO performs much better than competing MF approaches, it has two main
limitations which are demonstrated with a simple $1D$ example in Figure 1
where each of the two LF sources is more correlated with the HF function in
half of the domain, see Figure 1a. Firstly, MFBO excludes highly biased LF
sources from BO with the rationale that they can steer the search process in
the wrong direction. This exclusion is done before BO starts since the
decision is made based on the fidelity manifold of LMGP that is trained on the
initial data, see Figure 1b. In this manifold each data source is encoded with
a point and hence distances between these points correspond to global
correlations between the corresponding data sources. That is, Figure 1b
suggests that the HF source is barely correlated with the LF sources even
though they are close to the HF function in half of the domain (since
eliminating both sources convert the BO into a single-fidelity process in this
example, we assume that MFBO only excludes LF2 and hence LF1 and the HF
sources are sampled via MFBO in Figure 1a). This is obviously a sub-optimal
decision as it precludes the possibility of leveraging an LF source that is
valuable in a small portion of the search space which may include the global
optimum of the HF source. That is, in the example of Figure 1a MFBO should
ideally leverage LF2 but mostly sample from it in the $x>5$ region.

The second limitation of MFBO is that it assumes all sources are corrupted
with the same noise process (with unknown noise variance). However, MF
datasets typically have different levels of noise especially if some sources
represent deterministic computer simulations while others are physical
experiments [13, 14]. In such applications, MFBO overestimates the
uncertainties associated with the noise-free data sources which, in turn,
reduces the performance of MFBO.

To address these two limitations, we introduce MFBOUQ for multi-fidelity cost-
aware Bayesian optimization. MFBOUQ has the same AFs as MFBO and can leverage
an arbitrary number of LF sources in optimizing an HF source. Unlike MFBO,
MFBOUQ never discards an LF source (regardless of its bias with respect to the
HF source) and estimates a noise process for each data source. We argue that
MFBOUQ quantifies the uncertainties more accurately than MFBO and thus
achieves a higher performance in MFBO. Figure 1c schematically demonstrates
the advantages of MFBOUQ over MFBO in a $1D$ example where there are one HF
and two LF sources. As it can be observed the LF sources are mostly sampled
where they $(1)$ are well correlated with the HF source, or $(2)$ provide
attractive function values (e.g., very small values in minimization). The
advantages of MFBOUQ over MFBO in finding the optimum of HF ($y^{*}$) hold
over various initializations, see Figure 1d.

The rest of the paper is organized as follows. We provide the methodological
details in Section 2 and then evaluate the performance of MFBOUQ via multiple
ablation studies in Section 3. In Section 3, we also visualize how strategic
sampling in MFBOUQ, driven by accurate uncertainty quantification and
effective handling of biased data sources, results in its superior performance
compared to MFBO. This is demonstrated through two real-world high-dimensional
material design examples with noisy and highly biased data sources. We
conclude the paper in Section 4 by summarizing our contributions and providing
future research directions.

##  2 Methods

In this section, we first provide some background on LMGP and MF modeling with
LMGP in Section 2.1 and Section 2.2, respectively. We then propose our
efficient mechanism for inversely learning a noise process for each data
source in Section 2.3. Next, we introduce the cost-aware AF of MFBOUQ in
Section 2.4. Finally, in Section 2.5 we elaborate on our idea that improves
the uncertainty quantification (UQ) capabilities of LMGPs and, in turn,
benefits MFBO.

###  2.1 Latent Map Gaussian Process (LMGP)

Gaussian processes (GPs) are emulators which assume the responses or outputs
in the training data come from a multivariate normal distribution with
parametric mean and covariance functions that depend on the inputs. Based on
this assumption, the following equation can be written:

$\begin{split}y(\boldsymbol{x})=\beta+\xi(\boldsymbol{x})\end{split}$ (1)

where $\boldsymbol{x}=[x_{1},x_{2},\dots,x_{dx}]^{T}$ is the input vector,
$y(\boldsymbol{x})$ is the output, $\beta$ is an unknown coefficient, and
$\xi(\boldsymbol{x})$ is a zero-mean GP with the covariance function:

$\begin{split}cov(\ \xi(\boldsymbol{x}),\
\xi(\boldsymbol{x}^{\prime}))=c(\boldsymbol{x},\boldsymbol{x}^{\prime})=\
\sigma^{2}r(\boldsymbol{x},\boldsymbol{x}^{\prime})\end{split}$ (2)

where $\sigma^{2}$ is the variance of the process and $r(\cdot,\cdot)$ is the
parametric correlation function which measures the distance between any two
input vectors. In this paper we use the Gaussian correlation function defined
as:

$\begin{split}r(\boldsymbol{x},\boldsymbol{x}^{\prime})=\exp{\\{-\sum_{i=1}^{dx}{10}^{\omega_{i}}\
(x_{i}-\ x_{i}^{\prime})^{2}\\}}\end{split}$ (3)

where $\boldsymbol{\omega}=[\omega_{1},\omega_{2},\dots,\omega_{dx}]^{T}$ are
the scale parameters. The versatility of a GP emulator highly depends on its
correlation function in that, e.g., traditional GPs do not accommodate
categorical variables directly since kernels like the one in Equation 3 cannot
compare qualitative features.

To directly use GPs in MF modeling, we follow [15] who convert MF modeling to
a manifold learning problem via LMGPs which are extensions of GPs that can
handle categorical data [16] while providing a visualizable manifold that can
be used to interpret the correlation among data sources.

Denoting the categorical inputs by
$\boldsymbol{t}=[t_{1},t_{2},\dots,t_{dt}]^{T}$ where variable $t_{i}$ has
$l_{i}$ distinct levels, LMGP maps each combination of the categorical levels
to a point in a learned quantitative manifold. To this end, LMGP assigns a
unique vector to each combination of the categorical variables and then uses a
parametric function to map these unique vectors into a compact manifold with
dimensionality $dz$. Assuming a linear transformation is used in LMGP, the
mapping operation reads as:

$\begin{split}\boldsymbol{z}(\boldsymbol{t})=\boldsymbol{\zeta}(\boldsymbol{t})\boldsymbol{A}\end{split}$
(4)

where $\boldsymbol{t}$ denotes a specific combination of the categorical
variables, $\boldsymbol{z}(\boldsymbol{t})$ is the $1\times dz$ posterior
latent representation of $\boldsymbol{t}$,
$\boldsymbol{\zeta}(\boldsymbol{t})$ is a unique prior vector representation
of $\boldsymbol{t}$, and $\boldsymbol{A}$ is a rectangular matrix that maps
$\boldsymbol{\zeta}(\boldsymbol{t})$ to $\boldsymbol{z}(\boldsymbol{t})$. In
this paper, grouped one-hot encoding is used to generate the prior vectors and
hence the dimensionality of $\boldsymbol{\zeta}(\boldsymbol{t})$ and
$\boldsymbol{A}$ are $1\times\sum_{i=1}^{dt}l_{i}$ and
$\sum_{i=1}^{dt}l_{i}\times dz$, respectively. These mapped points can now be
directly embedded in the correlation function as:

$\begin{split}r(\boldsymbol{u},\boldsymbol{u}^{\prime})=\exp{\\{-\sum_{i=1}^{dx}{10}^{\omega_{i}}\
(x_{i}-\
x_{i}^{\prime})^{2}\\}}\sum\\{-\sum_{i=1}^{dz}(z_{i}(\boldsymbol{t})-\
z_{i}(\boldsymbol{t}^{\prime}))^{2}\\}\end{split}$ (5)

where $\boldsymbol{u}=[\boldsymbol{x};\ \boldsymbol{t}]$ and
$\boldsymbol{z}(\boldsymbol{t})=[z_{1}(\boldsymbol{t}),\
z_{2}(\boldsymbol{t}),\ \ldots,\ z_{dz}(\boldsymbol{t})]$ is the location in
the learned latent space corresponding to the specific combination of the
categorical variables denoted by $\boldsymbol{t}$.

LMGP estimates the hyperparameters $(\beta,\ \boldsymbol{A},\
\boldsymbol{\omega},\ \sigma^{2})$ via maximum a posteriori (MAP) which,
assuming $dz=2$, provides point-estimates for $dx+2\
\times\sum_{i=1}^{dt}l_{i}$ variables. Upon parameter estimation, LMGP uses
the conditional distribution formulas to predict the response distribution at
the arbitrary point $\boldsymbol{u}$ with the following mean and variance:

$\begin{split}\mathbb{E}[y(\boldsymbol{u})]=\mu(\boldsymbol{u})=\
\hat{\beta}+\ \boldsymbol{r}^{T}(\boldsymbol{u})\
\boldsymbol{R}^{-1}(\boldsymbol{y}-\boldsymbol{1}_{n\times
1}\hat{\beta})\end{split}$ (6)
$\begin{split}c(y(\boldsymbol{u}),y(\boldsymbol{u}))=\
\sigma^{2}(\boldsymbol{u})={\hat{\sigma}}^{2}(1-\boldsymbol{r}^{T}(\boldsymbol{u})\boldsymbol{R}^{-1}\boldsymbol{r}(\boldsymbol{u})+{(g(\boldsymbol{u}))}^{2}({\boldsymbol{1}_{1\times
n}\boldsymbol{R}}^{-1}\boldsymbol{1}_{n\times 1})^{-1})\end{split}$ (7)

where $n$ is number of training samples, $\mathbb{E}$ denotes expectation,
$\boldsymbol{1}_{a\times b}$ is an $a\times b$ matrix of ones,
$\boldsymbol{r}(\boldsymbol{u})$ is an $n\times 1$ vector with the $i^{th}$
element $r(\boldsymbol{u}^{i},\boldsymbol{u})$, $\boldsymbol{R}$ is an
$n\times n$ matrix with $R_{ij}=r(\boldsymbol{u}^{i},\boldsymbol{u}^{j})$, and
$g(\boldsymbol{u})=1-\mathbf{1}_{1\times
n}\boldsymbol{R}^{-1}\boldsymbol{r}(\boldsymbol{u})$.

###  2.2  Multi-fidelity Emulation via LMGP

The first step to MF emulation with LMGP is to augment the inputs with the
additional categorical variable $s$ that indicates the source of a sample,
i.e., $s=\\{^{\prime}1^{\prime},^{\prime}2^{\prime},\
\ldots,^{\prime}ds^{\prime}\\}$ where the $j^{th}$ element corresponds to
source $j$ for $j=1,\ldots,ds$. Subsequently, the training data from all
sources are concatenated and used in LMGP to build an MF emulator. Upon
training, to predict the objective value of a point $\boldsymbol{x}$ from
source $j$, $\boldsymbol{x}$ is concatenated with the categorical variable $s$
that corresponds to source $j$ and fed into the trained LMGP. We refer the
readers to [12] for more detail but note here that in case the input variables
already contain some categorical features (see Section 3.2 for an example), we
endow LMGP with two manifolds where one encodes the fidelity variable $s$
while the other manifold encodes the rest of the categorical variables. While
this choice does not noticeably affec the accuracy of LMGP during test time,
it increase interpretability. For instance, we use the learned manifold for
the categorical variables in Section 3.2 to show the trajectory of BO in the
design space.

It has been recently shown [16] that LMGPs have the following primary
advantages over other MF emulators: (1) they provide a more flexible and
accurate mechanism to build MF emulators since they learn the relations
between the sources in a nonlinear manifold, (2) they learn all the sources
quite accurately rather than just emulating the HF source, and (3) they
provide a visualizable global metric for comparing the relative
discrepancies/similarities among the data sources.

###  2.3 Source-dependent Noise Modeling

The presence of noise significantly affects the performance of BO and
incorrectly modeling it can cause over-exploration or under-exploration of the
search space. To mitigate the effects of noise in BO, we reformulate LMGPs to
independently model a noise process for each data source. This reformulation
can improve the accuracy of the model in noisy regions and, in turn, guide the
search toward the global optimum when the modeled is deployed in MFBO.

To model noise in GPs, the nugget or jitter parameter, $\delta$, is used [17]
to replace $\boldsymbol{R}$ with
$\boldsymbol{R}_{\delta}=\boldsymbol{R}+\delta\boldsymbol{I}$ where
$\boldsymbol{I}$ is an $n\times n$ identity matrix. With this approach, the
estimated stationary noise variance in the data is $\delta\sigma^{2}$ and the
mean and variance formulations in Eq. 6 and Eq. 7 are modified by using
$\boldsymbol{R}_{\delta}$ instead of $\boldsymbol{R}$.

Although incorporating this modification in the correlation matrix can enhance
the performance of the emulator and BO in single-fidelity (SF) problems, it
does not yield the same benefits in MF optimization. This is likely because of
the dissimilar nature of the data sources and their corresponding noises. When
dealing with multiple sources of data, each source may suffer from different
levels and types of noise. Consider a bi-fidelity dataset where the HF data
comes from an experimental setup and is subject to measurement noise, while
the LF data is generated by a deterministic computer code which has a
systematic bias due to missing physics. In this case, using only one nugget
parameter in LMGP for MF emulation is obviously not an optimum choice.

To address this issue effectively, we propose to use multiple nugget
parameters in the emulator. Specifically, we define the nugget vector
$\boldsymbol{\delta}=[\delta_{1},\delta_{2},\dots,\delta_{ds}]$ and update the
correlation matrix as follows:

$\begin{split}\boldsymbol{R}_{\delta}=\boldsymbol{R}+\boldsymbol{N}_{\delta}\end{split}$
(8)

where $\boldsymbol{N}_{\delta}$ denotes an $n\times n$ diagonal matrix whose
$(i,i)^{th}$ element is the nugget element corresponding to the data source of
the $i^{th}$ sample. For instance, suppose the $i^{th}$ sample
$(\boldsymbol{u}^{i})$ is generated by source $ds$. Then, $(i,i)^{th}$ element
of $\boldsymbol{N}_{\delta}$ is $\delta_{ds}$. Then, we use Eq. 8 to build the
correlation matrix of LMGP and jointly estimate all the parameters via MAP as:

$\begin{split}[\widehat{\beta},\hat{\sigma},\widehat{\boldsymbol{\omega}},\widehat{\boldsymbol{A}},\widehat{\boldsymbol{\delta}}]=\underset{\beta,\sigma^{2},{\boldsymbol{\omega}},\boldsymbol{A}}{\operatorname{argmin}}\frac{n}{2}\log(\sigma^{2})+\frac{1}{2}\log(|\boldsymbol{R}_{\delta}|)+\frac{1}{2\sigma^{2}}(\boldsymbol{y}-\boldsymbol{M}\beta)^{T}\boldsymbol{R}_{\delta}^{-1}(\boldsymbol{y}-\boldsymbol{M}\beta){+\log(\underset{\beta,\sigma^{2},{\boldsymbol{\omega}},\boldsymbol{A},\boldsymbol{\delta}}{\operatorname{P(\cdot)}})}\end{split}$
(9)

where $p(\cdot)$ is the prior of the hyperparameters. We define independent
priors for each parameter where $\omega^{i}\sim N(-3,3)$, $\beta\sim N(0,1)$,
$A^{ij}\sim N(0,3)$, $\sigma\sim LN(0,3)$111Log-Normal , and $\delta^{i}\sim
LHS(0,0.01)$222Log-Half-Horseshoe with zero lower bound and scale parameter
$0.01$. [18]. Our multi-noise approach increases the number of LMGP’s
hyperparameters to $dx+2\times\sum_{i=1}^{dt}l_{i}+ds$.

###  2.4 Multi-source Cost-aware Acquisition Function

The choice of AF is crucial in MFBO since it must consider the biases of LF
data and source-dependent sampling costs in addition to balancing exploration
and exploitation. To capture these goals, separate AFs are defined in [12] for
LF and HF sources with a focus on exploration and exploitation, respectively.

Following the idea of proposing an AF with a focus on exploration for the LF
sources, the AF of the $j^{th}$ LF source ($j\neq l$, $l$ denotes the HF
source) is defined as the exploration part of the expected improvement (EI) in
MFBO:

$\begin{split}\gamma_{LF}(\boldsymbol{u};j)=\sigma_{j}(\boldsymbol{u})\phi(\frac{y_{j}^{*}-\mu_{j}(\boldsymbol{u})}{\sigma_{j}(\boldsymbol{u})})\end{split}$
(10)

where $y_{j}^{*}$ is the best function value obtained so far from source $j$
and $\phi(\cdot)$ denotes the probability density function (PDF) of the
standard normal variable. $\sigma_{j}(\boldsymbol{u})$ and
$\mu_{j}(\boldsymbol{u})$ are the standard deviation and mean, respectively,
of point $\boldsymbol{u}$ from source $j$ which we estimate via :

$\begin{split}\mu(\boldsymbol{u})=\hat{\beta}+\boldsymbol{r}^{T}(\boldsymbol{u})\boldsymbol{R}_{\delta}^{-1}(\boldsymbol{y}-\mathbf{1}_{n\times
1}\hat{\beta})\end{split}$ (11)
$\begin{split}c(y(\boldsymbol{u}),y(\boldsymbol{u}))=&\sigma^{2}(\boldsymbol{u})=\hat{\sigma}^{2}(1-\boldsymbol{r}^{T}(\boldsymbol{u})\boldsymbol{R}_{\delta}^{-1}\boldsymbol{r}(\boldsymbol{u})+(g(\boldsymbol{u}))^{2}(\mathbf{1}_{1\times
n}\boldsymbol{R}_{\delta}^{-1}\mathbf{1}_{n\times
1})^{-1})+\hat{\delta}_{j}\end{split}$ (12)

where $g(\boldsymbol{u})=1-\mathbf{1}_{1\times
n}\boldsymbol{R}_{\delta}^{-1}\boldsymbol{r}(\boldsymbol{u})$ and
$\hat{\delta}_{j}$ is the estimated nugget parameter for source $j$.

MFBO utilizes improvement as the AF for the HF data source, since it is
computationally efficient and emphasizes exploitation. Accordingly, MFBOUQ
uses improvement for the HF source (source $l$) with the new mean calculated
based on the Eq. 11:

$\begin{split}\gamma_{HF}(\boldsymbol{u};l)=\mu_{l}(\boldsymbol{u})-y_{l}^{*}\end{split}$
(13)

In each iteration of BO, we first use the mentioned AFs to solve $ds$
auxiliary optimizations to find the candidate points with the highest
acquisition value from each source. We then scale these values by the
corresponding sampling costs to obtain the following composite AF :

$\begin{split}\gamma_{\texttt{MFBO\textsubscript{UQ}}}(\boldsymbol{u};j)=\begin{cases}\nicefrac{{\gamma_{LF}(\boldsymbol{u};j)}}{{O(j)}}&j=[1,\cdots,ds]\And
j\neq l\\\
\nicefrac{{\gamma_{HF}(\boldsymbol{u};l)}}{{O(l)}}&j=l\end{cases}\end{split}$
(14)

where $O(j)$ is the cost of acquiring one sample from source $j$. We determine
the final candidate point (and the source that it should be sampled from) at
for iteration $k+1$ via:

$\begin{split}[\boldsymbol{u}^{k+1},j^{k+1}]=\underset{\boldsymbol{u},j}{\operatorname{argmax}}\leavevmode\nobreak\
\gamma_{\texttt{MFBO\textsubscript{UQ}}}(\boldsymbol{u};j)\end{split}$ (15)

###  2.5 Emulation for Exploration

The composite AF in Eq.14 quantifies the information value of LF samples via
Eq.10 whose value scales with the prediction uncertainties, i.e.,
$\sigma(\boldsymbol{u})$. The source-dependent noise modeling of Section 2.3
improves LMGP’s ability in learning the uncertainty by introducing a few more
hyperparameters. However, the added hyperparameters may result in overfitting
and, in turn, deteriorate the predicted uncertainties [19, 20]. A related
issue is the effect of large local biases of LF sources which can inflate the
uncertainty quite substantially and, as a result, increase
$\gamma_{LF}(\boldsymbol{u};j)$. This increase causes MFBO to repeatedly
sample from the biased LF sources. Such repeated samplings reduce the
efficiency of MFBO and may cause numerical issues (due to ill-conditioning of
the covariance matrix) or even convergence to a suboptimal solution.

To address the above issues simultaneously, we argue that the training process
of the emulator should increase the importance of UQ which directly affects
the exploration part of MFBO. To this end, we leverage strictly proper scoring
rules while training LMGPs.

Scoring rules [21] evaluate a probabilistic prediction by assigning a
numerical score to it. The scoring rule of an emulator is (strictly) proper if
matching the predicted distribution with the underlying sample distribution
(uniquely) maximizes the expected score for any sample. The probabilistic
nature of LMGP’s prediction motivates us to use the negatively oriented
interval score (hereafter denoted by $IS$) to evaluate the UQ capabilities of
LMGPs. We choose $IS$ since it is robust to outliers, rewards narrow
prediction intervals, and is flexible in the choice of desired coverage levels
[22, 23].

$IS$ is a special case of quantile prediction that penalizes the model for
each observation that is not inside the $(1-v)\times 100\%$ prediction
interval. The lower $(\mathcal{L}^{i})$ and upper $\mathcal{U}^{i}$ endpoints
of this prediction interval for the $i^{th}$observation are their predictive
quantiles at levels $v/2$ and $1-v/2$, respectively. So, we calculate the $IS$
as:

$\begin{split}IS_{v}=\frac{1}{n}\sum_{i=1}^{n}(\mathcal{U}^{i}-\mathcal{L}^{i})+\frac{2}{v}(\mathcal{L}^{i}-y(\boldsymbol{u}^{i}))\mathbbm{1}\leavevmode\nobreak\
\\{y(\boldsymbol{u}^{i})<\mathcal{L}^{i}\\}+\frac{2}{v}(y(\boldsymbol{u}^{i})-\mathcal{U}^{i})\mathbbm{1}\\{y(\boldsymbol{u}^{i})>\mathcal{U}^{i}\\}\end{split}$
(16)

where $\mathbbm{1}\\{\cdot\\}$ is an indicator function which is $1$ if its
condition holds and zero otherwise [24, 25]. We use $v=0.05$ ($95\%$
prediction interval), so
$\mathcal{U}^{i}=\mu{(\boldsymbol{u}^{i})}+1.96\sigma{(\boldsymbol{u}^{i})}$
and
$\mathcal{L}^{i}=\mu{(\boldsymbol{u}^{i})}-1.96\sigma{(\boldsymbol{u}^{i})}$.

Having defined the $IS$, we now formulate the new objective function for
training LMGPs where $IS_{0.05}$ is used as a penalty term during
hyperparameter estimation to increase the focus on UQ. Since the effectiveness
of this penalization mechanism depends on the value of the posterior, we
introduce an adaptive coefficient whose magnitude depends on the posterior
value. With this penalty term, we estimate the hyperparameters of LMGP via:

$\begin{split}{[\hat{\beta},\widehat{\sigma},\widehat{\boldsymbol{\omega}},\widehat{\boldsymbol{A}},\widehat{\boldsymbol{\delta}}]=}\underset{\beta,\sigma,\boldsymbol{\omega},\boldsymbol{A},\boldsymbol{\delta}}{\operatorname{argmin}}L_{MAP}+\varepsilon|L_{MAP}|\times
IS_{0.05}\end{split}$ (17)

where $|\cdot|$ denotes the absolute function and $\varepsilon$ is a user-
defined scaling parameter. In this paper, we use $\varepsilon=0.08$ for all of
our examples.

##  3 Results and Discussion

We demonstrate the performance of MFBOUQ on two analytic examples (see Table 1
for details on functional forms, sampling costs, number of LF sources and
their accuracy with respect to the HF source, and size of initial data) and
two real-world problems. In each case, we compare the results against those of
MFBO and single-fidelity BO (SFBO). While SFBO uses EI as its AF, MFBOUQ and
MFBO use the AFs introduced in Section 2.4.

We assume that the cost of querying any of the data sources is much higher
than the computational costs of BO (i.e., fitting LMGP and solving the
auxiliary optimization problem). Therefore, we compare the methods based on
their capability to identify the global optimum of the HF source and the
overall data collection cost. By comparing these methods, we aim to
demonstrate: $(1)$ the advantages of estimating noise process for each data
source, $(2)$ that using IS improves the accuracy of LMGP and, in turn,
enhances the convergence of BO (since our defined AFs highly rely on the
quality of the prediction), and $(3)$ that deploying IS eliminates the need
for excluding highly biased fidelity sources from BO.

We use the same stop conditions across the three methods to clearly
demonstrate the benefits of our two contributions. In particular, the
optimization is stopped when either of the following happens: $(1)$ the
overall sampling cost exceeds a pre-determined maximum budget, or $(2)$ the
best HF sample does not change over $50$ iterations. The maximum budget for
the analytical examples is $40000$ units, while it is $1000$ and $1800$ for
the two real-world examples. These budgets are chosen based on the data
collection costs.

###  3.1 Analytical Examples

We consider two analytical examples, Wingand Borehole, whose input
dimensionality is $10$ and $8$, respectively. To challenge the convergence and
better illustrate the power of separate noise estimation, we only add noise to
the HF data (the noise variance is defined based on the range of each
function.). The added noise variance to the HF source of Wing and Borehole are
$9$ and $16$, respectively. Both examples are single response and details
regarding their formulation, initialization, and sampling cost is presented in
Table 1). To assess the robustness of the results and quantify the effect of
random initial data, we repeat the optimization process $20$ times for each
example with each of the three methods (all initial data are generated via
Sobol sequence).

In each example, the relative root mean squared error (RRMSE) is calculated
between LF sources and their corresponding HF source based on $1000$ samples
to show the relative accuracy of the LF sources (presented in Table 1). Based
on these ground truth numbers (which are not used in BO), in the case of
Boreholethe source ID, true fidelity level (based on the RRMSEs), and sampling
costs are not related (e.g., although the first LF source is the most
expensive one, it has the least accuracy compared to the HF source). In the
case of Wing, however, these numbers match (e.g., LF1 is the most accurate and
expensive LF source and is followed by LF2 and then LF3).

MFBO excludes the highly biased LF sources from BO before any new samples are
obtained (also, during BO, the initial samples from highly biased LF sources
are not used in emulation). This exclusion is done based on the latent map of
the LMGP model that is trained on the initial data. Figure 2 shows the latent
maps of Wing and Borehole examples. As shown in Figure 2, while all the
fidelity sources of Wing are beneficial (since the points encoding the LF
sources are very close to the HF point), the first two LF sources of Borehole
are not correlated enough with the HF (their latent positions are distant from
that of the HF) and hence are excluded in MFBO. However, MFBOUQ does not
require this exclusion because it leverages the biased LF sources merely in
the regions that they are correlated with the HF source. In this paper, we do
not exclude the biased sources in MFBO to better compare it with our proposed
method.

(a) Wing

(b) Borehole

Figure 2: Fidelity manifolds of analytic examples: The plots in (a) and (b)
are obtained by fitting an LMGP to the initial data in the Wing and Borehole
examples, respectively. Due to the consistency across the $20$ repetitions,
the plots are randomly chosen among them. In (b), the HF source is encoded far
from LF1 and LF2 which indicates that these two sources have large biases with
respect to the HF source. MFBO excludes these two sources from the BO while
MFBOUQ does not.

Figure 3 summarizes the convergence history of each example by depicting the
best HF sample ($y^{*}$) found by each method versus its accumulated sampling
cost. As we expect, MF methods (MFBOUQ and MFBO) outperform SFBO in Wing
(Figure 3a) by leveraging the inexpensive LF sources that are globally
correlated with the HF source. However, the superior performance of MFBOUQ is
more obvious in Borehole with biased data sources.

In Borehole (Figure 3b), all the thin red curves (MFBO) are straight lines,
except for two curves. This means that for $18$ repetitions, the optimization
process fails to improve. The reason behind this lack of improvement is that
MFBO cannot handle local correlation of the LF sources and samples points that
steer the optimization in the wrong direction. Consequently, MFBO cannot find
any efficient HF sample with large enough information value (that compensates
for its high sampling cost) which results in the lack of improvement.
Conversely, all the thin green curves (MFBOUQ) converge to a value very close
to the ground truth. In addition, MFBOUQ yields almost the same convergence
value as SFBO, but with lower overall computational cost. This instance
further demonstrates the effectiveness of our proposed AFs since SFBO is very
accurate (albeit more expensive) due to only sampling from HF source and not
dealing with the local biases of the LF sources.

(a) Wing

(b) Borehole

Figure 3: Convergence histories of analytic examples: The plots depict the
best HF sample found by each approach ($y^{*}$) versus their sampling costs
accumulated during the BO iterations (the cost of initial data is included).
(a) and (b) summarize the results for the Wing and Borehole examples,
respectively. The thin curves show the convergence history of each repetition
and the solid thick ones indicate the average behavior across the $20$
repetitions. In both examples, MFBOUQ outperforms MFBO in terms of both
convergence value and cost. In both examples, SFBO performs the worst. The
ground truth is represented by the black dashed line.

###  3.2 Real-world Datasets

In this section, we study two materials design problems where the aim is to
find the composition that best optimizes the property of interest. We do not
add noise to these two examples as they are inherently noisy. The design space
of both examples has categorical inputs (denoted by $\boldsymbol{t}$) and we
add one more categorical variable (denoted by $s$) to enable data fusion as
described in Section 2.2. We design our LMGP to map the categorical inputs
onto two $2D$ manifolds (one for $\boldsymbol{t}$ and the other for $s$) to
help with the visualization of the exploration-exploitation behavior of BO in
the design space. The HF and LF data are obtained via simulations (based on
the density functional theory or DFT) with different fidelity levels.

The first problem is bi-fidelity where the goal is to find the member of the
nanolaminate ternary alloy (NTA) family with the largest bulk modulus [26].
The HF and LF datasets have $224$ samples each and are $10$-dimensional ($7$
quantitative and $3$ categorical where the latter have $10$, $12$, and $2$
levels). The cost ratio between the HF and LF sources is $10/1$ and we
initialize the BO with $30$ HF and $30$ LF samples (the composition with the
largest bulk modulus is never in the initial data). To quantify the robustness
of the proposed method to the random initial data, we repeat this process $20$
times for each BO method.

Our second problem is on designing hybrid organic–inorganic perovskite (HOIP)
crystals where we aim to find the compound with the smallest inter-molecular
binding energy [27]. In this example, there are $3$ datasets ($1$ from HF and
$2$ from LF sources) which have the same dimensionality ($1$ output and $3$
categorical inputs with $10$, $3$, and $16$ levels) but different sizes. The
HF dataset has $480$ samples while the first and second LF datasets have $179$
and $240$ samples, respectively. The cost ratio between the three sources is
$15/10/5$ (where the HF and LF2 sources are the most expensive and cheapest,
respectively) and we initialize the BO with $(15,20,15)$ samples for the HF
and LF sources (the best compound is excluded from the initial data). We
repeat the BO process $20$ times to assess the sensitivity of the results to
the initial data.

As mentioned before, the first step in MFBO is to train an LMGP to the initial
data in each problem to exclude the highly biased sources. As NTA has
categorical variables, LMGP learns two manifolds. Based on Figure 4a, the
latent points of the fidelity sources of NTA are very close in the learned
fidelity manifold which indicates that there is a high correlation between the
corresponding two data sources. However, both latent points of LF sources in
HOIP are far from the HF one so they both should be excluded due to their
large global bias. By excluding both LF sources the MF problem in HOIPreduces
to an SF one so we do not exclude the biased LF sources from HOIP to be able
to compare the performance of MFBOUQ with MFBO.

(a) NTA

(b) HOIP

Figure 4: Fidelity manifolds of real-world examples: The plots in (a) and (b)
are obtained by fitting an LMGP to the initial data in the NTA and HOIP
examples, respectively. Due to the consistency across the $20$ repetitions,
the plots are randomly chosen among them. In (b), the HF source is encoded far
from LF sources which indicates that these two sources have large biases with
respect to the HF source and MFBO should exclude them. However, by excluding
these sources, the problem transforms into SF, rendering it incomparable to
MFBOUQ. To maintain comparability with MFBOUQ, we retain the LF sources in
MFBO.

A summary of the convergence history of NTA and HOIP is depicted in Figure 5
by showing the best HF sample ($y^{*}$) found by each method versus its
accumulated sampling cost. In Figure 5a, the LF sources are globally
correlated with the HF source and hence both MF methods perform better than
SFBO by using inexpensive and informative LF data. Additionally, the higher
prediction accuracy of the emulator of MFBOUQ results in a more efficient
sampling and faster convergence of BO in MFBOUQ compared to MFBO. Regarding
the spike in the convergence plot of MFBO in Figure 5a, we note that $18$
repetitions converge at costs below $500$. Consequently, the thick red line
(which is the average across the $20$ repetitions) becomes highly sensitive to
the convergence values after cost exceeds $500$ since it is an average of only
$2$ values. Specifically, in one of these two repetitions the maximum bulk
modulus found is $237$ for many iterations until the cost reaches $544$ when
MFBO suddenly converges to the ground truth (i.e., $255$). This sudden
convergence results in the spike in the corresponding history and, in turn,
the average behavior captured by the thick red line.

(a) NTA

(b) HOIP

Figure 5: Convergence histories of real-world examples: The plots depict the
best HF sample found by each approach ($y^{*}$) versus their sampling costs
accumulated during the BO iterations (the cost of initial data is included).
(a) and (b) summarize the results for the NTA and HOIP, respectively. The thin
curves show the convergence history of each repetition and the solid thick
ones indicate the average behavior across the $20$ repetitions. In (b), MFBO
fails to find the optimum due to it disability in handling biased LF sources.
In both examples, MFBOUQ outperforms other methods.

The superiority of MFBOUQ is more obvious in HOIP (see Figure 5b) which has
two highly biased LF sources. In the HOIP example, MFBO expectedly converges
to a sub-optimal compound since both LF sources are only locally correlated
with the HF source. So, the AFs fail to sample valuable points to improve the
optimization as they cannot find the region where the LF sources are
beneficial and informative. Additionally, each data source is obtained from a
distinct process so it suffers from different types and levels of noise.
Therefore, estimating a single noise for all the data sources in MFBO reduces
the emulation accuracy and further exacerbates the performance of AFs. MFBOUQ
overcomes these issues by focusing more on UQ and estimating separate noise
processes; resulting in a better performance compared to SFBO and especially
to MFBO.

The $2D$ manifolds in Figures 6 and 7 demonstrate the trajectory of BO in the
categorical design space of each data source in NTA and HOIP, respectively.
The top and bottom rows of these figures correspond to MFBO and MFBOUQ. In
these manifolds, each latent point indicates a compound and is color coded
based on the ground truth response value (i.e., the bulk modulus) from each
source. The marker shapes in these manifolds indicate whether a compound is
part of the initial data, sampled during BO, or never seen by LMGP. As
expected, most markers are triangles which indicates that most combinations
are never tested by either MFBO or MFBOUQ. The red arrows next to the legend
mark the range of response in the data sources which indicate that, unlike in
Figure 6 for NTA, the response ranges across the three sources are quite
different in the HOIP problem.

To benefit any MFBO approach, LF sources should be sampled in two primary
regions of their input space: $(1)$ the region that contains their own optima
since each data source is analyzed separately in the auxiliary optimization
problems (see Section 2.4 for details), and $(2)$ the region where the LF
sources are correlated with the HF source. These two regions may overlap with
each other (as is the case in NTA) or not (as is the case in HOIP or the $1D$
example in Figure 1b where MFBOUQ only samples LF2 once when $x<5$). We note
that exploring the correlation region (if it exists!) is crucial for capturing
the relationship between the LF and HF sources and as shown below the
effectiveness of this exploration highly depends on the accuracy of the
emulator in surrogating each source, estimating uncertainties, and identifying
the correlation patterns among different data sources.

(a) HF MFBO

(b) LF MFBO

(c) HF MFBOUQ

(d) LF MFBOUQ

Figure 6: BO sampling history in the encoded categorical design space of NTA:
The plots in the top and bottom row illustrate the exploration-exploitation
behavior of BO in MFBO and MFBOUQ, respectively. The left and right columns
correspond to the space of HF and LF sources, respectively. All latent points
are color-coded based on the ground truth bulk modulus from each source and
the marker shapes indicate whether the compound is part of the initial data,
sampled during BO, or never seen by LMGP. The red arrows next to the legend
indicate the range of response in the two data sources. This figure
effectively demonstrates how strategic sampling in MFBOUQ leads to faster
convergence compared to MFBO(see text for more detailed explanations).

As shown in Figure 6, for both MFBO and MFBOUQ manifolds with very similar
structures are learnt by LMGP for HF and LF data (this was expected per Figure
4 which indicates that the two sources are highly correlated). For instance,
for both LF and HF data, the optimum compound is located at the top-right
corner of the manifold and their values are also quite close ($255$ for HF and
$244$ for LF). This similarity indicates that MFBO and MFBOUQ are both able to
learn about the HF source by sampling the space of the LF source. However,
this sampling is more effective in the case of MFBOUQ since its emulator
quantifies the uncertainties more accurately. In particular, MFBOUQ correctly
samples compounds from the LF source that are mostly encoded in the top right
corner of the manifold (see Figure 6d) while MFBO tests compounds that explore
the entire design space (see Figure 6b).

(a) HF MFBO

(b) LF1 MFBO

(c) LF2 MFBO

(d) HF MFBOUQ

(e) LF1 MFBOUQ

(f) LF2 MFBOUQ

Figure 7: BO sampling history in the encoded categorical design space of HOIP:
The plots in the top and bottom row illustrate the exploration-exploitation
behavior of BO in MFBO and MFBOUQ, respectively. The left, middle and right
columns correspond to the space of HF, LF1 and LF2 sources, respectively. All
latent points are color-coded based on the ground truth binding energy from
each source and the marker shapes indicate whether the compound is part of the
initial data, sampled during BO, or never seen by LMGP. The red arrows next to
the legend indicate the range of responses in the data sources. This figure
demonstrates how the strategic sampling in MFBOUQ enables it to find the
optimum while MFBO fails.

As shown in Figure 7, for any of the sources and with either MFBO or MFBOUQ,
the compounds in the HOIP example are encoded by LMGP into two major clusters
where the smaller one contains the optimum design. By examining these two
clusters we observe that all the compounds in the smaller cluster have
Dimethylformamide (DMF) solvent. These observations are quite interesting in
that they provide engineers with insights into the most important design
variables that affect the materials properties (e.g., DMF solvent which
decreases the binding energy in this example).

The initial HF dataset used in either MFBO or MFBOUQ (see Figures 7(a) and
7(d)) is very small and does not have any compounds from the small cluster
that contains the optimum. However, there are some initial samples from LF1
and LF2 in this cluster and so we should expect BO to leverage these samples
(and the fact that they have some correlation with the unseen HF compounds) in
emulating the HF source and sampling compounds from it that belong to the
small cluster. While this expectation is met by MFBOUQ, MFBO fails to explore
the (encoded) design space that contains the optimum HF sample. This failure
is because $(1)$ both LF sources (especially LF1) provide smaller binding
energies than the HF source, and $(2)$ the emulator of MFBO overestimates the
uncertainties in LF sources. The combination of these two factors prevents
MFBO to find an HF sample that is valuable enough to be selected in Eq. 15.

##  4 Conclusion

In this paper, we develop a novel method to improve the performance of multi-
fidelity cost-aware BO techniques. Our method enhances the accuracy and
convergence rate of MFBO through two main contributions. Firstly, we enable
the emulator to estimate separate noise processes for each source of data.
This feature increases the accuracy of the trained model since different data
sources may exhibit different types and levels of noise. Secondly, we define a
new objective function penalized by strictly proper scoring rules to $(1)$
improve the prediction, $(2)$ increase the focus on UQ, and $(3)$ forgo the
need to exclude highly biased data sources from BO. Our BO method, MFBOUQ,
accommodates any number of data sources with any levels of noise, does not
require any prior knowledge about the relative accuracy of (or relation
between) these sources, and can handle both continuous and categorical
variables. In this paper, we illustrate these features via both analytic and
engineering problems.

In this work, we use two fixed AFs in each iteration. However, one can also
customize the choice of AFs for different iterations using adaptive
approaches. Additionally, the examples presented in this paper are limited to
single-objective problems and we do not aim to exclude the effect of noise in
the final solution (i.e., the best HF sample found is noisy). We intent to
study these direction in our future works.

## Acknowledgments

We appreciate the support from National Science Foundation (award number
$CMMI-2238038$), the Early Career Faculty grant from NASA’s Space Technology
Research Grants Program (award number $80NSSC21K1809$), and the UC National
Laboratory Fees Research Program of the University of California (Grant Number
$L22CR4520$).

## Appendices

##  Appendix A Table of Numerical Examples

Name |  Source ID | Formulation |  n | RRMSE |  Cost
---|---|---|---|---|---
Borehole |  HF | $\frac{2\pi T_{u}(H_{u}-H_{l})}{\ln(\frac{r}{r_{w}})(1+\frac{2LT_{u}}{\ln(\frac{r}{rw})r_{w}^{2}k_{w}}+\frac{T_{u}}{T_{l}})}$ |  $5$ | $-$ |  $1000$
LF1 | $\frac{2\pi T_{\mathrm{u}}(H_{\mathrm{u}}-0.8H_{l})}{\ln(\frac{r}{r_{w}})(1+\frac{1LT_{\mathrm{u}}}{\ln(\frac{r}{r_{w}})r_{w}^{2}k_{w}}+\frac{T_{u}}{T_{l}})}$ |  $5$ | $4.40$ |  $100$
LF2 | $\frac{2\pi T_{u}(H_{u}-H_{l})}{\ln(\frac{r}{r_{w}})(1+\frac{8LT_{u}}{\ln(\frac{r}{rw})r_{w}^{2}k_{w}}+0.75\frac{T_{u}}{T_{l}})}$ |  $50$ | $1.54$ |  $10$
LF3 | $\frac{2\pi T_{\mathrm{u}}(1.09H_{u}-H_{l})}{\ln(\frac{4r}{r_{w}})(1+\frac{3LT_{u}}{\ln(\frac{r}{r_{w}})r_{w}^{2}k_{w}}+\frac{T_{u}}{T_{l}})}$ |  $5$ | $1.30$ |  $100$
LF4 |  $\frac{2\pi T_{\mathrm{u}}(1.05H_{u}-H_{l})}{\ln(\frac{2r}{r_{w}})}$ $\times\left(1+\frac{3LT_{iu}}{\ln(\frac{r}{\tau w})r_{w}^{2}k_{\mathrm{W}}}+\frac{T_{\mathrm{u}}}{T_{l}}\right)$ |  $50$ | $1.3$ |  $10$
Wing |  HF |  $0.036s_{w}^{0.758}w_{fw}^{0.0035}(\frac{A}{\cos^{2}(\Lambda)})^{0.6}$ $\times\lambda^{0.04}(\frac{100t_{c}}{\cos(\Lambda)})^{-0.3}(N_{z}W_{dg})^{0.49}+s_{w}w_{p}$ |  $5$ | $-$ |  $1000$
LF1 |  $0.036s_{w}^{0.758}w_{fw}^{0.0035}(\frac{A}{\cos^{2}(\Lambda)})^{0.6}$ $\times\lambda^{0.04}(\frac{100t_{c}}{\cos(\Lambda)})^{-0.3}(N_{z}W_{dg})^{0.49}+w_{p}$ |  $5$ | $0.19$ |  $100$
LF2 |  $0.036s_{w}^{0.8}w_{fw}^{0.0035}(\frac{A}{\cos^{2}(\Lambda)})^{0.6}$ $\times\lambda^{0.04}(\frac{100t_{c}}{\cos(\Lambda)})^{-0.3}(N_{z}W_{dg})^{0.49}+w_{p}$ |  $10$ | $1.14$ |  $10$
LF3 |  $0.036s_{w}^{0.9}w_{fw}^{0.0035}(\frac{A}{\cos^{2}(\Lambda)})^{0.6}$ $\times\lambda^{0.04}(\frac{100t_{c}}{\cos(\Lambda)})^{-0.3}(N_{z}W_{dg})^{0.49}$ |  $50$ | $5.75$ |  $1$
Table 1: List of analytic functions: $n$ denotes the number of initial
samples. The relative root mean squared error (RRMSE) of an LF source is
calculated by comparing its output to that of the HF source at $10000$ random
points. The cost column is the cost of obtaining a sample from the
corresponding source.

## References

  * [1] Bobak Shahriari et al.  “Taking the human out of the loop: A review of Bayesian optimization”  In _Proceedings of the IEEE_ 104.1 IEEE, 2015, pp. 148–175
  * [2] Peter I Frazier  “A tutorial on Bayesian optimization”  In _arXiv preprint arXiv:1807.02811_ , 2018
  * [3] Shibo Li, Wei Xing, Robert Kirby and Shandian Zhe  “Multi-fidelity Bayesian optimization via deep neural networks”  In _Advances in Neural Information Processing Systems_ 33, 2020, pp. 8521–8531
  * [4] Jialin Song, Yuxin Chen and Yisong Yue  “A general framework for multi-fidelity bayesian optimization with gaussian processes”  In _The 22nd International Conference on Artificial Intelligence and Statistics_ , 2019, pp. 3158–3167  PMLR
  * [5] Shion Takeno et al.  “Multi-fidelity Bayesian optimization with max-value entropy search and its parallelization”  In _International Conference on Machine Learning_ , 2020, pp. 9334–9345  PMLR
  * [6] Shuhan Zhang et al.  “An efficient multi-fidelity bayesian optimization approach for analog circuit synthesis”  In _Proceedings of the 56th Annual Design Automation Conference 2019_ , 2019, pp. 1–6
  * [7] Leshi Shu, Ping Jiang and Yan Wang  “A multi-fidelity Bayesian optimization approach based on the expected further improvement”  In _Structural and Multidisciplinary Optimization_ 63 Springer, 2021, pp. 1709–1719
  * [8] Anh Tran, Tim Wildey and Scott McCann  “sMF-BO-2CoGP: A sequential multi-fidelity constrained Bayesian optimization framework for design applications”  In _Journal of Computing and Information Science in Engineering_ 20.3 American Society of Mechanical Engineers Digital Collection, 2020
  * [9] Manyu Xiao et al.  “Extended Co-Kriging interpolation method based on multi-fidelity data”  In _Applied Mathematics and Computation_ 323 Elsevier, 2018, pp. 120–131
  * [10] Marc C Kennedy and Anthony O’Hagan  “Bayesian calibration of computer models”  In _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ 63.3 Wiley Online Library, 2001, pp. 425–464
  * [11] Jacob Gardner et al.  “Gpytorch: Blackbox matrix-matrix gaussian process inference with gpu acceleration”  In _Advances in neural information processing systems_ 31, 2018
  * [12] Zahra Zanjani Foumani, Mehdi Shishehbor, Amin Yousefpour and Ramin Bostanabad  “Multi-fidelity cost-aware Bayesian optimization”  In _Computer Methods in Applied Mechanics and Engineering_ 407 Elsevier, 2023, pp. 115937
  * [13] P Jorge Escamilla-Ambrosio and Neil Mort  “Hybrid Kalman filter-fuzzy logic adaptive multisensor data fusion architectures”  In _42nd IEEE International Conference on Decision and Control (IEEE Cat. No. 03CH37475)_ 5, 2003, pp. 5215–5220  IEEE
  * [14] Ondrej Kreibich, Jan Neuzil and Radislav Smid  “Quality-based multiple-sensor fusion in an industrial wireless sensor network for MCM”  In _IEEE Transactions on Industrial Electronics_ 61.9 IEEE, 2013, pp. 4903–4911
  * [15] Jonathan Tammer Eweis-Labolle, Nicholas Oune and Ramin Bostanabad  “Data Fusion With Latent Map Gaussian Processes”  In _Journal of Mechanical Design_ 144.9 American Society of Mechanical Engineers, 2022, pp. 091703
  * [16] Nicholas Oune and Ramin Bostanabad  “Latent map Gaussian processes for mixed variable metamodeling”  In _Computer Methods in Applied Mechanics and Engineering_ 387 Elsevier, 2021, pp. 114128
  * [17] Ramin Bostanabad et al.  “Leveraging the nugget parameter for efficient Gaussian process modeling”  In _International journal for numerical methods in engineering_ 114.5 Wiley Online Library, 2018, pp. 501–516
  * [18] Carlos M Carvalho, Nicholas G Polson and James G Scott  “The horseshoe estimator for sparse signals”  In _Biometrika_ 97.2 Oxford University Press, 2010, pp. 465–480
  * [19] Yarin Gal, Mark Van Der Wilk and Carl Edward Rasmussen  “Distributed variational inference in sparse Gaussian process regression and latent variable models”  In _Advances in neural information processing systems_ 27, 2014
  * [20] Rekar O Mohammed and Gavin C Cawley  “Over-fitting in model selection with Gaussian process regression”  In _Machine Learning and Data Mining in Pattern Recognition: 13th International Conference, MLDM 2017, New York, NY, USA, July 15-20, 2017, Proceedings 13_ , 2017, pp. 192–205  Springer
  * [21] Dennis V Lindley  “Scoring rules and the inevitability of probability”  In _International Statistical Review/Revue Internationale de Statistique_ JSTOR, 1982, pp. 1–11
  * [22] Johannes Bracher, Evan L Ray, Tilmann Gneiting and Nicholas G Reich  “Evaluating epidemic forecasts in an interval format”  In _PLoS computational biology_ 17.2 Public Library of Science San Francisco, CA USA, 2021, pp. e1008618
  * [23] K Mitchell and CAT Ferro  “Proper scoring rules for interval probabilistic forecasts”  In _Quarterly Journal of the Royal Meteorological Society_ 143.704 Wiley Online Library, 2017, pp. 1597–1607
  * [24] Tilmann Gneiting and Adrian E Raftery  “Strictly proper scoring rules, prediction, and estimation”  In _Journal of the American statistical Association_ 102.477 Taylor & Francis, 2007, pp. 359–378
  * [25] Carlos Mora et al.  “Probabilistic Neural Data Fusion for Learning from an Arbitrary Number of Multi-fidelity Data Sets”  In _arXiv preprint arXiv:2301.13271_ , 2023
  * [26] Chan-Uk Yeom and Keun-Chang Kwak  “Performance evaluation of automobile fuel consumption using a fuzzy-based granular model with coverage and specificity”  In _Symmetry_ 11.12 MDPI, 2019, pp. 1480
  * [27] Henry C Herbol et al.  “Efficient search of compositional space for hybrid organic–inorganic perovskites via Bayesian optimization”  In _npj Computational Materials_ 4.1 Nature Publishing Group UK London, 2018, pp. 51


#  Ratio of baryon and electric-charge cumulants at second order with
acceptance corrections
J-PARC-TH-0272

Masakiyo Kitazawa  ShinIchi Esumi  Toshihiro Nonaka  Department of Physics,
Osaka University, Toyonaka, Osaka, 560-0043 Japan  J-PARC Branch, KEK Theory
Center, Institute of Particle and Nuclear Studies, KEK, 319-1106 Japan
Tomonaga Center for the History of the Universe, University of Tsukuba,
Tsukuba, Ibaraki, 305-8751, Japan

###### Abstract

We evaluate the ratio of baryon and electric-charge cumulants at second order
from the recent experimental results at $\sqrt{s_{{}_{NN}}}=200$ GeV by the
STAR Collaboration. The baryon number cumulant is reconstructed from the
proton number distribution, and effects of the finite acceptance on the
transverse momentum are corrected assuming the independent particle emission.
We show that the obtained ratio has a dependence on the rapidity window.
Comparison of the result with the hadron resonance gas model and lattice QCD
numerical simulations suggests that if the fluctuations are generated from a
thermal medium its temperature is significantly lower than the chemical
freezeout temperature.

###### keywords:

fluctuations of conserved charges , cumulants , event-by-event analysis

††journal: Physics Letters B

##  1 Introduction

In relativistic heavy-ion collisions (HIC), fluctuations are believed to be
useful observables for investigating phase transitions in the medium created
by the collisions [1, 2, 3]. Fluctuations are characterized by cumulants [1],
which are known to show anomalous behaviors in a thermal medium near the
boundaries of QCD phase transitions especially as the order becomes higher [4,
5, 6, 7, 8, 9, 10]. In the HIC, fluctuations are measurable by the event-by-
event analysis. Various cumulants have been analyzed up to sixth order [11,
12, 13, 14, 15, 16, 17, 18], and a suggestive non-monotonic behavior as
functions of $\sqrt{s_{{}_{NN}}}$ has been reported [14, 15]. The experimental
analyses will be refined further by the new data from the RHIC-BES-II [19] and
HIC in future facilities [20].

Among fluctuation observables, those of conserved charges have particularly
useful properties. First, the cumulants of conserved charges are calculable
unambiguously in thermal field theory, especially in lattice QCD numerical
simulations [21, 22, 23, 24, 25, 26, 27, 28, 29]. While the cumulants in a
thermal system are extensive variables and proportional to the spatial volume,
the volume dependence can be eliminated in their ratios [7, 30]. The
comparison between theoretical and experimental results in terms of the ratios
has been made in the literature [31, 32, 33, 34, 35]. Second, since the
evolution of conserved charges is achieved only by the diffusion it is
typically slow especially when the spatial volume is taken to be large [5, 6,
36]. In earlier studies [5, 6], it has been suggested to use this property of
conserved-charge fluctuations for investigating thermodynamics in the earlier
stage of the HIC. See Refs. [37, 38] for an extension of this idea for higher
order cumulants. Although fluctuations in the HIC are sometimes regarded to be
generated from a thermal medium around chemical freezeout, this assumption has
to be checked carefully because of this property [36].

In the present Letter, among the conserved-charge cumulants we focus on the
ratio of the second-order cumulants of net-baryon number and net-electric
charge, $\langle N_{\rm B}^{2}\rangle_{\rm c}$ and $\langle N_{\rm
Q}^{2}\rangle_{\rm c}$. An advantage to focus on this quantity is that this
would be the ratio between conserved-charge cumulants higher than the first
order that can be analyzed most reliably in the HIC. In the measurement of
cumulants in the HIC, effects of the imperfect performance of detectors, such
as inefficiencies and finite acceptance, have to be corrected [39, 40, 41,
42]. The statistical and systematic uncertainties arising from this procedure,
however, grow rapidly as the order becomes higher. Moreover, when comparing
experimental results with theoretical ones assuming thermodynamics, one has to
consider modifications of fluctuations arising from experimental environment
in the HIC, such as the volume fluctuations [43, 44, 45, 46, 17], global
charge conservation [47, 48, 49], collision pileups [50, 51, 52], dynamical
evolution of fluctuations [5, 6, 37, 38], and etc. The modification of the
cumulants due to these effects is more amplified in a non-trivial way as the
order becomes higher and makes the meaningful comparison more difficult. The
use of the second-order cumulants enables a stable comparison by suppressing
these effects.

Another advantage to focus on this ratio is that, as we will see later, in a
thermal medium the ratio $\langle N_{\rm B}^{2}\rangle_{\rm c}/\langle N_{\rm
Q}^{2}\rangle_{\rm c}$ is a monotonically increasing function of temperature
($T$) and behaves almost linearly as a function of $T$ around the pseudo
critical temperature $T_{\rm c}^{*}\simeq 155$ MeV. This linear $T$ dependence
is suitable for studying the nature of fluctuations in the HIC.

In the present study, we construct the values of $\langle N_{\rm
B}^{2}\rangle_{\rm c}$ and $\langle N_{\rm Q}^{2}\rangle_{\rm c}$ and their
ratio from the recent experimental data in Au+Au central collisions at
$\sqrt{s_{{}_{NN}}}=200$ GeV by the STAR Collaboration [11, 15]. In addition
to the reconstruction of baryon number cumulants [53, 39], we perform the
correction of the finite acceptance in the transverse momentum, $p_{T}$, space
assuming the independent particle emission. We show that the acceptance
correction has a large effect on the individual cumulants and the ratio,
suggesting the importance of the correction and the necessity to measure
fluctuations with wider acceptance.

We compare the value of $\langle N_{\rm B}^{2}\rangle_{\rm c}/\langle N_{\rm
Q}^{2}\rangle_{\rm c}$ obtained in this way with the ratio obtained in the
hadron resonance gas (HRG) model and the lattice QCD simulations. Provided
that the experimentally-observed fluctuations are emitted from a thermal
medium, the comparison shows the temperature $T\simeq 134-138$ MeV, which is
significantly lower than the chemical freezeout temperature $T_{\rm chem}$.
Our result also shows that the ratio $\langle N_{\rm B}^{2}\rangle_{\rm
c}/\langle N_{\rm Q}^{2}\rangle_{\rm c}$ has a dependence on the rapidity
window $\Delta y$ whereas it is independent of $\Delta y$ for thermal
fluctuations. These results suggest the violation of the assumption and
motivate further investigations on the nature of fluctuations in the HIC.

##  2 Experimental data and their correction

To obtain the cumulant ratio $\langle N_{\rm B}^{2}\rangle_{\rm c}/\langle
N_{\rm Q}^{2}\rangle_{\rm c}$ in the HIC, we use experimental results by the
STAR Collaboration in Refs. [11, 15]. We construct $\langle N_{\rm
B}^{2}\rangle_{\rm c}$ from the data on the proton number cumulants in Ref.
[15] according to the procedure in Refs. [53, 39], while we use the data in
Ref. [11] for $\langle N_{\rm Q}^{2}\rangle_{\rm c}$. Effects of the
detector’s efficiencies are corrected in these results. Throughout this study
we concentrate on the result for the most central ($0-5\%$) Au+Au collisions
at $\sqrt{s_{{}_{NN}}}=200$ GeV. Effects of the violation of the boost
invariance and the global charge conservation are most suppressed and the
method in Refs. [53, 39] is well justified at this $\sqrt{s_{{}_{NN}}}$. Since
the chemical potentials of conserved charges are small at this
$\sqrt{s_{{}_{NN}}}$, we neglect their effects [25, 26, 28] in the following.
Since the effects of the QCD critical point [4] would be suppressed at this
$\sqrt{s_{{}_{NN}}}$, the analysis is suitable for studying non-critical
behavior of fluctuations.

In Ref. [15] and Ref. [11], particles are observed in rapidity and pseudo-
rapidity spaces, respectively. To compare experimental results with
theoretical ones, the use of space-time rapidity is most desirable [1].
Rapidity and pseudo-rapidity are used for its proxy, while the former is
better in the Bjorken picture [54]. This is the reason why we use the data on
the proton number cumulant in Ref. [15]; for protons, rapidity and pseudo-
rapidity has a significant difference due to a heavy proton mass. On the other
hand, the difference is smaller in $\langle N_{\rm Q}^{2}\rangle_{\rm c}$
since electric charges are dominantly carried by pions whose mass is
comparable with the mean $p_{T}$. In the following, we thus regard the
analysis in Ref. [11] as that in the rapidity space.

Table 1:  Particle abundance in the $p_{T}$ acceptance $R_{p_{T}}$ obtained by the Blast-wave model in Ref. [55] for the most-central collisions at $\sqrt{s_{{}_{NN}}}=200$ GeV.  particle species ($p_{T}$ range) | $R_{p_{T}}$
---|---
pions ($0.4<p_{T}<1.6$ GeV) | 0.44
kaons ($0.4<p_{T}<1.6$ GeV) | 0.71
protons ($0.4<p_{T}<1.6$ GeV) | 0.71
$\pi$+$K$+$p$ ($0.4<p_{T}<1.6$ GeV) | 0.49
protons ($0.4<p_{T}<2.0$ GeV) | 0.82

The measurements in Refs. [11, 15] are performed within a finite
$p_{T}$-acceptance; $0.4<p_{T}<2.0$ GeV in Ref. [15] and $0.4<p_{T}<1.6$ GeV
in Ref. [11]. Due to the acceptance the particles in the final state are
observed only with imperfect probabilities

$\displaystyle
R_{p_{T}}=\frac{{\rm(particle~{}number~{}in}~{}p_{T}\mathrm{-acceptance)}}{\rm(total~{}particle~{}number)}.$
(1)

Using the Blast wave model with the parameters for $\sqrt{s_{{}_{NN}}}=200$
GeV in Ref. [55], the values of $R_{p_{T}}$ for individual particles are
obtained as shown in Table 1 for the $p_{T}$-acceptances in Refs. [11, 15].
“$\pi$+$K$+$p$” in the Table shows the weighted probability for the charged
particles with the particle abundances taken from Ref. [55].

Since the measurement in the finite $p_{T}$-acceptance modifies the particle-
number distributions, its effect on the cumulants has to be corrected before
comparing them with theoretical studies without such an acceptance cut. In the
present study, we perform this correction assuming that the individual
particles are emitted toward different $p_{T}$ with independent probabilities
according to a given $p_{T}$ distribution. In this case, the correction can be
carried out with the same procedure as the efficiency correction assuming the
binomial distribution [39, 1] with the probabilities $R_{p_{T}}$ in Table 1111
For $\langle N_{\rm Q}^{2}\rangle_{\rm c}$, we carry out this procedure with
the weighted probability, i.e. “$\pi$+$K$+$p$” in the Table. One may think
that the correction procedure for multi-particle species [41] with individual
probabilities for $\pi$, $K$, and $p$ should be used for this correction.
However, in Ref. [11] the electric charge $N_{\rm Q}$ is measured without
particle identification. One thus cannot employ the method in Ref. [41]
without a detailed knowledge on the detector’s response. In any case, as
discussed in Ref. [41] the systematic deviation due to the use of the weighted
$R_{p_{T}}$ is small at the second order and thus our results would be less
affected, while the effect is amplified for higher order cumulants. .

We note that another way to compensate the effect of the $p_{T}$-acceptance in
the comparison between experimental and theoretical analyses is to perform the
correction in theoretical calculations, if possible. Such comparisons have
been made in the studies employing the HRG model [56, 31, 33, 27, 57].
However, in general it is not possible to perform such a correction in
theoretical calculations. The correction of the experimental data enables
direct comparisons even for such cases in terms of the true values of
cumulants.

Figure 1:  Second-order cumulants $\langle N_{\rm B}^{2}\rangle_{\rm c}$,
$\langle N_{p}^{2}\rangle_{\rm c}$, and $\langle N_{\rm Q}^{2}\rangle_{\rm c}$
divided by $\Delta y$ obtained from the experimental results in Refs. [11, 15]
at $\sqrt{s_{{}_{NN}}}=200$ GeV. The circle and square symbols show the
results with and without the $p_{T}$-acceptance correction, respectively. The
triangles in the right panel show $\langle N_{p}^{2}\rangle_{\rm c}/\Delta y$.
The error bars show statistical errors. The dashed lines near the symbols are
the corresponding total particle numbers.

##  3 Cumulant ratio

In Fig. 1, we show the second-order cumulants $\langle N_{\rm
B}^{2}\rangle_{\rm c}$, $\langle N_{\rm Q}^{2}\rangle_{\rm c}$, as well as
that of the proton number $\langle N_{p}^{2}\rangle_{\rm c}$, divided by the
rapidity window $\Delta y$ as functions of $\Delta y$. These quantities are
constant if they are generated from a thermal system having a boost invariance
[1]. In the left panel, the triangles show $\langle N_{p}^{2}\rangle_{\rm
c}/\Delta y$ in Ref. [15]. The dashed line near the data shows the total
particle number $\langle N_{p}^{\rm(total)}\rangle/\Delta y$. The squares in
the same panel show $\langle N_{\rm B}^{2}\rangle_{\rm c}/\Delta y$ obtained
with the procedure in Refs. [53, 39], while the circles show $\langle N_{\rm
B}^{2}\rangle_{\rm c}/\Delta y$ for which the $p_{T}$-acceptance correction is
performed. The error bars show the statistical errors, which are negligibly
small in these results. The dashed lines near these results are the total
baryon number $\langle N_{\rm B}^{\rm(total)}\rangle/\Delta y=2\langle
N_{p}^{\rm(total)}\rangle/\Delta y$. One sees that the deviation of $\langle
N_{\rm B}^{2}\rangle_{\rm c}$ from $\langle N_{\rm B}^{\rm(total)}\rangle$ at
large $\Delta y$ is pronounced by the corrections. This result is reasonable
since the incomplete measurement tends to make the distribution close to the
Skellam distribution in which $\langle N_{\rm B}^{2}\rangle_{\rm c}=\langle
N_{\rm B}^{\rm(total)}\rangle$ [1].

Shown in the right panel of Fig. 1 are $\langle N_{\rm Q}^{2}\rangle_{\rm
c}/\Delta y$ with (circles) and without (squares) the $p_{T}$-acceptance
correction222 The total electric charge required for the $p_{T}$-acceptance
correction is provided by the private communication with Arghya Chatterjee. .
The meaning of the dashed lines is the same as the left panel. The panel shows
that the effect of the $p_{T}$-acceptance correction is more significant than
$\langle N_{\rm B}^{2}\rangle_{\rm c}$ because of the smaller $R_{p_{T}}$ for
the electric charge.

Figure 2:  Left panel shows the ratio $\langle\delta N_{\rm
B}^{2}\rangle/\langle\delta N_{\rm Q}^{2}\rangle$ obtained from the HIC as a
function of the rapidity window $\Delta y$. The right panel shows the ratio
obtained in the lattice QCD simulations [29] and the HRG model. The dotted
horizontal lines and shaded area show a comparison of the two results.

In the left panel of Fig. 2, we show the $p_{T}$-acceptance-corrected result
of $\langle N_{\rm B}^{2}\rangle_{\rm c}/\langle N_{\rm Q}^{2}\rangle_{\rm c}$
by the circles. The dashed lines show $\langle N_{\rm
B}^{\rm(total)}\rangle/\langle N_{\rm Q}^{\rm(total)}\rangle$. The shaded band
represents the systematic errors that account for the propagation from that of
$\langle N_{\rm B}^{2}\rangle_{\rm c}$ in Refs. [15]; we, however, note that
this error band should be regarded only as a guide since the estimate of the
systematic errors of $\langle N_{\rm B}^{2}\rangle_{\rm c}/\langle N_{\rm
Q}^{2}\rangle_{\rm c}$ needs a detailed knowledge on the experimental
analyses. In the panel, the result without the $p_{T}$-acceptance correction
is also shown by the squares as a reference. One sees that the correction
strongly modifies the ratio.

If fluctuations are emitted from a thermal system having a boost invariance,
the ratio $\langle N_{\rm B}^{2}\rangle_{\rm c}/\langle N_{\rm
Q}^{2}\rangle_{\rm c}$ is a constant as a function of $\Delta y$ [1]. While
$\langle N_{\rm B}^{2}\rangle_{\rm c}/\Delta y$ and $\langle N_{\rm
Q}^{2}\rangle_{\rm c}/\Delta y$ become decreasing functions when the effects
of the global charge conservation are taken into account [47, 48], these
$\Delta y$ dependence cancels out in the ratio $\langle N_{\rm
B}^{2}\rangle_{\rm c}/\langle N_{\rm Q}^{2}\rangle_{\rm c}$ for the thermal
case [48]. On the other hand, from the left panel of Fig. 2 one sees that the
acceptance-corrected ratio has a clear increasing trend as a function of
$\Delta y$. This result shows that the fluctuations in the HIC are not emitted
from a purely thermal system including the effect of global conservation.

##  4 Comparison with HRG model and lattice results

Assuming that the fluctuations observed in the HIC are those of a thermal
system, one can estimate the temperature of the system by comparing the ratio
of cumulants with the results obtained in lattice QCD simulations. Even when
the fluctuations are not thermal, such a comparison is useful for
investigating the nature of the fluctuations.

To perform the comparison using the ratio $\langle N_{\rm B}^{2}\rangle_{\rm
c}/\langle N_{\rm Q}^{2}\rangle_{\rm c}$, in the right panel of Fig. 2 we show
the $T$ dependence of $\langle N_{\rm B}^{2}\rangle_{\rm c}/\langle N_{\rm
Q}^{2}\rangle_{\rm c}$ obtained from a lattice QCD simulation [29] by the
solid line with an error band. Finite-volume effects of $\langle N_{\rm
Q}^{2}\rangle_{\rm c}$ are corrected according to Ref. [29]. The range of the
vertical axis is the same as the left panel. The lattice results on
thermodynamics are known to be well reproduced by the HRG model at low $T$. In
the panel, the ratio $\langle N_{\rm B}^{2}\rangle_{\rm c}/\langle N_{\rm
Q}^{2}\rangle_{\rm c}$ obtained in the HRG model is shown by the dashed line,
where we use the set of hadrons in “QMHRG2020” [29] for the HRG model. The
figure shows that the lattice result agrees well with the HRG model for
$T\lesssim 145$ MeV, which suggests the validity of the latter in this range
of $T$. From the panel one also finds that $\langle N_{\rm B}^{2}\rangle_{\rm
c}/\langle N_{\rm Q}^{2}\rangle_{\rm c}$ behaves almost linearly as a function
of $T$ in the range of $T$ shown in the panel. As discussed in Sec. 1, this is
an attractive feature of this ratio.

To compare the results in the left and right panels, in Fig. 2 we show the
dotted horizontal lines at the values of $\langle N_{\rm B}^{2}\rangle_{\rm
c}/\langle N_{\rm Q}^{2}\rangle_{\rm c}$ obtained from the HIC at $\Delta
y=0.2$ and $1.0$. By comparing these values with the ratio in the HRG model,
one finds that the temperature extracted from the naïve comparison gives
$T\simeq 134-138$ MeV depending on $\Delta y$ as shown by the shaded box in
the right panel. We note that this temperature is significantly smaller than
the chemical freezeout temperature $T_{\rm chem}\simeq 156$ MeV for
$\sqrt{s_{{}_{NN}}}=200$ GeV [55]. The value of $\langle N_{\rm
B}^{2}\rangle_{\rm c}/\langle N_{\rm Q}^{2}\rangle_{\rm c}$ in the HIC itself
is about twice smaller than the value in the HRG model at $T=T_{\rm chem}$.

##  5 Discussions

In the present study, we have investigated the cumulant ratio $\langle N_{\rm
B}^{2}\rangle_{\rm c}/\langle N_{\rm Q}^{2}\rangle_{\rm c}$ observed in the
HIC. Because this ratio consists only of the second-order conserved-charge
cumulants, various uncertainties in the experimental measurement that are
amplified for higher order cumulants are suppressed in its analysis. The ratio
in the HIC at $\sqrt{s_{{}_{NN}}}=200$ GeV is estimated from the experimental
results by the STAR Collaboration [11, 15]. In addition to the reconstruction
of the baryon number cumulant from those of protons, effects of the
$p_{T}$-acceptance are corrected assuming the independent particle emission.
Our result shows that this correction strongly modifies the resulting values
of the cumulants and their ratio and thus is crucial. Since the effect of the
correction becomes more significant for higher order cumulants [41], this
result also shows the importance of the correction in their analysis for the
search for the QCD critical point.

The naïve comparison of the obtained ratio with the HRG model suggests the
temperature $T\simeq 134-138$ MeV, which is significantly lower than $T_{\rm
chem}$. By taking this result seriously, it is suggested that the fluctuation
observables in the HIC are generated in the hadronic phase later than the
chemical freeze out in contrast to the earlier suggestions [5, 6].

However, we emphasize that this comparison is made assuming that the
fluctuations in the HIC are thermal. On the other hand, existence of the
$\Delta y$ dependence of $\langle N_{\rm B}^{2}\rangle_{\rm c}/\langle N_{\rm
Q}^{2}\rangle_{\rm c}$ shows the violation of this assumption in the HIC.
Therefore, to understand the experimental result correctly one needs further
investigations on the nature of fluctuations especially taking their dynamics
into account [5, 6, 37, 38]. The modifications of the cumulants due to the use
of (pseudo-)rapidity in place of space-time rapidity [54] and the resonance
decays after the chemical freezeout are other important effects to be
considered since they tend to make $\langle N_{\rm Q}^{2}\rangle_{\rm c}$
larger and suppress $\langle N_{\rm B}^{2}\rangle_{\rm c}/\langle N_{\rm
Q}^{2}\rangle_{\rm c}$. Because these effects are suppressed by extending
$\Delta y$ [1], the measurement of the fluctuations with larger $\Delta y$ is
an important experimental subject for resolving these issues.

Finally, we remark that the value of $\langle N_{\rm B}^{2}\rangle_{\rm
c}/\langle N_{\rm Q}^{2}\rangle_{\rm c}$ obtained in the present study from
Refs. [11, 15] would have a deviation from the true value. As discussed in
Sec. 2, while $\langle N_{\rm Q}^{2}\rangle_{\rm c}$ is measured in the
pseudo-rapidity space in Ref. [11], the measurement in the rapidity space is
more desirable. Although we have performed the $p_{T}$-acceptance correction
assuming the independent particle emission in the $p_{T}$ space, there is no a
priori justification of this assumption in the HIC. In addition to theoretical
studies on this point, the measurement of fluctuations with wider acceptance
is an important experimental subject for reducing uncertainties from the
violation of this assumption. As emphasized already, such measurements are
more important for the analyses of higher order cumulants. It thus is quite
interesting to realize such experiments at the future experiments at FAIR,
NICA and J-PARC-HI [20].

## Acknowledgments

The authors thank Arghya Chatterjee for providing unpublished data in Ref.
[11]. They also thank Krzysztof Redlich and Volodymyr Vovchenko for valuable
comments. This work was supported by JSPS KAKENHI (Grants No. JP19H05598,
JP20H01903, JP22K03619).

## References

  * [1] M. Asakawa, M. Kitazawa, Fluctuations of conserved charges in relativistic heavy ion collisions: An introduction, Prog. Part. Nucl. Phys. 90 (2016) 299–342.  arXiv:1512.05038, doi:10.1016/j.ppnp.2016.04.002.
  * [2] A. Bzdak, S. Esumi, V. Koch, J. Liao, M. Stephanov, N. Xu, Mapping the Phases of Quantum Chromodynamics with Beam Energy Scan, Phys. Rept. 853 (2020) 1–87.  arXiv:1906.00936, doi:10.1016/j.physrep.2020.01.005.
  * [3] M. Bluhm, et al., Dynamics of critical fluctuations: Theory – phenomenology – heavy-ion collisions, Nucl. Phys. A 1003 (2020) 122016\.  arXiv:2001.08831, doi:10.1016/j.nuclphysa.2020.122016.
  * [4] M. A. Stephanov, K. Rajagopal, E. V. Shuryak, Event-by-event fluctuations in heavy ion collisions and the QCD critical point, Phys. Rev. D 60 (1999) 114028\.  arXiv:hep-ph/9903292, doi:10.1103/PhysRevD.60.114028.
  * [5] M. Asakawa, U. W. Heinz, B. Muller, Fluctuation probes of quark deconfinement, Phys. Rev. Lett. 85 (2000) 2072–2075.  arXiv:hep-ph/0003169, doi:10.1103/PhysRevLett.85.2072.
  * [6] S. Jeon, V. Koch, Charged particle ratio fluctuation as a signal for QGP, Phys. Rev. Lett. 85 (2000) 2076–2079.  arXiv:hep-ph/0003168, doi:10.1103/PhysRevLett.85.2076.
  * [7] S. Ejiri, F. Karsch, K. Redlich, Hadronic fluctuations at the QCD phase transition, Phys. Lett. B 633 (2006) 275–282.  arXiv:hep-ph/0509051, doi:10.1016/j.physletb.2005.11.083.
  * [8] M. A. Stephanov, Non-Gaussian fluctuations near the QCD critical point, Phys. Rev. Lett. 102 (2009) 032301.  arXiv:0809.3450, doi:10.1103/PhysRevLett.102.032301.
  * [9] M. Asakawa, S. Ejiri, M. Kitazawa, Third moments of conserved charges as probes of QCD phase structure, Phys. Rev. Lett. 103 (2009) 262301.  arXiv:0904.2089, doi:10.1103/PhysRevLett.103.262301.
  * [10] B. Friman, F. Karsch, K. Redlich, V. Skokov, Fluctuations as probe of the QCD phase transition and freeze-out in heavy ion collisions at LHC and RHIC, Eur. Phys. J. C 71 (2011) 1694.  arXiv:1103.3511, doi:10.1140/epjc/s10052-011-1694-2.
  * [11] J. Adam, et al., Collision-energy dependence of second-order off-diagonal and diagonal cumulants of net-charge, net-proton, and net-kaon multiplicity distributions in Au + Au collisions, Phys. Rev. C 100 (1) (2019) 014902, [Erratum: Phys.Rev.C 105, 029901 (2022)].  arXiv:1903.05370, doi:10.1103/PhysRevC.100.014902.
  * [12] L. Adamczyk, et al., Collision Energy Dependence of Moments of Net-Kaon Multiplicity Distributions at RHIC (2017).  arXiv:1709.00773.
  * [13] L. Adamczyk, et al., Beam energy dependence of moments of the net-charge multiplicity distributions in Au+Au collisions at RHIC, Phys. Rev. Lett. 113 (2014) 092301.  arXiv:1402.1558, doi:10.1103/PhysRevLett.113.092301.
  * [14] J. Adam, et al., Net-proton number fluctuations and the Quantum Chromodynamics critical point (2020).  arXiv:2001.02852.
  * [15] M. Abdallah, et al., Cumulants and correlation functions of net-proton, proton, and antiproton multiplicity distributions in Au+Au collisions at energies available at the BNL Relativistic Heavy Ion Collider, Phys. Rev. C 104 (2) (2021) 024902.  arXiv:2101.12413, doi:10.1103/PhysRevC.104.024902.
  * [16] M. S. Abdallah, et al., Measurements of Proton High Order Cumulants in 3 GeV Au+Au Collisions and Implications for the QCD Critical Point (11 2021).  arXiv:2112.00240.
  * [17] J. Adamczewski-Musch, et al., Proton number fluctuations in $\sqrt{s_{NN}}$ = 2.4 GeV Au+Au collisions studied with HADES (2020).  arXiv:2002.08701.
  * [18] S. Acharya, et al., Global baryon number conservation encoded in net-proton fluctuations measured in Pb-Pb collisions at $\sqrt{s_{\rm NN}}$ = 2.76 TeV, Phys. Lett. B 807 (2020) 135564.  arXiv:1910.14396, doi:10.1016/j.physletb.2020.135564.
  * [19] A STAR white paper summarizing the current understanding and describing future plans (2014).
URL https://drupal.star.bnl.gov/STAR/system/files/BES_WPII_ver6.9_Cover.pdf

  * [20] T. Galatyuk, Future facilities for high $\mu_{B}$ physics, Nucl. Phys. A 982 (2019) 163–169.  doi:10.1016/j.nuclphysa.2018.11.025.
  * [21] S. Borsanyi, Z. Fodor, S. D. Katz, S. Krieg, C. Ratti, K. K. Szabo, Freeze-out parameters from electric charge and baryon number fluctuations: is there consistency?, Phys. Rev. Lett. 113 (2014) 052301.  arXiv:1403.4576, doi:10.1103/PhysRevLett.113.052301.
  * [22] H. T. Ding, S. Mukherjee, H. Ohno, P. Petreczky, H. P. Schadler, Diagonal and off-diagonal quark number susceptibilities at high temperatures, Phys. Rev. D 92 (7) (2015) 074043.  arXiv:1507.06637, doi:10.1103/PhysRevD.92.074043.
  * [23] M. D’Elia, G. Gagliardi, F. Sanfilippo, Higher order quark number fluctuations via imaginary chemical potentials in $N_{f}=2+1$ QCD, Phys. Rev. D 95 (9) (2017) 094503.  arXiv:1611.08285, doi:10.1103/PhysRevD.95.094503.
  * [24] A. Bazavov, et al., The QCD Equation of State to $\mathcal{O}(\mu_{B}^{6})$ from Lattice QCD, Phys. Rev. D 95 (5) (2017) 054504.  arXiv:1701.04325, doi:10.1103/PhysRevD.95.054504.
  * [25] A. Bazavov, et al., Skewness and kurtosis of net baryon-number distributions at small values of the baryon chemical potential, Phys. Rev. D 96 (7) (2017) 074510\.  arXiv:1708.04897, doi:10.1103/PhysRevD.96.074510.
  * [26] S. Borsanyi, Z. Fodor, J. N. Guenther, S. K. Katz, K. K. Szabo, A. Pasztor, I. Portillo, C. Ratti, Higher order fluctuations and correlations of conserved charges from lattice QCD, JHEP 10 (2018) 205.  arXiv:1805.04445, doi:10.1007/JHEP10(2018)205.
  * [27] R. Bellwied, S. Borsanyi, Z. Fodor, J. N. Guenther, J. Noronha-Hostler, P. Parotto, A. Pasztor, C. Ratti, J. M. Stafford, Off-diagonal correlators of conserved charges from lattice QCD and how to relate them to experiment, Phys. Rev. D 101 (3) (2020) 034506.  arXiv:1910.14592, doi:10.1103/PhysRevD.101.034506.
  * [28] A. Bazavov, et al., Skewness, kurtosis, and the fifth and sixth order cumulants of net baryon-number distributions from lattice QCD confront high-statistics STAR data, Phys. Rev. D 101 (7) (2020) 074502.  arXiv:2001.08530, doi:10.1103/PhysRevD.101.074502.
  * [29] D. Bollweg, J. Goswami, O. Kaczmarek, F. Karsch, S. Mukherjee, P. Petreczky, C. Schmidt, P. Scior, Second order cumulants of conserved charge fluctuations revisited: Vanishing chemical potentials, Phys. Rev. D 104 (7) (2021).  arXiv:2107.10011, doi:10.1103/PhysRevD.104.074512.
  * [30] F. Karsch, K. Redlich, Probing freeze-out conditions in heavy ion collisions with moments of charge fluctuations, Phys. Lett. B 695 (2011) 136–142.  arXiv:1007.2581, doi:10.1016/j.physletb.2010.10.046.
  * [31] P. Alba, R. Bellwied, M. Bluhm, V. Mantovani Sarti, M. Nahrgang, C. Ratti, Sensitivity of multiplicity fluctuations to freeze-out conditions in heavy ion collisions, Phys. Rev. C 92 (6) (2015) 064910.  arXiv:1504.03262, doi:10.1103/PhysRevC.92.064910.
  * [32] A. Chatterjee, S. Chatterjee, T. K. Nayak, N. R. Sahoo, Diagonal and off-diagonal susceptibilities of conserved quantities in relativistic heavy-ion collisions, J. Phys. G 43 (12) (2016) 125103.  arXiv:1606.09573, doi:10.1088/0954-3899/43/12/125103.
  * [33] M. Bluhm, M. Nahrgang, Freeze-out conditions from strangeness observables at RHIC, Eur. Phys. J. C 79 (2) (2019) 155.  arXiv:1806.04499, doi:10.1140/epjc/s10052-019-6661-3.
  * [34] P. Braun-Munzinger, B. Friman, K. Redlich, A. Rustamov, J. Stachel, Relativistic nuclear collisions: Establishing a non-critical baseline for fluctuation measurements, Nucl. Phys. A 1008 (2021) 122141.  arXiv:2007.02463, doi:10.1016/j.nuclphysa.2021.122141.
  * [35] S. Gupta, D. Mallick, D. K. Mishra, B. Mohanty, N. Xu, Limits of thermalization in relativistic heavy ion collisions, Phys. Lett. B 829 (2022) 137021.  doi:10.1016/j.physletb.2022.137021.
  * [36] M. Asakawa, M. Kitazawa, B. Müller, Issues with the search for critical point in QCD with relativistic heavy ion collisions, Phys. Rev. C 101 (3) (2020) 034913.  arXiv:1912.05840, doi:10.1103/PhysRevC.101.034913.
  * [37] M. Kitazawa, M. Asakawa, H. Ono, Non-equilibrium time evolution of higher order cumulants of conserved charges and event-by-event analysis, Phys. Lett. B 728 (2014) 386–392.  arXiv:1307.2978, doi:10.1016/j.physletb.2013.12.008.
  * [38] M. Kitazawa, Rapidity window dependences of higher order cumulants and diffusion master equation, Nucl. Phys. A 942 (2015) 65–96.  arXiv:1505.04349, doi:10.1016/j.nuclphysa.2015.07.008.
  * [39] M. Kitazawa, M. Asakawa, Relation between baryon number fluctuations and experimentally observed proton number fluctuations in relativistic heavy ion collisions, Phys. Rev. C 86 (2012) 024904, [Erratum: Phys.Rev.C 86, 069902 (2012)].  arXiv:1205.3292, doi:10.1103/PhysRevC.86.024904.
  * [40] A. Bzdak, V. Koch, Local Efficiency Corrections to Higher Order Cumulants, Phys. Rev. C 91 (2) (2015) 027901.  arXiv:1312.4574, doi:10.1103/PhysRevC.91.027901.
  * [41] T. Nonaka, M. Kitazawa, S. Esumi, More efficient formulas for efficiency correction of cumulants and effect of using averaged efficiency, Phys. Rev. C95 (6) (2017) 064912.  arXiv:1702.07106, doi:10.1103/PhysRevC.95.064912.
  * [42] T. Nonaka, M. Kitazawa, S. Esumi, A general procedure for detector–response correction of higher order cumulants, Nucl. Instrum. Meth. A 906 (2018) 10–17.  arXiv:1805.00279, doi:10.1016/j.nima.2018.08.013.
  * [43] M. I. Gorenstein, M. Gazdzicki, Strongly Intensive Quantities, Phys. Rev. C 84 (2011) 014904.  arXiv:1101.4865, doi:10.1103/PhysRevC.84.014904.
  * [44] V. Skokov, B. Friman, K. Redlich, Volume Fluctuations and Higher Order Cumulants of the Net Baryon Number, Phys. Rev. C 88 (2013) 034911.  arXiv:1205.4756, doi:10.1103/PhysRevC.88.034911.
  * [45] P. Braun-Munzinger, A. Rustamov, J. Stachel, Bridging the gap between event-by-event fluctuation measurements and theory predictions in relativistic nuclear collisions, Nucl. Phys. A960 (2017) 114–130.  arXiv:1612.00702, doi:10.1016/j.nuclphysa.2017.01.011.
  * [46] T. Sugiura, T. Nonaka, S. Esumi, Volume fluctuation and multiplicity correlation in higher-order cumulants, Phys. Rev. C100 (4) (2019) 044904.  arXiv:1903.02314, doi:10.1103/PhysRevC.100.044904.
  * [47] M. Sakaida, M. Asakawa, M. Kitazawa, Effects of global charge conservation on time evolution of cumulants of conserved charges in relativistic heavy ion collisions, Phys. Rev. C 90 (6) (2014) 064911.  arXiv:1409.6866, doi:10.1103/PhysRevC.90.064911.
  * [48] V. Vovchenko, R. V. Poberezhnyuk, V. Koch, Cumulants of multiple conserved charges and global conservation laws, JHEP 10 (2020) 089.  arXiv:2007.03850, doi:10.1007/JHEP10(2020)089.
  * [49] P. Braun-Munzinger, A. Rustamov, J. Stachel, The role of the local conservation laws in fluctuations of conserved charges (7 2019).  arXiv:1907.03032.
  * [50] P. Garg, D. Mishra, Higher moments of net-proton multiplicity distributions in a heavy-ion event pile-up scenario, Phys. Rev. C 96 (4) (2017) 044908.  arXiv:1705.01256, doi:10.1103/PhysRevC.96.044908.
  * [51] S. Sombun, J. Steinheimer, C. Herold, A. Limphirat, Y. Yan, M. Bleicher, Higher order net-proton number cumulants dependence on the centrality definition and other spurious effects, J. Phys. G45 (2) (2018) 025101.  arXiv:1709.00879, doi:10.1088/1361-6471/aa9c6c.
  * [52] T. Nonaka, et al., Pileup corrections on higher-order cumulants, Nucl. Instrum. Meth. A 984 (2020) 164632.  arXiv:2006.15809, doi:10.1016/j.nima.2020.164632.
  * [53] M. Kitazawa, M. Asakawa, Revealing baryon number fluctuations from proton number fluctuations in relativistic heavy ion collisions, Phys. Rev. C 85 (2012) 021901.  arXiv:1107.2755, doi:10.1103/PhysRevC.85.021901.
  * [54] Y. Ohnishi, M. Kitazawa, M. Asakawa, Thermal blurring of event-by-event fluctuations generated by rapidity conversion, Phys. Rev. C 94 (4) (2016) 044905\.  arXiv:1606.03827, doi:10.1103/PhysRevC.94.044905.
  * [55] B. I. Abelev, et al., Systematic Measurements of Identified Particle Spectra in $pp,d^{+}$ Au and Au+Au Collisions from STAR, Phys. Rev. C 79 (2009) 034909\.  arXiv:0808.2041, doi:10.1103/PhysRevC.79.034909.
  * [56] P. Alba, W. Alberico, R. Bellwied, M. Bluhm, V. Mantovani Sarti, M. Nahrgang, C. Ratti, Freeze-out conditions from net-proton and net-charge fluctuations at RHIC, Phys. Lett. B 738 (2014) 305–310.  arXiv:1403.4903, doi:10.1016/j.physletb.2014.09.052.
  * [57] V. Vovchenko, V. Koch, Thermodynamic approach to proton number fluctuations in baryon-rich heavy-ion matter created at moderate collision energies (3 2022).  arXiv:2204.00137.