Taxonomy of Machine Learning

ML can be classified in several ways. One way is to classify it by its function, i.e., into classification, regression, generation, and dimensionality reduction. The other way is to classify ML by the type of training data, i.e., into supervised learning, unsupervised learning, semi-supervised learning, self-supervised learning, active learning, and reinforcement learning. For example, supervised learning requires data to be labeled in such a way that the model can be trained to build a mapping between the input and the labels. Unsupervised learning does not need labeled data; it can learn patterns from data, assuming that the machine makes self-consistent predictions on data that are perturbed or slightly augmented. Semi-supervised learning requires a small amount of labeled data along with a large amount of unlabeled data. Self-supervised learning works with specific data such as natural language or images that are sequential. It allows the machine to predict one part of the sequence from the other part. Active learning is a type of semi-supervised learning that employs two pools of data: a small pool of labeled data and a large pool of unlabeled data. The machine is trained on the labeled data and validated on the unlabeled data. The performance of the simply trained machine differs for different samples from the unlabeled data pool. For example, the machine may be uncertain on one sample, predicting that the label of the sample is A with 51% probability and B with 49% probability. This sample is assumed to be more difficult and more important for the trained machine than simple samples for which the machine’s predictions are certain. For efficiency, this sample is labeled and moved from the unlabeled pool to the labeled pool for further training. Reinforcement learning uses data generated by interactions with the environment.

According to the previous description, the loss function for supervised learning in the regression task can be expressed as $l=\left|\right|y_{\text{pred}}-y_{\text{true}}\left|\right|\mathrm{,}$ (1) where y_pred=f(x, θ) is the function represented by ML models such as decision trees or deep neural networks (DNNs), x represents the input data, θ represents all trainable model parameters, and y_true is the label of the input data x. $\left|\right|\cdot \left|\right|$ usually represents the l₁ norm, which gives the mean absolute error, or the l₂ norm, which gives the root-mean-square error.

The cross-entropy loss is widely used for classification. It is defined as $l=-{\displaystyle \sum _{k=1}^K}{p}_k\log q_k$ (2) where K represents the number of possible categories of the input data x, pk = y_true is the true label (probability), and qk = f(x, θ) represents the network prediction. This loss is inspired by the Kullback–Leibler (KL) divergence, which quantifies the difference between two distributions p and q: $\text{KL}\left(p\left|\right|q\right)={\displaystyle \sum _{k=1}^K{p}_k}\log p_k{q}_k$ (3) $={\displaystyle \sum _k{p}_k}\log p_k-p_k\log q_k$ (4) $=-H(p)+H(p\mathrm{,}q),$ (5) where H(p) represents the entropy of the distribution p and the cross entropy H(p, q)=∑kpk log qk quantifies the average number of bits needed to encode the distribution p using the model q.

In binary classification, the cross entropy is reduced to $l=-\frac{1}{m}{\displaystyle \sum _{i=1}^m\left[p_i\log q_i+\left(1-p_i\right)\log \left(1-p_i\right)\right]}\mathrm{,}$ (6) where pi is the true label of the ith sample, whose value is 0 or 1. qi represents the network prediction obtained using the sigmoid activation function in the last layer to ensure 0 < qi < 1. m represents the number of samples in each minibatch. If the true label is pi = 0, only the 2nd part contributes to the loss function.

For multi-categorical classification, the loss function is the cross-entropy loss, with the activation function in the last layer replaced by the softmax function: $\text{Softmax}(z_i)=\frac{e^{z_i}}{{\displaystyle {\sum }_{k=1}^K{e}^{z_k}}}.$ (7) For unsupervised learning, the loss function can generally be expressed as $l=\Vert {\text{manipulate}}_1(x)-{\text{manipulate}}_2(x)\Vert ,$ (8) where ${\text{manipulate}}_{\mathrm{1,2}}$ represents two manipulations on the same data. For example, in clustering tasks, the manipulations on x are to compute the total distance s of samples to multiple centers. In image classification tasks, the manipulations are to compute the network prediction over two different augmentations of the same image, e.g., cropping or rotation. This loss is also called self-consistent loss.

For semi-supervised learning, the loss function is the combination of the supervised loss and unsupervised loss: $l=l_{\text{supervised}}+l_{\text{unsupervised}}.$ (9) For self-supervised learning, a widely used loss function is the reconstruction loss. For example, in computer vision, the reconstruction loss is defined as the difference between the original image and the image reconstructed by a neural network from a masked image. $l=\Vert x-f\left(\left(1-M\right)\cdot x\right)\Vert$ (10) Here, x represents the original image, M represents the binary mask used to remove M=0 pixels from the image, and f represents a neural network used to reconstruct the image. The same method can be used to reconstruct natural language, by predicting the next sentence or missing words in a sentence. The pretrained network can be used in many downstream tasks, such as classification, regression, or generation.

In active learning, the loss function is essentially the same as that in supervised learning. The difference is that the trained network ranks samples from the unsupervised pool for annotation. Thus, the key is to rank the samples. There are two main methods for this. One is to rank the samples according to the entropy of the predictions made by the pretrained network: $s=-{\displaystyle \sum _i{p}_i}\log p_i,$ (11) where pi represents the predicted probability that the sample is in class i. The other method is to rank the samples according to the diversity of the training dataset, by giving the highest rank to the sample that has the longest distance from the training data.

For reinforcement learning, the data are generated by subsequent interactions between the network policy and the environment. The network receives an observation ot from the environment at time t, makes a decision, and performs an action at on the environment. The environment returns a new observation o_t+1, an immediate reward r_t+1, and a done signal. The data are thus $\left\{o_t\mathrm{,}{a}_t\mathrm{,}{o}_{t+1}\mathrm{,}{r}_{t+1}\mathrm{,}\text{done}\right\}$ trajectories. The loss function of reinforcement learning is similar to that of supervised learning, with data ot and true labels at, r_t+1.

Optimization

The goal of ML is to minimize the loss for the prediction of new data not used for training. In gradient-based models, this is achieved simply via stochastic gradient descent (SGD) and its variants: $\theta =\theta -ϵ\frac{1}{m}{\displaystyle \sum _{i=1}^m\frac{\partial l_i}{\partial \theta }},$ (12) where θ represents all the trainable parameters of the ML model, ϵ represents a small positive number called the learning rate, and m represents the size of the mini-batch. Updating θ with the negative gradient $-ϵ\frac{1}{m}{\displaystyle {\sum }_{i=1}^m\frac{\partial l_i}{\partial \theta }}$ helps to gradually reduce the loss. This can be easily verified if there is only one trainable parameter θ and the loss is l = θ², whose negative gradient is -2θ.

The possible values of θ form a space called the parameter space. The initial value of θ is usually a random number. Updating θ using SGD is analogous to walking around the parameter space looking for the minimum value of the loss function. The loss function can be thought of as the potential surface whose negative gradients give the direction of acceleration $\overrightarrow{a}$. Thus, in simple SGD, the position θ in the parameter space is updated using the acceleration. Consequently, native SGD has two major drawbacks. First, if the gradient is 0, the optimization stops immediately. Second, the network update is far faster in the direction where the gradient is large r. These two drawbacks are partly solved using the momentum mechanism [16] and the adaptive learning rate [17].

In reinforcement learning, the goal is to maximize the accumulated rewards. The optimization method is stochastic gradient ascent. In the popular policy gradient method, the parameters of the policy network are updated as follows: $\theta =\theta +ϵG_t\nabla \ln \pi \left(a_t|o_t\mathrm{,}\theta \right),$ (13) where $G_t={\displaystyle {\sum }_{k=t+1}^T{\gamma }^{k-t-1}{r}_k}$ is the return representing the accumulated rewards in the future with a discounting factor γ < 1.

Automatic Differentiation

The number of trainable parameters in a DNN is large. To learn from the data, one must compute the negative gradients of loss with respect to each of the millions or trillions of model parameters $-\frac{\partial l}{\partial \theta }$. This is intractable using finite difference or analytic differentiation. Finite difference has both truncation and round-off errors that cannot be controlled. Analytical differentiation has exploding expressions for DNNs that are too complex to compute efficiently. In deep learning, the negative gradient is mainly computed using automatic differentiation (AD), which is computationally efficient meanwhile also has analytical precision.

AD has a forward mode and a backward mode. If the DNN is a $R^1\to R^n$ mapping, a forward pass gives derivatives of all output variables yi with respect to the input variable x. In contrast, if the network is an $R^n\to R^1$ mapping, each forward pass returns only the derivative of the output variable y on one of the input variables xi. In the SGD algorithm, the backward mode is far more efficient because the mapping from θ to the loss is an $R^n\to R^1$ mapping. In the following, the forward AD is briefly explained.

In the forward mode, AD is implemented by introducing a dual number for each variable: $x\to x+\dot{x}d,$ (14) $y\to y+\dot{y}d,$ (15) where x and y are two variables that require gradients, and $\dot{x}$ and $\dot{y}$ are the derivatives of x and y, respectively, with respect to some variable. As mentioned previously, with $\dot{x}=\mathrm{1,}\dot{y}=0$ gives $\partial l/\partial x$ in one pass of the forward mode, and setting $\dot{x}=\mathrm{0,}\dot{y}=1$ gives $\partial l/\partial y$ in another pass of the forward mode. $d$ is an infinitesimal symbol satisfying $d^2=0$, which is analogous to the imaginary symbol $I^2=-1$. With this definition, the traditional output z of each operator is a dual number $z+\dot{z}d$ whose coefficient $\dot{z}$ is the derivative of z, as follows: $x+\dot{x}d+y+\dot{y}d=\left(x+y\right)+\left(\dot{x}+\dot{y}\right)d,$ (16) $x+\dot{x}d-\left(y+\dot{y}d\right)=\left(x-y\right)+\left(\dot{x}-\dot{y}\right)d,$ (17) $\left(x+\dot{x}d\right)*\left(y+\dot{y}d\right)=x*y+\left(x\dot{y}+y\dot{x}\right)d,$ (18) $\left(x+\dot{x}d\right)/\left(y+\dot{y}d\right)=\frac{\left(x+\dot{x}d\right)\left(y-\dot{y}d\right)}{y^2-{\dot{y}}^2{d}^2}$ (19) $=\frac{x}{y}+\frac{y\dot{x}-x\dot{y}}{y^{2}d}$ (20) The calculations of dual numbers can easily be extended to polynomial functions: $P\left(x+\dot{x}d\right)=P(x)+P^{\prime }(x)d.$ (21) Using a computer, more complex functions such as sin(x), log x, and ex can be approximated by polynomial functions. In principle, AD works for these functions as well. In practice, these functions can be overloaded to produce outputs in the form of dual numbers, e.g., $\sin \left(x+\dot{x}d\right)\to \sin x+\cos x\dot{x}d$.

Because of the universal approximation capability of DNNs and the efficient and accurate auto-diff, DNNs are widely used to represent solutions of ordinary differential equations (ODEs) and partial differential equations (PDEs) that require gradients. Thus, many physical problems are translated into optimization problems. This method is commonly referred to as physics-informed neural networks (PINNs). Compared with traditional numerical solutions, PINNs are mesh-free, work for very high dimensions, and are easy to implement—particularly for multi-scale and multi-physics problems.

Convolutional Neural Networks

Convolutional neural networks (CNNs) are distinguished from other neural networks by their superior performance for image, speech, and audio signal inputs. A naive CNN consists of three main types of layers, i.e., convolutional layers, pooling layers, and fully connected layers, as shown in Fig. 1.

Fig. 1Full-screen View

(Color online) CNNs.

The convolutional layer is the core building block of a CNN. The term convolution refers to the convolution operation between the input features and the filters (or kernels). In the mathematical view, a convolution operation is a special type of linear operation where two functions are multiplied to produce a third function that expresses how the shape of one function is modified by the other. In the ML view, the convolutional layer uses the filters to extract the features from the input data and combines the extracted features as the output. In a well-trained convolutional layer, a filter is only sensitive to one specific type of feature. Usually, there are many filters in a convolutional layer, to satisfy the complex input features. After the convolution operation, a rectified linear unit (ReLU) activation function is typically used, which introduces nonlinearity into the neural network.

After the convolutional layer, a pooling layer is applied to reduce the number of parameters, which is also known as downsampling. There are two main types of pooling: max pooling and average pooling. Max pooling selects the maximum value to be the output, and average pooling uses the average of the pixels covered by the pooling kernel. The fully connected layer is used to map the features extracted by the previous layers to the final output.

The convolutional layers can be stacked to make the neural network deeper. Earlier layers break down the complex features from the input data into individual simple features. As the features pass through the subsequent convolutional layers, the filters begin to capture larger elements or shapes. Owing to its ability to extract complex features, the CNN architecture became a foundation of modern computer vision.

However, when neural networks are deep, the vanishing gradient problem is severe. To overcome this problem in CNN architectures, many complex neural networks have been developed, such as AlexNet, VGGNet, InceptionNet, GoogLeNet, and ResNet.

Recurrent Neural Networks

Recurrent neural networks (RNNs) are distinguished from other neural networks by their superior performance for sequence or time-series data.

Fig. 2 shows the structure of a basic RNN, where U denotes the weights for the connection of the input layer to the hidden layer, V denotes the weights for the connection of the hidden layer to the hidden layer, and W denotes the weights for the connection of the hidden layer to the output layer. Using self-connection with weights V, the RNN takes information from previous inputs to influence the current input and output. This feature, which is often referred to as “memory,” makes the RNN good at processing sequential data. The loss function $\mathcal{L}$ of all timesteps is defined according to the loss at each timestep, as follows: $\mathcal{L}\left(\widehat{Y}\mathrm{,}Y\right)={\displaystyle \sum _{t=1}^T\mathcal{L}}\left({\widehat{Y}}^t\mathrm{,}{Y}^t\right)\mathrm{.}$ (22) The RNN uses the backpropagation through time (BPTT) algorithm to determine the gradients. The error is backpropagated from the last timestep to the first timestep. At timestep T, the derivative of the loss 𝒧 with respect to the weight matrix W is expressed as follows: $\frac{\partial {\mathcal{L}}^{(T)}}{\partial W}={\displaystyle \sum _{t=1}^T\frac{\partial {\mathcal{L}}^{(t)}}{\partial W}}.$ (23) RNNs also suffer from the problems of gradient vanishing and exploding. To deal with the gradient problems, variant networks have been developed, such as long short-term memory (LSTM) networks and gated recurrent units (GRUs).

Fig. 2Full-screen View

(Color online) RNNs.

Point Cloud Network

The final-state particles from HICs form a point cloud in the momentum space. The data must be manipulated to use CNNs and RNNs, as these networks were originally designed for images and natural language. For example, to use a CNN, density estimation (histogram) is typically used to convert the particle cloud into images. However, this does not work well for a few particles in a three-dimensional (3D) space, because the particles are dilute and the resolution is poor. To use an RNN, the particle cloud must be sorted to one dimension, which can only keep the local information in one dimension. The point cloud network is designed to preserve the permutation symmetry of a set of particles.

Fig. 3 shows a simple demonstration of a point cloud network. The input to the network is a set of particles in the momentum space, including their 4-momenta, mass, and other quantum numbers. A fully connected neural network or multilayer perceptron (MLP) is applied to a particle to transform its m input features into 128 features in the high-dimensional latent space. The MLP is shared by all the particles in the cloud and is also called a 1DCNN. This step preserves the permutation symmetry of all the particles. Then, global max pooling (GMP) or global average pooling (GAP) is applied to these latent features of the particles to extract the global information of the particle cloud. The GMP and GAP extract the boundaries of the input particle cloud in the high-dimensional latent space, which learn the multi-particle correlation for the final decision. This extracted global information (128 features) is fed to another MLP for the final decision. The output neuron has a value in the range (0, 1) and uses 0.5 as the decision boundary.

Fig. 3Full-screen View

(Color online) Simple example of a point cloud network.

The network shown in Fig. 3 is used to classify nuclear phase transitions [18]. Some point cloud networks apply a Euclidean rotation to the point cloud to preserve rotational symmetry, i.e., the network should make self-consistent prediction if the point cloud is rotated globally [19]. Other variants use k-nearest neighbors in the spatial or momentum space to extract the high-dimensional latent features of each particle, for keeping more local correlation. The k-nearest neighbors of each particle can be calculated in the feature space to capture the long-range multiple particle correlation, because particles that are close in the feature space may be far apart in the spatial or momentum space. This technique is called dynamical edge convolution and was used to search for self-similarity between particles in the momentum space, which is associated with critical phenomena that may occur in HICs [20]. The dynamical edge convolutional neural network is a type of message-passing neural network that is also called a graph neural network.

Generative Modeling

In unsupervised learning, generative modeling is a class of techniques related to probability distribution learning. With regard to tasks, ML can generally be categorized into discriminative modeling and generative modeling. From probabilistic perspectives, discriminative modeling, such as pattern recognition, aims at learning a conditional probability p(y|x), which can be used to make predictions for a given input object (x) its associated properties or class identities (y), while the goal of generative modeling is to capture the joint distribution p(x,y), from which one can generate new data points following the same statistics as the training set. Generative modeling has achieved considerable success in numerous applications, including image synthesis, inpainting, super-resolution, text-to-image translation, speech generation, and chat robots. Many of the generative models were developed with profound influence from and on physics. Generative modeling also has numerous direct applications in science, e.g., computational fluid simulation, drug molecule design, anomaly detection, many-body physics, and lattice field configuration generation for QCD.

The central purpose of generative modeling is to sample data ($\tilde{x}$) from the same distribution of the training set pd(x). Most of the generative modes construct parametric (explicit or implicit) models pθ(x) to approach the desired data distribution. From information theory, the KL divergence (Eq.(3)), which measures the dissimilarity between the model and data distributions, provides an objective for this task. Per Jensen inequality, the KL divergence is non-negative and is zero only when the two distributions match exactly. The minimization of the KL divergence under given observational data for the system with the collected training set, i.e., $\mathcal{D}=\left\{x\right\}$, is equivalent to the minimization of the negative log-likelihood (NLL): $\mathcal{L}=\frac{1}{\left|\mathcal{D}\right|}{\displaystyle \sum _{x\in \mathcal{D}}\log }{p}_{\theta }(x);$ (24) thus, maximum likelihood estimation (MLE) is performed.

In the following, we briefly review several representative and popular deep generative models, including the variational autoencoder (VAE), generative adversarial networks (GANs), autoregressive modeling, and normalizing flows (NF).

VAE [21], introduces a latent variable z to facilitate the generation process; thus, it constructs a trainable conditional probability $p_{\theta }\left(x|z\right)$ (called the decoder or generator, usually modeled by a neural network). For generation convenience, the latent variable is assumed to follow an easy-to-sample prior distribution p(z), such as the multivariate Gaussian distribution. However, the introduction of the latent variable makes the data generation distribution (thus the likelihood) intractable, because the required marginalization is $p_{\theta }(x)={\displaystyle \int }{p}_{\theta }\left(x|z\right)p(z)dz$. Thus, the posterior distribution for the latent variable is intractable as well, because $p(z|x)=p_{\theta }\left(x|z\right)p(z)/p_{\theta }(x)$. The VAE employs a variational inference approach to approximately perform MLE on the training data. Specifically, an encoder model qϕ(z|x) (also a neural network) is introduced to approach the real posterior p(z|x), and the KL divergence ${\mathcal{D}}_{KL}(q_{\varphi }(z|x)\left|\right|p(z|x)$ provides the training objective, which is derived as a variational lower bound (also known as an evidence lower bound (ELBO), which is the cornerstone of the VAE) of the likelihood: $\begin{array}{l}\mathcal{L}={\mathbb{E}}_{q_{\varphi }(z|x)}[\log p_{\theta }\left(x|z\right)+\log p(z)\\ -\log q_{\varphi }(z|x)]\le \log p_{\theta }(x),\end{array}$ (25) $\theta ,\varphi =\underset{\theta ,\varphi }{\text{arg\hspace{0.17em}max\hspace{0.17em}}}\mathcal{L}.$ (26) The generative adversarial network, i.e., GAN, as another latent variable generative model, is developed to train the generator through an adversarial strategy. Intuitively, the GAN framework constructs two nonlinear differentiable functions (both represented by adaptive neural networks with dimensionality accordingly set). The first one—called the generator G(z)—maps latent variable z to the target data manifold $\tilde{x}=G(z)$, which gives an implicit synthesized data distribution pG(x) when the latent variable follows a prior latent distribution p(z), e.g., multivariate uniform or Gaussian, and the goal is training the generator to drive pG(x) approaching the target distribution ptrue(x). The other one—called the discriminator D(x)—maps the data manifold to a single scalar representing the fake-vs.-true distinguishing result of the discriminator for the input data. For a vanilla GAN, the discriminator is designed as a binary classifier; i.e., for real and generated data, it is trained to output D(x)=1 and $D(\tilde{x})=0$, respectively. The generator and discriminator are trained alternately to improve their abilities in competing against each other; this can be achieved by mimicking a two-player min-max game. Thus, the discriminator is trained to better distinguish the real data from the generated data, while the generator is trained to trick the discriminator into classifying generated data as “real” data.

It was proven mathematically that the adversarial training of a GAN is equivalent to minimizing the Jensen–Shannon divergence: $\text{JS}\left(p_{\text{real}}\left|\right|p_G\right)=\frac{1}{2}\left(\text{KL}\left(p_{\text{real}}\left|\right|p_{\text{mix}}\right)+\text{KL}\left(p_G\left|\right|p_{\text{mix}}\right)\right)\mathrm{,}$ (27) with $p_{\text{mix}}=\left(p_{\text{real}}+p_G\right)/2$. Thus, the GAN is an implicit MLE-based generative model. The optimally trained GAN converges in the Nash equilibrium state, where the generator excels in synthesizing samples that the discriminator cannot differentiate from the real data; thus, the generator-induced distribution indeed captures the real data distribution after training. This technique has been utilized in various scientific contexts, e.g., in condensed-matter physics [22, 23], particle physics [24, 25], cosmology [26, 27], and QFT study with lattice simulation [28, 29].

Principal Component Analysis

In ML, principal component analysis (PCA) is a statistical technique that involves transforming a set of correlated variables into independent variables through orthogonal transformations. The principal components, which are associated with the obtained main eigenvectors (or non-negligible singular values), reveal the most representative configurations of the data. As one of the unsupervised learning techniques, PCA implements singular value decomposition (SVD) on a real matrix [42]: $M=X\Sigma Z=VZ,$ (30) where M is a matrix of size N×m; X and Z are two orthogonal matrices of size N×N and m × m, respectively; and ∑ is a diagonal matrix with the singular values arranged in descending order. Then, the ith row of the matrix M⁽ⁱ⁾ can be expressed as $\begin{array}{c}{M}^{(i)}={\displaystyle \sum _{j=1}^m{x}_j^{(i)}{\sigma }_j{z}_j}={\displaystyle \sum _{j=1}^m{\tilde{v}}_j^{(i)}{z}_j}\\ \approx {\displaystyle \sum _{j=1}^k{\tilde{v}}_j^{(i)}{z}_j}(i)=\mathrm{1,...,}N\end{array}$ (31) where ${\tilde{v}}_j^{(i)}$ is the corresponding coefficient of zj for the ith row. In the last step, there is a cut on the indices, i.e., k, because PCA focuses n the most important components. Owing to its effectiveness for data mining, PCA has been widely used in various areas of physics research. For recent progress in HICs, please see Sect.7.

Impact parameter estimation

The impact parameter b describes the distance between the centers of the two colliding nuclei in the classical view, which is a crucial quantity determining the initial geometry of a collision. In experiments, the impact parameter is not directly measurable and is usually estimated from the multiplicity of final-state particles in track detectors or the energy deposited in calorimeters. ML approaches are proposed to determine the impact parameters from the final-state particles and exhibit better performance than conventional methods. Ref.[44] proposed the use of a DNN and CNN to reconstruct the impact parameters from the energy spectra of final-state charged hadrons of HICs at $\sqrt{s_{\text{NN}}}$ = 7.7 to 200 GeV, which were simulated with a multiphase transport (AMPT) model. Both the DNN and CNN can reconstruct the impact parameters with an MAE of approximately 0.4 fm. When the input feature is from a larger pseudorapidity window, the CNN has a higher prediction accuracy than the DNN. Ref.[45] reported the performance of a CNN and Light Gradient Boosting Machine (LightGBM) in reconstructing the impact parameter from HICs at $\sqrt{s_{\text{NN}}}$ = 0.2 to 1 GeV, which were simulated with the Ultra-relativistic Quantum Molecular Dynamics (UrQMD) model. The input features are constructed from the proton spectra for transverse momentum and rapidity. The average difference between the true impact parameter and the estimated one can be less than 0.1 fm. LightGBM has better performance than the CNN.

A model-independent Bayesian inference method for reconstructing the impact parameter distributions was proposed in Ref.[46]. The impact parameter distributions are inferred from model-independent data. This method is based on Bayes’ theorem: $P\left(b|X\right)=P\left(b\right)P\left(X|b\right)/P\left(X\right)\mathrm{,}$ (32) where P (X) represents the probability of the observable that can be measured in the experiment, and $P\left(X|b\right)$ represents the probability density distribution of X for a given impact parameter b. Fluctuation is taken into account by assuming $P\left(X|b\right)$ to be a Gaussian or gamma distribution, which can be determined by fitting the data with the formula $P\left(X\right)={\displaystyle \int }P\left(X|b\right)P\left(b\right)\text{d}b$. $P\left(X\right)$ can be multidimensional. In Ref.[46] two observables were used: $X=\left\{M\mathrm{,}{p}_{\text{t}}^{\text{tot}}\right\}$, where M represents the multiplicity of the charged particles and $p_{\text{t}}^{\text{tot}}$ represents the total transverse momentum of the light particles.

End-to-end centrality estimation for CBM

The compressed baryonic matter (CBM) detector is currently under construction for the Facility for Antiproton and Ion Research (FAIR) at Gesellschaft für Schwerionenforschung (GSI), which will study the properties of strongly compressed nuclear matter via HICs with beam energies ranging from 2 to 10 AGeV. A characteristic of the CBM experiment is its very high event rate and trigger rate, which will produce a large amount of raw data per second in real-time and pose a challenge for online event characterization and storage. To address the online event characterization, it is essential to be able to work on the direct output of the detector, which has an inherent point cloud structure—a collection of points as an unordered list with particles or tracks’ attributes recorded. One important property of the point cloud is that they as a whole should be invariant under permutation. The PointNet structure [47] was specially developed to respect this order invariance. Accordingly, for HICs, PointNet-based models can perform real-time physics analysis on the detector output directly.

Refs. [48, 19] proposed the use of PointNet-based models for event-by-event impact parameter determination for the CBM experiment using the direct output from the detector, where the trained model serves as an end-to-end centrality estimator. The supervised learning strategy is used for this regression task, where the training data are prepared from UrQMD followed by CBMRoot detector simulation to obtain the detector output, which are hits or tracks of the particles. A PointNet-based model is constructed and trained to capture the inverse mapping between the detector output and the impact parameter information. It was shown that PointNet-based models can perform accurate event-by-event impact parameter determination using hits of charged particles in different detector planes and/or the tracks reconstructed from these hits. With regard to both precision and accuracy, these models outperformed a baseline model using charged track multiplicity as the input inside a polynomial fit. While the baseline model had a similar resolution (relative precision) to the PointNet-based model in the semi-central collision region, it had a lower accuracy and more fluctuations in the accuracy for impact parameters ranging from 3 to 16 fm, as indicated by the mean of the prediction error for the impact parameter. This trend was more evident for a realistic event distribution (i.e., ~ bdb), as shown in Fig. 4 for the mean prediction error. Considering the natural parallelizability and high speed, the PointNet-based model paves the way for real-time end-to-end event characterization for HIC studies.

Fig. 4Full-screen View

(Color online) Taken from Ref. [48]. Mean error in predictions as a function of centrality. Dataset Test2 is used, in which peripheral events are more likely to occur than other centralities. The track multiplicity is used for the centrality binning. The points at 90% centrality are results from events with no tracks reconstructed. Therefore, the Polyfit and MS-Tracks models do not have a data point at 90% centrality.

Nuclear deformation estimation

The momentum distribution of final-state hadrons is sensitive to nuclear shape deformation. For example, owing to the different collision geometry, the elliptic flow as a function of charge multiplicity differs significantly between Pb+Pb and U+U collisions. As shown in Fig. 5, the ²⁰⁸Pb is a doubly magic nucleus with an almost perfectly spherical shape, whose collision patterns depend only on the impact parameter b. In contrast, the shape of ²³⁸U is similar to a watermelon, and the corresponding collision patterns are far more complex than those of Pb+Pb collisions. For example, the U+U collisions have body–body aligned, body–body crossed, tip–tip, and tip–body collisions. Different collision patterns correspond to different charge multiplicities and elliptic flows. Both the fully overlapped body–body aligned and central tip–tip collisions correspond to most central collisions with high charge multiplicity, but their elliptic flows differ significantly. This type of difference leads to a far larger variance in the elliptic flow for most central U+U collisions, compared with high-multiplicity Pb+Pb collisions. In principle, the complex collision patterns lead to many differences in the elliptic flow compared with the charged multiplicity diagram. Deep learning can be used to identify these differences and predict the nuclear shape deformation parameters using these patterns.

Fig. 5Full-screen View

(Color online) Collision geometries for Pb+Pb and U+U collisions.

It was demonstrated that by using nuclei with different deformation parameters β₂ and β₄, high-energy HICs can be simulated using the TRENTo Monte Carlo model to obtain the event-by-event total initial entropy (which is proportional to the final charged multiplicity) and the corresponding geometric eccentricity (which is approximately proportional to the elliptic flow). A deep residual neural network was trained to predict β₂ and β₄ using the two-dimensional (2D) images of total entropy vs. eccentricity [49]. The network accurately predicted the absolute values of β₂ and β₄ but failed to predict their signs using the information provided. Using the class activation map (CAM) method to map the last convolutional layer onto the input image, the authors found two regions in the image that are important for decision-making. One is the most central collision region, which is the most sensitive region to the variance of the elliptical flow.

Recently, Bayesian inference with a Gaussian process (GP) emulator was used for reconstructing the nuclear structure including deformation parameters based on HIC measurements [50]. As a first-step exploratory study, the collision observables (charged multiplicities Nch, elliptic flow v₂, triangular flow v₃, and mean transverse momentum $\langle p_T\rangle$) were estimated via Monte Carlo Glauber model calculation (total energy E, elliptic eccentricity ϵ₂, triangular eccentricity ϵ₃, and energy density $d_{\perp }$), can reasonably estimate the ratio of observables in isobaric collision systems owing to the cancellation of dynamics’ uncertainties [51]. Under this setup, nuclear structure reconstruction based on both single collision system and contrast isobaric collision system measurements were discussed. For single collision systems, it was found that the Woods–Saxon parameters of nuclei can be precisely inferred from final-state observables estimated with (P,ϵ₂,ϵ₃, $d_{\perp }$). For isobar collision systems, the simultaneous inference on the two set of nuclear structures fails with only the ratio of these final observables, while the further provision of the single collision system’s multiplicity distribution allows high-precision nuclear structure reconstruction. Additionally, the ratio of radial flow was found to be redundant in the presence of the ratio of elliptic flow and vise versa.

α-clustering structure

The clustering structure is an exotic phenomenon in nuclei, and it usually occurs in light nuclei [52]. In nuclear collisions between light clustering nuclei and heavy ions, the clustering structure can make the final-state particles anisotropically distributed [43, 53, 54]. It is crucial to extract the quantitative information about the clustering from the final observables.In the ¹²C / ¹⁶O + ¹⁹⁷Au collisions at relativistic energies, an ML method was used to obtain evidence of the cluster structures from the azimuthal angle and transverse momentum distributions of charged pions [55]. In this study, a Bayesian convolutional neural network (BCNN) was used. In addition to the input layer and output layer, there were hidden layers, consisting of four convolutional layers and three fully connected layers. The parameters of the three fully connected layers were sampled from distributions learned via Bayesian inference. A 2D histogram of azimuthal angle vs. transverse momentum was used as the input. Considering the detection efficiency in the experiments, charged pions with rapidity ranging from –1 to 1 and transverse momentum ranging from 0 to 2 GeV/c were selected. The dataset consisted of 1.6 × 10⁶ histograms with 64 × 64 bins (pixels), with different labels to indicate different configurations.

The typical spectra of 4000 merged events are shown in Fig. 6. Even with merging, the samples of different configurations are barely distinguishable to the naked eye. The number of merged events is denoted as NEvent, which is taken to be 1000, 2000, and 4000.

Fig. 6Full-screen View

(Color online) Taken from Ref. [55]. Two-dimensional azimuthal angle vs transverse momentum distributions of charged pions for non-clustered (Up) and clustered (Down) ¹²C from an AMPT-generated ¹²C+¹⁹⁷Au collision event at $\sqrt{S_{\text{NN}}}$ =200 GeV.

The learning curves are shown in Fig. 7. As more events were merged, the event-by-event fluctuations were reduced, and the network was able to learn the features of the final state for predicting the initial configuration. For ¹²C with NEvent = 4000 and ¹⁶O with NEvent = 2000, the validation accuracy reached 95% and 97%, respectively, and for ¹⁶O with NEvent = 4000, it reached 99%.

Fig. 7Full-screen View

(Color online) Taken from Ref. [55]. Validation accuracy during the training process for colliding systems ¹²C/¹⁶O+¹⁹⁷Au with NEvent = 1000, 2000, and 4000.

For the clustering phenomenon, it is extremely difficult to extract signals from the final particles, because fluctuations play such an important role in relativistic HICs. By averaging over multiple events, the BCNN model can learn the features with good performance.

Neutron skin estimation

The distribution of neutrons is important in determining the thickness of the neutron skin, the symmetry energy of the nucleus, the QCD equation of state (EoS) of dense nuclear matter, and astrophysical observables such as the mass–radius relationship of neutron stars and the gravitational wave emitted during neutron star mergers. However, extracting the distribution of neutrons inside the nucleus is extremely difficult. The distribution of neutrons inside the nucleus differs from the distribution of protons. The proton distribution is far easier to measure than the neutron distribution because the former is equivalent to the charge distribution, whereas the latter is associated with the weak charge distribution. The neutron skin, which is the difference between the root-mean-square radii of neutrons and protons, can be used to determine the neutron (weak charge) distribution in the nucleus. PREX2 measured the parity-violating asymmetry by scattering longitudinally polarized electrons on Pb208 to obtain a neutron skin thickness of approximately $R_{\text{n}}-R_{\text{p}}=0.283\pm 0.071$ fm [56]. The neutron skin is used as a constraint in the calculation of the positive and negative correlations between the symmetry energy and the slope parameter at the saturation density. With this constraint, the Bayesian analysis achieves a compromise between the “conflicting” data that lead to the famous “PREXII puzzle” and the “soft Tin puzzle” [57, 58].

There have been many attempts to determine the neutron skin thickness and the symmetry energy at low energy [59], e.g., by investigating the charge-exchange spin–dipole excitation [60], the supernova neutrinos [61], nuclear fragmentation reactions [62], and parity-violating electron scattering [56, 63].

For high-energy HICs, it was proposed that the isobar ratios of the charge multiplicities of the mean transverse momentum and the net charge multiplicities between ${}_{44}^{96}\text{Ru}{+}_{44}^{96}\text{Ru}$ and ${}_{40}^{96}\text{Zr}{+}_{40}^{96}\text{Zr}$ can be used to precisely determine the nucleon skin and the symmetry energy [64]. The authors claimed that the high-energy isobar collisions can significantly improve the result of the traditional low-energy method. In another paper, the yields of spectator protons and neutrons at the forward velocity of ultra-central collisions were proposed to be good probes of the neutron skin—sensitive to the neutron skin of ²⁰⁸Pb but insensitive to other parameters during the collision [65]. A more accurate method is to measure the free spectator neutron yield ratios between ${}_{44}^{96}\text{Ru}{+}_{44}^{96}\text{Ru}$ and ${}_{40}^{96}\text{Zr}{+}_{40}^{96}\text{Zr}$ in ultra-central collisions [66].

A large amount of data has already been collected from high-energy HICs. There may be a data-driven way to reuse these data to determine the neutron distribution and neutron skin thickness. It has been tested in [67] that nucleons sampled from nuclei with different neutron skin types can be classified with reasonable accuracy using deep CNNs and point cloud networks. However, once the nucleus is involved in HICs, it is almost impossible to distinguish the neutron skin types of the colliding nucleus using the momentum distribution of the final-state hadrons. For this task, the signal is weak in minimum bias collisions, and DNNs fail to solve the difficult inverse problem. A new ML method is needed to search for weak signals in data with large statistical fluctuations.

Shear and bulk viscosities

The shear and bulk viscosities are important properties that significantly affect the dynamical expansion of QGP and the momentum distribution of final-state hadrons, as indicated by relativistic fluid dynamics simulations [68-71]. In solving the inverse problem of HICs, it was found that the effects of viscosity are entangled with the initial thermalization time, the EoS of QGP, and the phase transition between QGP and HRG. Thus, determining the shear and bulk viscosities of hot nuclear matter is a notoriously difficult problem. Regarding the nucleonic degree of freedom, the shear viscosity has attracted considerable attention because it is related to the nuclear EoS, phase change, and strong interaction [72-75]. A similar feature to the QGP viscosity has been demonstrated for the behavior η/s(T). Bayesian analysis plays an important role in determining the temperature dependence of the ratio of shear viscosity over the entropy density ratio η/s(T) as well as the bulk viscosity over the entropy density ratio ζ/s(T) [76-78].

Suppose that all the parameters in the theoretical model of HICs form a set $\left\{\theta \right\}$ and all the experimental data from RHIC and LHC form another set $\left\{D\right\}$. Then, the posterior distribution of the model parameters is given by $P\left({\theta }_i|D\right)=\frac{P\left(D|{\theta }_i\right)P({\theta }_i)}{P(D)}=\frac{P\left(D|{\theta }_i\right)P({\theta }_i)}{{\displaystyle {\sum }_{j}P}\left(D|{\theta }_j\right)P({\theta }_j)},$ (33) where $P(D|{\theta }_i)$ represents the likelihood between the experimental data D and model output using parameter combinations θi; $P({\theta }_i)$ represents the a priori distribution of θi, which may be a belief based on past experience or physical considerations; and the denominator $P(D)={\displaystyle {\sum }_{j}P}\left(D|{\theta }_j\right)P({\theta }_j)$ is a normalization factor called evidence. Computing P(D) is too expensive because it requires the theoretical model to traverse the entire parameter space. Fortunately, in Bayesian analysis, the normalization factor is not needed, because the Markov chain Monte Carlo (MCMC) method can sample from the following un-normalized distributions: $P\left({\theta }_i|D\right)\propto P\left(D|{\theta }_i\right)P({\theta }_i),$ (34) The final output of the Bayesian analysis is a large number of different combinations of model parameters sampled from the above un-normalized posterior distribution function. Performing a density estimation for each parameter, e.g., the slope of η/s over Tc, gives a distribution (or histogram) of the slope parameter. The location of the maximum value in this distribution corresponds to the MAP estimate. Additionally, the distribution has a variance that corresponds to the uncertainty in the slope parameter, which comes from the experimental data, the prior distribution, and the likelihood function. Thus, it is clear that the extracted model parameters are well constrained when their posterior distribution has a narrow peak.

To estimate the temperature dependence of the shear and bulk viscosities, two parameterized functions based on physical a priori are required. In a Nature Physics paper [77], the shear and bulk viscosities were parameterized as follows: $\left(\eta /s\right)(T)={\left(\eta /s\right)}_{\min }+{\left(\eta /s\right)}_{\text{slope}}\left(T-T_{\text{c}}\right){\left(\frac{T}{T_{\text{c}}}\right)}^{{(\eta /s)}_{\text{crv}}},$ (35) $\left(\zeta /s\right)(T)=\frac{{\left(\zeta /s\right)}_{\max }}{1+{\left(\frac{T-{\left(\zeta /s\right)}_{T_{\text{peak}}}}{{\left(\zeta /s\right)}_{\text{width}}}\right)}^2},$ (36) where (η / s)_min and (ζ / s)_max represent the minimum and maximum shear and bulk viscosity values to be determined, respectively, Tc=154 MeV is the QCD transition temperature representing the location of the minimum in η/s(T), and ${(\zeta /s)}_{T_{\text{peak}}}$ denotes the location of the maximum bulk viscosity to be determined. Other parameters to be determined are the slope (η / s)_slope and the curvature of the shear viscosity (η / s)_crv and the width of the bulk viscosity peak (ζ / s)_width.

Without considering other parameters, these six parameters form a six-dimensional parameter space. The above Bayes formulae are used to traverse this space, with the trajectories forming a set of parameter combinations. This is equivalent to importance sampling using the posterior distribution of the six parameters. Density estimation indicates that the distribution of (η / s)_min is approximately normal, whose mean and variance give a quantitative estimate of ${(\eta /s)}_{\min }{=0.085}_{-0.025}^{+0.026}$. An anti-correlation is observed between (ζ / s)_max and (ζ / s)_width, indicating that it is the integral of (ζ / s)(T) that matters, not its specific form. The analysis also indicates that the experimental data used can not constrain the parameters (η / s)_crv and ${\left(\zeta /s\right)}_{T_{\text{peak}}}$, as there are no obvious peaks in the posterior distributions of these two parameters.

Crossover or first-order phase transition

In general, as mentioned in the Introduction, the challenge faced by high-energy nuclear collision studies can essentially be viewed as an inverse problem. Assuming that all related physical factors (e.g., initial condition/fluctuations, QGP bulk properties, transport coefficients, freeze-out parameter, hadronic interactions) are given, well-established theoretical models (e.g., relativistic viscous hydrodynamics with hadronic transport simulation) can be adopted to simulate the HIC process to give their final-state observables, and such a forward process is well understood. However, given instead only limited measurements of the final state of HICs, it is unclear how to disentangle those different influencing physical factors for decoding the corresponding early time dynamics. For high-energy HICs, there are two strategies for solving this inverse problem using statistical methods and ML: one is Bayesian inference with the task of parameter estimation for calibrating the chosen model (e.g., in Ref. [79]), and the other is supervised ML for directly capturing the inverse mapping from the final state to the corresponding physics of interest.

Ref. [80] proposed the use of a deep CNN to capture the direct inverse mapping from the final-state information to the types of QCD transition happened in early time. This is inspired by the success of image recognition in computer vision. Although the inverse mapping may be very implicit, DNNs can be used to decode it and represent it in the sense of Big Data in a supervised manner. The required training data can be prepared through well-established model simulation for HICs, e.g., using the state-of-the-art 3+1-dimensional viscous hydrodynamics [81-84], where diversity can be introduced by varying different physical factors (i.e., parameters in the simulation). As an exploratory study, a binary classification task was targeted, where the Deep CNN was trained to identify the QCD transition type embedded within the collision dynamics as crossover or first-order solely according to the final pion spectra $\rho \left(p_T\mathrm{,}\varphi \right)$, as shown in Fig. 8. The EoS of the hot and dense matter is a crucial ingredient in the hydrodynamic simulations. Embedded in it is the nature of the QCD transition (first-order or crossover), which can significantly affect the hydrodynamic evolution according to the shape of the pressure gradient. As the input to the deep CNN, the final charged pion’s spectra at mid-rapidity are obtained using the Cooper–Frye formula in each hydrodynamic simulation: $\rho \left(p_T\mathrm{,}\varphi \right)=\frac{d{N}_i}{\text{d}Y{p}_T\text{d}{p}_T\text{d}\varphi }=g_i{\displaystyle {\int }_{\sigma }{p}^{\mu }}\text{d}{\sigma }_{\mu }{f}_i\left(p\cdot u\right),$ (37) where Ni represents the particle number density, Y represents the rapidity, gi represents the degeneracy, $d{\sigma }_{\mu }$ represents the freeze-out hypersurface element, and fi represents the thermal distribution. The training dataset of $\rho (p_T\mathrm{,}\varphi )$ was generated from the event-by-event hydrodynamic package CLVisc [81] with fluctuating AMPT initial conditions, with which supervised learning using the CNN is performed for binary classification in identifying the QCD transition types.

Fig. 8Full-screen View

(Color online) Schematic of QCD transition classification with HIC final particle spectra.

Fig. 9 shows the space–time evolution histories of QGP expansion starting from the same initial condition model with different fluctuations, in relativistic hydrodynamic simulations using CLVisc. For EOSQ with a first-order phase transition, the pressure gradient is zero in the mixed phase. Multiple ridge structures are formed with a first-order phase transition in the EoS because the expansion of QGP is mainly driven by the pressure gradient and the acceleration is 0 in the mixed phase. However, the expansion histories are significantly different when the shear viscosity is not 0. Different evolution histories lead to different final-state particle spectra in the momentum space.

Fig. 9Full-screen View

Evolution history of QGP simulated using the relativistic hydrodynamic model CLVisc, starting from the same initial condition with four different parameter combinations. From top to bottom, each row presents four snapshots taken at different times, using different combinations of the EoS and shear viscosity over the entropy density ratio. Here, EOSL represents the lattice QCD EoS with a crossover transition between QGP and hadron resonance gas, and EOSQ represents an EoS with a first-order phase transition between QGP and hadron resonance gas.

To verify the robustness of the trained deep CNN in this QCD EoS recognition task, the test set was simulated from a different hydrodynamics package or with different initial fluctuating conditions (IP-Glasma or MC-Glauber) and different $\eta /s$ parameters. The conventional observables, such as the elliptic flows v₂ and integrated particle spectrum, were shown to be insufficient for distinguishing the two QCD transition classes for these test set, whereas the trained deep CNN achieved an average classification accuracy of 95%, indicating that it was robust against contamination from factors such as the initial fluctuations and shear viscosity. For comparison, the best classification accuracy among traditional ML algorithms such as decision tree, random forest, support vector machine, and gradient boosting was approximately 80%. The good performance of the trained deep CNN indicates that the imprint of the early time transition dynamics is not fully washed out by the collision evolution and is still embedded in the final-state information. Additionally, the inverse mapping from final-state observables to the QCD transition information can be well captured by the deep CNN from the supervised training strategy, providing a discriminative and traceable encoder for the dynamical information of QCD transition. Thus, the constructed deep CNN functions as an “EoS-meter" to efficiently bridge the HIC experiments to QCD bulk matter physics. This study paved a path to the success of experimental research on the QCD EoS and the search for the critical endpoint of the QCD phase diagram. In the study, the afterburner hadronic cascade effects were not considered; thus, the conclusion regarding the direct inverse mapping was drawn from the viewpoint of pure hydrodynamic evolution.

Later, this strategy was deepened in a series of studies for more realistic scenarios, e.g., to take into account the afterburner hadronic cascade by incorporating UrQMD following the hydrodynamics evolution [85, 86]; to consider non-equilibrium dynamics of the phase transition’s influence, e.g., spinodal decomposition [18, 87] or Langevin dynamics [88]; to include more realistic experimental detector effects through detector simulation with hits or tracks as the input [48, 89]; to perform unsupervised outlier detection for HICs [90]; and to determine the nuclear symmetry energy [91]. Specifically, in Ref. [89] it was shown that by using just the detector output directly, PointNet models can be employed to classify collision events simulated by an EoS associated with a first-order phase transition and those simulated by an EoS with a crossover transition. The PointNet models take the reconstructed tracks from the CBM detector (simulated with CBMRoot) followed by the hybrid UrQMD events. They achieved a binary classification accuracy of approximately 96% when trained on collision events for impact parameters ranging from 0 to 7 fm. When the model training set was shrunk to the mid-central region with b=0~3 fm, the model accuracy increased to approximately 99%. A combination of training sets from both peripheral and mid-central collisions resulted in a classifier being able to identify the phase transition type across different centralities, while not compromising the accuracy for the central region.

Active learning for QCD EoS

First-principles calculations using lattice QCD provide the EoS of hot nuclear matter at high temperatures and zero baryon chemical potential. Because of the fermionic sign problem, lattice QCD fails to compute the nuclear EoS at finite μB at present. Using Taylor expansion, it is possible to obtain the nuclear EoS at a small μB that is close to zero, approximately. The BEST collaboration formulated a nuclear EoS with a critical endpoint by mapping the 3D Ising model with the Tylor expansion result. However, the model contains four free parameters whose values determine the size and location of the critical endpoint. Some combinations of these parameters lead to an unphysical, e.g., acausal or unstable, EoS.

Supervised learning can help to map unphysical regions of parameter combinations. However, labeling is computationally expensive in this task. For thermodynamic stability, one must check the positivity of the energy density, pressure, entropy density, baryon density, second-order baryon susceptibility ${\chi }_2^B$, and heat capacity ${\left(\partial S/\partial T\right)}_{n_B}$, as well as the causality condition: $0\le c_s^2\le 1,$ (38) where cs represents the speed of sound in hot nuclear matter.

Active learning was used to find the most informative parameter combinations before labeling them [92]. In active learning, the network is first trained using a small amount of labeled data. Then, the trained network is employed to make predictions on all samples from a large unsupervised pool. If the network is uncertain about one parameter combination, e.g., it predicts that this group of parameter combinations will lead to an EoS that is unphysical with probability 51%, this sample lives on the decision boundary and should be informative and important for the network. Labeling this sample will improve the performance of the network more than labeling easy samples. The newly labeled sample is moved out of the pool and will be used in supervised learning later.

Accelerated relativistic hydrodynamic simulation via deep learning

Relativistic hydrodynamics is a powerful tool for simulating the QGP expansion and studying the flow observables in relativistic HICs at the RHIC and LHC energies [95-100]. For ideal hydrodynamics with zero net charge densities, it solves the transport equations of the energy momentum tensor: ${\partial }_{\mu }{T}^{\mu \nu }=0,$ (39) where $T^{\mu \nu }=\left(e+p\right)u^{\mu }{u}^{\nu }-p{g}^{\mu \nu }$, e represents the energy density, p represents the pressure, and $u^{\mu }$ represents the four-velocity. In traditional hydrodynamic simulations, these transport equations are numerically solved with an algorithm such as SHASTA or LCPFCT that transforms the initial conditions into final-state profiles through nonlinear evolutions [97, 100-102].

Recently [93, 94], a DNN called stacked U-net (sU-net) was designed and trained to learn the initial- and final-state mapping from the nonlinear hydrodynamic evolution. The constructed sU-net has an encoder–decoder architecture, which contains four U-net blocks with residual connections between them. For each U-net block, there are three convolutional and deconvolutional layers with Leaky ReLU and softplus activation functions employed for the inner and output layers, respectively. By concatenating the feature maps along the channel dimension, the output of the first two convolutional layers is fed to the last two deconvolution layers. For details, please refer to [93, 94].

The training and test data (the profiles of the initial and final energy momentum tensor $T^{\tau \tau }$, $T^{\tau x}$, $T^{\tau y}$) were generated from VISH2+1 hydrodynamics [103, 104] with zero viscosity, zero net baryon density, and longitudinal boost invariance. In more detail, sU-net was trained with 10000 initial and final profiles from VISH2+1 with MC-Glauber initial conditions [105, 106], and then its prediction accuracy was tested using the profiles of four different types of initial conditions: MC-Glauber [105, 106], MC-KLN [106, 107], AMPT [81, 108, 109], and TRENTo [110]. Fig. 10 presents the final energy density and flow velocity predicted by sU-net, together with a comparison with the hydrodynamic results. As shown, the trained sU-net captured the magnitudes and structures of both the energy density and the flow velocity. In particular, panels (b), (d), and (f) show that the network, which was trained with datasets generated with MC-Glauber initial conditions, was also capable of predicting the final profiles of other types of initial conditions. In Refs. [93, 94], the eccentricity coefficients, which indicate the deformation and inhomogeneity of a large number of the energy density profiles, were calculated, and the predictions from sU-net almost overlapped with the results from VISH2+1.

Fig. 10Full-screen View

(Color online) Energy density and flow velocity profiles predicted by sU-net and calculated from VISH2+1 for six test initial profiles of MC-Glauber, MC-KLN, AMPT, and TRENTo. Taken from Refs. [93, 94].

Compared with the 10∼20 minutes simulation time of VISH2+1 on a traditional CPU, sU-net took several seconds to directly generate the final profiles for different types of initial conditions on one P40 GPU, which significantly accelerated the traditional hydrodynamic simulations. However, the sU-net model designed and trained in Refs. [93, 94] mainly focuses on mimicking the 2+1-dimensional hydrodynamic evolution with a fixed evolution time. For more realistic implementation, it is important to explore the possibilities of mapping the initial profiles to the final profiles of the particles emitted on the freeze-out surface of the relativistic HICs.

Spectral function reconstruction

Accessing real-time properties of QCD (or a many-body system in general) remains a notoriously difficult problem, because the non-perturbative computations, such as lattice field simulations or functional methods, usually operate in Euclidean space–time (after a Wick rotation $t\to it\equiv \tau$) and thus can only provide Euclidean correlators (i.e., in imaginary time). Thus, the analytic continuation of these discrete noisy data is often ill-posed. Quantitatively understanding the real-time dynamics determined by the Minkowski correlator is important and interesting, e.g., for understanding scattering processes, transport, or non-equilibrium phenomena that occur in HICs. The Minkowski correlator is usually accessed from the Euclidean correlator via spectral reconstruction.

The associated ill-posed problem can be cast as a Fredholm equation of the first kind: $g(t)={\displaystyle {\int }_a^{b}K(t\mathrm{,}s)\rho (s)\text{d}s}\mathrm{,}$ (40) with the goal of retrieving the function ρ(s) given the kernel function K(t, s) but limited information about g(t). It has been well shown that the required inverse transform becomes ill-conditioned if only a finite set of data points with non-vanishing uncertainty are available for g(t). In the context of QFT, one can simply approach this problem via the Källén–Lehmann spectral representation of the correlators, taking the kernel function to be $K(t\mathrm{,}s)=s{\left(s^2+t^2\right)}^{-1}{\pi }^{-1}\mathrm{,}$ (41) where the ρ(s) functions involved are usually called spectral functions. The task of reconstructing the spectral function from the correlator measurements (from the lattice calculation) needs to be regularized to make sense of the inverse problem involved. Over the past few decades, many different regularization techniques have been explored for this ill-conditioned inverse problem, such as Tikhonov regularization, maximum entropy methods, and Bayesian inference techniques.

Recently, deep learning-based strategies have also been explored to tackle spectral reconstruction, which can be mainly categorized into two schemes: data-driven supervised learning approaches and unsupervised learning-based approaches. The first application of domain-knowledge-free deep-learning methods to this ill-conditioned spectral reconstruction (also called analytic continuation) was reported in Ref. [111] in the context of general quantum many-body physics. The results indicated the good performance of DNNs with supervised training in the cases of a Mott–Hubbard insulator and a metallic spectrum. In particular, a CNN was found to achieve better reconstruction than a fully connected network, with performance superior to that of the MEM—one of the most widely used conventional methods. In Ref. [112], the authors adopted a similar strategy but also introduced PCA to reduce the dimensionality of the QMC-simulated imaginary time correlation function of the position operator for a harmonic oscillator linearly coupled to an ideal heat bath.

The authors of Ref. [113] also adopted a data-driven perspective. They adopted a strategy similar to spectral function reconstruction in the QFT context and considered the Källen–Lehmann spectral representation as the accessible propagator, i.e., $G(p)={\displaystyle {\int }_0^{\infty }\frac{d\omega }{\pi }\frac{\omega \rho (\omega )}{{\omega }^2+p^2}}$, which takes the kernel in the Fredholm equation as the Källen–Lehmann kernel. For the dummy spectral functions, the superposition of Breit–Wigner peaks was used, according to the perturbative one-loop QFT-derived parameterization ${\rho }^{\text{BW}}(\omega )=4A\Gamma \omega /\left({\left(M^2+{\Gamma }^2-{\omega }^2\right)}^2+4{\Gamma }^2{\omega }^2\right)$. Two types of DNNs have been studied, both with a noisy propagator as the input, but with different outputs: one estimates the parameters (e.g., Γi and Mi for the collection of Breit–Wigner peaks) of the spectral function (denoted as PaNet), and the other attempts to directly reconstruct the discretized data points of the spectral function (denoted as PoNet).

As another type of non-parametric representation, GPs were used in the reconstruction of the 2+1 flavor QCD ghost and gluon spectral function in Ref. [115]. In general, the GP can define a probability distribution over families of functions, which is typically characterized by the chosen kernel function. In Ref. [115] the GP was assumed to describe the spectral function: $\rho (\omega )\sim \mathcal{G}\mathcal{P}(\mu (\omega )\mathrm{,}C(\omega \mathrm{,}{\omega }^{\prime }))\mathrm{,}$ (42) where the mean function $\mu (\omega )$ is often set to zero, and the covariance $C(\omega \mathrm{,}{\omega }^{\prime })$ is determined by the kernel function used, for which a common standard choice is the radial basis function (RBF) kernel $C\left(\omega \mathrm{,}{\omega }^{\prime }\right)={\sigma }_C{e}^{-\frac{{\left(\omega -{\omega }^{\prime }\right)}^2}{2l^2}}\mathrm{,}$ (43) with tunable hyperparameters σ_C for the overall magnitude and l for the length scale. The prior represented by this GP can be plugged into the Bayesian inference procedure with lattice data for the ghost dressing function and gluon propagator for evaluating the likelihood. In Ref. [115], the lattice data were specifically extended. The ghost dressing function was extended to the deep infrared range, and the low-frequency behavior was constrained by spectral DSE results [116]. The gluon dressing function was extended to the ultraviolet range with previous fRG computation results [117]. This reduced the variance in the solution space and enhanced the stability compared with the inference without such extensions. It was shown that while approximately fulfilling the Oehme–Zimmermann superconvengence (OZS) condition for gluons, the reconstruction with GP regression in this work accurately reproduced the lattice data within the uncertainties with deviations for a gluon propagator stronger in some regions than those for the ghost dressing function. For the spectral function, the reconstruction exhibited a similar peak structure to a previous fRG reconstruction of the Yang–Mills propagator [117].

In Refs. [114, 118, 119], the authors developed an unsupervised approach based on DNN representation for the spectral function together with automatic differentiation (AD) to reconstruct the spectral function, which does not need training data preparation for supervision (a similar DNN-based inverse problem solving strategy within the AD framework was used for reconstructing the neutron-star EoS from astrophysical observables [120, 121] and inferring the parton distribution function of pions in lattice QCD studies [122]). The introduced DNN representation can preserve the smoothness of the spectral function automatically, helping to regularize the degeneracy issue in this inverse problem. This is because, as analyzed in Ref. [119], the degeneracy is related to the null modes of the investigated kernel function, which usually induce oscillation for the reconstructed spectral function. Specifically, the DNN-represented spectral function, i.e., $\overrightarrow{\rho }=[{\rho }_1\mathrm{,}{\rho }_2\mathrm{,...,}{\rho }_{N_{\omega }}]$, can be converted into the propagator under the discretization scheme as $D(p)={\displaystyle {\sum }_i^{N_{\omega }}\left(p\mathrm{,}{\omega }_i\right)}K{\rho }_i\Delta \omega$. Then, the loss function over the propagator relative to lattice data, i.e., $\mathcal{L}={\displaystyle {\sum }_i^{N_p}{\left(D_i-D(p_i)\right)}^2/{\sigma }_i}$, can be evaluated and provide guidance for tuning over the DNN-represented spectral function. Taking gradient-based algorithms, the derivative of the loss with respect to network parameters can be derived as ${\nabla }_{\theta }\mathcal{L}={\displaystyle \sum _{j\mathrm{,}k}K}\left(p_j\mathrm{,}{\omega }_k\right)\frac{\partial \mathcal{L}}{\partial D(p_j)}{\nabla }_{\theta }{\rho }_k\mathrm{.}$ (44) ${\nabla }_{\theta }{\rho }_k$ is computed easily under standard backward propagation for the network.

For the DNN representation of the spectral function, two different schemes were investigated in this work: one uses the multiple outputs of an L-layer neural network to represent in list format the spectral function (denoted as NN), and the other directly uses a feedforward neural network for parameterization (denoted as NN-P2P) of the spectral function as a function of frequency, i.e., ρ(ω). For the training, the Adam optimizer is adopted, and the L₂ regularization is set in the warm-up beginning stage under an annealing strategy until the regularization strength value is sufficiently small (set as < 10^-8 in the calculation). This can relax the regularization to obtain hyperparameter-independent inference results. For the direct NN list representation, a quenched implementation of smoothness condition ${\lambda }_s{{\displaystyle {\sum }_{i=1}^{N_{\omega }}\left({\rho }_i-{\rho }_{i-1}\right)}}^2$ is also performed with λs reduced from 10^-2 to 0. This unsupervised spectral reconstruction method was validated with regard to the uniqueness of the solutions both analytically and numerically [119]. As shown in Fig. 11, for superposed Breit–Wigner peaks, this method outperformed the traditional MEM method—particularly for multi-peak spectra with large amounts of measurement noise.

Fig. 11Full-screen View

(Color online) Spectral functions reconstructed from MEM, NN, and NN-P2P under different amounts of Gaussian noises added to the propagator data with Np = 25, and $N_{\omega }=500$. Taken from Ref. [114].

In addition to Gaussian-like and Lorentzian-like spectral reconstruction tests, the newly devised framework presented in Refs. [114, 118] was validated through two physics-motivated tests. One was for non-positive definite spectral reconstruction, which is beyond the scope of classical MEM applicability but is often encountered for spectral functions related to confinement phenomenon of, e.g., gluons and ghosts, or thermal excitations with long-range correlation in strongly coupled systems. The other one was for the hadron spectral function encoded in the temperature-dependent thermal correlator with lattice QCD noise-level noises. For both of these physical cases, the proposed DNN and AD-based method with NN representation consistently works well, whereas traditional MEM based methods lose the peak information or fail to resolve the non-positiveness.

The spectral function can also be reconstructed from finite correlation data by implementing the radial basis function network (RBFN), which is an MLP model based on the RBF [127, 128]. The RBFN has been widely used in feature extraction, classification, regression, etc. [129-132]. In Ref. [123], the spectral function $\rho \left(\omega \right)$ was approximately described by a linear combination of RBFs: $\rho \left(\omega \right)={\displaystyle \sum _{j=1}^N{w}_j\varphi }\left(\omega -m_j\right)\mathrm{,}$ (45) where ϕ represents the active RBF with an adjustable weight wj and an adjustable center mj, which can take a Gaussian form $\varphi (r)=e^{-\frac{r^2}{2a^2}}$ or an MQ form $\varphi (r)=(r^2+a^2{)}^{\frac{1}{2}}$. Here, a is the shape parameter, which is adjustable and essential for the regularization. Then, the inverse mapping problem of constructing the spectral function is transformed into calculating the linear weights of the RBF, which allows smooth and continuous reconstruction.

For calculating these parameters in Eq.(45), in Ref. [123], a neutral network called the RBFN was constructed, which is a three-layer feedforward neural network with the active RBFs in the hidden layer. After discretization of the spectral function, Eq.(45) is converted to matrix form: $\left[\rho \right]=\left[\Phi \right]\left[W\right]$. Then, the correlation functions in the Euclidean space with the integral spectral representation $G\left(\tau \mathrm{,}T\right)={\displaystyle {\int }_0^{\infty }}\frac{d\omega }{2\pi }\rho \left(\omega \mathrm{,}T\right)K\left(\omega \mathrm{,}\tau \mathrm{,}T\right)$ are converted to matrix form: $G_i={\displaystyle \sum _{j=1}^M{\displaystyle \sum _{k=1}^N{K}_{ij}{\Phi }_{jk}{w}_k}}\equiv {\displaystyle \sum _{k=1}^M{\tilde{K}}_{ik}{w}_k}\mathrm{,}i=\mathrm{1...}\widehat{N}$ (46) where $\tilde{K}$ is a $\widehat{N}\times M$ matrix associated with the integration kernel, and $\widehat{N}$ represents the number of data points for the correlation function Gi. The spectral function $\rho ({\omega }_i)$ has been discretized into N parts with $m_i={\omega }_i\mathrm{,}i=\mathrm{1...}N$, and M is set to M=N=500. To obtain wj, one can use the truncated singular value decomposition (TSVD) method or a DNN [123]. Compared with other ML approaches based on supervised learning [133, 134, 113], this method allows faster training and is free from the overfitting problem.

Fig. 12 shows a comparison of the spectral functions reconstructed using RBFN, TSVD, Tikhonov, and MEM, using the correlation data generated by a mock SPF. The mock SPF was obtained by mixing two Breit–Wigner distributions: ${\rho }_{\text{Mock}}(\omega )={\rho }_{\text{BW}}\left(A_1\mathrm{,}{\Gamma }_1\mathrm{,}{M}_1\mathrm{,}\omega \right)+{\rho }_{\text{BW}}\left(A_2\mathrm{,}{\Gamma }_2\mathrm{,}{M}_2\mathrm{,}\omega \right)$ with ${\rho }_{\text{BW}}\left(A_i\mathrm{,}{\Gamma }_i\mathrm{,}{M}_i\mathrm{,}\omega \right)=\frac{4A_i{\Gamma }_i\omega }{{\left(M_i^2+{\Gamma }_i^2-{\omega }^2\right)}^2+4{\Gamma }_i^2{\omega }^2}$. The parameters for the mock SPF in Fig. 12 were set to $A_1=\mathrm{0.8,}{M}_1=\mathrm{2,}{\Gamma }_1=0.5;A_2=\mathrm{1,}{M}_2=\mathrm{5,}{\Gamma }_2=0.5$. Here, 30 discrete correlation data were generated using the Euclidean correlation functions of the mock SPF, with noise added, i.e., $G_{\text{noise}}({\tau }_i)=G({\tau }_i)+\text{noise}$.

Fig. 12Full-screen View

(Color online) Constructed spectral functions obtained from RBFN, TSVD, Tikhonov regularization, and MEM, using the correlation data generated by the mock SPF acquired by mixing two Breit–Wigner distributions. From left to right, different Gaussian noises are added to the correlation data with ϵ = 0.001, 0.0001, and 0.00001. Taken from Ref. [123].

Compared with the results of traditional methods, the RBFN provided a better description of the spectral functions—particularly for the low-frequency part. It almost reproduced the first peak of the mock SPF using the correlation data with a small amount of noise ϵ = 0.00001. In contrast, Tikhonov, TSVD, and MEM exhibited oscillation behavior at a low frequency. For such a task of extracting the transport coefficients from the Kubo relation, an improved reconstruction of the spectral functions at a low frequency is important. Although the RBFN failed to reconstruct the second peak of the mock SPF, it was the only method that reduced the oscillation at the low frequency, among the methods tested. In Ref. [123], the Gaussian and MQ RBFs used in the network were compared, and it was found that the Gaussian RBF provided better construction of the SPF, including the location and the width of the peak. Additionally, with mock data generated from the spectral function of the energy momentum tensor, it was demonstrated that the RBFN method allows precise and stable extraction of the transport coefficients.

In-medium heavy quark potential

As an important probe for the properties of the created QGP in HICs, heavy quarkonium (the bound state of a heavy quark and its anti-quark) has been intensively measured in experiments and analyzed in theoretical studies [135, 136], wherein the investigation and calculation require an understanding of the in-medium heavy quark interaction. The heavy quarkonium provides a calibrated QCD force, because in vacuum the simple Cornell potential can well reproduce the spectroscopy of heavy quarkonium, and when we put the bound state into the QCD medium, the color screening effects naturally occur and weaken the interactions between the heavy quarks, beyond which a non-vanishing imaginary part manifested as thermal width is argued to appear according to both one-loop hard thermal loop (HTL) perturbative QCD calculations [137, 138] and recent effective field theory (EFT) studies, e.g., those on PNRQCD [139, 140]. However, a non-perturbative treatment similar to that of lattice QCD is necessary because it is difficult to obtain a satisfactory description of the strong interaction dictated in-medium heavy quarkonium solely from perturbative calculations. These EFT studies suggested that a potential-based picture can provide a good approximation of the quarkonium, under which the Schrödinger equation can be employed to study the spectroscopy of the bound state. Recent lattice QCD studies involved quantification of the in-medium spectrum–mass shift and thermal widths of bottomonium ($b\overline{b}$) up to 3S and 2P states in QGP [125], where it was found cannot be reproduced by the one-loop HTL-motivated functional form of the heavy quark in-medium potential, i.e., $V_{\text{R}}(T\mathrm{,}r)$ and $V_{\text{I}}(T\mathrm{,}r)$. Note that the mass shift may affect the quarkonium production in HICs [141].

In Ref. [124], the authors developed a model-independent DNN-based method for reconstructing the temperature and inter-quark distance-dependent in-medium heavy quark potential according to the aforementioned lattice QCD results for bottomonium. Inspired by the universal approximation theorem, the authors introduced the DNN to parameterize the potential in an unbiased yet flexible manner (can be named as potential-DNN). The DNN-represented heavy quark potential is coupled to the Schrödinger equation solving process to be converted into complex valued energy eigenvalues En, which are related to the bound state in-medium mass and thermal width through Re[En]=mn-2mb and Im[En]=-Γ_n. Through comparison with the lattice QCD “measurements”, the corresponding χ₂ provide the loss function for optimizing the parameters of the potential-DNN: $\mathcal{L}=\frac{1}{2}{{\displaystyle \sum _{T\mathrm{,}n}\left(\frac{m_{T\mathrm{,}n}-m_{T\mathrm{,}n}^{\text{LQCD}}}{\delta m_{T\mathrm{,}n}^{\text{LQCD}}}\right)}}^2+{\left(\frac{{\Gamma }_{T\mathrm{,}n}-{\Gamma }_{T\mathrm{,}n}^{\text{LQCD}}}{\delta {\Gamma }_{T\mathrm{,}n}^{\text{LQCD}}}\right)}^2\mathrm{,}$ (47) with $T\in \left\{\mathrm{0,}\mathrm{151,}\mathrm{173,}\mathrm{199,}\mathrm{251,}334\right\}$ MeV and $n\in \left\{1S\mathrm{,}2S\mathrm{,}3S\mathrm{,}1P\mathrm{,}2P\right\}$ according to the lattice QCD evaluation conditions. Gradient descent with backpropagation can be applied for the DNN optimization here, where the gradient is estimated efficiently from perturbative analysis based on the Schrödinger equation with respect to the perturbative change of the potential and just arrived at the Hellman–Feynman theorem. Furthermore, the uncertainty of the reconstructed potential is quantified via Bayesian inference; thus, the posterior distribution of the DNN parameters is evaluated. With the outlined approach, in Ref. [124], good agreement with the lattice QCD results for the masses and thermal widths of bottomonium was achieved simultaneously; see the left and middle panels of Fig. 13. Additionally, the temperature- and distance-dependent heavy quark potential was obtained, as shown in the right panel of Fig. 13. Clearly, the color screening effect emerged for the reconstruction with a flatter structure appearing in $V_{\text{R}}(T\mathrm{,}r)$ with the increasing temperature at a long distance, but the temperature dependence was mild compared with the perturbative analysis-based results in the same temperature range. In contrast, the imaginary part, i.e., $V_{\text{I}}(T\mathrm{,}r)$, exhibited significant growth with respect to both temperature and distance and also exhibited a larger magnitude than the one-loop HTL-motivated results.

Fig. 13Full-screen View

(Color online) Picture taken from Ref. [124]. Left and middle: In-medium mass shifts with respect to the vacuum mass (left) and the thermal widths (right) of different bottomonium states obtained from fits to LQCD results of Ref. [125] (lines and shaded bands) using weak-coupling-motivated functional forms [126] (open symbols) and DNN-based optimization (solid symbols). The points are shifted horizontally for better visualization. Υ(1S), ${\chi }_{b_0}(1P)$, Υ(2S), ${\chi }_{b_0}(2P)$, and Υ(3S) states are represented by red circles, orange pluses, green squares, blue crosses, and purple diamonds, respectively. Right: The DNN-reconstructed real (top) and imaginary (bottom) parts of the heavy quark potential at temperatures of T = 0 (black), 151 (purple), 173 (blue), 199 (green), 251 (orange), and 334 MeV (red). The uncertainty bands represent the 68%(1σ) confidence region.

Deep learning for quasi-particle mass

The EoS of hadron resonance gas in the QCD phase diagram can be calculated using a simple statistical formula with the following partition function: $\ln Z(T)={\displaystyle \sum _i\ln Z_{hi}}(T),$ (48) where $Z_{hi}(T)$ is the partition function for one of the several hundred hadrons in HRG, assuming that there is no interaction between different hadrons. The obtained EoS agrees with lattice QCD calculations. It is impossible to obtain the lattice QCD EoS for QGP using the same formula, as quarks and gluons interact with each other and form a many-body quantum system. However, if one assumes that the quarks and gluons are non-interacting quasi-particles whose masses depend on the local temperature, the lattice QCD EoS can be reproduced using the following simple statistical formula: $\begin{array}{l}\ln Z(T)=\ln Z_{\text{g}}(T)+{\displaystyle \sum _i\ln Z_{q_i}}(T)\\ \ln Z_g(T)=-\frac{d_{\text{g}}V}{2{\pi }^2}{\displaystyle {\int }_0^{\infty }}{p}^2\text{d}p\\ \ln \left[1-\exp \left(-\frac{1}{T}\sqrt{p^2+m_{\text{g}}^2(T)}\right)\right]\\ \ln Z_{q_i}(T)=+\frac{d_{q_i}V}{2{\pi }^2}{\displaystyle {\int }_0^{\infty }}{p}^2\text{d}p\\ \ln \left[1+\exp \left(-\frac{1}{T}\sqrt{p^2+m_{q_i}^2(T)}\right)\right],\end{array}$ (49) where Zg represents the partition function of quasi-gluons; $Z_{q_i}$ represents the partition function of quasi-quarks; dg and $d_{q_i}$ represent the spin and color degeneracy for gluons and quarks, respectively; p represents the magnitude of momentum; and T represents the local temperature. Gluons, along with up, down, and strange quarks, are considered in this calculation. It is assumed that the temperature quasi-particle masses $m_{\text{u/d}}(T)$ are the same for up and down quarks but different for gluons $m_{\text{g}}(T)$ and strange quarks $m_{\text{s}}(T)$. Thus, there are three variational functions whose forms are unknown and must be determined by mapping the following EoS to the lattice QCD EoS: $\begin{array}{l}P(T)=T{\left(\frac{\partial \ln Z(T)}{\partial V}\right)}_T\\ ϵ(T)=T^{2}V{\left(\frac{\partial \ln Z(T)}{\partial T}\right)}_V\end{array}$ (50) Several deep residual neural networks were constructed to represent the variational functions m_u/d(T), m_s(T), and mg(T). The mass functions of these quasi-partons are used in Eq. 49 to compute the partition function. The resulting partition function is used in Eq. 50 to compute the pressure and energy density as a function of the temperature. This procedure involves both numerical integration and differentiation. The integration is implemented using Gaussian quadrature with the TensorFlow library, while the differentiation is given by auto-differentiation. The loss function is designed as $\text{loss}={\left|s_{\text{dnn}}-s_{\text{lattice}}\right|}^2+{\left|{\Delta }_{\text{dnn}}-{\Delta }_{\text{lattice}}\right|}^2+L_{\text{constrain}},$ (51) where $s=(ϵ+P)/T$ represents the entropy density, and $\Delta =(ϵ-3P)/T^4$ represents the trace anomaly. The L_constrain contains physical constraints in the high-temperature region whose theoretical function form is given by HTL calculations. The learned quasi-partons reproduce the lattice QCD EoS. Using these mass functions, the authors calculated η/s(T) and found that its minimum was located at approximately 1.25 Tc [142].

PCA for flow analysis

In relativistic HICs, the collective flow provides important information about the properties of the QGP and its initial state fluctuations [95-100]. The flow observables are generally defined by a Fourier decomposition of the produced particle distribution in the momentum space, such as $\frac{\text{d}N}{\text{d}\phi }=\frac{1}{2\pi }{\displaystyle \sum _{-\infty }^{\infty }{\overrightarrow{V}}_n{e}^{-in\phi }}=\frac{1}{2\pi }\left(1+2{\displaystyle \sum _{n=1}^{\infty }{v}_n{e}^{-in\left(\phi -{\Psi }_n\right)}}\right),$ (54) where ${\overrightarrow{V}}_n=v_n{e}^{in{\Psi }_n}$ is the flow vector of order n, vn represents flow harmonics of order n, and Ψn represents the corresponding event plane angle. Additionally, the flow coefficients can be obtained from the two-particle correlations associated with a Fourier decomposition: $\langle \frac{d{N}_{\text{pairs}}}{d{p}_{1}d{p}_2}\rangle \propto 1+2{\displaystyle \sum _{n=1}^{\infty }{V}_{n\Delta }}\left(p_{\text{T1}}\mathrm{,}{p}_{\text{T2}}\right)\cos \left(n\Delta \varphi \right),$ (55) where $V_{n\Delta }(p_{\text{T1}}\mathrm{,}{p}_{\text{T2}})$ is a symmetric covariance matrix and $\Delta \varphi ={\varphi }^a-{\varphi }^b$ represents the relative azimuthal angle between two emitted particles. Under the assumption of flow factorization, $V_{n\Delta }(p_{\text{T1}}\mathrm{,}{p}_{\text{T2}})$ is related to the flow harmonics $v_n(p_{\text{T}})$ as follows: $V_{n\Delta }(p_{\text{T1}}\mathrm{,}{p}_{\text{T2}})\approx v_n(p_{\text{T1}})v_n(p_{\text{T2}})$ [161] (For other flow methods and flow measurements, see [97, 100, 162-164]).

Recently, an ML technique called PCA based on SVD has been used to study the collective flow in relativistic HICs. For the two-particle correlations with the Fourier expansion [166-169], the event-by-event flow fluctuations have been investigated via PCA, revealing the substructures of the flow fluctuations [166-168]. Using PCA, V_nΔ(P_T1,P_t2) can be expressed as [167] $V_{n\Delta }\left(p_{\text{T1}}\mathrm{,}{p}_{\text{T2}}\right)={\displaystyle \sum _{\alpha }{v}_n^{(\alpha )}}\left(p_{\text{T1}}\right)v_n^{(\alpha )}\left(p_{\text{T2}}\right),$ (56) $\text{with }{\displaystyle \int }d{p}_{\text{T}}{w}^2\left(p_{\text{T}}\right)v_n^{(\alpha )}\left(p_{\text{T}}\right)v_n^{(\beta )}\left(p_{\text{T}}\right)={\lambda }_{\alpha }{\delta }_{\alpha \beta }$ (57) where $v_n^{(\alpha )}(p_{\text{T}})$ are the eigenvectors of the two-particle covariance matrix, and $w(p_{\text{T}})$ is the weight for the particle. α=1 denotes the leading modes, α=2 denotes the subleading modes, α=3 denotes the subsubleading modes, and so on. It was found that the leading modes correspond to the traditional flow harmonics and that the subleading modes lead to the breakdown of the flow factorization. In Ref. [167, 168], a linear relationship $V_n^{(\alpha )}\propto {\mathcal{E}}_n^{(\alpha )}$ was demonstrated for the leading, subleading, and subsubleading modes via hydrodynamic simulations. In Ref. [169], PCA was used to study the mode coupling between flow harmonics, which revealed hidden mode-mixing patterns that had not been previously discovered. Recently, the CMS collaboration extracted the subleading flow modes for Pb+Pb and p+Pb collisions at the LHC, reporting qualitative agreement between experimental measurements and theoretical calculations [170]. Using AMPT and HIJING simulations, Ref. [171] showed that the PCA modes depend on the choice of the p_T range and the particle weight w. In addition, the leading modes are influenced by non-flow effects, and the mixing between the non-flow and leading flow modes leads to fake subleading modes. Therefore, it is important to carefully handle the non-flow effects and the choice of the weight and phase space when implementing PCA to extract the subleading flow modes in both experimental and theoretical studies.

The aforementioned PCA studies on collective flow [166-171] were all based on the correlation data obtained with a Fourier expansion. Recently, PCA has been applied directly to single particle distributions dN/dφ without prior treatment with a Fourier transform, for exploring whether it can be used to directly discover flow without the guidance from humans [165]. Specifically, with PCA matrix multiplication, the i_th row of a particle distribution matrix with N events generated from VISH2+1 hydrodynamics can be expressed as $\begin{array}{c}dN/d{\phi }^{(i)}={\displaystyle \sum _{j=1}^m{x}_j^{(i)}{\sigma }_j{z}_j}={\displaystyle \sum _{j=1}^m{\tilde{v}}_j^{(i)}{z}_j}\\ \approx {\displaystyle \sum _{j=1}^k{\tilde{v}}_j^{(i)}{z}_j}(i)=\mathrm{1,...,}N\end{array}$ (58) Here, (i)=1,2,..., N is the index of the event. j is the index for the azimuthal angle, where the total azimuthal angle [-π,π] is divided into m bins to count the particles in each bin. After SVD, $dN/d{\phi }^{(i)}$ can be expressed by a linear combination of the eigenvectors zj with the corresponding coefficient ${\tilde{v}}_j^{(i)}$ (where j=1,2,...,m), and σ_i represents the diagonal elements (singular values) of the particle distribution matrix, which are arranged in descending order. In the spirit of PCA, in the last step, a cut is made at the indices k to focus only on the most important components.

Figs. 14 and 15 show the first 12 eigenvectors zj and the first 20 singular values σ_j of the PCA in descending order for the final-state matrix constructed from 2000 dN/dφ distributions with the azimuthal angle [-π,π] equally divided into 50 bins. Such dN/dφ distributions are generated from the VISH2+1 hydrodynamics with event-by-event fluctuating TRENTo initial conditions for 2.76 A TeV Pb+Pb collisions at 10%–20% centrality. Fig. 14 shows that the PCA eigenvectors are similar to the traditional Fourier bases. For example, the 1st and 2nd eigenvectors are close to sin(2φ) and cos(2φ), and the 3rd and 4th eigenvectors are close to sin(3φ) and cos(3φ). The corresponding singular values in Fig. 15 are arranged in pairs, which correspond to the real and imaginary parts of the anisotropic flow. It was found that for n≤6, the values of these PCA flow harmonics were very close to those of the traditional event-averaged flow harmonics obtained from the Fourier expansion but not exactly the same. Fig. 16 presents a comparison of the event-by-event flow harmonics obtained from PCA and from the traditional Fourier expansion. As shown, the elliptic flow with n=2 and the triangular flow with n=3 from the two methods agreed well. However, for higher flow harmonics with $n\ge 4$, the PCA and Fourier expansion results differed significantly owing to the mode-mixing effects. In Ref. [165], with these PCA flow harmonics v’n, the symmetric cumulants SC’(m,n) were calculated. Except for SC’(2,3), these PCA symmetric cumulants were significantly reduced compared with the traditional Fourier ones, because of the significantly increased linearity between the PCA flow harmonics and the initial eccentricities. These results indicated that PCA could define the collective flow on its own basis. Compared with the traditional ones obtained from the Fourier decomposition, the PCA method reduces the mode coupling effects between different flow harmonics [165].

Fig. 14Full-screen View

(Color online) PCA eigenvectors zj for the final-state matrix of particle distributions, generated from VISH2+1 hydrodynamics in 2.76 A TeV Pb+Pb collisions at 10%–20% centrality [165].

Fig. 15Full-screen View

(Color online) Singular values of PCA for the final-state matrix of particle distributions in Pb+Pb collisions at 10%–20% centrality [165].

Fig. 16Full-screen View

(Color online) Comparison between the event-by-event flow harmonics vn’ from PCA and vn from the Fourier expansion in Pb+Pb collisions at 10%-20% centrality [165].

CME detection

In the presence of a magnetic field, the chiral magnetic effect (CME) can occur when the system has a chiral imbalance, i.e., the numbers of left- and right-handed particles differ. Essentially, a current of electric charge (known as chiral magnetic current) can be induced to flow along the direction of the magnetic field. The use of the CME to reveal the vacuum structure of QCD has been proposed. In HICs, a strong magnetic field can be created by the motion of the colliding ions, and it is predicted that in the formed hot and dense QGP, the topological fluctuations of gluon fields may cause chiral imbalance for quarks. Accordingly, the CME may occur, which can manifest as a separation of electric charge along the magnetic-field direction. However, several challenges hinder the detection of the CME in HICs, among which the chief difficulty is disentangling the CME signal from other possible sources of charge separation (CS), e.g., elliptic flow, the global polarization, and other background noises, although multiple observables are proposed.

Despite the challenges, there is long-term and continuing interest in the search for the CME in HICs because of its general importance to QCD. Recently, Ref. [172] proposed the use of deep learning to construct an end-to-end CME-meter that can efficiently analyze the final-state hadronic spectrum as a whole in the sense of Big Data with a deep CNN to reveal the fingerprints of the CME. For supervised learning, the training set was prepared from the string melting AMPT model with CME implemented under a global CS scheme. Essentially, the CME events are generated by switching the y-components of momenta of a fraction of a downward moving light quark and its corresponding anti-quarks with upward moving direction. The fraction defines the CS fraction f, which separates the events into the “no CS” (label as “0”) class for those with f=0% and the “CS” class (label as “1”) for those with f>0%. Each event is represented as 2D transverse momentum and azimuthal angle spectra of charged pions in the final state, i.e., ${\rho }_{\pi }(p_T\mathrm{,}\varphi )$. Then, the deep CNN is trained to perform binary classification on the labeled events with the spectra to be the input. Fig. 17 shows the architecture of the developed deep CNN for CME-meter construction.

Fig. 17Full-screen View

(Color online) Taken from [172]. The CNN architecture with π⁺ and π^- spectra ${\rho }^{\pm }(p_T\mathrm{,}\varphi )$ as inputs.

As shown in Fig. 17, the output of the network has two nodes, each of which is naturally interpreted as the probability resulting from the network decision in recognizing any given input spectrum as CME (P₁) or non-CME (P₀=1-P₁) events. The training set contains multiple collision beam energies and centralities for diversity consideration. The pion spectrum is obtained by averaging over 100 events with the same collision condition to reduce the fluctuations, which also reduces the backgrounds and thus should be considered a prerequisite for realistic application in experiments. For the training, different levels of the CS fraction are used, and it is found that the classification validation accuracy is lower for a smaller number of CS fraction training events. This indicates that a larger CS fraction can be identified more easily, which is expected. Despite the different induced discernibility, the trained deep CNNs all exhibited robust performance against varying collision centrality and energy. One can conclude that at least under the AMPT modeling level, the CS signals can survive into the final state of the collision dynamics at different collision conditions, which can be recognized by the deep CNN-based CME-meter.

Note that the network was trained only on Au+Au collision systems, while the extrapolation to other collision system was validated. Specifically, the obtained CME-meter was applied to isobaric collisions of ${}_{40}^{96}$Zr+ ${}_{40}^{96}$Zr and ${}_{44}^{96}$Ru+ ${}_{44}^{96}$Ru, which were proposed for the search of the CME. Because Ru contains more protons than Zr, which induce a stronger magnetic field, there is expected to be a stronger CS signal in Ru+Ru collisions. To reveal this difference and also the distinguishable difference of the CME-meter for the two isobaric collision systems from $P_1^{\text{Ru}}>P_1^{\text{Zr}}$, the Riso was evaluated, and the results validated the developed CNN-based CME-meter: $R_{\text{iso}}=2\times \frac{\langle \text{logit}\left(P_1^{\text{Ru}}\right)\rangle -\langle \text{logit}\left(P_1^{\text{Zr}}\right)\rangle }{\langle \text{logit}\left(P_1^{\text{Ru}}\right)\rangle +\langle \text{logit}\left(P_1^{\text{Zr}}\right)\rangle }\mathrm{,}$ (59) where the function $\text{logit}(x)=\text{log}\left[x/\left(1-x\right)\right]$ is used to restore the derivative in the saturation region of the activation in the last layer of the NN, i.e., softmax. From Tab. 1 on Riso, the trained CME-meter was well validated beyond the training collision system, indicating its robust capture of the general CME signal in the collisions.

The CME-meter was also validated through a different model simulation, i.e., anomalous-viscous fluid dynamics (AVFD). P₁ exhibited a consistent positive correlation with N_5/S, which controlled the CME strength, while the contamination from local charge conservation (LCC) up to 30% did not augment the performance of the CME-meter on the testing events from AVFD. In Ref. [172], to reveal the underlying account for the trained CME-meter, the network output P₁ and γ-correlator were compared. The γ-correlator—a conventional CME probe—can measure the event-by-event two-particle azimuthal correlation of charged hadrons. It was shown that for averaged events, both the CME signal and the background from δγ (difference between correlations within particles of the same charge and within particles of opposite charge) are suppressed. Being differently, the CME-meter output P₁ works well in classifying CS and no-CS classes on the averaged events.

The direct implementation of this trained CME-meter in real experiments would require reconstructing the reaction plane of each collision event to form the averaged events as input for the meter. In general, the reaction plane reconstruction can be achieved by measuring correlations of final-state particles, and it inevitably contains finite resolution and background effects. It was shown that even in restricted event plane reconstruction, the trained CME-meter can recognize the CS signals. For the deployment of the trained CME-meter on single event measurements, Ref. [172] proposed a hypothesis test perspective.

Another way to interpret the trained deep-learning algorithm is the DeepDream method, which was used in Ref. [172] to reconstruct the network most responding input pion spectrum, manifesting the “CME pattern” that the CNN-based CME-meter essentially captured for its further CME signal recognition. The key idea is to perform variational tuning on the input pion spectrum with the trained and frozen network to maximize its output (i.e., pushing $P_1\to 1$), driven by the gradient $\delta P_1\left({\rho }_{\pi }\left(p_T\mathrm{,}\varphi \right)\right)/\delta {\rho }_{\pi }\left({\partial }_T\mathrm{,}\varphi \right)$. The resultant “CME pattern” from the trained network is displayed in Fig. 18, where the charge conservation and a clear dipole structure appear, both being CME-related features.

Fig. 18Full-screen View

(Color online) DeepDream map for the (0%+10%) model [172].

Summary

As a modern computational paradigm, AI—particularly machine- and deep-learning techniques—has introduced a wealth of applications and new possibilities in scientific research. Owing to its special ability to recognize patterns and structures hidden in complex data, these learning-based strategies make physics exploration with a Big Data or smart computation mindset feasible. In the context of HENP revolving around HIC programs to understand nuclear matter properties under different conditions, various research fields have benefited from the incorporation of these techniques.

In this mini-review, we presented the recent progress in the field of HICs, including initial state physics inference, QCD matter transport and bulk properties, thermal medium modifications for partons or hadrons, and recognition of physical observables in HICs.

We first reviewed different loss functions l used in supervised learning, unsupervised learning, semi-supervised learning, self-supervised learning, and active learning. During training, the negative gradients $-\partial l\partial \theta$ are used to optimize the network in SGD-like algorithms. AD is employed to efficiently compute the derivatives of the loss with regard to model parameters θ. As auto-diff has analytical precision for the variational function represented by the neural network, it has been widely used in physics-informed neural networks to solve ODEs and PDEs. We then introduced the widely used neural-network architectures, such as the MLP, CNN, RNN, and point cloud network. Next, we explained the generative models, such as the autoencoder, GAN, flow model, and diffusion models, in detail. These models are widely used in lattice QCD to generate field configurations.

For the initial condition, ML has been widely used to determine the centrality classes and impact parameters using the final-state hadrons in the momentum space, for extracting the initial nuclear structures, such as the nuclear deformation, the α clustering, and the neutron skin. In general, it is easier to extract the nuclear deformation than the α clustering and neutron skin according to the current literature.

For bulk matter, Bayesian parameter estimation has been successfully used to determine the temperature-dependent shear and bulk viscosities of QGP. An unsupervised autoencoder was used to reconstruct the charged multiplicity distributions, which helps to determine the source temperature and the temperature of the nuclear liquid gas phase transition. Deep CNNs, point cloud networks, and event-averaging techniques are employed to classify the crossover and first-order phase transition regions in the QCD phase diagram, using data generated by relativistic hydrodynamic models and hadronic transport models. Active learning is used to map out thermodynamically unstable regions near the critical endpoint in the QCD phase diagram. For hydrodynamic evolution, a well-designed network called sU-net can capture the nonlinear mapping between the initial and final profiles with sufficient precision, which is also far faster than the traditional hydrodynamic simulations.

For QGP in-medium effects, we first reported some of the recently proposed ML-based methods for spectral function reconstruction, which is a notorious ill-posed inverse problem. Both supervised and unsupervised methods have been discussed for the inference of spectral out-of-Euclidean correlator measurements from Monte Carlo simulations (e.g., lattice study). Then, in-medium heavy quark interaction inference based on in-medium heavy quarkonium spectroscopy was introduced. A novel DNN representation integrated inside the forward problem-solving pipeline with AD was proposed. This strategy is also used for in-medium quasi-particle effective model construction from the lattice QCD EoS.

For hard probes, Bayesian analysis is widely used to extract the temperature-dependent jet (or heavy quark) transport coefficient $\widehat{q}(T)$ and the jet energy loss distributions. Recently, deep learning-assisted jet tomography was developed to locate the initial jet production positions. This is important for the study of jet substructures and the medium response. Using this technique, it was observed that the signal of jet-induced Mach cones is amplified by selecting jet events.

For the observables, PCA has been implemented to study the collective flow in relativistic HICs. This revealed the substructures of the flow fluctuations, which can potentially be used to extract the subleading modes of flow with efforts from both the experimental and theoretical sides. When applied directly to the single particle distributions, PCA can directly discover flow with a basis similar to the Fourier expansion ones, which significantly reduces the mode coupling between different flow harmonics.

Outlook

Despite the impressive progress, the interplay between HENP and ML is still inducing hectic evolution. Many questions and challenges remain and deserve further exploration. In addition to the aforementioned applications of ML in the field of HICs, several other topics can be explored with ML, e.g., critical endpoint searching for eRHIC and the Electron Ion Collider (EIC) regime [173], spin polarization study, the upcoming FAIR program, Nuclotron-based Ion Collider facility (NICA) experiments, nuclear structure inference, and High Intensity Heavy-ion Accelerator Facility (HIAF) experiments in China. Regarding the future prospects of applying ML techniques for HIC physics research, because this field is rapidly evolving, we present questions that we consider worthy of future investigation:

• Can ML provide more efficient "observables” to pin down the desired physics?

• Can the algorithms provide new physical knowledge to advance our understanding of nuclear matter?

• How can we make the ML algorithms to be confronted with realistic experiments? Is on-line analysis possible? How can experimental raw data be accessed to test neural networks pretrained with model simulations?

• Is it possible to accelerate HIC dynamical simulations for performing high statistics measurements or Bayesian inference?

• How can Bayesian inference be combined with ML to advance our field and better connect experiment to theory?

• How can symmetries be fully incorporated into the analysis using ML, e.g., Lorentz Group Equivariant Autoencoders [174]? How can dimensionality analysis (constraints) be incorporated into the ML methods properly and consistently?

It is also important to consider how we can adopt potentially useful approaches from other fields, e.g., particle physics, condensed-matter physics, and astrophysics, and how the community can better organize with joint efforts, e.g., for maximizing the potential of these novel computational techniques to advance the field of HENP.