marginparsep has been altered.
topmargin has been altered.
marginparpush has been altered.
The page layout violates the ICML style. Please do not change the page layout,
or include packages like geometry, savetrees, or fullpage, which change it for
you. We’re not able to reliably undo arbitrary changes to the style. Please
remove the offending package(s), or layout-changing commands and try again.

Neural Basis Functions for Accelerating Solutions
to High Mach Euler Equations

Anonymous Authors1

††footnotetext:  1Anonymous Institution, Anonymous City, Anonymous Region,
Anonymous Country. Correspondence to: Anonymous Author
<EMAIL_ADDRESS>
Preliminary work. Under review by the 2nd AI4Science Workshop at ICML 2022. Do
not distribute.

###### Abstract

We propose an approach to solving partial differential equations (PDEs) using
a set of neural networks which we call Neural Basis Functions (NBF). This NBF
framework is a novel variation of the POD DeepONet operator learning approach
where we regress a set of neural networks onto a reduced order Proper
Orthogonal Decomposition (POD) basis. These networks are then used in
combination with a branch network that ingests the parameters of the
prescribed PDE to compute a reduced order approximation to the PDE. This
approach is applied to the steady state Euler equations for high speed flow
conditions (Mach 10-30) where we consider the 2D flow around a cylinder which
develops a shock condition. We then use the NBF predictions as initial
conditions to a high fidelity Computational Fluid Dynamics (CFD) solver
(CFD++) to show faster convergence. Lessons learned for training and
implementing this algorithm will be presented as well.

##  1 Introduction

The past five years have seen significant interest in the use of neural
networks to understand the behavior of systems governed by differential
equations Han et al. (2018); Willard et al. (2020); Karniadakis et al. (2021).
For example, methods like physics-informed neural networks Raissi et al.
(2019); Lu et al. (2019b) can solve a given instance of a differential
equation or infer its properties given some amount of data and have been
applied in domains like electromagnetism Khan & Lowther (2022) and blood cell
mechanics Yazdani et al. (2021).

More recently still, the subfield of operator-learning has attained prominence
Lu et al. (2019a); Li et al. (2020); Bhattacharya et al. (2021); Mao et al.
(2021). The objective of operator learning is to learn a data-driven
approximation to an operator that defines a differential equation; a trained
operator can then produce a solution to this equation at a greatly reduced
computational cost. Despite the field’s novelty, applications have already
been considered in areas like climate science Kashinath et al. (2021).

Prior to some of these deep network architectures, efforts to reduce the
computational cost often attempted to solve a reduced form of the equations of
interest. Approaches like the Proper Orthogonal Decomposition (POD) Willcox &
Peraire (2002); Gunzburger et al. (2007); Witman et al. (2017); Quarteroni et
al. (2015); Berkooz et al. (1993) use an eigen-decomposition to extract the
primary modes of a given solution space. These modes could then be re-
projected in the reduced space to significantly reduce the degrees of freedom
required to solve the equations.

The approach being proposed here borrows concepts from the more traditional
Reduced Order Modeling (ROM) community and fuses it with the newer deep
network architectures to build a Neural Basis Function (NBF) framework for
solving differential equations.

One of the motivating interests of this work will be the focus on how well
this approach can resolve solutions with complicated phenomena. The specific
phenomenon we will focus on will be solutions with large gradients or
discontinuities, as these problems have posed difficulties for both
traditional techniques and the more recent Neural Network based approaches.
Existing work that considers similar classes of problems includes Mao et al.
(2020) for physics-informed neural networks and Mao et al. (2021) for
operator-learning. We first propose the network architecture and training
methodology that will be used to learn the underlying basis functions and
associated unknown vectors.

In this paper, we introduce the formalism of the neural basis function
(Section 2.1) and describe how to apply it to data sampled from a given PDE
(Section 2.2). In Section 2.3, we demonstrate how to apply the NBF approach to
the steady-state Euler equations. We evaluate our model in Section 3.1 with a
comparison to vanilla DeepONet in Section 3.2 and then show, in Section 3.3,
that predicted solutions from the NBF can be used to accelerate the process of
obtaining high-fidelity data from a CFD solver.

##  2 Methods

###  2.1 The Neural Basis Function Formulation

The neural basis function framework at its core takes a traditional reduced
basis Finite Element Method (FEM) Quarteroni et al. (2015) approach and uses a
set of approximating neural networks to learn the underlying bases and
explicit unknowns of the governing equations. Our basis-approximation networks
are similar to the trunk component of a Deep Operator Net (DeepONet), and our
unknowns-approximation networks are similar to the branch component Lu et al.
(2019a), specifically as found in the recent POD DeepONet extension Lu et al.
(2022).

We consider a subclass of PDEs defined for a vector field
$\mathbf{w}_{\psi}:\Omega\to R^{L}$, where $\Omega\subseteq R^{M}$, that
satisfies the following equations:

$\mathcal{S}_{\psi}\mathbf{w}_{\psi}=\mathbf{0},\,\,\,\,\mathbf{x}\in\Omega$
(1)
$\mathcal{B}_{\psi}\mathbf{w}_{\psi}=\mathbf{0},\,\,\,\,\mathbf{x}\in\partial\Omega,$
(2)

where $\mathcal{S}_{\Psi}$ is a (potentially nonlinear) operator involving
derivatives with respect to space, $\partial\Omega$ is the boundary of
$\Omega$, $\mathcal{B}_{\psi}$ is a boundary condition operator, and
$\psi\in\Psi$ is a set of parameters that fully specifies the PDE and
determines its solution. Section 2.3 illustrates this formulation for a
specific problem. In this paper, we consider steady-state problems with no
time-dependence, but, similar to physics-informed neural networks Raissi et
al. (2019), the NBF approach can be extended to time-dependent problems.

In the operator-learning problem Li et al. (2020); Lu et al. (2019a);
Bhattacharya et al. (2021), we learn an operator $\mathcal{F}$ that, for a
given novel parameter $\psi$, produces an approximate solution
$\hat{\mathbf{w}}=\mathcal{F}(\psi)$ satisfying
$||\hat{\mathbf{w}}-\mathbf{w}_{\psi}||$ under some appropriate norm. This
requires a set of training samples of solutions $\mathbf{w}_{\psi}$ for
varying values of $\psi$.

To define this vector solution $\hat{\mathbf{w}}$ we introduce a linear
combination of reduced basis functions
($\mathbf{\phi}_{j}(\mathbf{x}):\Omega\to R^{M}$) and a set of unknowns
($\mathbf{C}_{j}(\psi):\Psi\to R^{M}$), where $j=1,\ldots,n_{BF}$ indexes the
set of bases used to define $\hat{\mathbf{w}}$. Specifically, Equation 3 shows
how the $i$th component of $\hat{\mathbf{w}}$ is calculated as a linear
combination of the basis functions $\mathbf{\phi}_{j}(\mathbf{x})$ and the
unknowns functions $\mathbf{C}_{j}(\psi)$. The discrete form of this approach
has been used successfully as a way to approximate solutions, within the
parameters space $\psi$, of PDEs at a fraction of the computational cost
Willcox & Peraire (2002); Gunzburger et al. (2007); Witman et al. (2017).

$\hat{w}_{i}(\mathbf{x},\psi)=\sum_{j}^{n_{BF}}C_{ij}(\psi)\phi_{ij}(\mathbf{x})$
(3)

Basis functions are trained by solving the data-driven minimization problem:

$\min_{\phi_{i,j}}\sum_{l}|\phi_{i,j}(\mathbf{x}_{l})-U_{i,j}(\mathbf{x}_{l})|^{2},$
(4)

where the ground truth basis data points $U_{i,j}(\mathbf{x}_{l})$, with
$l=1,...,n$ representing the index of all the spatial data-points. We obtain
the truth basis data via a process described in Section 2.2.

Unknowns functions are trained by solving the following physics-informed
minimization problem:

$\min_{C_{i,j}}\sum_{\psi\sim\mathbb{P}_{\psi}}\sum_{\mathbf{x}\sim\mathbb{P}_{\mathbf{x}}}||\mathcal{S}_{\psi}\hat{\mathbf{w}}(\mathbf{x},\psi)]||^{2}+||\mathcal{B}_{\psi}\hat{\mathbf{w}}(\mathbf{x},\psi)]||^{2},$
(5)

where $\mathbb{P}_{\psi}$ and $\mathbb{P}_{\mathbf{x}}$ are probability
distributions defined for the PDE parameters $\psi$ and the spatial inputs
$\mathbf{x}$, respectively.

Due to the presence of $\mathcal{S}_{\psi}$ in Equation 5, evaluating the loss
function will require taking derivatives of $\hat{\mathbf{w}}$ with respect to
$\mathbf{x}$. This can be done using modern automatic differentiation
techniques Baydin et al. (2017). In practice, we implement our networks in
PyTorch Paszke et al. (2019) and use the functorch He & Zou (2021) library to
calculate loss function derivatives.

###  2.2 Data-driven basis construction

Training the basis functions Equation 4 requires data, which we obtain by
using a high-fidelity simulator like CFD++ to solve instantiations of the PDE.
We sweep a subset of our parameter space $\Psi$ and generate field values
$\mathbf{w}$ for each sampled $\psi\in\Psi$. These samples are divided into
training, validation, and test sets as normal in ML.

For training, we will refer to the total number of training instances as $D$,
and we let $\psi_{d}$ refer to the parameters for the $d$th training instance.
We assume that all of our snapshots
$\\{(\mathbf{w}_{d},\psi_{d}):d=1,\ldots,D\\}$ are evaluated at the same space
points. We index these points as $\\{\mathbf{x}_{l}\\}_{l=1}^{n}$. Thus, the
training set is specified by
$\\{(\\{\mathbf{w}_{d}(\mathbf{x}_{l})\\}_{l=1}^{n},\psi_{d})\\}_{d=1}^{D}$

Once the snapshot data is collected, we follow the Snapshot Proper Orthogonal
Decomposition method Sirovich (1987) and perform singular value decompositions
(SVDs) on the data to extract their bases. For each component $w_{i}$ of
$\mathbf{w}$, we construct a matrix $\mathcal{W}_{i}$ containing the observed
data points, and then we decompose it as
$\mathcal{W}_{i}=U_{i}\Sigma_{i}V_{i}^{T}$.

$U_{i}\Sigma_{i}V_{i}^{T}=\mathcal{W}_{i}=\begin{bmatrix}w_{i}(\mathbf{x}_{1}|\psi_{1})&\ldots&w_{i}(\mathbf{x}_{1}|\psi_{D})\\\
w_{i}(\mathbf{x}_{2}|\psi_{1})&\ldots&w_{i}(\mathbf{x}_{2}|\psi_{D})\\\
\vdots&\ddots&\vdots\\\
w_{i}(\mathbf{x}_{n}|\psi_{1})&\ldots&w_{i}(\mathbf{x}_{n}|\psi_{D})\end{bmatrix}$
(6)

The $U_{i}$ matrix captures the row-space of the snapshot set and represents
the primary space modes of the set of solutions.

Once we have generated the set of decomposed spatial modes describing the
parametric solution space, we can generate a bound on the expected
reconstructed error levels. A well known POD error bound Witman et al. (2017)
can be derived from the singular values in the decomposition via the equation
$\sum_{d=n_{BF}+1}^{D}\sigma_{d}^{2}$ where $D$ here represents the total
number of parameter snapshots available. This error bound allows use to
a-priori choose a suitable error tolerance before initiating the full training
process. For any of the results presented here, we consider the total number
of basis functions as equal to the number of training scenarios.

###  2.3 Euler Equations

Figure 1: The NBF formulation ingests spatial points $\mathbf{x}$ and PDE
parameters $\psi$ into separate sets of networks $\mathbf{\phi}$ and
$\mathbf{C}$, and then the predicted state $\hat{\mathbf{w}}$ is a linear
combination of these terms.

We consider the steady-state compressible Euler equations Certik (2017),
defined for a state vector $\mathbf{w}_{\psi}=[\rho,u,v,E]^{T}$ on a domain
$\Omega\subseteq R^{2}$:

$\mathcal{S}_{\psi}\mathbf{w}_{\psi}:=\partial_{x}\mathbf{F}_{1}(\mathbf{w}_{\psi})+\partial_{y}\mathbf{F}_{2}(\mathbf{w}_{\psi})=0$
(7)

where

$\mathbf{F}_{1}=\begin{bmatrix}\rho u\\\ \rho u^{2}+p\\\ \rho vu\\\
(E+p)u\end{bmatrix},\,\,\,\,\mathbf{F}_{2}=\begin{bmatrix}\rho v\\\ \rho uv\\\
\rho v^{2}+p\\\ (E+p)v\end{bmatrix}.$ (8)

The components of the state vector $\mathbf{w}_{\psi}$ include the fluid
density $w_{0}=\rho$, the $x$ and $y$ velocities $w_{1}=u$ and $w_{2}=v$, and
the total energy $w_{3}=E$.

The Euler equations Equation 7 also include the system pressure $P$, which is
determined by a closure relation in terms of the other state variables:

$\displaystyle P$ $\displaystyle=$
$\displaystyle(\gamma-1)[E-\frac{1}{2}\rho(u^{2}+v^{2})]$ $\displaystyle=$
$\displaystyle(\gamma-1)[w_{3}-\frac{1}{2}w_{0}(w_{1}^{2}+w_{2}^{2})]$

where $\gamma$ is the gas constant and is set to 1.4, using the perfect gas
assumption.

Ultimately, we are interested in using solutions to these equations as initial
conditions for a high fidelity CFD solvers to accelerate convergence. In
additional term the state variable terms we will also need a relation for the
temperature, $T$, which can be expressed as:

$T=\frac{P}{R\rho}=\frac{P}{Rw_{0}}$ (9)

where $R$ represents the ratio of specific heats and for this problem we set
it to a constant $287.058$, also under the perfect gas assumption. The final
term that will be referenced in this example is the speed of the flow and
accounts for the magnitude of the $x$ and $y$ components of velocity
($|\mathbf{u}|=\sqrt{u^{2}+v^{2}}$). With these variables we can write the
equation for the Mach number:

$M=\frac{|\mathbf{u}|}{\sqrt{\gamma RT}}$ (10)

which will be used in the inlet boundary condition described in Section 3.1.
For further expansion of the Euler equations in terms of our state variables
please see Appendix A.

###  2.4 Related work

The operator-learning problem is frequently viewed as learning mappings
between function-spaces Lu et al. (2019a); Li et al. (2020); Bhattacharya et
al. (2021), where the input function might be an initial condition or forcing
function. We consider a more restricted setting, in that we require our
operators $\mathcal{S}_{\phi}$ Equation 1 and $\mathcal{B}_{\psi}$ Equation 2
to be parametric functions of a finite-dimensional vector $\psi$. In both
cases, however, the output is a function Equation 3 defined over an entire
domain.

Of operator-learning approaches, the NBF is most similar to the DeepONet Lu et
al. (2019a). Both use multiple types of networks that ingest different types
of data (called “branch” and “trunk” networks in Lu et al. (2019a)). In
particular, the NBF can be viewed as a extension of the recently-proposed POD
DeepONet Lu et al. (2022). The POD DeepONet retains the full discrete set of
bases Equation 6, but the NBF regresses these bases onto a set of neural
networks $\phi_{i,j}$. This has a few benefits: It 1) can interpolate the
modes to new regions of space when data are only available on sparse or
unstructured grids; 2) is useful for lower memory applications where storing
the potentially large full reduced basis is not feasible; and 3) allows for
the ability to compute spatial and/or temporal derivatives of modeled
quantities like fluid velocity in a batch process via GPU or another
distributed architecture.

The NBF formulation, like the DeepONet formulation, does not impose any
structure to the domain $\Omega$. This is in contrast to some methods like the
Fourier Neural Operator Li et al. (2020). Our use case for the Euler equations
in Section 3.1 relies on this flexibility. See also Lu et al. (2022) for a
discussion of how different operator-learning approaches deal with non-uniform
domains.

##  3 Results

###  3.1 Evaluation of NBF for Euler Equation

For this work, we will make use of a rectangular domain with a quarter circle
cutout located at $x=[0,0]$ with radius 1 (see Figure 3). There is an inflow
condition on the left side of the domain which is controlled by the Mach
number, which for this study is the sole parameter in $\psi$. The domain
extends to $-4$ in the $x$-direction and $7$ in the $y$-direction. This domain
was chosen to allow sufficient room for the shock condition based on our
predefined parametric Mach number range. This 2D cylinder in cross-flow
problem is often used as a benchmark in the CFD space for testing different
algorithms and their ability to capture shockwaves.

We then generated a steady-state dataset of high-fidelity solutions ranging
from Mach 10 to Mach 30 in increments of 1. This data was generated using
CFD++ in a time-dependent configuration and allowed to evolve until the time
derivative portion of Euler equation sufficiently decreased. CFD++ is an
unstructured finite volume flow solver that is maintained by Metacomp
Technologies Inc. It has the capability of solving the full Navier-Stokes
equations with high-order spatial and temporal accuracy. To resolve the
inviscid Euler equations, the solver uses a pseudo time-marching approach
where the initial condition is integrated over time until convergence is
reached. The convergence can be accelerated by using spatially varying time
steps. CFD++ version 20.1 was utilized to perform all the simulations.

Once we have the data generated, the next step is to extract the basis from
the set of snapshots. In general, we hope to be able to approximate solutions
within the parameter space, so, for analysis of this approach, we hold out a
set of Mach cases as our test set. An ablation study of various number of
training data points considered will be presented at the end of this section.
But unless specifically noted, the results presented will use a dataset of 20
Mach numbers with just a single held-out solution. This is done to illustrate
the best case scenarios for how well this approach applies.

Given the training data considered, we extract the state variables of interest
and subtract the mean state from each of the solutions. Empirically, we have
found that this tends to improve the speed of convergence when regressing the
basis networks. The basis networks we consider are general fully-connected
networks with five layers and 40 nodes per layer. We use leaky rectifed linear
units Maas et al. (2013) as activation functions (with a negative slope of
0.01) on the hidden layers, which tends to help resolve the small-scale
features in the basis functions. For example in Figure 2, the shock is very
narrow in space and changes rather rapidly. We use the Adam optimizer Kingma &
Ba (2015) and 250 total training epochs with an initial learning rate of
$10^{-3}$ while applying an exponential learning rate decay factor (0.9) after
70 training epochs.

For PDEs with complicated phenomena (shockwaves), we have found accurately
representing the high-frequency basis functions can be challenging. To resolve
this, we use an importance sampling technique during training that prioritizes
diverse data points in the training routine. We can describe the probability
$p_{i,j,l}$ of selecting a single training data point when minimizing Equation
4, $(\mathbf{x}_{l}$, $U_{i,j}(\mathbf{x}_{l}))$, as:

$p_{i,j,l}=\frac{p^{*}_{i,j,l}}{\sum_{l^{\prime}}p^{*}_{i,j,l^{\prime}}},$
(11)

where $p^{*}_{i,j,l}=|U_{i,j}(\mathbf{x}_{l})-\bar{U}_{i,j}(\mathbf{x}_{l})|$,
and $\bar{U}$ is the mean basis value for all spatial points in the domain
$\mathbf{x}$. In practice, we have found that overusing this sampling approach
leads to poor generalization. However, applying this methodology every so
often (specifically every 4 epochs) improves the resolution of nuanced
features.

Figure 2: To demonstrate the complex spatial features inherent in higher-
frequency basis functions, we visualize the predicted fourth basis for the
density $\rho$ as a function of spatial coordinates $\mathbf{x}$. This is a
zoomed in version around the shock region.

Once the basis networks have been sufficiently trained, the next step is to
learn the set of unknown networks. For unknown network configurations, we use
a seven layer feed forward network with 120 nodes per layer with leaky ReLU
activation functions, similar to the basis networks. In practice, we have
found that Equation 5 is non-trivial to learn with standard weight
initialization techniques. To aid convergence of this loss term, we use a pre-
training approach that attempts to regress the $C_{i,j}$ terms to known
solution data to get the unknown network weights in a sufficient neighborhood
before applying the operator loss in Equation 5. This pre-training
minimization loss can be written as:

$\min_{\tilde{C}_{i,j}}\sum_{\psi}\sum_{l}|C_{i,j}(\mathbf{x}_{l},\psi)-\tilde{C}_{i,j}(\mathbf{x}_{l},\psi)|^{2},$
(12)

where $\tilde{C}_{i,j}$ can be calculated through standard minimization
according to the loss function:

$\min_{\tilde{C}_{i,j}}\sum_{\psi}\sum_{l}\left|\sum_{j}\tilde{C}_{i,j}(\mathbf{x}_{l},\psi)\phi(\mathbf{x}_{l})-\mathbf{w}(\mathbf{x}_{l},\psi)\right|^{2}$
(13)

this allows the unknown network to start off with a reasonable estimate of the
desired states. It is important to note that, during this learning process, we
place no dependencies on the operator portion,
$\mathcal{S}_{\psi}\hat{\mathbf{w}}(\mathbf{x},\psi)$, in the loss function.

After pre-training, the final step is to optimize the unknown network with
respect to the governing equations, boundary conditions and any other
constraints. Minimization of the governing equations and boundary conditions
follows the definition of the Euler equations and the loss function defined in
Equation 5. For this problem we add two additional constraints to minimize the
mean square error between predicted and CFD Pressures and Temperatures.
Unfortunately incorporating this constraint restricts our loss function to
only use $\psi$ parameters that exist in our training set. However, we have
found that without incorporating this constraint leads to negative values of
pressure/temperature which for this problem are non-physical.

Figure 3: The predicted NBF approximation of speed $\sqrt{u^{2}+v^{2}}$ at the
test point, Mach 25. The majority of the domain (solid color) has high values
of speed, then, at the shockwave (bottom right), the speed distribution
rapidly changes. Figure 4: We evaluate our NBF prediction of speed for the
Mach 25 test case with the log absolute error compared to the ground truth.
The errors are highest near the boundary of the shockwave.

Finally with the basis and unknown networks trained, we are able to infer new
solutions within the defined parameter space. We seek to generalize
performance across the parameter space $\Psi$ of Mach numbers. Figure 3 shows
the NBF predicted solution for the speed within the domain while considering a
Mach 25 inflow condition. Figure 4 shows the log absolute error of this
predicted solution with respect to the CFD++ data. It is worth noting that
this particular scenario is held out of our training data set and still shows
good predictive performance.

Figure 5: Relative $L_{2}$ error for variables of interest with respect to
Mach number, for the test data (Mach 25) and training data (remainder).

It is also useful to understand how well the NBF approximation works across
the parameter space. Figure 5 shows the relative $L_{2}$ error as a function
of the Mach number. The errors appear to be relatively stable across the
parameter space with slightly better performance for the higher Mach numbers.

###  3.2 Comparison with Vanilla DeepONet

In order to understand how well the NBF approach compares against existing
approaches, we will benchmark our solution against the Vanilla DeepONet
operator approach. The equation for the vanilla DeepONet is given by Lu et al.
(2022):

$\mathcal{G}(v)(\xi)=\sum_{j=1}^{q}b_{j}(v)t_{j}(\xi)+b_{0}$ (14)

where $v$ and $\xi$ are the parameter (in our case $\psi$) and spatial values
(in our case $\mathbf{x}$). Lu et al. (2022) focuses on the case where $v$ is
a function that parameterizes the PDE – for example, an initial condition. For
comparison purposes, we will consider it a parameterization of $\Psi$.
Additionally, $b_{j}(v)$ and $t_{j}(\xi)$ represent the branch and trunk
networks respectively, with $b_{0}$ acting as a bias imposed on the outputs.

Since we are interested in the two dimensional steady state Euler equations,
we require a DeepONet representation of multiple state variables defined in
Equation 8. Making use of recommendations from Lu et al. (2022), we implement
a set of feed-forward DeepONets for each state variable considered. The
network architectures that we found generated the best performance for our use
case comprised of: 5 layers of 120 nodes with hyperbolic tangent activation
functions for the branch networks and 6 layers of 120 nodes with rectified
linear units for the trunk networks. These networks were then connected in a
shared layer of 64 nodes combined with a bias ($b_{0}$) to produce predictions
for each state variable. During training, we used the same Adam optimizer,
learning rate decay and number of epochs as was used for the NBF approach.

We added a few components that appeared to aid the training process for our
problem specifically: adding prioritized spatial sampling and zero-mean unit-
variance (ZMUV) normalization functions for spatial/parameter inputs as well
as the output state variables. For spatial sampling prioritization we used the
Equation 11 to more often select training examples with more deviation from
the mean solution. Additionally, adding the ZMUV normalization to input/output
parameters appeared to speed convergence and prevent collapsing the solution
to the mean flow.

Figure 6 shows the comparison for the vanilla DeepONet vs the Neural Basis
Function approach considering the total number of training solutions as
inputs. We are able to show better performance in the approximation across the
number of training examples. Specifically we are able to produce significantly
better results for smaller training data sets, indicating better
generalization for sparse data. We believe we are able to show this better
performance due to the pre-conditioned nature of basis functions included in
the NBF estimates as well as application of the governing equation as
regularization during training.

It is worth noting, that although the training process lines up reasonably
closely between the two approaches, there is some overhead incurred by the NBF
approach in: calculating the SVD, then training the unknown network.
Additionally, the DeepONet architecture is significantly smaller than the NBF
architecture as it is described here. We believe there are architecture
optimizations that could reduce the size and will consider this for future
work. Thus, it is not a completely direct comparison but is as close as we
could get to provide a baseline for the NBF approach.

Figure 6: As an ablation study, we vary the amount of training data and
evaluate $L_{2}$ relative error for the test set. The dots represent the
median error with the black tails indicating inner 50th percentile. Blue show
NBF relative error compared to the orange of the vanilla DeepONet

We found that for many of the cases within the ablation study, the Energy
($E$) state variable was the most difficult to get to converge. This is likely
due to the magnitude of values compared to the other state variables ($\rho$,
$u$ and $v$). For DeepONet, applying the ZMUV transformation on the output
state improved the convergence but not significantly so.

###  3.3 Acceleration of CFD solutions

In the conventional CFD approach an initial condition is employed with
iterative updates until steady state is reached. If one was to use an initial
condition that is close to the converged solution, fewer iterations would be
required to reach the high-fidelity steady-state. The NBF solution can be
leveraged to accelerate the CFD solver by providing this nearly converged
initial condition. In the following section, we will compare the performance
of the conventional CFD approach with NBF informed CFD approach. The two
approaches will be evaluated by comparing their average residual history. The
time-dependent form of the Euler equation includes the the term:

$\frac{\partial\mathbf{w}}{\partial t}$ (15)

and adds to Equation 8. Additional details on the setup and implementation of
the CFD++ configuration can be found in Appendix C

The convergence is determined by the evolution of the time-dependent residual.
Once the residual levels off, the simulation is deemed to be converged and the
solution has reached steady-state. In the CFD++ finite volume approach, the
governing equations are integrated over the cell volumes. For the case of the
Euler equations, the right hand side becomes the sum of the scalar product of
the inviscid fluxes and the cell face areas. For a given equation, the
residual is defined as the absolute value of the residual divided by the cell
volume averaged over all the cells. Since we are working with a system of
equations, we monitor the average residual of our system equations, i.e. the
sum of the residual of each equation divided by the number of equations.
Additional detail of all the numerical procedures can be found in the CFD++
user-manual111https://www.metacomptech.com/index.php/features/icfd.

Figure 7: Steady state convergence acceleration of a Mach 15 scenario. Using
the NBF prediction as the initialization results in a speedup of roughly 750
iterations compared to the freestream initialization. Figure 8: Steady state
convergence acceleration of a Mach 25 scenario. Using the NBF prediction as
the initialization results in a speedup of roughly 1000 iterations, compared
to the freestream initialization.

Figure 7 and Figure 8 show the time-dependent residual as a function of the
number of iterations for Mach 15 and 25 respectively. As we can see from these
plots, the NBF solution provides a noticeable speed up in the convergence of
the high fidelity solver.

##  4 Discussion

This work has shown the potential for using an operator learning based
approach, like the Neural Basis Function, as an acceleration mechanism for
high-fidelity CFD solvers. We found that regressing a set of networks to POD
basis data, then using an additional unknown network in linear combination can
provide sufficient estimates of state variables within a configurable
parameter space. High Mach flow cases often induce complicated phenomena, like
discontinuous shockwaves, which are difficult to represent numerically, even
with standard techniques. Using the 2D steady state Euler equations as an
example, we were able to produce accurate estimates to the flow field for high
Mach (10-30) scenarios. These solutions were then used as initial conditions
to a high-fidelity CFD solver (CFD++) to improve the convergence, resulting in
speed ups up to 1000 iterations.

As future work on the application side, we intend to extend these results to
the full Navier-Stokes equations, while transitioning to a three-dimensional
domain. Along the way improvements will be made to the methodology to better
understand network architecture performance as well as improving training
techniques for the basis and unknown networks. We also hope to develop a
methodology using NBF (or other operator approaches) to handle problems where
the spatial domain might vary across snapshots. This would allow for the
ability to do shape optimization or analysis of multiple geometric
configurations.

On the learning side, we hope to factor in an active learning feedback loop
Settles (2009) that will be able to identify regions of the $\psi$ parameter
space with poor predictive performance, in order to execute more high-fidelity
scenarios in these regions to improve overall accuracy. Finally, some
combination of DeepOnet and NBF seems like a reasonable path forward
especially for speeding convergence of these techniques. For example, one
could create a set of trunk networks that are regressed from the orthogonal
SVD basis, then freeze their weights and biases. Then adding a vanilla
DeepONet in combination with these bases could greatly speed the convergence
and make it easier to implement an active learning type of approach for data
generation.

## References

  * Baydin et al. (2017) Baydin, A. G., Pearlmutter, B. A., Radul, A. A., and Siskind, J. M.  Automatic differentiation in machine learning: A survey.  _J. Mach. Learn. Res._ , 18(1):5595–5637, jan 2017.  ISSN 1532-4435.
  * Berkooz et al. (1993) Berkooz, G., Holmes, P., and Lumley, J. L.  The proper orthogonal decomposition in the analysis of turbulent flows.  _Annual review of fluid mechanics_ , 25(1):539–575, 1993.
  * Bhattacharya et al. (2021) Bhattacharya, K., Hosseini, B., Kovachki, N. B., and Stuart, A. M.  Model Reduction And Neural Networks For Parametric PDEs.  _The SMAI journal of computational mathematics_ , 7:121–157, 2021.  doi: 10.5802/smai-jcm.74.
  * Certik (2017) Certik, O.  Theoretical physics reference.  2017\.
  * Gunzburger et al. (2007) Gunzburger, M. D., Peterson, J. S., and Shadid, J. N.  Reduced-order modeling of time-dependent pdes with multiple parameters in the boundary data.  _Computer methods in applied mechanics and engineering_ , 196(4-6):1030–1047, 2007.
  * Han et al. (2018) Han, J., Jentzen, A., and E, W.  Solving high-dimensional partial differential equations using deep learning.  _Proceedings of the National Academy of Sciences_ , 115(34):8505–8510, 2018.  ISSN 0027-8424.  doi: 10.1073/pnas.1718942115.
  * He & Zou (2021) He, H. and Zou, R.  functorch: Jax-like composable function transforms for pytorch.  https://github.com/pytorch/functorch, 2021.
  * Karniadakis et al. (2021) Karniadakis, G. E., Kevrekidis, I. G., Lu, L., Perdikaris, P., Wang, S., and Yang, L.  Physics-informed machine learning.  _Nature Reviews Physics_ , 3(6):422–440, 2021\.
  * Kashinath et al. (2021) Kashinath, K., Mustafa, M., Albert, A., Wu, J.-L., Jiang, C., et al.  Physics-informed machine learning: case studies for weather and climate modelling.  _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 379(2194):20200093, 2021.  doi: 10.1098/rsta.2020.0093.
  * Khan & Lowther (2022) Khan, A. and Lowther, D. A.  Physics informed neural networks for electromagnetic analysis.  _IEEE Transactions on Magnetics_ , pp. 1–1, 2022.  doi: 10.1109/TMAG.2022.3161814.
  * Kingma & Ba (2015) Kingma, D. P. and Ba, J.  Adam: A method for stochastic optimization.  In Bengio, Y. and LeCun, Y. (eds.), _3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings_ , 2015.
  * Li et al. (2020) Li, Z., Kovachki, N., Azizzadenesheli, K., Liu, B., Bhattacharya, K., Stuart, A., and Anandkumar, A.  Fourier neural operator for parametric partial differential equations.  _arXiv preprint arXiv:2010.08895_ , 2020.
  * Lu et al. (2019a) Lu, L., Jin, P., and Karniadakis, G. E.  Deeponet: Learning nonlinear operators for identifying differential equations based on the universal approximation theorem of operators.  _arXiv preprint arXiv:1910.03193_ , 2019a.
  * Lu et al. (2019b) Lu, L., Meng, X., Mao, Z., and Karniadakis, G. E.  Deepxde: A deep learning library for solving differential equations.  _CoRR_ , abs/1907.04502, 2019b.
  * Lu et al. (2022) Lu, L., Meng, X., Cai, S., Mao, Z., Goswami, S., Zhang, Z., and Karniadakis, G. E.  A comprehensive and fair comparison of two neural operators (with practical extensions) based on fair data.  _Computer Methods in Applied Mechanics and Engineering_ , 393:114778, 2022.
  * Maas et al. (2013) Maas, A. L., Hannun, A. Y., and Ng, A. Y.  Rectifier nonlinearities improve neural network acoustic models.  In _in ICML Workshop on Deep Learning for Audio, Speech and Language Processing_ , 2013.
  * Mao et al. (2020) Mao, Z., Jagtap, A. D., and Karniadakis, G. E.  Physics-informed neural networks for high-speed flows.  _Computer Methods in Applied Mechanics and Engineering_ , 360:112789, 2020.  ISSN 0045-7825.  doi: https://doi.org/10.1016/j.cma.2019.112789.
  * Mao et al. (2021) Mao, Z., Lu, L., Marxen, O., Zaki, T. A., and Karniadakis, G. E.  Deepm&mnet for hypersonics: Predicting the coupled flow and finite-rate chemistry behind a normal shock using neural-network approximation of operators.  _Journal of Computational Physics_ , 447:110698, 2021.  ISSN 0021-9991.  doi: https://doi.org/10.1016/j.jcp.2021.110698.
  * Paszke et al. (2019) Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., et al.  Pytorch: An imperative style, high-performance deep learning library.  In Wallach, H., Larochelle, H., Beygelzimer, A., d'Alché-Buc, F., Fox, E., and Garnett, R. (eds.), _Advances in Neural Information Processing Systems 32_ , pp. 8024–8035. Curran Associates, Inc., 2019.
  * Quarteroni et al. (2015) Quarteroni, A., Manzoni, A., and Negri, F.  _Reduced basis methods for partial differential equations: an introduction_ , volume 92.  Springer, 2015.
  * Raissi et al. (2019) Raissi, M., Perdikaris, P., and Karniadakis, G. E.  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\.
  * Settles (2009) Settles, B.  Active learning literature survey.  2009\.
  * Sirovich (1987) Sirovich, L.  Turbulence and the dynamics of coherent structures. i. coherent structures.  _Quarterly of applied mathematics_ , 45(3):561–571, 1987.
  * Willard et al. (2020) Willard, J., Jia, X., Xu, S., Steinbach, M., and Kumar, V.  Integrating scientific knowledge with machine learning for engineering and environmental systems, 2020.  URL https://arxiv.org/abs/2003.04919.
  * Willcox & Peraire (2002) Willcox, K. and Peraire, J.  Balanced model reduction via the proper orthogonal decomposition.  _AIAA journal_ , 40(11):2323–2330, 2002.
  * Witman et al. (2017) Witman, D. R., Gunzburger, M., and Peterson, J.  Reduced-order modeling for nonlocal diffusion problems.  _International Journal for Numerical Methods in Fluids_ , 83(3):307–327, 2017.
  * Yazdani et al. (2021) Yazdani, A., Deng, Y., Li, H., Javadi, E., Li, Z., Jamali, S., Lin, C., Humphrey, J. D., Mantzoros, C. S., and Em Karniadakis, G.  Integrating blood cell mechanics, platelet adhesive dynamics and coagulation cascade for modelling thrombus formation in normal and diabetic blood.  _Journal of The Royal Society Interface_ , 18(175):20200834, 2021.  doi: 10.1098/rsif.2020.0834.

##  Appendix A Euler Equation Expansion

We can expand the Euler equations in terms of our state variables
$\mathbf{w}$:

$\displaystyle\mathbf{F}_{1}$ $\displaystyle=$
$\displaystyle[f_{11},f_{12},f_{13},f_{14}]^{T}$ $\displaystyle f_{11}$
$\displaystyle=$ $\displaystyle\rho u$ $\displaystyle=$ $\displaystyle
w_{0}w_{1}$ $\displaystyle f_{12}$ $\displaystyle=$ $\displaystyle\rho
u^{2}+p$ $\displaystyle=$ $\displaystyle
w_{0}w_{1}^{2}+(\gamma-1)[w_{3}-\frac{1}{2}w_{0}(w_{1}^{2}+w_{2}^{2})]$
$\displaystyle f_{13}$ $\displaystyle=$ $\displaystyle\rho uv$
$\displaystyle=$ $\displaystyle w_{0}w_{1}w_{2}$ $\displaystyle f_{14}$
$\displaystyle=$ $\displaystyle(E+p)u$ $\displaystyle=$
$\displaystyle(w_{3}+(\gamma-1)[w_{3}-\frac{1}{2}w_{0}(w_{1}^{2}+w_{2}^{2})])w_{1}$
$\displaystyle\mathbf{F}_{2}$ $\displaystyle=$
$\displaystyle[f_{21},f_{22},f_{23},f_{24}]^{T}$ $\displaystyle f_{21}$
$\displaystyle=$ $\displaystyle\rho v$ $\displaystyle=$ $\displaystyle
w_{0}w_{2}$ $\displaystyle f_{22}$ $\displaystyle=$ $\displaystyle\rho uv$
$\displaystyle=$ $\displaystyle w_{0}w_{1}w_{2}$ $\displaystyle f_{23}$
$\displaystyle=$ $\displaystyle\rho v^{2}+p$ $\displaystyle=$ $\displaystyle
w_{0}w_{2}^{2}+(\gamma-1)[w_{3}-\frac{1}{2}w_{0}(w_{1}^{2}+w_{2}^{2})]$
$\displaystyle f_{24}$ $\displaystyle=$ $\displaystyle(E+p)v$ $\displaystyle=$
$\displaystyle(w_{3}+(\gamma-1)[w_{3}-\frac{1}{2}w_{0}(w_{1}^{2}+w_{2}^{2})])w_{2}$

##  Appendix B Other Variable Field Visuals

The following set of plots provide a visual for the other variables of
interest using the NBF approximation to the Mach 25 scenario.

Figure 9: Density ($\rho$) Figure 10: $x$-velocity component ($u$) Figure 11:
$y$-velocity component ($v$) Figure 12: Energy ($E$) Figure 13: Pressure ($P$)
Figure 14: Temperature ($T$)

##  Appendix C CFD++ Details

In our investigation, CFD++ was run with a second order spatial discretization
scheme. The inviscid fluxes were calculated based on the minmod TVD limiter.
The fluxes on the cell faces were reconstructed based on centroidal
polynomials. Spatial scheme blending was leveraged for the sake of numerical
stability. During the initial 250 iterations, first order fluxes were used
exclusively. Then a linear blending of the first and second order
discretization scheme was used until iteration 750, beyond which the second
order fluxes were used exclusively.

The time-residual from iteration to iteration depends on the integration
scheme. In this study, an implicit Gauss-Seidel relaxation time integration
scheme was utilized with no additional convergence acceleration methods. The
time step sizes were based on the local (spatially varying) CFL (Courant
Friedrichs-Lewy) number. The CFL number is a dimensionless number that is
equal to the time step size scaled by the time it takes for the fasts moving
characteristics to travel across a cell
$CFL=\frac{dt}{\frac{dx}{|\mathbf{u}|+a}}$, recalling $|\mathbf{u}|$ is the
scaler speed field. $a$ is the speed of sound ($a=\sqrt{T*R*\gamma}$) In each
simulations, the CFL number was linearly increased from 0.01 to 10 over the
first 1000 iterations beyond which the CFL number was kept constant.


# Numerical Radius of Products of Special Matrices

Mohammad Alakhrass  Department of Mathematics, University of Sharjah , Sharjah
27272, UAE<EMAIL_ADDRESS>
###### Abstract.

The purpose of this note is to present upper bounds estimations for the
numerical radius of a products and Hadamard products of special matrices,
including sectorial and accretive-dissipative matrices.

###### Key words and phrases:

Numerical radius, Sectorial matrices, accretive-dissipative matrices, Matrix
inequalities

###### 2010 Mathematics Subject Classification:

15A45, 15A60, 47A12, 47A30.

##  1\. Introduction

Let $\mathbb{M}_{n}$ be the algebra of all $n\times n$ complex matrices. If
$X=[x_{i,j}],Y=[y_{i,j}]\in\mathbb{M}_{n}$, then their Hadamard product
$X\circ Y$ is the matrix $[x_{i,j}y_{i,j}].$ The cartesian decomposition of
$X\in\mathbb{M}_{n}$ is presented as

$X=A+iB,$ (1)

where $A$ and $B$ are the Hermitian matrices $A=Re(X)=\frac{X+X^{*}}{2}$ and
$B=Im(X)=\frac{X-X^{*}}{2i}$. A matrix $X$ is said to be accretive (resp.
dissipative) if in its cartesian decomposition (1) the matrix $A$ (resp. $B$)
is positive definite. If both $A$ and $B$, in the decomposition (1), are
positive definite, $X$ is called accretive-dissipative.

The numerical range of $X\in\mathbb{M}_{n}$ is the compact convex subset of
the complex plane defined as follows:

$W(X)=\\{\langle{Xx,x}\rangle:x\in\mathbb{C}^{n},||x||=1\\},$

where $\langle{\cdot,\cdot}\rangle$ is the standard inner product on
$\mathbb{C}^{n}$ and $||\cdot||$ is the Euclidean norm on $\mathbb{C}^{n}$. A
very important result is that

$\sigma(X)\subset W(X),$

where $\sigma(X)$ is the spectrum of $X$.

For $\alpha\in[0,\pi/2)$, let $S_{\alpha}$ be the sector defined in the
complex plane by

$S_{\alpha}=\\{z\in\mathbb{C}:Re(z)>0,|Im(z)|\leq\tan(\alpha)Re(z)\\}.$

A matrix $X$ is called sectorial if $W(zX)\subset S_{\alpha}$ for some complex
number $z$ with $|z|=1$. The smallest possible such $\alpha$ is called the
index of sectoriality.

For $\alpha\in[0,\pi/2)$, let $\mathbb{M}^{s}_{n,\alpha}$ be the class of all
$n\times n$ matrices $X$ with $W(zX)\subset S_{\alpha}$ for some complex
number $z$ with $|z|=1$.

It is clear that $X$ is accretive-dissipative if and only if $W(e^{-\pi
i/4}X)\subset S_{\pi/4}$, and hence $X\in\mathbb{M}^{s}_{n,\alpha}$ with
$\alpha=\frac{\pi}{4}$. For more study of sectorial matrices see [1, 2, 3, 4,
9, 10, 11, 12] and the references therein.

A norm $N$ on $\mathbb{M}_{n}$ is said to be unitarily invariant if it
satisfies the property $N(UXV)$ for all $X\in\mathbb{M}_{n}$ and all unitaries
$U,V\in\mathbb{M}_{n},$ and it is said to be multiplicative if $N(AB)\leq
N(A)N(B)$ for all $A,B\in\mathbb{M}_{n}$. Examples of such unitarily invariant
multiplicative norms are the Schatten $p$-norm defined by
$||X||_{p}=\left(\sum_{j=1}^{n}s^{p}_{j}(X)\right)^{p},\;p\geq 1.$ When
$p=\infty$, this last norm is just the usual operator norm defined by
$\|X\|=\sup_{\|x\|=1}\|Xx\|.$

Associated with numerical range, the numerical radius of $X$ is defined by

$\omega(X)=\sup\\{|z|:z\in W(X)\\}.$

It is well known that $w(\cdot)$ defines a norm on $\mathbb{M}_{n}$ which is
equivalent to the usual operator norm $\|\cdot\|$. In fact we have

$\frac{1}{2}\|X\|\leq\omega(X)\leq\|X\|;\quad\forall X\in\mathbb{M}_{n}.$ (2)

Moreover, if $X\in\mathbb{M}_{n}$ is normal then $\omega(X)=\|X\|.$ Therefore,
the inequalities in (2) are sharp.

Obviously, $\omega(\cdot)$ defines a weakly unitarily invariant norm on
$\mathbb{M}_{n}$; that is it satisfies the property $\omega(UXU)=\omega(X)$
for all $X\in\mathbb{M}_{n}$ and all unitary $U\in\mathbb{M}_{n}.$

######  Example 1.1.

Let $X=\left(\begin{array}[]{cc}0&2\\\ 0&0\\\ \end{array}\right)$ and
$Y=\left(\begin{array}[]{cc}0&0\\\ 2&0\\\ \end{array}\right).$ Then
$\omega(X)=\omega(Y)=1$ and $\omega(XY)=4$

Example 1.1 shows that $\omega(\cdot)$ is not a multiplicative norm. However,
the inequalities in (2) implies that $4\omega(\cdot)$ is a multiplicative
norm. That is for all $X,Y\in\mathbb{M}_{n}$,
$4\omega(XY)\leq\left(4\omega(X)\right)\left(4\omega(Y)\right)$. Equivalently,

$\omega(XY)\leq 4\omega(X)\omega(Y);\quad\forall X,Y\in\mathbb{M}_{n}.$ (3)

Obviously, Example 1.1 shows that the inequality (3) is sharp and the constant
4 is the best possible in (3).

The Hadamard product version of (3) can be written as

$\omega(X\circ Y)\leq 2\omega(X)\omega(Y).$ (4)

That is, $2\omega(\cdot)$ is a multiplicative norm over the Hadamard product.
See [6, p. 73].

The following example shows that the constant 2 is the best possible in (4).

######  Example 1.2.

Let $X=Y=\left(\begin{array}[]{cc}0&2\\\ 0&0\\\ \end{array}\right)$ Then
$\omega(X)=\omega(Y)=1$ and $\omega(X\circ Y)=2$.

By considering special matrices $X$ and $Y$, it is possible to obtain better
estimations than those in (3) and (4).

If $XY=YX$, then

$\omega(XY)\leq 2\omega(X)\omega(Y).$ (5)

See [7, Theorem 2.5-2].

If $X$ or $Y$ is normal such that $XY=YX$, then

$\text{}\quad\omega(XY)\leq\omega(X)\omega(Y).$ (6)

See [7, Corollary 2.5-6].

If $X$ or $Y$ is normal. Then

$\omega(X\circ Y)\leq\omega(X)\omega(Y).$ (7)

See [6, Corollary 4.2.17].

If $A,B\in\mathbb{M}_{n}$ and $A=[a_{ij}]$ is positive semidefinite, then

$\omega(A\circ B)\leq\left(\max_{j}a_{jj}\right)\,\,\omega(B).$ (8)

See [5, Corollary 4] and [8, Proposition 4.1].

The purpose of this short note is to add more inequalities to the above list.
More precisely, we give estimations of the numerical radius of products or
Hadamard products of sectorial matrices and related matrices such as accretive
and dissipative matrices.

##  2\. Main Results

We start this section by the following two observations.

######  Lemma 2.1.

Let $X\in\mathbb{M}_{n}$. If $W(X)\subset S_{\alpha}$, then

$\|X\|\leq\sec(\alpha)\|Re(X)\|.$

The above Lemma can be found in [1], [3] , [2] and [12].

######  Remark 2.1.

  1. (1)

We recall that $\omega(\cdot)$ defines a self-adjoint norm on
$\mathbb{M}_{n}$, that is it is a norm satisfies the properties
$w(X^{*})=w(X)$ for all $X\in\mathbb{M}_{n}$. Therefore,

$\displaystyle\omega\left(Re(X)\right)$
$\displaystyle=\omega\left(\frac{X+X^{*}}{2}\right)$
$\displaystyle\leq\frac{1}{2}\left(\omega\left(X+X^{*}\right)\right)$
$\displaystyle\leq\frac{1}{2}\left(\omega\left(X\right)+\omega\left(X^{*}\right)\right)$
$\displaystyle\leq\frac{1}{2}\left(\omega\left(X\right)+\omega\left(X\right)\right)$
$\displaystyle\leq\omega\left(X\right).$

  2. (2)

Let $X\in\mathbb{M}_{n}$ and let $X=A+iB$ be its cartesian decomposition. If
$W(X)\subset S_{\alpha}$, then for any $v\in\mathbb{C}^{n}$ with $||v||=1$
$\langle Av\;,\;v\rangle,\langle Bx\;,\;x\rangle\in\mathbb{R}$ and

$\langle Xv\;,\;v\rangle=\langle Av\;,\;v\rangle+i\;\langle Bv\;,\;v\rangle\in
S_{\alpha}.$

Therefore,

$|\langle Bv\;,\;v\rangle|\leq\tan(\alpha)\langle Av\;,\;v\rangle.$

Taking the supremum over all such v’s gives

$\omega(B)\leq\tan(\alpha)\omega(A).$

Now we are ready to state the first result.

######  Theorem 2.1.

Let $X\in\mathbb{M}^{s}_{n,\alpha_{1}}$ and
$Y\in\mathbb{M}^{s}_{n,\alpha_{2}}$, where
$\alpha_{2},\alpha_{2}\in[0,\pi/2)$. Then

$\omega(XY)\leq\sec(\alpha_{1})\sec(\alpha_{2})\omega(X)\omega(Y).$

###### Proof.

Since $X\in\mathbb{M}^{s}_{n,\alpha_{1}}$ and
$Y\in\mathbb{M}^{s}_{n,\alpha_{2}}$, there are two complex numbers
$z,w\in\mathbb{C}$ with $|z|=|w|=1$ such that $W(zX)\subset S_{\alpha_{1}}$
and $W(wY)\subset S_{\alpha_{2}}$.

Notice that

$\displaystyle\omega(XY)$ $\displaystyle\leq||XY||\quad\text{( by \eqref{A}
))}$ $\displaystyle\leq||X||\,\,||Y||$ $\displaystyle=||zX||\,\,||wY||$
$\displaystyle\leq
sec(\alpha_{1})\sec(\alpha_{2})||Re\left(zX\right)|\,\,|||Re\left(wY\right)||\quad\text{(by
Lemma \ref{L1})}$
$\displaystyle=sec(\alpha_{1})\sec(\alpha_{2})\omega(Re\left(zX\right))\,\,\omega\left(Re\left(wY\right)\right)$
(since $Re\left(zX\right))$ and $Re\left(wY\right)$ are Hermitian )
$\displaystyle\leq
sec(\alpha_{1})\sec(\alpha_{2})\omega\left(zX\right)\,\,\omega(wY)\quad\text{(by
Remark \ref{R1}) }$
$\displaystyle=sec(\alpha_{1})\sec(\alpha_{2})\omega\left(X\right)\,\,\omega(Y).$

∎

We remark that if $X,Y\in\mathbb{M}^{s}_{n,\alpha}$, Theorem 2.1 implies that

$\omega(XY)\leq\sec^{2}(\alpha)\omega(X)\omega(Y).$ (9)

The inequality (9) presents a refinement of the inequality (3) when
$0\leq\alpha\leq\frac{\pi}{3}$. A particular case is when $X$ and $Y$ are
accretive-dissipative as in the following result.

######  Corollary 2.1.

If $X,Y\in\mathbb{M}_{n}$ are accretive-dissipative, then

$\omega(XY)\leq 2\omega(X)\omega(Y).$

###### Proof.

The result follows from Theorem 2.1 and the fact that if
$X,Y\in\mathbb{M}_{n}$ are accretive-dissipative then
$X,Y\in\mathbb{M}^{s}_{n,\alpha}$ with $\alpha=\pi/4.$ ∎

######  Remark 2.2.

  1. (1)

If $X\in\mathbb{M}_{n}$ is accretive, i.e. $Re(X)>0$, then $X$ is sectorial
with sectorial index

$\alpha_{X}=\tan^{-1}(|\lambda_{1}\left(\left(ReX\right)^{-1}Im(X))|\right)=||\left(ReX\right)^{-1/2}\left(ImX\right)\left(ReX\right)^{-1/2}||.$
(10)

See [4]. Therefor $X\in\mathbb{M}^{s}_{n,\alpha}$ with $\alpha=\alpha_{X}$.

  2. (2)

If $X\in\mathbb{M}_{n}$ is dissipative, i.e. $Im(X)>0$, then $iX$ is accretive
and hence $X\in\mathbb{M}^{s}_{n,\alpha}$ with $\alpha=\alpha_{X}$.

######  Corollary 2.2.

If $X,Y\in\mathbb{M}_{n}$ are accretive (or dissipative), then

$\omega(XY)\leq(1+a^{2})\omega(X)\omega(Y),$

where

$a=\max\\{|\lambda_{1}\left(\left(ReX\right)^{-1}ImX\right)|,|\lambda_{1}\left(\left(ReY\right)^{-1}ImY\right)|\\}.$

###### Proof.

Since $X,Y\in\mathbb{M}_{n}$ are accretive (or dissipative), we have
$X\in\mathbb{M}^{s}_{n,\alpha_{X}}$ and $Y\in\mathbb{M}^{s}_{n,\alpha_{Y}}$,
where $\alpha_{X}$ and $\alpha_{Y}$ are as in (10). By Theorem 2.1, we have

$\omega(XY)\leq\sec(\alpha_{X})\sec(\alpha_{Y})\omega(X)\omega(Y).$

Now the result follows by Noting that

$\sec(\alpha_{X})=\sec(\tan^{-1}(|\lambda_{1}\left(\left(Re(X)\right)^{-1}Im(X))|\right))=\sqrt{1+\left(\lambda_{1}\left(\left(ReX\right)^{-1}ImX\right)\right)^{2}},$

and

$\sec(\alpha_{Y})=\sec(\tan^{-1}(|\lambda_{1}\left(\left(ReY\right)^{-1}ImY)|\right))=\sqrt{1+\left(\lambda_{1}\left(\left(ReY\right)^{-1}ImY\right)\right)^{2}}.$

∎

Same proof method used to prove Theorem 2.1 can be used to prove the following
more general result.

######  Theorem 2.2.

Let $X_{j}\in\mathbb{M}^{s}_{n,\alpha_{j}},j=1,2,3,...,m$. Then

$\omega\left(\prod_{j=1}^{m}X_{j}\right)\leq\prod_{j=1}^{m}\sec(\alpha_{j})\omega(X_{j}).$

Consequently,

######  Corollary 2.3.

If $X_{1},X_{2},...,X_{m}\in\mathbb{M}_{n}$ are accretive-dissipative, then

$\omega\left(\prod_{j=1}^{m}X_{j}\right)\leq
2^{m/2}\prod_{j=1}^{m}\omega(X_{j}).$

######  Corollary 2.4.

If $X_{1},X_{2},...,X_{m}\in\mathbb{M}_{n}$ are accretive (or dissipative),
then

$\omega\left(\prod_{j=1}^{m}X_{j}\right)\leq(1+a^{2})^{m/2}\prod_{j=1}^{m}\omega(X_{j}),$

where

$a=\max\\{|\lambda_{1}\left(\left(ReX_{j}\right)^{-1}ImX_{j}\right)|,j=1,2,...,m\\}.$

In what follows, we present upper bounds for the numerical ranges of hadamard
products of special matrices. The following Lemma is important in our
analysis. It can be found in [1, 2, 3].

######  Lemma 2.2.

Let $T\in\mathbb{M}_{n}$. If $W(T)\subset S_{\alpha}$ for some
$\alpha\in[0,\pi/2)$. Then

$\left(\begin{array}[]{cc}\sec(\alpha)\Re(T)&T\\\ T^{*}&\sec(\alpha)\Re(T)\\\
\end{array}\right)\geq 0.$

Now we estimate the numerical range for a Hadamard product of sectorial
matrices.

######  Theorem 2.3.

Let $X\in\mathbb{M}^{s}_{n,\alpha_{1}}$ and
$Y\in\mathbb{M}^{s}_{n,\alpha_{2}}$. Then

$\omega(X\circ Y)\leq\sec(\alpha_{1})\sec(\alpha_{1})\omega(X)\omega(Y).$

###### Proof.

Since $X\in\mathbb{M}^{s}_{n,\alpha_{1}}$ and
$Y\in\mathbb{M}^{s}_{n,\alpha_{2}}$, there are two complex numbers
$z,w\in\mathbb{C}$ with $|z|=|y|=1$ such that $W(zX)\subset S_{\alpha_{1}}$
and $W(wY)\subset S_{\alpha_{2}}$. Therefore, by Lemma 2.2, the following two
block matrices

$\left(\begin{array}[]{cc}\sec(\alpha_{1})Re(zX)&zX\\\
\overline{z}X^{*}&\sec(\alpha_{1})Re(zX)\\\
\end{array}\right),\,\,\,\left(\begin{array}[]{cc}\sec(\alpha_{2})Re(wY)&wY\\\
\overline{w}Y^{*}&\sec(\alpha_{2})Re(wY)\\\ \end{array}\right)$

are positive semidefinite. Hence

$\left(\begin{array}[]{cc}\sec(\alpha_{1})\sec(\alpha_{2})Re(zX)\circ
Re(wY)&zw\left(X\circ Y\right)\\\ \overline{zw}\left(X\circ
Y\right)^{*}&\sec(\alpha_{1})\sec(\alpha_{2})Re(zX)\circ Re(wY)\\\
\end{array}\right)$

is also positive semidefinite. Since $||\cdot||$ is a Lieb Function, we have

$||X\circ Y||=||zw\left(X\circ
Y\right)||\leq\sec(\alpha_{1})\sec(\alpha_{2})||Re(zX)\circ Re(wY)||.$

Now, observe that

$\displaystyle\omega(X\circ Y)$ $\displaystyle\leq||X\circ Y||$
$\displaystyle\leq\sec(\alpha_{1})\sec(\alpha_{2})||Re(zX)\circ Re(wY)||$
$\displaystyle=\sec(\alpha_{1})\sec(\alpha_{2})\omega(Re(zX)\circ Re(wY))$
$\displaystyle\leq\sec(\alpha_{1})\sec(\alpha_{2})\omega(Re(zX))\omega(Re(wY))\quad(\text{by
\eqref{H}})$ (11)
$\displaystyle\leq\sec(\alpha_{1})\sec(\alpha_{2})\omega(zX)\omega(wY)\quad(\text{by
Remark \ref{R1}})$
$\displaystyle=\sec(\alpha_{1})\sec(\alpha_{2})|z|\omega(X)|w|\omega(Y)$
$\displaystyle=\sec(\alpha_{1})\sec(\alpha_{2})\omega(X)\omega(Y).$

∎

######  Remark 2.3.

  1. (1)

If $X,Y\in\mathbb{M}^{s}_{n,\alpha}$, then

$\omega(X\circ Y)\leq\sec^{2}(\alpha)\omega(X)\omega(Y).$ (12)

It is clear that inequality (12) presents a refinement of inequality (4) for
$0\leq\alpha<\frac{\pi}{4}$. In particular if $X$ and $Y$ are accretive-
dissipative, then both (12) and (4) give the same estimation.

  2. (2)

If $X,Y\in\mathbb{M}_{n}$ are accretive (or dissipative), then (12) implies
that

$\omega(X\circ Y)\leq(1+a^{2})\omega(X)\omega(Y),$

where

$a=\max\\{|\lambda_{1}\left(\left(ReX\right)^{-1}ImX\right)|,|\lambda_{1}\left(\left(ReY\right)^{-1}ImY\right)|\\}.$

The argument used to prove Theorem 2.3 can be easily modified to prove the
following more general result.

######  Theorem 2.4.

Let $X_{j}\in\mathbb{M}^{s}_{n,\alpha_{1}},j=1,2,...,m$. Then

$\omega(X_{1}\circ...\circ
X_{m})\leq\prod_{j=1}^{m}\sec(\alpha_{j})\omega(X_{j}).$

Consequently,

######  Corollary 2.5.

If $X_{1},X_{2},...,X_{m}\in\mathbb{M}_{n}$ are accretive-dissipative, then

$\omega\left(X_{1}\circ...\circ X_{m}\right)\leq
2^{m/2}\prod_{j=1}^{m}\omega(X_{j}).$

######  Corollary 2.6.

If $X_{1},X_{2},...,X_{m}\in\mathbb{M}_{n}$ are accretive (or dissipative),
then

$\omega\left(X_{1}\circ...\circ
X_{m}\right)\leq(1+a^{2})^{m/2}\prod_{j=1}^{m}\omega(X_{j}),$

where

$a=\max\\{|\lambda_{1}\left(\left(ReX_{j}\right)^{-1}ImX_{j}\right)|,j=1,2,...,m\\}.$

In the following result we give an estimation for $\omega(X\circ Y)$ in terms
of the diagonal entries of $X$ and $Y$.

######  Theorem 2.5.

Let $X=[x_{ij}]\in\mathbb{M}^{s}_{n,\alpha_{1}}$ and
$Y=[y_{ij}]\in\mathbb{M}^{s}_{n,\alpha_{2}}$. Then

$\omega(X\circ
Y)\leq\sec(\alpha_{1})\sec(\alpha_{2})\min\left\\{\max_{j}|x_{jj}|\,\,\omega(Y),\max_{j}|y_{jj}|\,\,\omega(X)\right\\}.$

###### Proof.

Since $X\in\mathbb{M}^{s}_{n,\alpha_{1}}$ and
$Y\in\mathbb{M}^{s}_{n,\alpha_{2}}$, the inequality in (11) implies that

$\omega(X\circ Y)\leq sec(\alpha_{1})\sec(\alpha_{2})\omega(Re(zX)\circ
Re(wY)),$

for some $z,w\in\mathbb{C}$ with $|z|=|w|=1$. Since $Re(zX)$ is positive
semidefinite, (8) implies that

$\omega(Re(zX)\circ
Re(wY))\leq\max_{j}\,Re(zx_{jj})\,\omega\left(Re(wY)\right).$

Now, we have

$\displaystyle\omega(X\circ Y)$
$\displaystyle=sec(\alpha_{1})\sec(\alpha_{2})\omega(Re(zX)\circ Re(wY))$
$\displaystyle\leq
sec(\alpha_{1})\sec(\alpha_{2})\max_{j}\,Re(zx_{jj})\,\omega\left(Re(wY)\right)$
$\displaystyle\leq
sec(\alpha_{1})\sec(\alpha_{2})\max_{j}\,|zx_{jj}|\,\omega\left(wY\right)$
$\displaystyle=sec(\alpha_{1})\sec(\alpha_{2})\max_{j}\,|x_{jj}|\,\omega\left(Y\right).$

Hence,

$\omega(X\circ Y)\leq
sec(\alpha_{1})\sec(\alpha_{2})\max_{j}\,|x_{jj}|\,\omega\left(Y\right).$ (13)

Similarly, we have

$\omega(X\circ Y)\leq
sec(\alpha_{1})\sec(\alpha_{2})\max_{j}\,|y_{jj}|\,\omega\left(X\right).$ (14)

The result follows by combining (13) and (14). ∎

######  Corollary 2.7.

If $X=[x_{ij}],Y=[y_{ij}]\in\mathbb{M}_{n}$ are accretive-dissipative, then

$\omega(X\circ Y)\leq
2\min\\{\max_{j}|x_{jj}|\,\,\omega(Y),\max_{j}|y_{jj}|\,\,\omega(X)\\}.$

######  Corollary 2.8.

If $X,Y\in\mathbb{M}_{n}$ are accretive (or dissipative), then

$\omega(X\circ
Y)\leq(1+a^{2})\min\left\\{\max_{j}|x_{jj}|\,\,\omega(Y),\max_{j}|y_{jj}|\,\,\omega(X)\right\\},$

where

$a=\max\\{|\lambda_{1}\left(\left(ReX\right)^{-1}ImX\right)|,|\lambda_{1}\left(\left(ReY\right)^{-1}ImY\right)|\\}.$

Another upper bound for the numerical radius of a Hadamard product of two
sectorial matrices can be obtained as follows.

######  Theorem 2.6.

Let $X\in\mathbb{M}^{s}_{n,\alpha_{1}}$ and
$Y\in\mathbb{M}^{s}_{n,\alpha_{2}}$. Then

$\omega(X\circ
Y)\leq\min\left\\{(1+\tan\alpha_{1})\omega(Re\,X)\omega(Y),(1+\tan\alpha_{2})\omega(X)\omega(ReY)\right\\}.$

Consequently,

$\omega(X\circ Y)\leq(1+\tan\alpha)\omega(X)\omega(Y),$

where $\alpha=\max\\{\alpha_{1},\alpha_{1}\\}$.

###### Proof.

The second inequality follows from the the first one and fact that
$\omega(ReX)\leq\omega(X)\,\forall X\in\mathbb{M}_{n}$. To prove the first
inequality, let $X=A+iB$ be the cartesian decomposition of $X$. Then

$\displaystyle\omega(X\circ Y)$ $\displaystyle=\omega((A+iB)\circ Y)$
$\displaystyle=\omega(A\circ Y+iB\circ Y)$ $\displaystyle\leq\omega(A\circ
Y)+\omega(B\circ Y)$
$\displaystyle\leq\omega(A)\omega(Y)+\omega(B)\omega(Y)\quad\text{by
\eqref{H}}$
$\displaystyle\leq\omega(A)\omega(Y)+\tan(\alpha_{1})\,\omega(A)\omega(Y)\quad\text{(by
part 2 of Remark \ref{R1})}$
$\displaystyle\leq(1+\tan\alpha_{1})\omega(A)\omega(Y)$
$\displaystyle=(1+\tan\alpha_{1})\omega(Re(X))\omega(Y).$

Hence,

$\omega(X\circ Y)\leq(1+\tan\alpha_{1})\omega(ReX)\omega(Y).$ (15)

Similarly, one can show that

$\omega(X\circ Y)\leq(1+\tan\alpha_{2})\omega(ReY)\omega(X).$ (16)

The result follows by combining (15) and (16). ∎

######  Corollary 2.9.

Let $X,Y\in\mathbb{M}^{s}_{n,\alpha}$. Then

$\omega(X\circ
Y)\leq(1+\tan\alpha)\min\left\\{(\omega(Re\,X)\omega(Y),\omega(X)\omega(ReY)\right\\}.$

Consequently,

$\omega(X\circ Y)\leq(1+\tan\alpha)\omega(X)\omega(Y).$ (17)

Finally, we remark that inequality (17) presents an improvement to (4) for
$0\leq\alpha<\frac{\pi}{4}$.

### Conflict of interest

The authors declare that they have no conflict of interest.

## References

  * [1] M. Alakhrass, On sectorial matrices and their inequalities, Linear Algebra and its Applications, 617, 179-189 (2021).
  * [2] M. Alakhrass, M. Sababheh, Lieb functions and sectorial matrices. Linear Algebra Appl. 586, 308–324 (2020)
  * [3] M. Alakhrass, A note on sectorial matrices, Linear Multilinear Algebra, Vol. 68, no. 11, pp. 2228–2238, 2020.
  * [4] Yu.M. Arlinski, A.B. Popov, On sectorial matrices, Linear Algebra Appl. 2003;370:133-146.
  * [5] T. Ando, K. Okubo, Induced norms of the Schur multiplier operator, Linear Algebra Appl. 147 (1991) 181–199.
  * [6] Horn, R. A. and C. A. Johnson (1991). Topics in Matrix Analysis. Cambridge, England: Cambridge University Press.
  * [7] K.E. Gustafson, D.K.M. Rao, Numerical Range, Springer, New York, 1997.
  * [8] H.-L. Gau, P.Y. Wu, Numerical radius of Hadamard product of matrices, Linear Algebra Appl. 504 (2016) 292–308.
  * [9] S. Drury, M. Lin, Singular value inequalities for matrices with numerical ranges in a sector,Oper. Matrices, 8(2014) 1143-1148.
  * [10] M. Lin, Some inequalities for sector matrices, Operators and Matrices, 10(2016) 915-921.
  * [11] M. Lin, F. Sun, A property of geometric mean of accretive operators, Linear Multilinear Algebra 65(3), 433-437 (2017).
  * [12] F. Zhang, A matrix decomposition and its applications, Linear Multilinear Algebra 10 (2015) 2033-2042.


# Applying Deep Learning Technique to Chiral Magnetic Wave Search

Yuan-Sheng Zhao  Physics Department and Center for Particle Physics and Field
Theory, Fudan University, Shanghai 200433, China  Xu-Guang Huang
<EMAIL_ADDRESS>Physics Department and Center for Particle Physics
and Field Theory, Fudan University, Shanghai 200433, China  Key Laboratory of
Nuclear Physics and Ion-beam Application (MOE), Fudan University, Shanghai
200433, China  Shanghai Research Center for Theoretical Nuclear Physics,
National Natural Science Foundation of China and Fudan University, Shanghai
200438, China

###### Abstract

The chiral magnetic wave (CMW) is a collective mode in quark-gluon plasma
originated from the chiral magnetic effect (CME) and chiral separation effect.
Its detection in heavy-ion collisions is challenging due to significant
background contamination. In Ref. Zhao  _et al._ (2022), we have constructed a
neural network which can accurately identify the CME-related signal from the
final-state pion spectra. In this paper, we generalize such a neural network
to the case of CMW search. We show that, after a updated training, the neural
network can effectively recognize the CMW-related signal. Additionally, we
assess the performance of the neural network compared to other known methods
for CMW search.

Keywords: Deep learning, Chiral magnetic wave, relativistic heavy-ion
collisions

##  I Introduction

The interplay between the chiral anomaly and external electromagnetic or
vortical fields can lead to intriguing anomalous transport phenomena in many-
body systems with chiral fermions. A notable example is the chiral magnetic
effect (CME) Kharzeev  _et al._ (2008); Fukushima  _et al._ (2008), which
induces an electric current aligned with an external magnetic field. In heavy-
ion collisions, the CME may cause charge separation relative to the reaction
plane, which can potentially be observed by analyzing the azimuthal-angle
distribution of charged hadrons using specific observables Voloshin (2004);
Abelev  _et al._ (2009). Other notable anomalous transports include the chiral
separation effect (CSE) Son and Zhitnitsky (2004); Metlitski and Zhitnitsky
(2005), the chiral vortical effect Erdmenger  _et al._ (2009); Banerjee  _et
al._ (2011); Son and Surowka (2009); Landsteiner  _et al._ (2011), and the
chiral electric separation effect (CESE) Huang and Liao (2013); Jiang  _et
al._ (2015). For reviews, see Refs Huang (2016); Kharzeev  _et al._ (2016);
Liu and Huang (2020); Kharzeev and Liao (2021); Hattori  _et al._ (2022).

In the presence of an external magnetic field, the coupled evolution of CME
and CSE gives rise to a gapless collective mode known as the chiral magnetic
wave (CMW) Kharzeev and Yee (2011). The CMW can transfer both chirality and
electric charge, potentially resulting in distinct charge and chirality
distributions. In heavy-ion collisions, given that the fireball contains a
small amount of positive charges inherited from the colliding nuclei,
theoretical studies have suggested that the CMW can induce a charge quadrupole
in the fireball, with an accumulation of positive charges at the tips and
negative charges around the equator. As the fireball expands, this quadrupole
leads to an imbalance in the elliptic flow of charged pions, specifically,
$v_{2}(\pi^{-})>v_{2}(\pi^{+})$ Burnier  _et al._ (2011). Owing to event-by-
event fluctuation of charges, some events could have net negative charges in
the fireball thus leading to $v_{2}(\pi^{-})<v_{2}(\pi^{+})$. This
characteristic feature of CMW provides a method to detect it in heavy-ion
collisions, and a series of experiments have found signals of charged pion
elliptic flow consistent with CMW expectations Adamczyk  _et al._ (2015); Adam
_et al._ (2016); Sirunyan  _et al._ (2019); Abdulhamid  _et al._ (2023);
Acharya  _et al._ (2023). However, like CME, CMW in heavy-ion collisions faces
strong background noise Deng and Huang (2012); Stephanov and Yee (2013);
Dunlop  _et al._ (2011); Bzdak and Bozek (2013); Hatta  _et al._ (2016); Xu
_et al._ (2012); Hattori and Huang (2017), which significantly obscures the
observables designed for CMW detection.

In Ref. Zhao  _et al._ (2022), we developed a CME-meter based on convolutional
neural network (CNNs) (For reviews of deep learning techniques in nuclear
physics, see Refs. Boehnlein  _et al._ (2022); He  _et al._ (2023a, b); Zhou
_et al._ (2024)). After training this CME-meter with AMPT-generated data
simulating CME (by introducing initial charge separation into AMPT model Lin
_et al._ (2005)) for Au + Au collisions at 200 GeV, the CME-meter demonstrated
exceptional robustness in distinguishing events with CME from those without.
Additionally, the CME-meter maintained strong performance across different
charge separation fractions, collision energies, and collision systems. This
success suggests the potential for creating a similar CMW-meter. As an
extension of our earlier work, we aim to increase the upper limit of salience
at the cost of some generalization capability. This approach could pave the
way for future studies on CMW physics and its detection.

In this paper, we report on the construction and performance of such a CMW-
meter. Section II details the training process, including the generation of
training samples using AMPT, the structure of the neural network, and the
training procedure. Section III examines the analysis of the trained model,
including its basic properties, comparisons to flows and observables, and a
hypothesis test. Section IV provides a summary of our findings.

##  II Construction and training of the CMW-meter

In this section, we introduce the deep learning model, data set preparation,
and training strategies employed in constructing the CMW-meter. The pion
spectra of heavy-ion collision final states serve as the input of this deep
learning model. Pions carry most of the electric charges in the final state,
making them an appropriate representation of charge distribution. A
convolutional neural network (CNN) is utilized, trained within a supervised
learning scheme to find out the CMW signals. The training data is generated
from the string-melting AMPT model Lin  _et al._ (2005), a transport model
which is widely used to simulate the evolution of both partonic and hadronic
matter in heavy-ion collisions.

In order to incorporate the CMW effect into the AMPT model, we adopt a global
charge quadrupole scheme introduced in Ref. Ma (2014). For an AMPT event with
$A_{\rm ch}>-0.01$, we interchange the positions of certain $u$ (or $\bar{d}$)
quarks in the initial state with those of $\bar{u}$ (or $d$) quarks if the
former are relatively farther from the reaction plane(RP), while the opposite
is done for events with $A_{\rm ch}<-0.01$. Here, $A_{\rm ch}$ stands for the
asymmetry of charged particle number, given by $A_{\rm
ch}=(N^{+}-N^{-})/(N^{+}+N^{-})$, where $N^{+}$ denotes the number of
positively charged particles measured in a given event, and $N^{-}$ denotes
negatively charged particles. The RP of all events is set in the $zOx-$plane.
The fraction of particles that are interchanged is represented by a relative
percentage with respect to the total number of quarks,

$f=\frac{\text{\\# Exchanged particles}}{\text{\\# All particles}}.$ (1)

According to previous study Ma (2014), switching $f=2-3\%$ of quarks generates
a CMW signal comparable to experimental observables. For training and
validation purposes, we thus choose events with a $f=2\%$ switching fraction.
The transition point, $A_{\rm ch}=-0.01$ in this scheme is based on STAR
experimental results Adamczyk  _et al._ (2015), where one may find more
details. Events at $\sqrt{s_{NN}}=200\,$GeV and different centrality are
generated for training and validation.

There are two primary reasons for training a model that results in bias and
overfitting at 200 GeV. Firstly, the pivotal issue pertains to the occurrence
of CMW in heavy-ion collisions, rather than the magnitude of the signal.
Consequently, any technique that can distinctly distinguish CMW signals from
background noise is considered valuable, irrespective of the $\sqrt{s_{NN}}$
or event centrality. Secondly, our research on the application of neural
networks for CME detection Zhao  _et al._ (2022) confirmed the robustness of
the trained network against variations in collision energy and event
centrality. The training was successful on the most comprehensive dataset,
demonstrating high accuracy levels. This means only small variance in the
network’s detection feasibility was made by them. Therefore, a model trained
on a single energy is capable of enhancing the signal detection in certain
events, while still maintaining a considerable degree of generalization.
However, further examinations involving various energies and centralities have
also been conducted to provide a more nuanced analysis.

Figure 1: A VGG-like network Simonyan and Zisserman (2015) with 4 hidden
layers is chosen in this study. Batch normalization(BN) is applied after each
hidden layer. The convolutional layers are modified to satisfy the periodic
boundary condition of the input data, and each followed by an average pooling.
10% dropout is set for the second to last dense layers.

The structure of CNN used in this work is shown in Fig. 1. There are three
2D-convolutional layers and two dense layer that contain parameters to be fit.
Some pooling layers are also applied here, providing proper reduction to data
and keeping the network from being too complicated. To encode “knowledge”
about CMW in the model, samples with and without CMW, labeled ‘1’ and ‘0’
separately, are fed to it during training, and the model is set to classify
these samples. The last activation function of the network is SoftMax, which
returns a pair of numbers $(P_{0},\,P_{1})$ for this binary classification
problem. Samples with $P_{1}>0.5$ are divided into class ‘1’, the else into
class ‘0’. Therefore, $P_{1}$ can be interpreted as the probability that a
specific sample is recognized by the neural network as containing the CMW
signal.

Data pre-processing involves several steps aimed at converting events into
analyzable samples. Initially, the spectra for mid-rapidity ($|\eta|<1$)
pions, denoted as $\rho^{\pm}(p_{T},\phi)$, are calculated, with the symbol
$\pm$ representing either $\pi^{+}$ or $\pi^{-}$, while $p_{T}$ denotes the
transverse momentum in the range of $0-2$ GeV and $\phi$ indicates the
azimuthal angle. These spectra are then segmented into histograms consisting
of 20 by 24 bins. Second, every spectrum is normalized itself so that the sum
of all bins is 1. Subsequently, a random selection of events is made, and for
each type of pion, their spectra are averaged bin by bin. These resulting
normalized and averaged pion spectra serve as the datasets for the neural
network’s training, validation, and testing phases. Unless otherwise stated,
the number of events we take in the last step is 100 in the rest of this
article. The model’s training encompasses 250 epochs, within which each epoch
contains 64 batches, and each batch comprises 100 samples. In total, 1.6
million samples are generated for the training.

##  III Performance of the CMW-meter

Accuracy, robustness and extrapolations.— As mentioned above, the model is
trained (and also validated) on samples generated at $\sqrt{s_{NN}}=200\,$GeV
that mimic final-state CMW behavior. It turns out that the model reaches high
accuracy on most events with signal at different $\sqrt{s_{NN}}$ and
centrality, as is shown in Fig. 2, which indicates preferable generalization
of the trained model. Reduction of accuracy can be found at low collision
energy and large centrality, of which the reason can be various. Different
pattern of CMW, weaker signals, stronger backgrounds or just overfitting, all
of them can account for the reduction of accuracy.

Figure 2: Accuracy of the trained model on samples (a) with CMW, and (b)
without CMW. Samples from various $\sqrt{s_{NN}}$ and centralities are
considered. The accuracy is remarkably high if the signal is surely encoded in
the samples, while those without can be mistaken as containing it, especially
at lower energy and more peripheral cases.

One of the ways to detect CMW in experiments is based on the dependency of
charge distribution and flow analysis. Specifically, the linear order
dependency of $A_{\rm ch}$ of the difference in charged-particle elliptic
flow,

$\Delta v_{2}\equiv v_{2}^{-}-v_{2}^{+}\simeq rA_{\rm ch},$ (2)

gives a measure for the CMW signal. Here $v_{2}^{-}$ and $v_{2}^{+}$ are
elliptic flows of negative- and positive-charge particles separately, and the
slope $r$ is related to the strength of the signal. Experimental results from
the STAR experiment for Au + Au collisions Abdulhamid  _et al._ (2023)
indicate that the uncertainty of the $\pi$ slope $r$ increases at lower
energies and higher centralities. Although the neural network is trained on
pions from a larger kinematic window than used in experimental analyses,
suggesting improved completeness and distinguishing capability, its
performance aligns with traditional statistical analysis trends ($\Delta
v_{2}(\pi)$). The decrease in accuracy is likely due to strong backgrounds in
those scenarios comparing to the signal strength. However, a new model can be
trained using low-energy samples or combined with high-energy samples to
create a more comprehensive training set, enhancing robustness. This approach
will be for a future, more detailed study. Overall, the trained model’s
accuracy in decoding the CMW signal is sufficient across all tested
$\sqrt{s_{NN}}$ levels, especially when focusing on high-energy samples.

As a potential detector for CMW, a measure of performance is the prediction on
non-labeled samples, where accuracy cannot be defined and the sample by sample
output becomes important. The two components of the model output are
identified as probabilities, i.e., $P_{1}$ stands for probability that the
neural network regards the input spectrum to be with CMW, and $P_{0}$ is the
probability for the other class so the two components satisfy $P_{0}+P_{1}=1$.
It is obvious that a positive correlation of $P_{1}$ with CMW signal is
expected. Events with different initial charge quadrupole fraction $f$ are
simulated, and prepared into samples as mentioned above. Tests on these
samples keep the high true-positive accuracy, yet the returned $P_{1}$ for all
$f$ are close to 1 so that little difference is made among them. To make a
clear comparison, we enlarge the difference of their output by introducing an
additional logit function,

$\text{logit}(x)=\text{log}\frac{x}{1-x}\,.$ (3)

The logit function is the inverse function of SoftMax, thus acting logit to
$P_{1}$ reveals feature space information that is encoded in the neural
network one layer before output. Besides, the logit function is monotonically
increasing, so logit($P_{1}$) keeps the correlation of $P_{1}$ and $f$
qualitatively.

Figure 3: Distribution of logit($P_{1}$) on events @ $\sqrt{s_{NN}}=200$ GeV,
centrality 30-40%. Tests are done on events with different initial charge
quadrupole($f=1,2,3,4\%$). The distribution are normalized to 1.

In Fig. 3, we present the outcomes with varying initial charge quadrupoles. As
the $f$ increases, the peak of the logit($P_{1}$) distribution shifts to the
right. In cases where $f$ equals 4%, $P_{1}$ approaches 1 so closely that the
logit function becomes numerically unstable with single precision
calculations. However, the pattern of the $f=4\%$ distribution is still in
line with the general trend. Additionally, the width of the peak remains
essentially unaffected by $f$, indicating that the model introduces minimal
error and reliably extracts the expected CMW signal. The width of the peak is
due to the event-to-event initial-state fluctuations and the method of
implementing the initial charge quadrupole. The reasonable extrapolation of
$P_{1}$ for various $f$ values suggests that the CMW strength for $f$ has been
correctly aligned to $P_{1}$ by the neural network. Consequently, it is also
indicative of the CMW signal intensity.

The model is also validated in some other tests. In test on no-CMW events
generated by UrQMD, it classifies most events correctly into ‘0’ class. To see
whether CME signal effects this CMW detector, a test set consists of AMPT
events with CME is prepared, and the trained model gives negative prediction
mostly.

Comparison with observables.— In above, we demonstrate that the trained model
efficiently decodes CMW information from $\rho^{\pm}(p_{T},\phi)$, providing a
potential measure of CMW in heavy-ion collisions. However, further comparisons
with experimental observables are necessary before constructing a measurement
based on the model. In Fig.3, we show that logit($P_{1}$) is correlated with
$f$, which in turn has a positive correlation with the slope $r$. In addition
to the slope, the following covariance between $v_{n}$ and $q_{3}$, which is
essentially a three-particle correlator, is another noteworthy observable Adam
_et al._ (2016),

$\lambda_{n}\equiv\langle v_{n}q_{3}\rangle-\langle q_{3}\rangle\langle
v_{n}\rangle,$ (4)

where $v_{n}$ is the n-th harmonic flow of the event, $q_{3}$ the charge of
the third particle, and $\langle\cdots\rangle$ denotes event average. The
differential three-particle correlator, which measures the correlation between
the flow at a particular kinematic region, and the charge of the third
particle at another particular coordinate, is more convenient when comparing
across experiments as no correction for efficiency is needed. In the following
we set $n=2$ for correlation with the elliptic flow. Using
$v_{2}^{\mp}\propto\bar{v}_{2}\pm rA_{\rm ch}/2$ and $A_{\rm ch}\sim\langle
q_{3}\rangle$, one notice that $\lambda_{2}\approx\pm r(\langle A_{\rm
ch}^{2}\rangle-\langle A_{\rm ch}\rangle^{2})/2$ for positive-charge/negative-
charge cases. In the following $\lambda_{2}$ is obtaiend by making half the
difference between the positive-charge and negative-charge cases.

Figure 4: Distribution of logit($P_{1}$) on events for Au + Au at
$\sqrt{s_{NN}}=200$ GeV and centrality 30-40%. Events are devided into
logit($P_{1}$) bins, and their $\lambda_{2}$ are averaged separately. Events
are all embed with initial charge quadrupole. A range of logit($P_{1}$) is
chosen where most events are included to avoid statistical minority. The
three-particle correlator demonstrate clearly a positive correlation with
logit($P_{1}$).

Figure 4 shows the results of the comparison between logit($P_{1}$) and
$\lambda_{2}$. The average $\lambda_{2}$ of events increase gently as the
response of the model gets stronger. Knowing logit($P_{1}$) is positvely
correlated to CMW signal, it indicates a reasonable trend in $\lambda_{2}$
when signal gets stronger. This agrees with early studies on $\lambda_{2}$
Adam  _et al._ (2016), and also shows the model prediction is qualified for
measurement.

Performance under backgrounds is necessary before we go further. There are
several mechanisms that may cause final-state $\Delta v_{2}$ and $A_{\rm ch}$
dependency as discussed in Refs. Deng and Huang (2012); Stephanov and Yee
(2013); Dunlop  _et al._ (2011); Bzdak and Bozek (2013); Hatta  _et al._
(2016); Xu  _et al._ (2012); Hattori and Huang (2017). To see how the trained
neural network works under such backgrounds, the prediction of the neural
network with different $\Delta v_{2}$ range (either with or without initial
charge quadrupoles) is studied, and the results are shown in Fig. 5. For input
sample without CMW, the prediction increases when events with larger absolute
$\Delta v_{2}$ are chosen. This shows that the neural network tends to regard
events with larger $\Delta v_{2}$ as events containing CMW signals although
they do not actually include CMW signals. However, it should be emphasized
that even in this situation, $P_{1}$ is still less than $0.5$, meaning that
the neural network still correctly classify them as events without CMW. For
samples with CMW, the model shows strong robustness against the background,
and the model makes classification correctly for all the samples.

Figure 5: Distribution of $P_{1}$ on events at $\sqrt{s_{NN}}=200$ GeV against
$\Delta v_{2}$. Events are devided in to 10 $\Delta v_{2}$ bins, and their
$P_{1}$ are averaged separately. As the magnitude of $\Delta v_{2}$ increases,
the tendency of the model to make a false positive classification also
increases. Nevertheless, in events involving an initial quadrupole, the model
consistently maintains a high level of accuracy.

Hypothesis test.— As previously discussed, the neural network demonstrates
good accuracy in predicting the CMW signal and exhibits robustness across
different collision energies, centralities, and background effects after
training. This makes it feasible to create a CMW-meter based on the neural
network. However, the need to average events poses a challenge for deploying
this measurement experimentally, as it is not possible to know the charge
quadrupole pattern in advance or align events according to their charge
distribution patterns.

However, from a hypothesis test perspective, the CMW-meter also holds
experimental feasibility. For a fixed finite number, $M$, of events, one can
assume the presence of a sufficiently large residual quadrupole that can be
detected through our meter if CMW is believed to exist in these events.
Conversely, if no CMW is observed in experiments, the neural network’s
predictions will consistently fall within the ‘0’ class. As demonstrated in
Fig. 3, the intensity of the CMW signal significantly alters the distribution
of logit$(P_{1})$ or $P_{1}$, thus influencing the distribution of $P_{1}$
itself (denoted by $\text{P}(P_{1})$). This distribution responds differently
depending on the presence or absence of CMW in the data set. If CMW exists in
the heavy-ion collisions, the neural network model’s prediction regarding the
residual quadrupole of a sample will align with $\text{P}(P_{1})$ in the
$f\neq 0$ case. In contrast, without any CMW, the distribution will match the
$f=0$ scenario. To establish a reasonable estimation of $\text{P}(P_{1})$ for
testing $M$ events, we treat $f$ as a latent variable representing CMW in a
single event, as defined by the initial charge quadrupole fraction used in
this work. Given event-by-event fluctuations, we model $f$ as a random
variable following a Gaussian distribution, $f\sim N(\mu,\sigma^{2})$, where
$\mu$ is the mean of the latent variable $f$, which is expected to be around
0. The variance $\sigma$ is estimated based on Ma (2014), where the average of
$|f|$ is approximately $2\%$,

$\displaystyle 2\%=\langle|f|\rangle$
$\displaystyle=\int|f|N^{>}(|f|;\sigma^{2})\;d|f|,$ (5)

here $N^{>}(|f|;\sigma^{2})$ is the half normal distribution, and $|f|$ is a
positive-definite variable because the model prediction is independent to the
sign of $f$. Solving (5) gives $\sigma\simeq 0.025$. Because we adopt averaged
events in preparing the CMW-meter, the way to compose a $\rho(f_{\text{eff}})$
from single events $\\{\rho(f_{i})\\}$ becomes crucial, where $f_{\text{eff}}$
is the effective charge quadrupole rate of averaged events. One can choose the
arithmetic mean as,

$\frac{1}{M}\sum_{i}^{M}\rho(f_{i})=\rho(\frac{1}{M}\sum_{i}^{M}f_{i})=\rho(f_{\text{eff}}).$
(6)

Therefore, the distribution of $|f_{\text{eff}}|$ can be achieved as

$F_{\text{eff}}\sim
N({\mu}/{M},M\,{\sigma^{2}}/{M^{2}})=N(0,{\sigma^{2}}/{M}),$ (7)

with $F_{\text{eff}}\equiv|f_{\text{eff}}|\sim N^{>}(\sigma^{2}/M)$. The
conditional probability $\text{P}(P_{1}|\;F_{\text{eff}})$ can be approximated
as a Beta distribution,

$\text{Be}(x;\alpha,\beta)=\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}x^{\alpha-1}(1-x)^{\beta-1},$
(8)

with $\alpha$ and $\beta$ the parameters of the Beta distribution, and
$\Gamma$ is Gamma function. To describe $\text{P}(P_{1}|\;F_{\text{eff}})$ at
any $F_{\text{eff}}$, we assume $\alpha$ and $\beta$ are functions of
$F_{\text{eff}}$, and fit several sets of $(\alpha,\beta)$ from the fitted
beta distribution with polynomial(for $\alpha$) and Softplus (for $\beta$, to
reach proper asymptotic behavior around $\text{P}_{1}=1$).

After parameterizing $\text{P}(P_{1}|\;F_{\text{eff}})$, $\text{P}(P_{1})$ is
derived as

$\displaystyle\text{P}(P_{1})$
$\displaystyle=\int\text{P}(P_{1}|\,F_{\text{eff}})\text{P}(F_{\text{eff}})\,dF_{\text{eff}}$
(9)
$\displaystyle=\int_{0}^{\infty}\text{Be}(P_{1};\alpha(F_{\text{eff}}),\beta(F_{\text{eff}}))N^{>}(\frac{\sigma^{2}}{M})\,dF_{\text{eff}}.$

The numerical results are shown in Fig. 6. $\text{P}(P_{1})$ of the “existing
CMW” has an obvious rise around $P_{1}=1$ compared to the “no CMW” case, which
suggests a non-zero probability of composing a large residual quadrupole. With
a smaller $M$, the width of $f_{\text{eff}}$ becomes larger, which allows one
to get a visible $P_{1}$. In Fig. 6 we also present results of randomly mixing
events of both charge quadrupole patterns generated by AMPT, where $M=$ 25\.
For both large- and small-$\text{P}_{1}$ area, it is consistent with our
hypothesis test analysis qualitatively, which indicates that the trained
neural network is capable of recognizing charge quadrupole with less averaged
events.

Figure 6: $\text{P}(P_{1})$ of AMPT events with $M=25$. The peak of this
distribution near $P_{1}=0$ suggests that the charge quadrupole in most events
negate one another, and the peak near $P_{1}=1$ is because the strong response
to some large enough residual quadrupole.

##  IV Summary

In this paper, we propose a deep convolutional neural network model for CMW
detection. Building upon our prior work Zhao  _et al._ (2022) focused on deep-
learning-based CME detection, this model expands its application to CMW
detection. We train the neural network using data generated from the AMPT
model for Au + Au collisions at 200 GeV, with CMW-like initial charge
quadrupole encoded. The trained model exhibits a robust capability to discern
events with CMW from those without, and it can quantitatively measure the
fraction or strength of the initial charge quadrupole, effectively functioning
as a CMW-meter. Furthermore, we validate the model’s performance across a
broad collision energies and centralities demonstrating its resilience. We
have also checked that the trained model can be well qualified even for other
collision systems like Zr + Zr and Ru + Ru collisions. Comparative analysis
against three-particle correlators and $\Delta v_{2}$ illustrates the model’s
effectiveness even in the presence of strong backgrounds. By employing a
hypothesis test, an experimentally viable analysis based on the model can be
established, wherein the distribution of model predictions serves as an
indicator of CMW occurrence in the data.

One drawback of the model is that its brightness is achieved at the cost of
generalization ability, as the training data is confined to a narrow range of
collision energies. In the future, it would be interesting to enhance the
model’s generalization capabilities and transform it into an end-to-end CMW
meter.

Acknowledgement — We acknowledge the useful discussions with L.-X. Wang and K.
Zhou. This work is supported by the Natural Science Foundation of Shanghai
(Grant No. 23JC1400200), the National Natural Science Foundation of China
(Grant No. 12225502, No. 12075061, and No. 12147101), and the National Key
Research and Development Program of China (Grant No. 2022YFA1604900).

## References

  * Zhao  _et al._ (2022) Y.-S. Zhao, L. Wang, K. Zhou, and X.-G. Huang, Phys. Rev. C 106, L051901 (2022), arXiv:2105.13761 [hep-ph] .
  * Kharzeev  _et al._ (2008) D. E. Kharzeev, L. D. McLerran, and H. J. Warringa, Nucl. Phys. A 803, 227 (2008), arXiv:0711.0950 [hep-ph] .
  * Fukushima  _et al._ (2008) K. Fukushima, D. E. Kharzeev, and H. J. Warringa, Phys. Rev. D 78, 074033 (2008), arXiv:0808.3382 [hep-ph] .
  * Voloshin (2004) S. A. Voloshin, Phys. Rev. C 70, 057901 (2004), arXiv:hep-ph/0406311 .
  * Abelev  _et al._ (2009) B. I. Abelev _et al._ (STAR), Phys. Rev. Lett. 103, 251601 (2009), arXiv:0909.1739 [nucl-ex] .
  * Son and Zhitnitsky (2004) D. T. Son and A. R. Zhitnitsky, Phys. Rev. D 70, 074018 (2004), arXiv:hep-ph/0405216 .
  * Metlitski and Zhitnitsky (2005) M. A. Metlitski and A. R. Zhitnitsky, Phys. Rev. D 72, 045011 (2005), arXiv:hep-ph/0505072 .
  * Erdmenger  _et al._ (2009) J. Erdmenger, M. Haack, M. Kaminski, and A. Yarom, JHEP 01, 055 (2009), arXiv:0809.2488 [hep-th] .
  * Banerjee  _et al._ (2011) N. Banerjee, J. Bhattacharya, S. Bhattacharyya, S. Dutta, R. Loganayagam, and P. Surowka, JHEP 01, 094 (2011), arXiv:0809.2596 [hep-th] .
  * Son and Surowka (2009) D. T. Son and P. Surowka, Phys. Rev. Lett. 103, 191601 (2009), arXiv:0906.5044 [hep-th] .
  * Landsteiner  _et al._ (2011) K. Landsteiner, E. Megias, and F. Pena-Benitez, Phys. Rev. Lett. 107, 021601 (2011), arXiv:1103.5006 [hep-ph] .
  * Huang and Liao (2013) X.-G. Huang and J. Liao, Phys. Rev. Lett. 110, 232302 (2013), arXiv:1303.7192 [nucl-th] .
  * Jiang  _et al._ (2015) Y. Jiang, X.-G. Huang, and J. Liao, Phys. Rev. D 91, 045001 (2015), arXiv:1409.6395 [nucl-th] .
  * Huang (2016) X.-G. Huang, Rept. Prog. Phys. 79, 076302 (2016), arXiv:1509.04073 [nucl-th] .
  * Kharzeev  _et al._ (2016) D. E. Kharzeev, J. Liao, S. A. Voloshin, and G. Wang, Prog. Part. Nucl. Phys. 88, 1 (2016), arXiv:1511.04050 [hep-ph] .
  * Liu and Huang (2020) Y.-C. Liu and X.-G. Huang, Nucl. Sci. Tech. 31, 56 (2020), arXiv:2003.12482 [nucl-th] .
  * Kharzeev and Liao (2021) D. E. Kharzeev and J. Liao, Nature Rev. Phys. 3, 55 (2021), arXiv:2102.06623 [hep-ph] .
  * Hattori  _et al._ (2022) K. Hattori, M. Hongo, and X.-G. Huang, Symmetry 14, 1851 (2022), arXiv:2207.12794 [hep-th] .
  * Kharzeev and Yee (2011) D. E. Kharzeev and H.-U. Yee, Phys. Rev. D 83, 085007 (2011), arXiv:1012.6026 [hep-th] .
  * Burnier  _et al._ (2011) Y. Burnier, D. E. Kharzeev, J. Liao, and H.-U. Yee, Phys. Rev. Lett. 107, 052303 (2011), arXiv:1103.1307 [hep-ph] .
  * Adamczyk  _et al._ (2015) L. Adamczyk _et al._ (STAR), Phys. Rev. Lett. 114, 252302 (2015), arXiv:1504.02175 [nucl-ex] .
  * Adam  _et al._ (2016) J. Adam _et al._ (ALICE), Phys. Rev. C 93, 044903 (2016), arXiv:1512.05739 [nucl-ex] .
  * Sirunyan  _et al._ (2019) A. M. Sirunyan _et al._ (CMS), Phys. Rev. C 100, 064908 (2019), arXiv:1708.08901 [nucl-ex] .
  * Abdulhamid  _et al._ (2023) M. I. Abdulhamid _et al._ (STAR), Phys. Rev. C 108, 014908 (2023), arXiv:2210.14027 [nucl-ex] .
  * Acharya  _et al._ (2023) S. Acharya _et al._ (ALICE), JHEP 12, 067 (2023), arXiv:2308.16123 [nucl-ex] .
  * Deng and Huang (2012) W.-T. Deng and X.-G. Huang, Phys. Rev. C 85, 044907 (2012), arXiv:1201.5108 [nucl-th] .
  * Stephanov and Yee (2013) M. Stephanov and H.-U. Yee, Phys. Rev. C 88, 014908 (2013), arXiv:1304.6410 [nucl-th] .
  * Dunlop  _et al._ (2011) J. C. Dunlop, M. A. Lisa, and P. Sorensen, Phys. Rev. C 84, 044914 (2011), arXiv:1107.3078 [hep-ph] .
  * Bzdak and Bozek (2013) A. Bzdak and P. Bozek, Phys. Lett. B 726, 239 (2013), arXiv:1303.1138 [nucl-th] .
  * Hatta  _et al._ (2016) Y. Hatta, A. Monnai, and B.-W. Xiao, Nucl. Phys. A 947, 155 (2016), arXiv:1507.04690 [hep-ph] .
  * Xu  _et al._ (2012) J. Xu, L.-W. Chen, C. M. Ko, and Z.-W. Lin, Phys. Rev. C 85, 041901 (2012), arXiv:1201.3391 [nucl-th] .
  * Hattori and Huang (2017) K. Hattori and X.-G. Huang, Nucl. Sci. Tech. 28, 26 (2017), arXiv:1609.00747 [nucl-th] .
  * Boehnlein  _et al._ (2022) A. Boehnlein _et al._ , Rev. Mod. Phys. 94, 031003 (2022), arXiv:2112.02309 [nucl-th] .
  * He  _et al._ (2023a) W. He, Q. Li, Y. Ma, Z. Niu, J. Pei, and Y. Zhang, Sci. China Phys. Mech. Astron. 66, 282001 (2023a), arXiv:2301.06396 [nucl-th] .
  * He  _et al._ (2023b) W.-B. He, Y.-G. Ma, L.-G. Pang, H.-C. Song, and K. Zhou, Nucl. Sci. Tech. 34, 88 (2023b), arXiv:2303.06752 [hep-ph] .
  * Zhou  _et al._ (2024) K. Zhou, L. Wang, L.-G. Pang, and S. Shi, Prog. Part. Nucl. Phys. 135, 104084 (2024), arXiv:2303.15136 [hep-ph] .
  * Lin  _et al._ (2005) Z.-W. Lin, C. M. Ko, B.-A. Li, B. Zhang, and S. Pal, Phys. Rev. C 72, 064901 (2005), arXiv:nucl-th/0411110 .
  * Ma (2014) G.-L. Ma, Phys. Lett. B 735, 383 (2014), arXiv:1401.6502 [nucl-th] .
  * Simonyan and Zisserman (2015) K. Simonyan and A. Zisserman, “Very deep convolutional networks for large-scale image recognition,” (2015), arXiv:1409.1556 [cs.CV] .