# Model fusion with physics-guided machine learning

Suraj Pawar
School of Mechanical & Aerospace Engineering,
Oklahoma State University,
Stillwater, Oklahoma - 74078, USA.
<EMAIL_ADDRESS>
&Omer San
School of Mechanical & Aerospace Engineering,
Oklahoma State University,
Stillwater, Oklahoma - 74078, USA.
<EMAIL_ADDRESS>
&Aditya Nair
Department of Mechanical Engineering,
University of Nevada,
Reno, NV 89557, USA &Adil Rasheed
Department of Engineering Cybernetics,
Norwegian University of Science and Technology,
7465 Trondheim, Norway.
Department of Mathematics and Cybernetics,
SINTEF Digital,
7034 Trondheim, Norway.
Trond Kvamsdal
Department of Mathematics and Cybernetics,
SINTEF Digital,
7034 Trondheim, Norway.
Department of Mathematical Sciences,
Norwegian University of Science and Technology,
7491 Trondheim, Norway

###### Abstract

The unprecedented amount of data generated from experiments, field
observations, and large-scale numerical simulations at a wide range of spatio-
temporal scales have enabled the rapid advancement of data-driven and
especially deep learning models in the field of fluid mechanics. Although
these methods are proven successful for many applications, there is a grand
challenge of improving their _generalizability_. This is particularly
essential when data-driven models are employed within outer-loop applications
like optimization. In this work, we put forth a physics-guided machine
learning (PGML) framework that leverages the interpretable physics-based model
with a deep learning model. The PGML framework is capable of enhancing the
generalizability of data-driven models and effectively protect against or
inform about the inaccurate predictions resulting from extrapolation. We apply
the PGML framework as a novel model fusion approach combining the physics-
based Galerkin projection model and long-short term memory (LSTM) network for
parametric model order reduction of fluid flows. We demonstrate the improved
generalizability of the PGML framework against a purely data-driven approach
through the injection of physics features into intermediate LSTM layers. Our
quantitative analysis shows that the overall model uncertainty can be reduced
through the PGML approach especially for test data coming from a distribution
different than the training data. Moreover, we demonstrate that our approach
can be used as an inverse diagnostic tool providing a confidence score
associated with models and observations. The proposed framework also allows
for multi-fidelity computing by making use of low-fidelity models in the
online deployment of quantified data-driven models.

_K_ eywords Physics-guided machine learning, Reduced order modeling, Multi-
fidelity modeling

##  1 Introduction

Data-driven approaches are emerging as the new paradigm in computational
modeling of various problems in different branches of science and engineering,
including fluid dynamics [1, 2]. The universal approximation capability of
neural networks [3] makes them the powerful algorithm for complicated problems
like turbulence closure modeling [4], spatio-temporal super-resolution [5],
state-estimation [6], and nonlinear model order reduction [7]. One of the key
issues with deep learning is that they exhibit poor _generalizability_ , i.e.,
they produce inaccurate prediction when the test data is from a distribution
far from the training data. This adversely affects the trustworthiness of
neural networks for scientific applications and embedding the known physics to
enhance the generalizability is a challenge and opportunity for data-driven
algorithm developments [8, 9, 10, 11].

While many techniques have sought to enforce physics into data-driven models,
such as regularizing the neural network with the governing equations or
statistical constraints [12, 13, 14], imposing conservation of physical
quantities directly into neural network [15, 16], embedding invariance
property into neural network architecture [17, 18], and incorporating
inductive biases to respect exact conservation laws [19], they offer many
possibilities to fuse domain knowledge to improve the generalizability and
explainability of these models [11]. In this work, we propose the physics-
guided machine learning (PGML) framework that can answer the question of how
to hybridize data-driven (i.e., non-intrusive) and first-principles (i.e.,
intrusive) approaches that are often utilized in fluid dynamics applications.
The PGML framework is based on a composite modeling strategy where the
physics-based model is integrated within the neural network architecture
leading to a more generalizable learning engine.

The high-dimensional and multiscale nature of fluid flows makes the numerical
discretization methods computationally infeasible for many practical
applications. There are several techniques that aim at at constructing the
lower-order representation of high-dimensional systems that can capture the
essential features and are computationally order of magnitudes faster than
full order model [20, 21, 22]. These reduced order models (ROMs) are
particularly appealing for outer-loop applications, such as data assimilation
[23, 24, 25], uncertainty quantification [26], model predictive control [27,
28] that require multiple evaluations of a model for multiple inputs. Proper
orthogonal decomposition (POD) [29] and dynamic mode decomposition [30]
enabled ROMs are some of the most widely used methods among the fluid dynamics
community. POD is often combined with Galerkin projection (GP) to model the
dynamics of the ROM [31, 32]. One of the limitations of intrusive approaches
like GP is that it requires the complete knowledge about the equations
governing the system’s dynamics. However, for many physical systems, the exact
governing equations are not available due to imperfect knowledge about the
system. For example, in many geophysical flows we might not have accurate
models for processes such as wind forcing, bottom friction, stratification,
and other parameterization [33].

Lately, equation-free or non-intrusive reduced order modeling (NIROM) has
drawn great attention in the fluid mechanics’ community due to its flexibility
and efficiency for systems with incomplete dynamics and is undergoing rapid
development with various algorithms emerging from different perspectives [34].
In most of the NIROMs, the main idea is to employ deep learning methods to
construct nonlinear manifold [35, 7, 36] and to model the temporal dynamics of
ROMs [37, 38, 39, 40]. Despite the success of NIROMs for many complex
nonlinear problems, the naive deployment of NIROMs in multi-query applications
is limited because they lack connection with the physics-based model and might
suffer from poor performance in extrapolation regime [41]. Therefore, the
hybrid modeling approach that can overcome these drawbacks is warranted and to
this end, we apply the PGML framework that exploits the Galerkin projection
model for the known physics and long-short term memory (LSTM) neural network
to model the unknown physics. We remark here that the ideas from other
approaches like physics-reinforced neural network (PRNN) [42] and physics-
embedded convolutional autoencoder (PhyCAE) [16] can be easily integrated
within the proposed PGML framework.

In this study, we consider a dynamical system whose evolution is governed by
partial differential equations such as the Navier-Stokes equations as follows

$\dot{\boldsymbol{u}}=\boldsymbol{f}(\boldsymbol{u};\boldsymbol{x},t;\boldsymbol{\mu}),$
(1)

where $\boldsymbol{u}$ is the prognostic variable that depends on a set of
parameters, $\boldsymbol{f}$ is the nonlinear function that models the known
processes and conservation laws, and $\boldsymbol{\mu}\in\mathbb{R}^{N_{p}}$
is the parameterization such as Reynolds or Rayleigh numbers. However, there
could be unknown physics that might not be covered within the physics-based
model $\boldsymbol{f}$. The augmented representation can then be written at an
abstract level as

$\dot{\boldsymbol{u}}=\boldsymbol{f}(\boldsymbol{u};\boldsymbol{x},t;\boldsymbol{\mu})+\boldsymbol{\pi}(\boldsymbol{u};\boldsymbol{x},t;\boldsymbol{\mu}),$
(2)

where $\boldsymbol{\pi}$ encompasses all the unknown processes and
parameterizations. Here, we can define $\delta\boldsymbol{u}$ as a natural way
to quantify the unknown physics $\boldsymbol{\pi}$ as follows

$\delta\boldsymbol{u}=\sqrt{(||\boldsymbol{u}_{1}||-||\boldsymbol{u}_{2}||)^{2}}/\sqrt{||\boldsymbol{u}_{1}||^{2}+||\boldsymbol{u}_{2}||^{2}},$
(3)

where $\boldsymbol{u}_{1}$ is the solution of Eq. 2 without
$\boldsymbol{\pi}$, $\boldsymbol{u}_{2}$ is the solution of Eq. 2 with
$\boldsymbol{\pi}$, and $||\cdot||$ is the $L_{2}$ norm. One can easily
quantify the level of unknown processes using the definition given by Eq. 3,
i.e. $\delta\boldsymbol{u}\rightarrow 0$ when $\boldsymbol{\pi}\rightarrow 0$.
On the other hand, a higher value of $\delta\boldsymbol{u}$, i.e. close to
1.0, means that our known physical model $\boldsymbol{f}$ is poor.

The chief motivation behind the PGML framework is to exploit domain knowledge
to produce generalizable data-driven models. We demonstrate the application of
the PGML framework to fuse intrusive and non-intrusive parametric projection-
based ROMs for incomplete dynamical systems that can be represented using Eq.
2. Apart from improving the extrapolation performance of data-driven models,
the PGML framework can also be considered as an enabler for multi-fidelity
computing through a synergistic combination of low-cost/low-fidelity models
with the high-fidelity data. The field of multi-fidelity computing is gaining
attention in computational science and engineering, as it provides a way to
utilize computationally cheap approximate models with high-fidelity
observations [43, 44]. This is particularly relevant to fluid mechanics
applications where there is a wealth of simplified models like potential flow
theory, Blasius boundary layer model, and continuum Darcy/pipe flow model. One
of the major challenges in multi-fidelity computing is determining the
correlation between low and high fidelity models, and the neural networks
employed within the PGML framework is capable of discovering this nonlinear
relation. This first step in the assessment of the new PGML model shows that
the proposed hybrid reduced order modeling framework could represent a viable
tool for emerging digital twin applications [10].

The paper is structured as follows. We introduce the parametric ROM framework
in Section 2. We then discuss the a specific formulation of the PGML framework
in Section 3 for combining the domain knowledge with data. The numerical
results for the two-dimensional vorticity transport equation are detailed in
Section 4. Finally, Section 5 will present conclusions and ideas for the
future work.

##  2 Parametric ROM framework

The first step in parametric ROMs is the dimensionality reduction. We choose
proper orthogonal decomposition (POD) to compute the linear orthogonal basis
functions. POD is one of the most popular technique for dimensionality
reduction in the fluid mechanics community [45, 46]. POD provides the linear
basis functions that minimizes the error between the true data and its
projection in the $L_{2}$ sense compared to all linear basis functions of the
same dimension. Given the $N_{s}$ snapshots of the data for a variable of
interest in the physical system, we can form the matrix $\boldsymbol{A}$ as
follows

$\boldsymbol{A}=\bigg{[}\boldsymbol{u}^{(1)},\boldsymbol{u}^{(2)},\dots,\boldsymbol{u}^{(N_{s})}\bigg{]}\in\mathbb{R}^{N\times
N_{s}}$ (4)

where $\boldsymbol{u}^{(k)}\in\mathbb{R}^{N}$ corresponds to individual
snapshot of the solution in time. Then, we perform the singular value
decomposition (SVD) as follows

$\boldsymbol{A}=\boldsymbol{W}\boldsymbol{\Sigma}\boldsymbol{V}^{T}=\sum_{k=1}^{N_{s}}\sigma_{k}\boldsymbol{w}_{k}\boldsymbol{v}_{k}^{T},$
(5)

where $\boldsymbol{W}\in\mathbb{R}^{N\times N_{s}}$ is a matrix with
orthonormal columns which represent the left-singular vectors,
$\boldsymbol{V}\in\mathbb{R}^{N_{s}\times N_{s}}$ is matrix with orthonormal
columns representing the right-singular vectors, and
$\boldsymbol{\Sigma}\in\mathbb{R}^{N_{s}\times N_{s}}$ is a matrix with
positive diagonal entries, called singular values, arranged such that
$\sigma_{1}\geq\sigma_{2}\geq\dots\geq\sigma_{N_{s}}\geq 0$. In dimensionality
reduction, only first $N_{r}$ columns of $\boldsymbol{W}$ and $\boldsymbol{V}$
are retained along with the upper-left $N_{r}\times N_{r}$ sub-matrix of
$\boldsymbol{\Sigma}$. The singular values ${\sigma_{i}}$ give a measure of
the quality of information that is retailed in $N_{r}$-order approximation of
the matrix $\boldsymbol{A}$. The $N_{r}$ columns of $\boldsymbol{W}$ are
referred to as the POD modes or basis functions, denoted as
$\boldsymbol{\Phi}_{\boldsymbol{\mu}}=\\{\phi_{k}\\}_{k=1}^{N_{r}}$. The
amount of the energy retained by POD modes can be computed using the quantity
called the relative information content (RIC) as follows

$\text{RIC}(N_{r})=\frac{\sum_{j=1}^{N_{r}}\sigma_{j}^{2}}{\sum_{j=1}^{N_{s}}\sigma_{j}^{2}}$
(6)

We utilize the Grassmann manifold interpolation [47, 48] to compute the basis
functions for a new set of parameters from available POD basis function
computed in the offline phase. In the offline phase, we calculate a set of POD
basis functions $\\{\boldsymbol{\Phi}_{i}\\}_{i=1}^{N_{c}}$ corresponding to a
set of parameters $\\{\boldsymbol{\mu}_{i}\\}_{i=1}^{N_{c}}$. The Grassman
manifold interpolation consists of choosing a reference point $S_{0}$
corresponding to the basis set $\mathbf{\Phi}_{0}$, and then mapping each
point $S_{i}$ to a matrix $\mathbf{\Gamma}_{i}$ which represents a point on
the tangent space at $S_{0}$ using logarithmic map $\text{Log}_{S_{0}}$ as
follows

$(\boldsymbol{\Phi}_{i}-\boldsymbol{\Phi}_{0}\boldsymbol{\Phi}_{0}^{T}\boldsymbol{\Phi}_{i})(\boldsymbol{\Phi}_{0}^{T}\boldsymbol{\Phi}_{i})^{-1}=\boldsymbol{W}_{i}\boldsymbol{\Sigma}_{i}\boldsymbol{V}_{i}^{T},$
(7)
$\boldsymbol{\Gamma}_{i}=\boldsymbol{W}_{i}\text{tan}^{-1}(\boldsymbol{\Sigma}_{i})\boldsymbol{V}_{i}^{T}.$
(8)

The matrix $\boldsymbol{\Gamma}_{t}$ corresponding to target operating point
$\boldsymbol{\mu}_{t}$ is computed by interpolating the corresponding entries
of matrices $\\{\mathbf{\Gamma}_{i}\\}_{i=1}^{N_{c}}$. When $N_{p}=1$,
typically Lagrange-type interpolation method is applied as follows

$\mathbf{\Gamma}_{t}=\sum_{i=1}^{N_{c}}\bigg{(}\prod_{\begin{subarray}{c}j=1\\\
j\neq
i\end{subarray}}^{N_{c}}\frac{\mu_{t}-\mu_{j}}{\mu_{i}-\mu_{j}}\bigg{)}\mathbf{\Gamma}_{i}.$
(9)

The POD basis functions $\boldsymbol{\Phi}_{t}$ corresponding to the test
parameter $\boldsymbol{\mu}_{t}$ is computed using the exponential mapping as
follows

$\boldsymbol{\Gamma}_{t}=\boldsymbol{W}_{t}\boldsymbol{\Sigma}_{t}\boldsymbol{V}_{t}^{T},$
(10)
$\boldsymbol{\Phi}_{t}=[\boldsymbol{\Phi}_{0}\boldsymbol{V}_{t}\text{cos}(\boldsymbol{\Sigma}_{t})+\boldsymbol{W}_{t}\text{sin}(\boldsymbol{\Sigma}_{t})]\boldsymbol{V}_{t}^{T},$
(11)

where the trigonometric operators apply only to diagonal elements.

Once the basis functions corresponding to the test operating point are
computed, the reduced order representation of the high-dimensional state
$\boldsymbol{u}(\boldsymbol{x},t;\boldsymbol{\mu})$ is

$\boldsymbol{u}(\boldsymbol{x},t;\boldsymbol{\mu})\approx\boldsymbol{u}_{r}(\boldsymbol{x},t;\boldsymbol{\mu})=\boldsymbol{\Phi}_{\boldsymbol{\mu}}\boldsymbol{a}(t;\boldsymbol{\mu}),$
(12)

where $\boldsymbol{a}(t;\boldsymbol{\mu})\in\mathbb{R}^{N_{r}}$ is the reduced
state on the basis functions and is also called as the modal coefficients and
$\boldsymbol{\Phi}_{\boldsymbol{\mu}}$ is the POD basis functions for a set of
parameters $\boldsymbol{\mu}$. The ROM is obtained by substituting the low-
rank approximation given in Eq. 12 into the full order model defined in Eq. 2
and then taking the inner product with test basis functions to yield a system
of $N_{r}$ ordinary differential equations (ODEs). For many fluid mechanics
problems, the mechanistic description of some variables or processes is not
available or insufficient for the desired task (i.e., the term
$\boldsymbol{\pi}(\cdot)$ is unknown). Therefore, the physics-based intrusive
approaches like Galerkin projection ROM (GROM) provide us

$\dot{\boldsymbol{a}}=\boldsymbol{\Phi}_{\boldsymbol{\mu}}^{T}\boldsymbol{f}(\boldsymbol{\Phi}_{\boldsymbol{\mu}}\boldsymbol{a};\boldsymbol{x},t;\boldsymbol{\mu}).$
(13)

We note here that in GROM, the test basis functions are the same as the trial
basis which allows us to make use of the orthonormality, i.e.,
$\boldsymbol{\Phi}_{\boldsymbol{\mu}}^{T}\boldsymbol{\Phi}_{\boldsymbol{\mu}}=\boldsymbol{I}_{N_{r}}$.

One of the limitation with intrusive ROMs like GROM is that the governing
equations of the system have to be known exactly. If the governing equations
are not known exactly, the effect of unknown dynamics is truncated in GROMs
and accurate prediction is not achieved. This issue is mitigated in NIROMs
that exploit machine learning (ML) algorithms to learn the reduced order
dynamics from the observed data [34]. The training data for the ML algorithm
can be generated by using the same data used for computing the POD basis
functions as follows

$\boldsymbol{\alpha}(t;\boldsymbol{\mu})=\langle\boldsymbol{u}(\boldsymbol{x},t;\boldsymbol{\mu}),\boldsymbol{\Phi}_{\boldsymbol{\mu}}\rangle,$
(14)

where $\boldsymbol{\alpha}$ is the data projection coefficient, and the angle-
parentheses refer to the Euclidean inner product defined as
$\langle\boldsymbol{p},\boldsymbol{q}\rangle=\boldsymbol{p}^{T}\boldsymbol{q}$.
The neural network is one of the most successful algorithms for learning the
dynamics on reduced order basis functions. However, as the complexity of the
problems increases, the number of trainable parameters quickly explodes and
neural networks suffer from high epistemic uncertainty associated with limited
training data. We address this problem by hybridizing intrusive and non-
intrusive approaches through our PGML framework.

##  3 Physics-guided machine learning (PGML) framework

Figure 1: Proposed physics-guided machine learning (PGML) framework for
reduced order modeling, where the physics-based features are directly embedded
into hidden layers of the neural network along with learned latent variables.

In this Section, we detail different components of the proposed PGML framework
for modeling the reduced order representation of the high-dimensional
nonlinear dynamical system. Different components of the PGML framework are
illustrated in Figure 1. First, we present the formulation of the LSTM neural
network for time series modeling and then discuss its application within the
PGML framework. In a typical ML algorithm for time series modeling, the model
is trained on a time series of observable
$\\{\mathbf{o}^{(1)},\dots,\mathbf{o}^{(T)}\\}$, where
$\mathbf{o}\in\mathbb{R}^{d_{o}}$, sampled at a fixed time interval. The LSTM
neural network is one of the most popular ML algorithms for time-series
modeling due to its ability to model long-term temporal dependencies without
suffering from the vanishing gradient problem of recurrent neural network
[49]. Recently, LSTM has also been very successful in modeling high-
dimensional spatio-temporal chaotic systems and for ROMs [50, 51, 52, 53, 54].
The general functional form of the models used for time series forecasting can
be written as

$\displaystyle\mathbf{h}^{(t)}$
$\displaystyle=f_{h}^{h}(\mathbf{o}^{(t)},\mathbf{h}^{(t-1)}),$ (15)
$\displaystyle\tilde{\mathbf{o}}^{(t+1)}$
$\displaystyle=f_{h}^{o}(\mathbf{h}^{(t)}),$ (16)

where $\mathbf{o}^{(t)}$ is the current observable and
$\tilde{\mathbf{o}}^{(t+1)}$ is the forecast for the observable at the next
time step $\mathbf{o}^{(t+1)}$ and $\mathbf{h}^{(t)}\in\mathbb{R}^{d_{h}}$ is
the hidden state. Here. $f_{h}^{h}$ is the hidden-to-hidden mapping and
$f_{h}^{o}$ is the hidden-to-output mapping.

The LSTM mitigates the issue with vanishing (or exploding) gradient by
employing the gating mechanism that allows information to be forgotten. The
equations that implicitly define the hidden-to-hidden mapping from hidden
state from the previous time step (i.e., $\mathbf{h}^{(t-1)}$) and input
vector at the current time step (i.e., $\mathbf{o}^{(t)}$) to the forecast
hidden state (i.e., $\mathbf{h}^{(t)}$) can be written as

$\displaystyle\mathbf{g}^{(t)}_{f}$
$\displaystyle=\sigma_{f}(\mathbf{W}_{f}[\mathbf{h}^{(t-1)},\mathbf{o}^{(t)}]+\mathbf{b}_{f}),$
(17) $\displaystyle\mathbf{g}^{(t)}_{i}$
$\displaystyle=\sigma_{i}(\mathbf{W}_{i}[\mathbf{h}^{(t-1)},\mathbf{o}^{(t)}]+\mathbf{b}_{i}),$
(18) $\displaystyle\tilde{\mathbf{c}}^{(t)}$
$\displaystyle=\text{tanh}(\mathbf{W}_{c}[\mathbf{h}^{(t-1)},\mathbf{o}^{(t)}]+\mathbf{b}_{c}),$
(19) $\displaystyle\mathbf{c}^{(t)}$
$\displaystyle=\mathbf{g}^{(t)}_{f}\odot\mathbf{c}^{(t-1)}+\mathbf{g}^{(t)}_{i}\odot\tilde{\mathbf{c}}^{(t)},$
(20) $\displaystyle\mathbf{g}^{(t)}_{o}$
$\displaystyle=\sigma_{o}(\mathbf{W}_{o}[\mathbf{h}^{(t-1)},\mathbf{o}^{(t)}]+\mathbf{b}_{o}),$
(21) $\displaystyle\mathbf{h}^{(t)}$
$\displaystyle=\mathbf{g}^{(t)}_{o}\odot\text{tanh}(\mathbf{c}^{(t)}),$ (22)

where $\mathbf{g}^{(t)}_{f},\leavevmode\nobreak\
\mathbf{g}^{(t)}_{i},\leavevmode\nobreak\
\mathbf{g}^{(t)}_{o}\in\mathbb{R}^{d_{h}}$ are the forget gate, input gate,
and output gate, respectively. The $\mathbf{o}^{(t)}\in\mathbb{R}^{d_{o}}$ is
the input vector at time $t$, $\mathbf{h}^{(t)}\in\mathbb{R}^{d_{h}}$ is the
hidden state, $\mathbf{c}^{(t)}\in\mathbb{R}^{d_{h}}$ is the cell state,
$\mathbf{W}_{f},\leavevmode\nobreak\ \mathbf{W}_{i},\leavevmode\nobreak\
\mathbf{W}_{c},\leavevmode\nobreak\
\mathbf{W}_{o}\in\mathbb{R}^{d_{h}\times(d_{h}+d_{o})}$ are the weight
matrices, and $\mathbf{b}_{f},\leavevmode\nobreak\
\mathbf{b}_{i},\leavevmode\nobreak\ \mathbf{b}_{c},\leavevmode\nobreak\
\mathbf{b}_{o}\in\mathbb{R}^{d_{h}}$ are the bias vectors. The symbol $\odot$
denotes the element-wise multiplication, and $\sigma$ is the sigmoid
activation function. The activation functions
$\sigma_{f},\sigma_{i},\sigma_{o}$ are sigmoids. The hidden-to-output mapping
is given by a fully connected layer with a linear activation function as
follows

$\tilde{\mathbf{o}}^{(t+1)}=\mathbf{W}_{o}\mathbf{h}^{t},$ (23)

where $\mathbf{W}_{o}\in\mathbb{R}^{d_{o}\times d_{h}}$.

When we apply the LSTM neural network for building NIROMs, the reduced order
state of the system at a future time step, i.e.,
$\boldsymbol{\alpha}^{(t+1)}$, is learned as the function of a short history
of $d$ past temporally consecutive reduced order states as follows

$\boldsymbol{\alpha}^{(t+1)}=\mathcal{F}(\underbrace{\boldsymbol{z}^{(t)},\boldsymbol{z}^{(t-1)},\dots,\boldsymbol{z}^{(t-d+1)}}_{\boldsymbol{z}^{(t):(t-d+1)}};\boldsymbol{\theta}),$
(24)

where $\mathcal{F}(\cdot\leavevmode\nobreak\ ;\boldsymbol{\theta})$ is the
nonlinear function parameterized by a set of parameters $\boldsymbol{\theta}$,
and $\boldsymbol{z}$ refers to input features consisting of the POD modal
coefficients and a set of parameters governing the system, i.e.,
$\boldsymbol{z}\in\mathbb{R}^{N_{r}+N_{p}}$. The LSTM is trained using the
backpropagation through time (BPTT) algorithm and the parameters are optimized
as

$\boldsymbol{\theta}^{*}=\operatorname{arg\,min}_{\theta}\frac{1}{N-d+1}\sum_{n=d}^{N}||\mathcal{F}(\underbrace{\boldsymbol{z}^{(t):(t-d+1)}}_{\zeta};\boldsymbol{\theta})-\boldsymbol{\alpha}^{(t+1)}||_{2}^{2}.$
(25)

We employ $(N_{l}-1)$ LSTM layers and a single dense layer with a linear
activation function as the last layer. Therefore, the ML model can be written
as

$\mathcal{F}_{\text{ML}}(\boldsymbol{\zeta};\boldsymbol{\theta})=\boldsymbol{h}_{N_{l}}(\cdot;\boldsymbol{\Theta}_{N_{l}})\circ\cdots\circ\boldsymbol{h}_{2}(\cdot;\boldsymbol{\Theta}_{2})\circ\boldsymbol{h}_{1}(\boldsymbol{\zeta};\boldsymbol{\Theta}_{1}),$
(26)

where the output of each LSTM layer ($i=1,\dots,N_{l}-1$) is
$\boldsymbol{h}_{i}(\cdot;\boldsymbol{\Theta}_{i})\in\mathbb{R}^{d\times
N_{h}}$ and the last dense layer maps the final hidden state to the output,
i.e.,
$\boldsymbol{h}_{N_{l}}(\cdot;\boldsymbol{\Theta}_{N_{l}}):\mathbb{R}^{N_{h}}\rightarrow\mathbb{R}^{N_{r}}$.
Here, $N_{h}$ is the number of hidden units in the LSTM cell and we use the
constant number of hidden units across all LSTM layers.

In the PGML framework, the features extracted from the physics-based model are
embedded into the $i$th intermediate hidden layer along with the latent
variables as follows

$\mathcal{F}_{\text{PGML}}(\boldsymbol{\zeta};\boldsymbol{\theta})=\boldsymbol{h}_{N_{l}}(\cdot;\boldsymbol{\Theta}_{N_{l}})\circ\cdots\circ\underbrace{\mathcal{C}(\boldsymbol{h}_{i}(\cdot;\boldsymbol{\Theta}_{i}),\boldsymbol{a}^{(t):(t-d+1)})}_{\text{Physics
injection}}\circ\cdots\circ\boldsymbol{h}_{1}(\boldsymbol{\zeta};\boldsymbol{\Theta}_{1}),$
(27)

where $\mathcal{C(\cdot,\cdot)}$ represents the concatenation operation and
$\mathbf{a}$ are the physics-based features (i.e., the GROM modal
coefficients). Therefore, the output of $i$th layer will be in
$\mathbb{R}^{d\times(N_{r}+N_{h})}$. We highlight here that the physics-based
features are embedded only at a single layer of the neural network. Therefore,
at which layer shall the physics-based features be embedded becomes another
hyperparameter of the PGML framework. While some guidelines or rule of thumb
can be established based on the interplay between known and unknown physics, a
thorough systematic investigation is needed. Methods like layer-wise relevance
propagation can be adopted to bring interpretability to complex neural
networks in the PGML framework and we consider it as part of our future work.
The physics-based features assist the LSTM network in constraining the output
to a manifold of the physically realizable solution and leads to improved
generalizability for the extrapolation regime. Even though we demonstrate the
PGML framework for POD-ROM, we emphasize here that the PGML framework is
highly modular and can also be applied to nonlinear, cluster-based, and
network-based ROMs. For example, we do not need to restrict to only linear
basis construction methods like POD and the PGML framework can be easily
extended for nonlinear manifold generation methods like convolutional
autoencoders [35, 7, 55].

##  4 Vortex merging experiments

We apply the PGML framework to the two-dimensional vorticity transport
equation for the vortex merger problem as a prototypical test case. This
problem involves the study of the merging of two co-rotating vortices. The
two-dimensional vorticity transport equation can be written as

$\dfrac{\partial\omega}{\partial
t}+J(\omega,\psi)=\dfrac{1}{\text{Re}}\nabla^{2}\omega+\boldsymbol{\pi}(\boldsymbol{x},t;\boldsymbol{\mu}),$
(28)

$\displaystyle\nabla^{2}\psi=\omega,$ (29)

where $\omega$ is the vorticity defined as $\omega=\nabla\times\mathbf{u}$,
$\mathbf{u}=[u,v]^{T}$ is the velocity vector, $\psi$ is the streamfunction,
Re is the Reynolds number parameterizing the system and
$\boldsymbol{\pi}(\cdot)$ represent the unknown physics. We can rewrite Eq. 28
as

$\dfrac{\partial\omega}{\partial
t}=\dfrac{1}{\text{Re}}\nabla^{2}\omega-J(\omega,\psi)+\boldsymbol{\pi}(\boldsymbol{x},t;\boldsymbol{\mu}).$
(30)

The first two terms on the right hand side of Eq. 30, i.e., the linear
Laplacian term and the nonlinear Jacobian term represent the known physics of
the system and are defined as follows

$\displaystyle\nabla^{2}\omega$
$\displaystyle=\dfrac{\partial^{2}\omega}{\partial
x^{2}}+\dfrac{\partial^{2}\omega}{\partial y^{2}},$ (31) $\displaystyle
J(\omega,\psi)$ $\displaystyle=\dfrac{\partial\psi}{\partial
y}\dfrac{\partial\omega}{\partial x}-\dfrac{\partial\psi}{\partial
x}\dfrac{\partial\omega}{\partial y}.$ (32)

As discussed in Section 2, the reduced order approximation of the vorticity
and streamfunction field can be written as

$\displaystyle\omega(x,y,t)\approx\sum_{k=1}^{N_{r}}a_{k}(t)\phi_{k}^{\omega}(x,y),$
(33)
$\displaystyle\psi(x,y,t)\approx\sum_{k=1}^{N_{r}}a_{k}(t)\phi_{k}^{\psi}(x,y),$
(34)

We have omitted the parameter $\boldsymbol{\mu}$ in Eq. 33-34 for simplified
notation. The vorticity and streamfunction share the same time dependent modal
coefficients as they are both related by the kinematic relationship given in
Eq. 29. Furthermore, the POD basis functions of the streamfunction can be
obtained by solving the below Poisson equation

$\nabla^{2}\phi_{k}^{\psi}=\phi_{k}^{\omega}.$ (35)

The GROM for the vorticity transport equations considering only the known
physics can be written as

$\boldsymbol{\dot{a}}={\mathfrak{L}}\boldsymbol{a}+\boldsymbol{a}^{T}{\mathfrak{N}}\boldsymbol{a},$
(36)

or more explicitly,

$\dfrac{\text{d}a_{k}}{\text{d}t}={\sum_{i=1}^{N_{r}}\mathfrak{L}^{i}_{k}a_{i}+\sum_{i=1}^{N_{r}}\sum_{j=1}^{N_{r}}\mathfrak{N}^{ij}_{k}a_{i}a_{j}},$
(37)

where $\mathfrak{L}$ and $\mathfrak{N}$ are the linear and nonlinear
operators. The linear and nonlinear operators for the vorticity transport
equation are as follows

$\displaystyle\mathfrak{L}_{ik}$
$\displaystyle=\bigg{\langle}\dfrac{1}{\text{Re}}\nabla^{2}\phi_{i}^{\omega},\phi_{k}^{\omega}\bigg{\rangle},$
(38) $\displaystyle\mathfrak{N}_{ijk}$
$\displaystyle=\bigg{\langle}-J(\phi_{i}^{\omega},\phi_{j}^{\psi}),\phi_{k}^{\omega}\bigg{\rangle},$
(39)

where the angle-parentheses refer to the Euclidean inner product defined as
$\langle\mathbf{x},\mathbf{y}\rangle=\mathbf{x}^{T}\mathbf{y}=\sum_{i=1}^{N}x_{i}y_{i}$.
We apply the third-order Adams-Bashforth (AB3) numerical scheme to solve the
system of ODEs given in Eq. 37. In discrete sense, the update formula can be
written as

${a}_{k}^{(n+1)}={a}_{k}^{(n)}+{\Delta
t}\sum_{q=0}^{s}\beta_{q}G({a}_{k}^{(n-q)}),$ (40)

where $s$ and $\beta_{q}$ are the constants corresponding to AB3 scheme, which
are $s=2,\leavevmode\nobreak\ \beta_{0}=23/12,\leavevmode\nobreak\
\beta_{1}=-16/12,$ and $\beta_{2}=5/12$. Here, the operator $G({a}_{k})$ is
the right hand side of Eq. 37 as follows

$G(a_{k})={\sum_{i=1}^{N_{r}}\mathfrak{L}^{i}_{k}a_{i}+\sum_{i=1}^{N_{r}}\sum_{j=1}^{N_{r}}\mathfrak{N}^{ij}_{k}a_{i}a_{j}}.$
(41)

We start with an initial vorticity field of two Gaussian-distributed vortices,

$\displaystyle\omega(x,y,0)$
$\displaystyle=\exp\left(-\pi\left[(x-x_{1})^{2}+(y-y_{1})^{2}\right]\right)$
$\displaystyle+\exp{\left(-\pi\left[(x-x_{2})^{2}+(y-y_{2})^{2}\right]\right)},$
(42)

where their centers are initially located at $(x_{1},y_{1})=(3\pi/4,\pi)$ and
$(x_{2},y_{2})=(5\pi/4,\pi)$. We utilize an arbitrary array of decaying
Taylor-Green vortices as the unknown processes parameterized by
${\mu}=\\{\text{Re}\\}$ i.e., we have

$\boldsymbol{\pi}(x,y,t;\text{Re})=-\gamma\leavevmode\nobreak\
e^{-t/\text{Re}}\text{cos}(3x)\text{cos}(3y).$ (43)

The parameter $\gamma$ controls the strength of unknown physics compared to
the known physics. The synthetic data for building the ROM is generated by
solving Eq. 28 on the computational domain $(x,y)\in[0,2\pi]$ with periodic
boundary conditions and $256^{2}$ spatial discretization. The system is
evolved for 20 time units with a time-step size of $\Delta t=1\times 10^{-3}$.
The data snapshots for the training are obtained for
$\text{Re}=\\{200,400,600,800\\}$ and the trained models are evaluated for
$\text{Re}=\\{1000,1500\\}$. Specifically, the physical variable of interest
is the vorticity field, i.e., $\boldsymbol{u}=\\{\omega\\}$ and the unknown
source term $\boldsymbol{\pi}(\cdot)$ is given in Eq. 43. We retain $N_{r}=8$
POD modes for all Reynolds numbers and these POD modes captures more than 99%
of the energy of the system as illustrated in Fig. 2. Instantaneous snapshots
of the vorticity field for different Reynolds numbers at different time
instances are displayed in Fig. 3. We note here that the vorticity field shown
in Fig. 3 is obtained by solving the vorticity transport equation assuming
that there is no unknown physics, i.e., setting $\gamma=0.0$.

Figure 2: Amount of the kinetic energy captured by each POD modes.

Figure 3: Temporal evolution of the vorticity field at different time
instances for several Reynolds number investigated in this study. We note here
that the vorticity field is presented for $\gamma=0.0$.

The LSTM network in the ML model consists of $N_{l}=5$ layers with the first
four layers as the LSTM layers and the dense layer with a linear activation
function as the output layer. Each LSTM layer has $N_{h}=80$ hidden units and
tanh activation function is applied. From our experimentation with the
hyperparameters, we found that the results were not improving by employing
deeper network architecture and in some cases, deeper networks led to poor
performance due to overfitting. We utilize $d=3$, i.e., the state of the
system for three previous time steps are used to predict the future state.
Since we are using the information of only $d$ past temporally consecutive
states as the input, the LSTM can capture dependencies up to $d$ previous time
steps. During the online deployment, the initial condition for the first $d$
time steps is provided. This information is used to predict the forecast state
at $(d+1)$th time step. Then the state of the system from $2-(d+1)$ is used to
predict the forecast state at $(d+2)$th time step. This procedure is continued
until the final time step. The PGML framework also uses the same LSTM network
architecture and the features from the Galerkin projection model (i.e., GROM
modal coefficients) are embedded into the third LSTM layer. We reiterate here
that the layer at which the physics based features should be embedded is
another hyperparameter of PGML framework and methods like neural architecture
search can be used for complex problems [56]. The online deployment of the
PGML framework is similar to the ML framework. We characterize the
generalizability of the ML and PGML model through uncertainty estimate
computed using the deep ensembles method where ensemble of neural networks are
trained with different initialization of parameters [57, 58, 59]. In deep
ensembles method, the prediction from the ensemble of neural networks is
interpreted as the epistemic uncertainty. Despite the simplicity of this
method, it is appealing due to its scalability and strong empirical results
showing that the uncertainty estimate is as good as the Bayesian neural
networks [59]. The quantitative performance for ML and PGML framework is
measured through the statistics of the weighted root mean squared error (RMSE)
defined as

$\epsilon(i)=\sum_{k=1}^{N_{r}}\lambda_{k}\bigg{(}{\frac{1}{N_{t}}\sum_{n=1}^{N_{t}}(\alpha_{k}^{(n)}-\tilde{\alpha}_{k}^{(n)}(i))^{2}}\bigg{)}^{1/2},$
(44)

where $\lambda_{k}={\sigma_{k}^{2}}/({\sum_{j=1}^{N_{r}}\sigma_{j}^{2}})$
represents the energy in each POD mode and $\tilde{\alpha}_{k}(i)$ is the
prediction of the modal coefficients for $i$th ensemble. Each neural network
in the ensemble set is initiated by using a different initial seed. The
statistics of the weighted RMSE, i.e., the mean and the standard deviation are
computed as follows

$\bar{\epsilon}=\frac{1}{N_{e}}\sum_{i=1}^{N_{e}}\epsilon(i),\quad
s_{\epsilon}=\bigg{(}{\frac{\sum_{i=1}^{N_{e}}(\epsilon(i)-\bar{\epsilon})^{2}}{N_{e}-1}}\bigg{)}^{1/2},$
(45)

where $N_{e}$ is the size of ensemble of the neural network. We set $N_{e}=30$
for both ML and PGML models in our study.

Fig. 4 shows the temporal evolution of selected modal coefficients predicted
by true projection, GP, ML, and PGML models at $\text{Re}=1500$ for
$\gamma=0.01$. The ensemble-averaged modal coefficients predicted by the PGML
model matches very accurately with the true projection compared to the GP and
ML model. We also observe that the ensemble-average prediction by the ML model
is also improved compared to the GP model. However, the prediction of the ML
model is marred by high uncertainty compared to the PGML model prediction.
From Fig. 3, we can notice that the flow is topologically very similar at
different Reynolds numbers and the main difference is in the amplitude of the
vorticity field. Therefore, the task of extrapolation for this problem is not
very challenging for neural networks. However, in this study, we demonstrate
the proof of concept to show the benefit of PGML over a purely data-driven ML
model, and evaluating the PGML framework for topologically different data sets
is a part of our future work. In Fig. 5, the temporal evolution of the last
mode (i.e., $a_{8}$) predicted by true projection, GP, ML, and PGML model is
displayed for different magnitudes of the unknown source term. Overall, we can
observe that the PGML framework can predict the modal coefficients with high
confidence up to $\gamma=0.04$ and the prediction for $\gamma=0.1$ is less
certain. One of the reasons behind inaccurate prediction at $\gamma=0.1$ can
be due to a very poor correlation between the known physics in the GP model
and the actual physics of the system.

We also evaluate the capability of GP, ML, and PGML framework in
reconstructing the vorticity field using the root mean squared error (RMSE)
between the true vorticity field and the predicted vorticity field. The RMSE
is defined as

$\text{RMSE}(t)=\bigg{(}\frac{1}{N}\sum_{i=1}^{N}\bigg{(}\omega_{T}(\boldsymbol{x}_{i},t)-\omega_{R}(\boldsymbol{x}_{i},t)\bigg{)}^{2}\bigg{)}^{1/2},$
(46)

where $N$ is the spatial resolution (i.e., $N=N_{x}\times N_{y}$),
$\omega_{T}$ is the true projection of the vorticity field, and $\omega_{R}$
is the vorticity field predicted by ROM (computed using Eq. 33). The
$\omega_{R}$ is computed using the ensemble-averaged modal coefficients
predicted by ML and PGML model. Fig. 6 shows the time evolution of RMSE for
different magnitudes of the source term at Re = 1000 and Re = 1500. We can
observe that when the complete physics of the vortex-merging process is known
(i.e., $\gamma=0.0$), the RMSE for both ML and PGML framework is less than the
GP model at Re = 1000. However, at Re = 1500, the RMSE for the ML model is
higher than the GP model. This is due to the poor extrapolation capability of
a purely data-driven ML model. On the other hand, the PGML model can
successfully use the physics-based information from the GP model leading to
less RMSE compared to both GP and ML models. The vorticity field predicted by
the GP model quickly diverges from the true vorticity field in the presence of
unknown physics and overall the PGML prediction is better than the ML model as
illustrated in Fig. 6. The difference between the vorticity field at the final
time $t=20$ predicted by true projection modal coefficients and GP, ML, and
PGML models is depicted in Fig. 7.

Figure 4: Temporal evolution of selected modal coefficients for the vortex-
merger test case at Re = 1500 and $\gamma=0.01$. The ML (left) and PGML
(right) represent the average of the modal coefficients predicted by all
neural networks trained using MSE loss function with different seeds for
initialization of the parameters. Although the ensemble-averaged ML model
provides more accurate predictions than the GP model, there is a large
standard deviation (SD) band over all ML model predictions.

Figure 5: Temporal evolution of the last modal coefficient for the vortex-
merger test case at Re = 1500 and for different magnitudes of the source term.
The dashed red curve represent the average of the modal coefficients predicted
by ML model (left), the PGML model (right).

Figure 6: Temporal evolution of the RMSE (defined in Eq. 46) for the vortex-
merger test case at Re = 1000 (left) and Re = 1500 (right) for different
magnitudes of the source term. The dashed red curve represent the average of
the modal coefficients predicted by ML model (left), the PGML model (right). .

Figure 7: The error in the prediction of the vorticity field at the final time
$t=20$ for the GP model (top), ML model (middle), and PGML model (bottom).

Table 1 reports the quantitative analysis of GP, ML, and PGML models at
$\text{Re}=\\{1000,1500\\}$ for different magnitudes of the unknown source
term. The weighted RMSE mentioned in Table 1 are calculated using Eq. 44 and
its statistics is computed with Eq. 45. Here, $\gamma=0.0$ corresponds to the
special case where the physics is completely available to us. In that case,
there is no need to make a stochastic learning machine, so a physics-based
model (e.g., GP) will provide an accurate prediction. From Table 1 we see that
the weighted RMSE for the ML model is considerably higher for Re = 1500 at
$\gamma=0.0$ due to the poor extrapolation ability of pure data-driven models.
The weighted RMSE for the PGML model is slightly higher than the GP model even
though the reconstructed vorticity field is more accurate as illustrated in
Fig. 6. The weighted RMSE assigns more weightage to the first few modes and
therefore the inaccurate prediction for dominant modes leads to a higher value
of weighted RMSE for the PGML model. For higher values of $\gamma$, the
dynamics of the system is not known exactly and it gives us an opportunity to
introduce a stochastic neural network model. As demonstrated in Table 1, the
trustworthiness of the pure data-driven ML model is less due to the high value
of the uncertainty in the prediction (i.e., high standard deviation as shown
in parenthesis). However, once we incorporate physics-based features from the
GP model in the PGML framework, the uncertainty in the prediction is reduced
by almost one order magnitude. We also consider the hypothetical situation
where $\delta\boldsymbol{u}=0.6$ (so the physics is almost completely unknown
and as shown in Figure 7, unknown physics severely dominates the vortex
merging processes), then it is not surprising that the PGML framework would
not be in favor of the standard ML because of the injection of the
unrepresentative physics. Considering a false-negative analogy, this further
verifies the success of the PGML methodology and will automatically indicate
that the injected physics is not representative of the underlying processes.
This feature of the PGML framework could act as a diagnostic approach in many
digital twin applications to physical systems. There are numerous simulation
packages specifically designed for modeling and computing $\boldsymbol{f}$ in
these systems. On the other hand, exhaustive data sets become more and more
available to generate versatile machine learning based predictive models. To
this end, a confidence score $\delta\boldsymbol{u}$ defined in Eq. 3 can be
utilized to assess the relative modeling strength of principled and data-
driven models, where lower $\delta\boldsymbol{u}$ might indicate that the
principled model is a reliable tool and the physics-based features from the
principled model should be incorporated into data-driven models.

$\gamma$ |  $\delta\boldsymbol{u}$ | $\text{Re}=1000$ |  | $\text{Re}=1500$
---|---|---|---|---
|  |  GP |  ML |  PGML |  |  GP |  ML |  PGML
$0.00$ |  0.00 |  0.075 |  0.262 (0.190) |  0.064 (0.012) |  |  0.138 |  0.959 (0.632) |  0.192 (0.024)
$0.01$ |  0.03 |  1.564 |  0.258 (0.194) |  0.110 (0.012) |  |  1.563 |  0.852 (0.521) |  0.300 (0.029)
$0.02$ |  0.09 |  6.033 |  0.316 (0.113) |  0.181 (0.026) |  |  6.061 |  0.998 (0.458) |  0.498 (0.040)
$0.03$ |  0.17 |  12.809 |  0.465 (0.387) |  0.259 (0.040) |  |  12.905 |  1.352 (0.582) |  0.784 (0.092)
$0.04$ |  0.26 |  20.843 |  0.726 (0.591) |  0.312 (0.047) |  |  21.082 |  1.828 (0.992) |  1.116 (0.149)
$0.10$ |  0.60 |  70.935 |  1.910 (0.804) |  1.424 (0.356) |  |  71.923 |  5.252 (2.387) |  5.842 (6.227)
Table 1: The mean of weighted root mean square error (and its standard
deviation) in predicting the modal coefficients for higher Reynolds number
flows. The training has been performed using data for
$\text{Re}=\\{200,400,600,800\\}$. Note that the control parameter $\gamma$
adjusts the relative strength between known-physics and unknown-physics (i.e.,
$\gamma=0$ refers to the special case when physics is fully known), and
$\delta\boldsymbol{u}$ quantifies the level of unknown processes.

##  5 Discussion and conclusion

Although advanced machine learning (ML) models like deep learning networks are
powerful tools for finding patterns in complicated datasets, the number of
trainable parameters (weights) quickly explodes as these neural network models
become complex, adversely affecting their interpretability and hence their
trustworthiness. Our chief aim in this study is to illustrate how physics can
be injected into these neural networks to improve the trustworthiness of the
overall system. With this in mind, we introduce a physics-guided ML (PGML)
framework that fuses first-principles with neural-network models, and explore
complementary physics versus statistics issues that arise in the development
of the PGML framework. The PGML framework puts particular emphasis on the
physics-based features and embeds them at an intermediate layer of the neural
network. The robustness of the PGML framework is successfully demonstrated for
the vortex-merging process in the presence of an array of Taylor-Green
vortices as the unknown source term. Our results indicate that the PGML
framework can give considerably accurate prediction in the extrapolation
regime compared to its counterpart, i.e., pure data-driven ML model. The
physics-based features from the GP model also ensures that the model
uncertainty in the neural network prediction is reduced. Finally, we also
emphasize that the PGML framework is highly modular, and it naturally allows
for the multi-fidelity model fusion in different branches of science and
engineering.

The natural extension of the present study is to assess the PGML framework for
much more complex flows, where there is a significant variation between the
train and test data. Other unanswered questions are understanding how neural
network assigns importance to physics-based features, at which layer shall the
physics-based features be embedded, and methods from machine learning
community can be utilized for this [60, 61]. Another direction for future work
will be to use other uncertainty quantification methods like Monte Carlo
dropouts, Bayesian neural networks to quantify model uncertainty.

## Acknowledgements

This material is based upon work supported by the U.S. Department of Energy,
Office of Science, Office of Advanced Scientific Computing Research under
Award Number DE-SC0019290. O.S. gratefully acknowledges their support.
Disclaimer: This report was prepared as an account of work sponsored by an
agency of the United States Government. Neither the United States Government
nor any agency thereof, nor any of their employees, makes any warranty,
express or implied, or assumes any legal liability or responsibility for the
accuracy, completeness, or usefulness of any information, apparatus, product,
or process disclosed, or represents that its use would not infringe privately
owned rights. Reference herein to any specific commercial product, process, or
service by trade name, trademark, manufacturer, or otherwise does not
necessarily constitute or imply its endorsement, recommendation, or favoring
by the United States Government or any agency thereof. The views and opinions
of authors expressed herein do not necessarily state or reflect those of the
United States Government or any agency thereof.

## References

  * [1] Steven L Brunton, Bernd R Noack, and Petros Koumoutsakos.  Machine learning for fluid mechanics.  Annual Review of Fluid Mechanics, 52:477–508, 2020.
  * [2] MP Brenner, JD Eldredge, and JB Freund.  Perspective on machine learning for advancing fluid mechanics.  Physical Review Fluids, 4(10):100501, 2019.
  * [3] Kurt Hornik, Maxwell Stinchcombe, and Halbert White.  Multilayer feedforward networks are universal approximators.  Neural Networks, 2(5):359–366, 1989.
  * [4] Karthik Duraisamy, Gianluca Iaccarino, and Heng Xiao.  Turbulence modeling in the age of data.  Annual Review of Fluid Mechanics, 51:357–377, 2019.
  * [5] Kai Fukami, Koji Fukagata, and Kunihiko Taira.  Machine-learning-based spatio-temporal super resolution reconstruction of turbulent flows.  Journal of Fluid Mechanics, 909:A9, 2021.
  * [6] Nirmal J Nair and Andres Goza.  Leveraging reduced-order models for state estimation using deep learning.  Journal of Fluid Mechanics, 897:R1, 2020.
  * [7] Kookjin Lee and Kevin T Carlberg.  Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders.  Journal of Computational Physics, 404:108973, 2020.
  * [8] Anuj Karpatne, Gowtham Atluri, James H Faghmous, Michael Steinbach, Arindam Banerjee, Auroop Ganguly, Shashi Shekhar, Nagiza Samatova, and Vipin Kumar.  Theory-guided data science: A new paradigm for scientific discovery from data.  IEEE Transactions on knowledge and data engineering, 29(10):2318–2331, 2017.
  * [9] Markus Reichstein, Gustau Camps-Valls, Bjorn Stevens, Martin Jung, Joachim Denzler, Nuno Carvalhais, et al.  Deep learning and process understanding for data-driven Earth system science.  Nature, 566(7743):195–204, 2019.
  * [10] Adil Rasheed, Omer San, and Trond Kvamsdal.  Digital twin: Values, challenges and enablers from a modeling perspective.  IEEE Access, 8:21980–22012, 2020.
  * [11] K Kashinath, M Mustafa, A Albert, JL Wu, C Jiang, S Esmaeilzadeh, K Azizzadenesheli, R Wang, A Chattopadhyay, A Singh, et al.  Physics-informed machine learning: case studies for weather and climate modelling.  Philosophical Transactions of the Royal Society A, 379(2194):20200093, 2021.
  * [12] Maziar Raissi, Paris Perdikaris, and George E Karniadakis.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  Journal of Computational Physics, 378:686–707, 2019.
  * [13] Jin-Long Wu, Karthik Kashinath, Adrian Albert, Dragos Chirila, Heng Xiao, et al.  Enforcing statistical constraints in generative adversarial networks for modeling chaotic dynamical systems.  Journal of Computational Physics, 406:109209, 2020.
  * [14] Nicholas Geneva and Nicholas Zabaras.  Modeling the dynamics of pde systems with physics-constrained deep auto-regressive networks.  Journal of Computational Physics, 403:109056, 2020.
  * [15] Tom Beucler, Stephan Rasp, Michael Pritchard, and Pierre Gentine.  Achieving conservation of energy in neural network emulators for climate modeling.  arXiv preprint arXiv:1906.06622, 2019.
  * [16] Arvind T Mohan, Nicholas Lubbers, Daniel Livescu, and Michael Chertkov.  Embedding hard physical constraints in neural network coarse-graining of 3d turbulence.  arXiv preprint arXiv:2002.00021, 2020.
  * [17] Julia Ling, Andrew Kurzawski, and Jeremy Templeton.  Reynolds averaged turbulence modelling using deep neural networks with embedded invariance.  Journal of Fluid Mechanics, 807:155–166, 2016.
  * [18] Laure Zanna and Thomas Bolton.  Data-driven equation discovery of ocean mesoscale closures.  Geophysical Research Letters, 47(17):e2020GL088376, 2020.
  * [19] Sam Greydanus, Misko Dzamba, and Jason Yosinski.  Hamiltonian neural networks.  arXiv preprint arXiv:1906.01563, 2019.
  * [20] Clarence W Rowley and Scott TM Dawson.  Model reduction for flow analysis and control.  Annual Review of Fluid Mechanics, 49:387–417, 2017.
  * [21] Peter Benner, Serkan Gugercin, and Karen Willcox.  A survey of projection-based model reduction methods for parametric dynamical systems.  SIAM Review, 57(4):483–531, 2015.
  * [22] Kunihiko Taira, Maziar S Hemati, Steven L Brunton, Yiyang Sun, Karthik Duraisamy, Shervin Bagheri, Scott TM Dawson, and Chi-An Yeh.  Modal analysis of fluid flows: Applications and outlook.  AIAA Journal, page DOI:10.2514/1.J058462, 2019.
  * [23] Yanhua Cao, Jiang Zhu, I Michael Navon, and Zhendong Luo.  A reduced-order approach to four-dimensional variational data assimilation using proper orthogonal decomposition.  International Journal for Numerical Methods in Fluids, 53(10):1571–1583, 2007.
  * [24] Răzvan Ştefănescu, Adrian Sandu, and Ionel Michael Navon.  POD/DEIM reduced-order strategies for efficient four dimensional variational data assimilation.  Journal of Computational Physics, 295:569–595, 2015.
  * [25] Marcello Meldi and Alexandre Poux.  A reduced order model based on kalman filtering for sequential data assimilation of turbulent flows.  Journal of Computational Physics, 347:207–234, 2017.
  * [26] Sofia Guzzetti, LA Mansilla Alvarez, PJ Blanco, Kevin Thomas Carlberg, and A Veneziani.  Propagating uncertainties in large-scale hemodynamics models via network uncertainty quantification and reduced-order modeling.  Computer Methods in Applied Mechanics and Engineering, 358:112626, 2020.
  * [27] Carlos E Garcia, David M Prett, and Manfred Morari.  Model predictive control: Theory and practice—a survey.  Automatica, 25(3):335–348, 1989.
  * [28] Svein Hovland, Jan Tommy Gravdahl, and Karen E Willcox.  Explicit model predictive control for large-scale systems via model reduction.  Journal of guidance, control, and dynamics, 31(4):918–926, 2008\.
  * [29] Philip Holmes, John L Lumley, and Gal Berkooz.  Turbulence, coherent structures, dynamical systems and symmetry.  Cambridge University Press, New York, 1998.
  * [30] J Nathan Kutz, Steven L Brunton, Bingni W Brunton, and Joshua L Proctor.  Dynamic mode decomposition: data-driven modeling of complex systems.  SIAM, 2016.
  * [31] D Rempfer.  On low-dimensional Galerkin models for fluid flow.  Theoretical and Computational Fluid Dynamics, 14(2):75–88, 2000\.
  * [32] Clarence W Rowley, Tim Colonius, and Richard M Murray.  Model reduction for compressible flows using POD and Galerkin projection.  Physica D: Nonlinear Phenomena, 189(1-2):115–129, 2004.
  * [33] Vladimir M Krasnopolsky and Michael S Fox-Rabinovitz.  Complex hybrid models combining deterministic and machine learning components for numerical climate modeling and weather prediction.  Neural Networks, 19(2):122–134, 2006.
  * [34] Jian Yu, Chao Yan, and Mengwu Guo.  Non-intrusive reduced-order modeling for fluid problems: A brief review.  Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering, 233(16):5896–5912, 2019.
  * [35] Francisco J Gonzalez and Maciej Balajewicz.  Deep convolutional recurrent autoencoders for learning low-dimensional feature dynamics of fluid systems.  arXiv preprint arXiv:1808.01346, 2018.
  * [36] Romit Maulik, Bethany Lusch, and Prasanna Balaprakash.  Reduced-order modeling of advection-dominated systems with recurrent neural networks and convolutional autoencoders.  Physics of Fluids, 33(3):037106, 2021.
  * [37] Jiayang Xu and Karthik Duraisamy.  Multi-level convolutional autoencoder networks for parametric prediction of spatio-temporal dynamics.  Computer Methods in Applied Mechanics and Engineering, 372:113379, 2020.
  * [38] Sk M Rahman, Suraj Pawar, Omer San, Adil Rasheed, and Traian Iliescu.  Nonintrusive reduced order modeling framework for quasigeostrophic turbulence.  Physical Review E, 100(5):053306, 2019.
  * [39] Zheng Wang, Dunhui Xiao, Fangxin Fang, Rajesh Govindan, Christopher C Pain, and Yike Guo.  Model identification of reduced order fluid dynamics systems using deep learning.  International Journal for Numerical Methods in Fluids, 86(4):255–268, 2018.
  * [40] Mengwu Guo and Jan S Hesthaven.  Data-driven reduced order modeling for time-dependent problems.  Computer methods in applied mechanics and engineering, 345:75–99, 2019.
  * [41] Hugo F. S. Lui and William R. Wolf.  Construction of reduced-order models for fluid flows using deep feedforward neural networks.  Journal of Fluid Mechanics, 872:963–994, 2019.
  * [42] Wenqian Chen, Qian Wang, J. S. Hesthaven, and C. Zhang.  Physics-informed machine learning for reduced-order modeling of nonlinear problems.  2020\.
  * [43] Benjamin Peherstorfer, Karen Willcox, and Max Gunzburger.  Survey of multifidelity methods in uncertainty propagation, inference, and optimization.  SIAM Review, 60(3):550–591, 2018.
  * [44] Paris Perdikaris, Maziar Raissi, Andreas Damianou, Neil D Lawrence, and George Em Karniadakis.  Nonlinear information fusion algorithms for data-efficient multi-fidelity modelling.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 473(2198):20160751, 2017.
  * [45] Lawrence Sirovich.  Turbulence and the dynamics of coherent structures. i. coherent structures.  Quarterly of applied mathematics, 45(3):561–571, 1987.
  * [46] Gal Berkooz, Philip Holmes, and John L Lumley.  The proper orthogonal decomposition in the analysis of turbulent flows.  Annual review of fluid mechanics, 25(1):539–575, 1993.
  * [47] David Amsallem and Charbel Farhat.  Interpolation method for adapting reduced-order models and application to aeroelasticity.  AIAA Journal, 46(7):1803–1813, 2008.
  * [48] Ralf Zimmermann, Benjamin Peherstorfer, and Karen Willcox.  Geometric subspace updates with applications to online adaptive nonlinear model reduction.  SIAM Journal on Matrix Analysis and Applications, 39(1):234–261, 2018.
  * [49] Sepp Hochreiter and Jürgen Schmidhuber.  Long short-term memory.  Neural computation, 9(8):1735–1780, 1997.
  * [50] Jaideep Pathak, Brian Hunt, Michelle Girvan, Zhixin Lu, and Edward Ott.  Model-free prediction of large spatiotemporally chaotic systems from data: A reservoir computing approach.  Physical Review Letters, 120(2):024102, 2018.
  * [51] Pantelis R Vlachas, Wonmin Byeon, Zhong Y Wan, Themistoklis P Sapsis, and Petros Koumoutsakos.  Data-driven forecasting of high-dimensional chaotic systems with long short-term memory networks.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 474(2213):20170844, 2018.
  * [52] Romit Maulik, Bethany Lusch, and Prasanna Balaprakash.  Non-autoregressive time-series methods for stable parametric reduced-order models.  Physics of Fluids, 32(8):087115, 2020.
  * [53] Pantelis R Vlachas, Jaideep Pathak, Brian R Hunt, Themistoklis P Sapsis, Michelle Girvan, Edward Ott, and Petros Koumoutsakos.  Backpropagation algorithms and reservoir computing in recurrent neural networks for the forecasting of complex spatiotemporal dynamics.  Neural Networks, 2020.
  * [54] Kazuto Hasegawa, Kai Fukami, Takaaki Murata, and Koji Fukagata.  Cnn-lstm based reduced order modeling of two-dimensional unsteady flows around a circular cylinder at different reynolds numbers.  Fluid Dynamics Research, 52(6):065501, 2020.
  * [55] Hamidreza Eivazi, Hadi Veisi, Mohammad Hossein Naderi, and Vahid Esfahanian.  Deep neural networks for nonlinear model order reduction of unsteady flows.  Physics of Fluids, 32(10):105104, 2020.
  * [56] Romit Maulik, Romain Egele, Bethany Lusch, and Prasanna Balaprakash.  Recurrent neural network architecture search for geophysical emulation.  arXiv preprint arXiv:2004.10928, 2020.
  * [57] Robert Tibshirani.  A comparison of some error estimates for neural network models.  Neural Computation, 8(1):152–163, 1996.
  * [58] Tom Heskes.  Practical confidence and prediction intervals.  In Advances in neural information processing systems, pages 176–182, 1997.
  * [59] Balaji Lakshminarayanan, Alexander Pritzel, and Charles Blundell.  Simple and scalable predictive uncertainty estimation using deep ensembles.  In Proceedings of the 31st International Conference on Neural Information Processing Systems, NIPS’17, page 6405–6416, Red Hook, NY, USA, 2017. Curran Associates Inc.
  * [60] Haakon Robinson, Adil Rasheed, and Omer San.  Dissecting deep neural networks.  arXiv preprint arXiv:1910.03879, 2019.
  * [61] Grégoire Montavon, Alexander Binder, Sebastian Lapuschkin, Wojciech Samek, and Klaus-Robert Müller.  Layer-wise relevance propagation: an overview.  Explainable AI: interpreting, explaining and visualizing deep learning, pages 193–209, 2019.


[a]Shinji Ejiri

# Particle density probability distribution function and center symmetry
breaking in finite density lattice gauge theories

###### Abstract

We study the nature of the phase transition at high temperature and high
density in lattice gauge theories by focusing on the probability distribution
function, which represents the probability that a certain density will be
realized in a heat bath. The probability distribution function is obtained by
constructing a canonical partition function by fixing the number of particles
from the grand partition function. However, if the $Z_{3}$ center symmetry,
which is important for understanding the finite temperature phase transition
of $SU(3)$ lattice gauge theory, is maintained on a finite lattice, the
probability distribution function is always zero, except when the number of
particles is a multiple of 3. For $U(1)$ gauge theory, this problem is more
extreme. The probability distribution becomes zero when the particle number is
not zero. In this study, we find a solution to this problem and propose a
method of avoiding the sign problem, which is an important problem at finite
density, using the center symmetry. This problem is essentially the same as
the problem that the expectation value of the Polyakov loop is always zero
when calculating with finite volume. In the case of $U(1)$ lattice gauge
theory with heavy fermions, numerical simulations are actually performed, and
we demonstrate that the probability distribution function at a finite density
can be calculated by the method proposed in this study.

##  1 Introduction

Theoretical study of the probability distribution of particle density is
important for investigating the aspect of density fluctuation in the heat bath
generated by heavy ion collision experiments and finding the critical point at
finite density. We develop a method for calculating the probability
distribution function by the first-principle calculation of lattice QCD. The
probability distribution function is obtained by constructing a canonical
partition function by fixing the number of particles from the grand partition
function. The relation between the grand partition function $Z_{GC}(T,\mu)$
with the chemical potential $\mu$ and the canonical partition function
$Z_{C}(T,N)$ with the particle number $N$ is given by the fugacity expansion,

$\displaystyle Z_{GC}(T,\mu)=\sum_{N}Z_{C}(T,N)e^{N\mu/T}.$ (1)

The left-hand side of this equation $Z_{GC}(T,\mu)$ is the normalization
factor of the Boltzmann weight, and is classified by $N$ in the right-hand
side. Hence, $Z_{C}(T,N)e^{N\mu/T}$ can be regarded as a weight factor for
each $N$, and the probability distribution $W(N)$ is in proportion to
$Z_{C}(T,N)e^{N\mu/T}$. Thus, the effective potential of $N$ can be defined as
$-\ln W(N)=-\ln Z_{C}(T,N)-N\mu/T$.

However, if the $Z_{3}$ center symmetry, which spontaneously breaks in the
deconfinement phase transition of $SU(3)$ lattice gauge theory, is strictly
maintained on a finite lattice, the canonical partition function of QCD
$Z_{C}(T,N)$ is zero when the number of particles is not a multiple of 3.
Moreover, in the case of $U(1)$ lattice gauge theory, the situation is
extreme. The partition function $Z_{C}(T,N)$ is exactly zero for $N\neq 0$.
This means that the existence probability of charged particles is zero, which
is unacceptable. To solve this problem, we focus on $U(1)$ lattice gauge
theory. In the calculation of $Z_{C}$, it is essential to break the center
symmetry adding a small external field. Then, we find that the central
symmetry can be used to avoid sign problems.

In the next section, we explain the canonical partition function when the
theory has the center symmetry. We calculate the canonical partition function
for $U(1)$ gauge theory with heavy dynamical fermions in Sec. 3. We also
propose a method to avoid the sign problem using the $U(1)$ center symmetry.
Our conclusions are given in Sec. 4

##  2 Center symmetry breaking in $U(1)$ gauge theory

Figure 1: Distribution of Polyakov loops in the complex plane generated by
simulations of $U(1)$ lattice gauge theory at $\beta=0.90$ (left) and 1.10
(right). The temporal lattice size $N_{t}$ is 4. Purple, blue, green and red
symbols are the results of the spatial lattice size $N_{s}=12$, 16, 24 and 32,
respectively.

The grand partition function is

$\displaystyle Z_{GC}(T,\mu)=\int{\cal D}U(\det M(\kappa,\mu))^{N_{\rm f}}\
e^{-S_{g}}$ (2)

for the $N_{\rm f}$ degenerate flavor case. Here, $S_{g}$ is the gauge action
and $M(\kappa,\mu)$ is the fermion kernel. The fermion determinant $\det M$
can be expressed as the sum of Wilson loops. We perform a hopping parameter
expansion of $\ln\det M$. Then, the $n^{\rm th}$-order Taylor expansion
coefficients are given by the sum of $n$-step connected Wilson loops [1, 2, 3,
4, 5]. When the hopping parameter $\kappa$ is small, the leading order
contribution consists of the plaquette $P$, and the Polyakov loop $\Omega$ [2,
5]:

$\displaystyle\ln\det M(\kappa,\mu)=96N_{\rm
site}\kappa^{4}P+2^{N_{t}+1}N_{s}^{3}\left(\kappa^{N_{t}}e^{\mu/T}\Omega+\kappa^{N_{t}}e^{\mu/T}\Omega^{*}\right)+\cdots$
(3)

for the standard Wilson fermion on a lattice with the spatial size $N_{s}$ and
temporal size $N_{t}$. At finite density with the chemical potential $\mu$,
the $\mu$-dependence of each Wilson loop term, which is wound $m$ times by the
anti-periodic boundary condition in the time direction, appears as a factor
$e^{m\mu/T}$ because the hopping term in the time direction is
$(1-\gamma_{4})U_{x,4}e^{\mu a}$ or $(1+\gamma_{4})U_{x,4}^{\dagger}e^{-\mu
a}$ and $(e^{\mu a})^{N_{t}}=e^{\mu/T}$, where $U_{x,4}$ is a link variable
and $a$ is the lattice spacing. We classify these expansion terms by the
winding number $m$. If we reconstruct the expansion of $\ln\det M$ into the
expansion of $\det M$, $Z_{GC}$ is expressed in the form of the fugacity
expansion, Eq. (1). Therefore, the fugacity expansion is a winding number
expansion.

Under the $U(1)$ center transformation: the time components of all link
variables in one time slice are changed as $U_{(\vec{x},t),4}\to
e^{i\theta}U_{(\vec{x},t),4}$, $Z_{C}(T,N)$ changes as

$\displaystyle Z_{C}(T,N)\to e^{iN\theta}Z_{C}(N,T),$ (4)

since $Z_{C}$ is multiplied by $e^{i\theta}$ $N$ times at the time slice in
total. Because $S_{g}$ and the integral measure are invariant,
$Z_{C}(T,N)=e^{iN\theta}Z_{C}(N,T)$. Thus, the canonical partition function is
zero,

$\displaystyle
Z_{C}(T,N)=\frac{1}{2\pi}\int_{0}^{2\pi}e^{iN\theta}Z_{C}(N,T)d\theta=0,$ (5)

except for $N=0$. This means that there cannot exist particles that interact
with the gauge field. In the confined phase, the probability of existence of
charged particles is zero because of the confinement. However, charged
particles must be exist in the deconfinement phase. $Z_{C}(T,N)=0$ is
unacceptable. In such a case, the grand partition function $Z_{GC}(T,\mu)$ is
equal to $Z_{C}(T,0)$, and $Z_{GC}$ does not depend on the chemical potential.

This problem is a common problem when calculating the expectation value of the
order parameter for investigating spontaneous symmetry breaking. As seen in
Eq. (3), the first nontrivial term of $\ln\det M$ is proportional to $\Omega\
e^{\mu/T}$ when $\kappa$ is small. Thus, the canonical partition of $N=1$ with
a small $\kappa$ is proportional to $\langle\Omega\rangle$, where the
expectation value is computed in a quenched simulation. We perform quenched
QCD simulations. Figure 1 shows the distribution of the Polyakov loop in the
complex plane. The left panel shows the result of the confinement phase
(symmetric phase) at an inverse gauge coupling $\beta=1/g^{2}=0.90$, and the
right panel shows the result of the deconfinement phase (broken phase) at
$\beta=1.10$. The temporal lattice size $N_{t}$ is fixed to be 4 and four
spatial lattice sizes $N_{s}$ are adopted. Purple, blue, green and red lines
are $N_{s}=12$, 16, 24 and 32, respectively.

The Polyakov loop is an order parameter of the deconfinement phase transition.
In the deconfinement phase, the center symmetry is spontaneous broken and
$\langle\Omega\rangle$ should be nonzero. However, because of the center
symmetry, the probability distribution is symmetric under $U(1)$
transformation: $\Omega\to e^{i\theta}\Omega$ for an arbitrary real number
$\theta$, as shown in Fig. 1. Therefore, the expectation value of $\Omega$ is
always zero even in the broken phase. The symmetry does not break in an actual
simulation with finite volume. Therefore, to discuss spontaneous symmetry
breaking, it is required to break the center symmetry adding an explicit
breaking term in the action. Then, the breaking term dependence and spatial
volume dependence are investigated. If $\langle\Omega\rangle$ is non-zero in
the double limit of zero breaking term and infinite volume, we identify that
spontaneous symmetry breaking has occurred.

Figure 2: Expectation value of the Polyakov loop as a function of the hopping
parameter $\kappa$ measured at $\beta=0.90$ (left) and 1.10 (right) on
$N_{s}^{3}\times 4$ lattices with $N_{s}=12$ (purple), 16 (blue), 24 (green)
and 32 (red).

To break the center symmetry, we add a heavy fermion having an infinitesimal
$\kappa$. We approximate the fermion determinant with the leading term of the
hopping parameter expansion given in Eq. (3) . The plaquette term can be
absorbed into the gauge action by shifting $\beta\to\beta^{*}=\beta+16N_{\rm
f}\kappa^{4}$. The expectation value is computed by the reweighting method
where the reweighting factor is the Polyakov loop term,

$\displaystyle\langle{\rm Re}\Omega\rangle\approx\frac{1}{Z_{GC}}\int{\cal
D}U{\rm Re}\ \Omega\ e^{\epsilon V{\rm
Re}\Omega}e^{-S_{g}}=\frac{\left\langle{\rm Re}\Omega\ e^{\epsilon V{\rm
Re}\Omega}\right\rangle}{\left\langle e^{\epsilon V{\rm
Re}\Omega}\right\rangle},$ (6)

where $\epsilon=4\times 2^{N_{t}}\kappa^{N_{t}}$, $V=N_{s}^{3}$ and
$\langle\cdots\rangle$ means the average over quenched configurations [1].
Using the data of $N_{t}=4$ in Fig. 1, the expectation value of ${\rm
Re}\Omega$ is computed. The number of Monte-Caro updates for each $\beta$ is
200,000 for $N_{s}=12,16$ and 1,000,000 for $N_{s}=24,32$. We plot the results
in Fig. 2 as a function of $\kappa^{N_{t}}$. The left figure is the result of
the confinement phase at $\beta=0.90$. Although the results of $N_{s}=12$
(purple), 16 (blue), 24 (green) and 32 (red) are plotted, no volume dependence
is seen. For this case, it is not necessary to take the volume infinity limit,
thus $\langle{\rm Re}\Omega\rangle\sim\kappa^{N_{t}}$ and $\langle{\rm
Re}\Omega\rangle$ vanishes in the $\kappa\to 0$ limit in the confinement phase
(symmetric phase). On the other hand, the right figure is the result of the
deconfinement phase at $\beta=1.10$. In the deconfinement phase (broken
phase), Polyakov loop behaves as $\langle{\rm Re}\Omega\rangle\sim
V\kappa^{N_{t}}$. Of course, if $V$ is fixed, $\langle{\rm Re}\Omega\rangle$
becomes always zero in the limit of $\kappa=0$. However, in the thermodynamic
limit $V\to\infty$ and $\kappa\to 0$ limit, this figure suggests that
$\langle{\rm Re}\Omega\rangle$ becomes a non-trivial finite value in the
broken phase.

Furthermore, when $\kappa$ is sufficiently small, the Polyakov loop can be
evaluated by a $\kappa^{N_{t}}$ expansion. Using the distribution function of
the Polyakov loop in the complex plane, the expectation value can be
calculated as follows,

$\displaystyle\langle{\rm Re}\Omega^{n}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{Z_{GC}}\int{\cal D}U\ {\rm Re}\Omega^{n}\ e^{\epsilon
V{\rm Re}\Omega}e^{-S_{g}}=\int|\Omega|^{N}\cos(N\phi)\ e^{\epsilon
V|\Omega|\cos\phi}W(|\Omega|)\ d\phi d|\Omega|$ (7) $\displaystyle=$
$\displaystyle\frac{2\pi(\epsilon
V)^{n}}{2^{n}n!}\int|\Omega|^{2n}\,W(|\Omega|)\,d|\Omega|+\cdots,$

where $W(|\Omega|)$ is the Polyakov loop distribution function that is
independent of the complex phase and is a function of the absolute value
$|\Omega|$. $\phi$ is the complex phase of $\Omega$. $\langle{\rm
Re}\Omega^{n}\rangle$ becomes zero in the limit of $\epsilon\to 0$ for finite
$V$ due to the cancelation of the phase $\phi$. However, after integrating
$\phi$, the complex phase of $\Omega$ have been already removed. Although, in
the double limit, the value of $\epsilon V$ cannot be determined, the leading
term of a ratio $\langle{\rm Re}\Omega^{4}\rangle/\langle{\rm
Re}\Omega^{2}\rangle^{2}$ does not depend on $\epsilon V$, for example. The
volume dependence of such a quantity is expected to be small. The explicit
breaking term is also important for the calculation of the canonical partition
function. As discussed below, this method can be applied to solve the sign
problem in the calculation of the derivative of the canonical partition
function.

##  3 Canonical partition function with a saddle point approximation

When the probability distribution function is maximum, the chemical potential
satisfies $-d\ln W/d\rho=-d\ln Z_{C}/d\rho-V\mu/T=0$, where $V=N_{s}^{3}$ and
$\rho=N/V$. Thus, the chemical potential with the maximum probability at a
density $\rho$ is $\mu/T=-(1/V)d\ln Z_{C}/d\rho$. This derivative is given by
the ratio $(d\ln Z_{C}/d\rho)/Z_{C}$. Both numerator and denominator are
exactly zero for $N\neq 0$ due to the $U(1)$ center symmetry. Thus, $d\ln
Z_{C}/d\rho$ cannot be computed without breaking the center symmetry in $U(1)$
lattice gauge theory. To compute $Z_{C}$, we apply the method we introduced in
the calculation of $\langle{\rm Re}\Omega^{n}\rangle$ in Eq. (7).

For simplicity, we calculate the canonical partition function by a saddle
point approximation proposed in Ref. [6]. The derivative of $Z_{C}$ is
obtained by the following equation,

$\displaystyle-\frac{1}{V}\frac{\partial\ln Z_{\rm C}(T,\rho
V)}{\partial\rho}\approx\frac{\left\langle z_{0}\
\exp\left[V\left(D(z_{0})-\rho
z_{0}\right)\right]e^{-i\alpha/2}\sqrt{\frac{1}{V|D^{\prime\prime}(z_{0})|}}\right\rangle}{\left\langle\exp\left[V\left(D(z_{0})-\rho
z_{0}\right)\right]e^{-i\alpha/2}\sqrt{\frac{1}{V|D^{\prime\prime}(z_{0})|}}\right\rangle}.$
(8)

Here, $D(z)=\ln[\det M(\kappa,\mu)/\det M(\kappa_{0},0)]/V$, and
$D^{\prime\prime}=d^{2}D/dz^{2}=|D^{\prime\prime}|e^{i\alpha}$. $z_{0}$ is the
saddle point which satisfies $D^{\prime}(z_{0})=dD/dz(z_{0})=\rho$. This
bracket $\langle\cdots\rangle$ means the expected value calculated at
$\kappa_{0}$ and $\mu=0$. The saddle point approximation is valid when the
spatial volume is large. The numerical computation of this equation is
possible by the ordinary Monte Carlo method, though this calculation has the
sign problem and the overlap problem.

We calculate the derivative of $Z_{C}$ in $U(1)$ lattice gauge theory for the
case that dynamical fermions are heavy, i.e. $\kappa$ is small, and show that
this calculation is possible avoiding the problem of the center symmetry and
the sign problem. To investigate the region with a small $\kappa$, we evaluate
$D(z)$ by the hopping parameter expansion on each configuration,

$\displaystyle D(z)=96N_{t}N_{\rm f}\kappa^{4}P+2\times 2^{N_{t}}N_{\rm
f}\kappa^{N_{t}}\left[e^{z}\Omega+e^{-z}\Omega^{*}\right]+\cdots,$ (9)
$\displaystyle D^{\prime}(z)=2\times 2^{N_{t}}N_{\rm
f}\kappa^{N_{t}}\left[e^{z}\Omega-e^{-z}\Omega^{*}\right]+\cdots,\hskip
14.22636ptD^{\prime\prime}(z)=2\times 2^{N_{t}}N_{\rm
f}\kappa^{N_{t}}\left[e^{z}\Omega+e^{-z}\Omega^{*}\right]+\cdots,$

and we perform quenched simulations with $\kappa_{0}=0$. The first term
proportional to $P$ can be absorbed by the shift of $\beta$ in the quenched
simulations, and the shift is very small. We thus omit this term. The saddle
point $z_{0}=x_{0}+iy_{0}$ is given by

$\displaystyle 2\times 2^{N_{t}}N_{\rm
f}\kappa^{N_{t}}\left(e^{z_{0}}\Omega-e^{-z_{0}}\Omega^{*}\right)=\rho$ (10)

on each configuration. Because $\rho$ is a real number, the complex phase of
$(e^{x_{0}+iy_{0}}\Omega-e^{-x_{0}-iy_{0}}\Omega^{*})$ is zero. The imaginary
part of saddle point $y_{0}$ is determined as

$\displaystyle
y_{0}=-\arctan\left(\frac{\mathrm{Im}\Omega}{\mathrm{Re}\Omega}\right)=-\mathrm{Arg}\Omega.$
(11)

Then, the real part $x_{0}$ satisfies $2\times 2^{N_{t}}N_{\rm
f}\kappa^{N_{t}}\left[e^{x_{0}}|\Omega|-e^{-x_{0}}\ |\Omega|\right]=\rho$, and
thus

$\displaystyle x_{0}={\rm arcsinh}\left(\frac{\rho}{4\times 2^{N_{t}}N_{\rm
f}\kappa^{N_{t}}|\Omega|}\right).$ (12)

Using the saddle point approximation, the derivative of $Z_{C}$ is obtained by

$\displaystyle-\frac{1}{V}\frac{\partial\mathrm{ln}Z_{C}(T,V\rho)}{\partial\rho}\approx\frac{\left\langle
z_{0}\exp(F+i\theta)\right\rangle}{\left\langle\exp(F+i\theta)\right\rangle},$
(13)

where $F$ and $\theta$ are real numbers and are given as

$\displaystyle F=V(D(z_{0})-\rho
x_{0})-\frac{1}{2}\ln[VD^{\prime\prime}(z_{0})],\hskip
14.22636pt\theta=-V\rho\ \mathrm{Arg}\Omega=-V\rho\
\arctan\left(\frac{\Omega_{I}}{\Omega_{R}}\right),$ (14)

with

$\displaystyle D(z_{0})=96N_{t}N_{\rm f}\kappa^{4}P+4\times 2^{N_{t}}N_{\rm
f}\kappa^{N_{t}}|\Omega|\cosh x_{0},\hskip
14.22636ptD^{\prime\prime}(z_{0})=4\times 2^{N_{t}}N_{\rm
f}\kappa^{N_{t}}|\Omega|\cosh x_{0}.$ (15)

We note that, omitting the plaquette term, $x_{0}$ and $F$ are functions of
$|\Omega|$ and are independent of $\mathrm{Arg}\Omega$, whereas $\theta$ is a
function of $\mathrm{Arg}\Omega$. In simulations of $U(1)$ gauge theory
without dynamical fermions, the probability distribution of $\Omega$ is $U(1)$
symmetric due to the center symmetry. Thus, the numerator and denominator of
Eq. (13) are exactly zero because of the phase factor $e^{i\theta}$, and
Equation (13) cannot be calculated. This is, so to say, the ultimate sign
problem.

Figure 3: Left: Real part of the saddle point $x_{0}$ as a function of the
absolute value of Polyakov loop for $N=10$, 50, 100 and 150. Right: Effective
potential of the Polyakov loop defined in Eq. (17) at $\beta=1.00$.

This problem is solved by adding an additional fermion with a large mass, as
in the case of discussing the expectation value of Polyakov loop. As we
calculated Eq. (7), we consider the $U(1)$ symmetric distribution function of
the Polyakov loop and integrate $\cos(N\phi)$ with respect to the complex
phase $\phi={\rm Arg}\Omega$. Then, $\int\cos(N\phi)\,e^{\epsilon
V|\Omega|\cos\phi}d\phi=2\pi\,(\epsilon V)^{N}|\Omega|^{N}/(2^{N}N!)+\cdots$.
Thus, for a small $\epsilon$, the derivative of $Z_{C}$ becomes

$\displaystyle-\frac{1}{V}\frac{\partial\ln
Z_{C}(T,V\rho)}{\partial\rho}\approx\frac{\int
x_{0}\,e^{F}W(|\Omega|)\frac{2\pi}{2^{N}N!}(\epsilon
V)^{N}|\Omega|^{N}d|\Omega|}{\int
e^{F}W(|\Omega|)\frac{2\pi}{2^{N}N!}(\epsilon V)^{N}|\Omega|^{N}\
d|\Omega|}=\frac{\int x_{0}\,e^{-V_{\rm eff}}\ d|\Omega|}{\int e^{-V_{\rm
eff}}\ d|\Omega|},$ (16)

where the effective potential $V_{\rm eff}(|\Omega|)$ is defined as

$\displaystyle V_{\rm eff}(|\Omega|)=-\ln W(|\Omega|)-F-N\ln|\Omega|.$ (17)

The factor $(\epsilon V)^{N}$ cancels in Eq. (16). Here, a term of $\int
y_{0}\sin(N\phi)d\phi$ is omitted, since this term depend on the upper bound
and lower bound of the integral of $\phi$ and vanishes if we impose an
appropriate dumping factor in the limit of $\phi\to\pm\infty$.

We perform simulations of $U(1)$ lattice gauge theory with the standard Wilson
gauge action at several inverse gauge couplings $\beta=1/g^{2}$ near the
deconfining transition point $\beta_{c}$. The lattice size is $N_{s}^{3}\times
N_{t}=24^{3}\times 6$. Using a pseudo heat bath algorithm, the configurations
are generated at thirteen $\beta$ values in the range from $\beta=1.0000$ to
$1.0240$. The data are taken until 1,000,000 heat-bath sweeps at each $\beta$.
The multipoint reweighting method are used to combine the data generated at
different $\beta$. The phase transition at $\beta\approx 1.0096$ is very weak
first order transition, where two phases coexist. We adopt $N_{\rm f}=2$ and
$\kappa^{6}=1.92\times 10^{-10}$. The real part of saddle point $x_{0}$ for
each $|\Omega|$ and $N=\rho V$ is given by Eq. (12), which is plotted in the
left panel of Fig. 3.

To compute the Polyakov loop distribution function $W(|\Omega|)$, we use a
Gaussian approximation for the delta function:
$\delta(x)\approx\exp[-(x/\Delta)^{2}]/(\Delta\sqrt{\pi})$ with
$\Delta=0.0025$. Also, $F$ of Eq. (14) is calculated from only $|\Omega|$ and
$\rho$. The result of the effective potential $V_{\rm eff}(|\Omega|)$ in Eq.
(16) is plotted in the right panel of Fig. 3 for $N=10,50,100,$ and $150$ when
$\beta=1.00$. When the spatial volume is sufficiently large, Equation (16)
means that $-(1/V)d\ln Z_{C}/d\rho$ is the value of $x_{0}$ at $|\Omega|$
where the effective potential $V_{\rm eff}(|\Omega|)$ is the minimum. The
minimum point increases as $N$ increases. On the other hand, $x_{0}$ increases
as $N$ increases and decreases as $|\Omega|$ increases.

Figure 4: Chemical potential $\mu/T$ for which the number of particles with
the maximum generation probability is $N$ for each $\beta$ in $U(1)$ lattice
gauge theory.

We then calculate the derivative of the canonical partition function. In the
actual calculation, we remove $W(|\Omega|)$ as $W(|\Omega|)d|\Omega|\to{\cal
D}U$ so as not to use the approximate delta function. $-(1/V)d\ln Z_{C}/d\rho$
is $\mu/T$ where the density with the maximum generation probability is
$\rho$, and is computable without the sign problem. We plot the result in Fig.
4 for $\beta=0.996$ – $1.012$. When $\beta$ is in the deconfined phase at zero
density, the chemical potential is monotonically increasing. However, in the
case of the confined phase at zero density, as the density increases, the
chemical potential drops once and increases again. This means that there are
multiple $N$’s for a given $\mu/T$. This is a sign when crossing the first-
order phase transition.

##  4 Conclusions

We studied the probability distribution function of particle density. The
probability distribution function is obtained by constructing the canonical
partition function by fixing the number of particles from the grand partition
function. However, if the system has the center symmetry on a finite lattice,
the canonical partition function is zero when the number of particles is not a
multiple of 3 for $SU(3)$ gauge theory and when the number of particles is not
zero for $U(1)$ gauge theory. Thus, the probability distribution function is
zero in these cases. This situation is natural in the confined phase, but is
unacceptable in the deconfinement phase because there should be states of
various particle numbers.

This problem is essentially the same as the problem that the expectation value
of the Polyakov loop is always zero when calculating with finite volume, as
long as the center symmetry is not broken. To solve this problem, it is
necessary to add an infinitesimal external field to break the symmetry and
take the limit of infinite volume. Moreover, in the case of $U(1)$ gauge
theory, the sign problem can be solved using the $U(1)$ center symmetry at the
same time.

We performed numerical simulations of $U(1)$ lattice gauge theory near the
deconfinement phase transition point. When the dynamical fermions are heavy,
we actually demonstrated that the calculation of the probability distribution
function at finite density is possible using the method proposed in this
study. We calculated the derivative of the canonical partition function using
a saddle point approximation [6], and found that our method to avoid the sign
problem works well. The derivative is equal to $\mu/T$ for which the density
of particles with the maximum generation probability is $\rho$. From the
results of $\rho$ vs. $\mu/T$, changes in the nature of the phase transition
can be investigated. The application of this method to QCD ($SU(3)$ gauge
theory) is a future work.

### Acknowledgments

The author thanks the members of the WHOT-QCD Collaboration for useful
discussions. This work was supported by JSPS KAKENHI Grant Numbers JP21K03550,
JP20H01903, JP19H05146, and the Uchida Energy Science Promotion Foundation.

## References

  * [1] H. Saito, S. Ejiri, S. Aoki, T. Hatsuda, K. Kanaya, Y. Maezawa, H. Ohno, T. Umeda [WHOT-QCD], “Phase structure of finite temperature QCD in the heavy quark region,” Phys. Rev. D 84, 054502 (2011).
  * [2] H. Saito, S. Ejiri, S. Aoki, K. Kanaya, Y. Nakagawa, H. Ohno, K. Okuno and T. Umeda, “Histograms in heavy-quark QCD at finite temperature and density,” Phys. Rev. D 89, 034507 (2014).
  * [3] S. Ejiri, S. Itagaki, R. Iwami, K. Kanaya, M. Kitazawa, A. Kiyohara, M. Shirogane, T. Umeda [WHOT-QCD], “End point of the first-order phase transition of QCD in the heavy quark region by reweighting from quenched QCD,” Phys. Rev. D 101, 054505 (2020).
  * [4] A. Kiyohara, M. Kitazawa, S. Ejiri and K. Kanaya, “Finite-size scaling around the critical point in the heavy quark region of QCD,” arXiv:2108.00118 [hep-lat].
  * [5] N. Wakabayashi, S. Ejiri, K. Kanaya and M. Kitazawa, “Scope and convergence of the hopping parameter expansion in finite temperature QCD with heavy quarks around the critical point,” arXiv:2112.06340 [hep-lat].
  * [6] S. Ejiri, “Canonical partition function and finite density phase transition in lattice QCD,” Phys. Rev. D 78, 074507 (2008).