# Performance Analysis of IRS-Assisted Cell-Free Communication

Diluka Loku Galappaththige, Dhanushka Kudathanthirige, and Gayan Amarasuriya
School of Electrical, Computer, and Biomedical Engineering, Southern Illinois
University, Carbondale, IL, USA 62901
Email: {diluka.lg, dhanushka.kudathanthirige<EMAIL_ADDRESS>
###### Abstract

In this paper, the feasibility of adopting an intelligent reflective surface
(IRS) in a cell-free wireless communication system is studied. The received
signal-to-noise ratio (SNR) for this IRS-enabled cell-free set-up is optimized
by adjusting phase-shifts of the passive reflective elements. Then, tight
approximations for the probability density function and the cumulative
distribution function for this optimal SNR are derived for Rayleigh fading. To
investigate the performance of this system model, tight bounds/approximations
for the achievable rate and outage probability are derived in closed form. The
impact of discrete phase-shifts is modeled, and the corresponding detrimental
effects are investigated by deriving an upper bound for the achievable rate in
the presence of phase-shift quantization errors. Monte-Carlo simulations are
used to validate our statistical characterization of the optimal SNR, and the
corresponding analysis is used to investigate the performance gains of the
proposed system model. We reveal that IRS-assisted communications can boost
the performance of cell-free wireless architectures.

##  I Introduction

Recently, wireless architectures based on the notion of cell-free have gained
much interest [1, 2]. In a cell-free system set-up, the cell-boundaries can be
relaxed, and thus, a vast number of access-points (APs) can be spatially
distributed to serve all users with a uniformly better quality-of-service
(QoS) over a much larger geographical region [1, 2]. Moreover, cell-free set-
ups may render spectral/energy efficiency gains, mitigate impediments caused
by spatial-correlated fading in compact/co-located antenna arrays, and
circumvent shadow fading impairments [1, 2]. Thus, cell-free architecture is a
foundation for practically realizing extremely large antenna arrays for next-
generation wireless standards.

An intelligence reflective surface (IRS) consists of a large number of passive
reflectors, whose reflective coefficients can be adjusted to attain desired
propagation effects for the impinging electromagnetic (EM) waves [3, 4]. The
feature of intelligently adjustable phase-shifts at an IRS can be used to
boost the signal-to-noise ratio (SNR) and to mitigate co-channel interference
at an intended destination through constructive and destructive signal
combining, respectively [5]. This leads to the notion of recycling of EM waves
within a propagation medium, and thereby, spectral/energy efficiency gains and
implementation cost reduction can be realized as IRSs are made out of low-cost
meta-atoms without active radio-frequency (RF) chains/amplifiers [4].

###  I-A Our motivation

In this paper, we aim to investigate the feasibility of embedding an IRS
within a cell-free set-up. Specifically, our objective is to investigate the
performance of an IRS-assisted cell-free set-up, and thereby, we explore the
feasibility of jointly reaping the aforementioned benefits of cell-free
architectures and IRS-assisted wireless channels. Moreover, to the best of the
authors knowledge, the fundamental performance metrics for an IRS-assisted
cell-free set-up have not yet been reported in open literature. To this end,
we aim to fill this important gap in IRS literature by presenting a
performance analysis for an IRS-assisted cell-free set-up.

###  I-B A literature survey for cell-free architecture and performance
analysis of IRS-assisted channels

In [1, 2], the basic concept of cell-free architectures is investigated, and
thereby, the performance metrics are compared against those of the co-located
antenna arrays. The analyses in [1, 2, 6] reveal that the cell-free set-ups
can outperform the co-located counterparts by serving users with a uniformly
better QoS, minimizing the impediments of spatial-correlation, and shortening
the end-to-end transmission distances to boost the overall energy/spectral
efficiency [1, 2]. Reference [7] proposes max-min power optimization
algorithms for cell-free massive multiple-input multiple-output (MIMO). In
[8], the performance of cell-free massive MIMO with underlay spectrum sharing
is investigated.

References [3, 4] present core architectural design principles of IRSs for
wireless communications. Ray-tracing techniques are used in [9] to generate a
novel path-loss model for IRS-assisted wireless channels. In [10], joint
optimization of precoder at the base-station (BS) and phase-shifts at the IRS
is studied through semi-definite relaxation and alternative optimization
techniques. Reference [5] studies the fundamental performance limits of
distributed IRS-assisted end-to-end channels with Nakagami-$m$ fading
channels. In [11], by using the statistical channel state information (CSI),
an optimal phase-shift design framework is developed to maximize the
achievable rates of IRS-assisted wireless channels. In [12], joint beamforming
and reflecting coefficient designs are investigated for IRSs to provision
physical layer security. Reference [13] proposes a practical IRS phase-shift
adjustment model, and thereby, the achievable rate is maximized through
jointly optimizing the transmit power and the BS beamformer by using
alternative optimization techniques.

###  I-C Our contribution

In above-referred prior research [10, 5, 11, 13, 12] for IRS-assisted
communications, a BS with either a single-antenna or a co-located antenna
array is used. Having been inspired by this gap in IRS/cell-free literature,
in this paper, we investigate an IRS-assisted wireless channel embedded within
a cell-free set-up over Rayleigh fading, and thereby, we present fundamental
performance metrics. To this end, first, we invoke the central limit theorem
(CLT) to tightly approximate the end-to-end optimal SNR to facilitate a
mathematically tractable probabilistic characterization. Then, we derive the
probability density function (PDF) and the cumulative density function (CDF)
of this approximated optimal SNR in closed-form. Thereby, we present a tight
approximation to the outage probability. Moreover, we derive tight upper/lower
bounds for the achievable rate. In particular, we investigate the impediments
of discrete phase-shifts in the presence of phase-shift quantization errors.
Finally, we present a set of rigorous numerical results to explore the
performance gains of the proposed system, and we validate the accuracy of our
analysis through Monte-Carlo simulations. From our numerical results, we
observe that by using an IRS with controllable phase-shift adjustments, the
performance of cell-free wireless set-ups can be enhanced.

Notation: The transpose of vector $\mathbf{y}$ is denoted as $\mathbf{y}^{T}$.
The expectation and variance of a random variable $Y$ are represented by
$\mathbb{E}\\!\left[{Y}\right]$ and
$\mathbb{V}\mathrm{ar}\\!\left[{Y}\right]$, respectively.
$Y\sim\mathcal{CN}\left(\mu_{Y},\sigma_{Y}^{2}\right)$ denotes that $Y$ is
complex-valued circularly symmetric Gaussian distributed with $\mu_{Y}$ mean
and $\sigma_{Y}^{2}$ variance. Moreover, $C_{n}=\\{0,1,\cdots,n\\}$ and
$C_{n}^{\prime}=C_{n}/\\{0\\}$.

$D$$\mathrm{IRS}$$\mathrm{AP}_{1}$$\mathrm{AP}_{m}$${\mathbf{g}}$${\mathbf{h}_{m}}$${{u_{m}}}$$\mathrm{AP}_{M}$$\mathrm{CPU}$
Figure 1: System model - IRS-aided cell-free communication set-up

##  II System, Channel and Signal Models

###  II-A System and channel model

We consider a cell-free communication set-up consisting of $M$ single-antenna
APs ($\mathrm{AP}_{m}$ for $m=1,\cdots,M$) and a single-antenna destination
$(D)$. An IRS having $N$ passive reflective elements is embedded within this
cell-free set-up as shown in Fig. 1. For the sake of exposition, we denote the
set of APs as $\mathcal{M}=\\{1,\cdots,M\\}$ and the set of reflective
elements at the IRS as $\mathcal{N}=\\{1,\cdots,N\\}$.

The direct link between the $m$th AP and $D$ is represented by $u_{m}$, while
$h_{mn}$ denotes the channel between the $m$th AP and the $n$th reflective
element of the IRS. Moreover, $g_{n}$ is used to represent the channel between
the $n$th reflective element of the IRS and $D$. We model the envelops of all
aforementioned channels to be independent Rayleigh distributed [14], and the
corresponding polar-form of these channels is given by

$\displaystyle v=\lambda_{v}\mathrm{e}^{j\theta_{v}},$ (1)

where $v\in\\{u_{m},h_{mn},g_{n}\\}$ for $m\in\mathcal{M}$ and
$n\in\mathcal{N}$. In (1), the envelop and the phase of $v$ are given by
$\lambda_{v}$ and $\theta_{v}$, respectively. The PDF of $\lambda_{v}$ is
given by [15]

$\displaystyle
f_{\lambda_{v}}(x)=\left({x}/{\xi_{v}}\right)\mathrm{exp}\left({-x^{2}}/{\left(2\xi_{v}\right)}\right),$
(2)

where $\xi_{v}=\zeta_{v}/2$ is the Rayleigh parameter, and $\zeta_{v}$
captures the large-scale fading/path-loss of the channel $v$. Since all
reflective elements are co-located within the IRS, it is assumed that all
large-scale fading parameters are the same.

###  II-B Signal model

The signal transmitted by the $m$th AP reaches $D$ through the direct and IRS-
assisted reflected channels. Thus, we can write the signal received at $D$ as

$\displaystyle
r=\sqrt{P}\sum\nolimits_{m\in{\mathcal{M}}}\left(u_{m}+\mathbf{g}^{T}\mathbf{\Theta}\mathbf{h}_{m}\right)x+w,$
(3)

where $x$ is the transmit signal from $S$ satisfying
$\mathbb{E}\\!\left[{|x|^{2}}\right]=1$, $P$ is the transmit power at each AP,
and $w$ is an additive white Gaussian noise (AWGN) at $D$ with zero mean and
variance of $\sigma_{w}^{2}$ such that $w\sim\mathcal{CN}(0,\sigma_{w}^{2})$.
In (3),
$\mathbf{h}_{m}=[h_{m1},\cdots,h_{mn},\cdots,h_{mN}]^{T}\in\mathbb{C}^{N\times
1}$ is the channel vector between the $m$th AP and the IRS. Moreover,
$\mathbf{g}^{T}=[g_{1},\cdots,g_{n},\cdots,g_{N}]\in\mathbb{C}^{1\times N}$
denotes the channel vector between the IRS and $D$. The diagonal matrix,
$\mathbf{\Theta}=\mathrm{diag}\left(\beta_{1}\mathrm{e}^{j\theta_{1}},\cdots,\beta_{n}\mathrm{e}^{j\theta_{n}},\cdots,\beta_{N}\mathrm{e}^{j\theta_{N}}\right)\in\mathbb{C}^{N\times
N}$, captures the reflective properties of the IRS through complex-valued
reflection coefficients $\beta_{n}\mathrm{e}^{j\theta_{n}}$ for
$n\in\mathcal{N}$, where $\beta_{n}$ and $\theta_{n}$ are the magnitude of
attenuation and phase-shift of the $n$th reflective element of the IRS,
respectively. Thus, we can rewrite the received signal at $D$ in (3) as

$\displaystyle
r=\sqrt{P}\sum\nolimits_{m\in{\mathcal{M}}}\left(u_{m}+\sum\nolimits_{n\in{\mathcal{N}}}\beta_{n}{g}_{n}{h}_{mn}\mathrm{e}^{j\theta_{n}}\right)x+w.$
(4)

Thereby, we derive the SNR at $D$ from (4) as

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\gamma$
$\displaystyle=$
$\displaystyle\bar{\gamma}\left|\sum\nolimits_{m\in{\mathcal{M}}}\left(u_{m}+\sum\nolimits_{n\in{\mathcal{N}}}\beta_{n}{g}_{n}{h}_{mn}\mathrm{e}^{j\theta_{n}}\right)\right|^{2}$
(5) $\displaystyle=$
$\displaystyle\bar{\gamma}\left|\sum\nolimits_{m\in{\mathcal{M}}}\\!u_{m}\\!+\\!\sum\nolimits_{n\in{\mathcal{N}}}\beta_{n}{g}_{n}\\!\left(\sum\nolimits_{m\in{\mathcal{M}}}{h}_{mn}\right)\mathrm{e}^{j\theta_{n}}\right|^{2}\\!\\!,$

where the average transmit SNR is denoted by $\bar{\gamma}=P/\sigma_{w}^{2}$.
Then, we define $u=\sum_{m\in{\mathcal{M}}}u_{m}$ and
$h_{n}=\sum_{m\in{\mathcal{M}}}h_{mn}$. Since $u_{m}$ and $h_{mn}$ are
independent complex Gaussian distributed for $m\in\mathcal{M}$ and
$n\in\mathcal{N}$, the polar-form of $u$ and $h_{n}$ can be also expressed
similar to (1), where $\lambda_{u}$ and $\lambda_{{h}_{n}}$ are the envelops
of $u$ and $h_{n}$, respectively. Thus, $\lambda_{u}$ and $\lambda_{{h}_{n}}$
are independent Rayleigh distributed with parameters
$\xi_{u}=\sum_{m\in{\mathcal{M}}}\zeta_{u_{m}}/2$ and
$\xi_{h_{n}}=\sum_{m\in{\mathcal{M}}}\zeta_{h_{mn}}/2$, respectively. From
(1), we can rewrite the SNR in (5) in terms of the channel phases as

$\displaystyle\gamma=\bar{\gamma}\left|\lambda_{u}\mathrm{e}^{j\theta_{u}}+\sum\nolimits_{n\in{\mathcal{N}}}\beta_{n}\lambda_{{g}_{n}}\lambda_{{h}_{n}}\mathrm{e}^{j\left(\theta_{n}+\theta_{{g}_{n}}+\theta_{{h}_{n}}\right)}\right|^{2}.$
(6)

It can be seen from (6) that the received SNR at $D$ can be maximized by
smartly adjusting the phase-shifts at each IRS reflecting elements
$(\theta_{n})$. Thus, it enables a constructive addition of the received
signals through the direct channels and IRS-aided reflected channels [10, 16].
To this end, the optimal choice of $\theta_{n}$ is given by
$\theta_{n}^{*}=\underset{-\pi\leq\theta_{n}\leq\pi}{\mathrm{argmax}}\;{\gamma}=\theta_{u}-\left(\theta_{g_{n}}+\theta_{h_{n}}\right)$.
Then, we can derive the optimal SNR at $D$ as

$\displaystyle\gamma^{*}=\bar{\gamma}\left|\lambda_{u}+\sum\nolimits_{n\in{\mathcal{N}}}\beta_{n}\lambda_{{g}_{n}}\lambda_{{h}_{n}}\right|^{2}.$
(7)

##  III Preliminaries

In this section, we present a probabilistic characterization of the optimal
received SNR at $D$ in (7). First, we denote the weighted sum of the product
of random variables in (7) by
$Y=\sum_{n\in{\mathcal{N}}}\beta_{n}\lambda_{{g}_{n}}\lambda_{{h}_{n}}$. Then,
we use the fact that $\lambda_{{g}_{n}}$ and $\lambda_{{h}_{n}}$ for
$n\in\mathcal{N}$ are independently distributed Rayleigh random variables to
tightly approximate $Y$ through an one-sided Gaussian distributed random
variable $(\tilde{Y})$ by invoking the CLT [15] as [5]

$\displaystyle\\!\\!\\!\\!\\!f_{Y}(y)\approx
f_{\tilde{Y}}(y)=\frac{\psi}{\sqrt{2\pi\sigma_{Y}^{2}}}\mathrm{exp}\left(\\!\frac{-(y-\mu_{Y})^{2}}{2\sigma_{Y}^{2}}\\!\right),\,\text{for}\,\,y\geq
0,$ (8)

where $\psi\triangleq 1/\mathcal{Q}\left(-\mu_{Y}/\sigma_{Y}\right)$ is a
normalization factor, which is used to ensure that
$\int_{-\infty}^{\infty}f_{\tilde{Y}}(x)dx=1$, and $\mathcal{Q}(\cdot)$ is the
Gaussian-$\mathcal{Q}$ function [15]. In (8), $\mu_{Y}$ and $\sigma_{Y}^{2}$
are given by

$\displaystyle\mu_{Y}$ $\displaystyle=$
$\displaystyle\sum\nolimits_{n\in{\mathcal{N}}}\pi\beta_{n}\left(\xi_{g_{n}}\xi_{h_{n}}\right)^{1/2}/2,$
(9a) $\displaystyle\sigma_{Y}^{2}$ $\displaystyle=$
$\displaystyle\sum\nolimits_{n\in{\mathcal{N}}}\beta_{n}^{2}\xi_{g_{n}}\xi_{h_{n}}\left(16-\pi^{2}\right)/4.$
(9b)

Next, we derive a tight approximation for the PDF of $R=\lambda_{u}+Y$ as (see
Appendix A)

$\displaystyle\\!\\!\\!\\!f_{R}(x)\\!$ $\displaystyle\approx$
$\displaystyle\\!f_{\tilde{R}}(x)\\!=\\!\sqrt{\pi}\rho\left(\frac{x-\mu_{Y}}{2\sigma_{Y}^{2}\sqrt{a}}\right)\mathrm{exp}\left(-\Delta\left(\frac{x-\mu_{Y}}{2\sigma_{Y}^{2}\sqrt{a}}\right)^{2}\right)$
(10)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\times\left(\mathrm{erf}\left(\frac{x-\mu_{Y}}{2\sigma_{Y}^{2}\sqrt{a}}\right)+1\right)+\rho\mathrm{exp}\left(-\left(\frac{x-\mu_{Y}}{2\sigma_{Y}^{2}}\right)^{2}\right)\\!,$

where $\mathrm{erf}\left(x\right)=2/\sqrt{\pi}\int_{0}^{x}\mathrm{e}^{-t}dt$
is the error function [17, Eqn. 8.250.1]. Here, $a$, $\rho$, and $\Delta$ are
given by

$\displaystyle\\!\\!\\!\\!\\!\\!a$ $\displaystyle=$
$\displaystyle{1}/{2\xi_{u}}+{1}/{2\sigma_{Y}^{2}},\qquad\rho={\psi}\Big{/}\left({2a\xi_{u}\sqrt{2\pi\sigma_{Y}^{2}}}\right),$
(11a) $\displaystyle\\!\\!\\!\\!\\!\\!\Delta$ $\displaystyle=$
$\displaystyle\left(1-{1}/{2\sigma_{Y}^{2}}\right)2\sigma_{Y}^{2}a.$ (11b)

In particular, (10) serves as the exact PDF of
$\tilde{R}=\lambda_{u}+\tilde{Y}$, where $\tilde{Y}$ is the one-sided Gaussian
approximated random variable for $Y$ in (7). Then, we derive an approximated
PDF for $\gamma^{*}=\bar{\gamma}R^{2}$ as

$\displaystyle f_{\gamma^{*}}(y)$ $\displaystyle\approx$ $\displaystyle
f_{\tilde{R}}\left(\sqrt{{y}/{\bar{\gamma}}}\right)\times{1}\big{/}{2\sqrt{\bar{\gamma}y}}.$
(12)

Specifically, (12) serves as the the exact PDF of
$\gamma^{*}\approx\tilde{\gamma}^{*}=\bar{\gamma}\tilde{R}^{2}$. From (10), we
derive the CDF of $\tilde{R}$ as (see Appendix B)

$\displaystyle F_{\tilde{R}}(x)$ $\displaystyle=$ $\displaystyle
1-\int_{x}^{\infty}f_{\tilde{R}}(u)du=1-\left(I_{a}+I_{b}\right),$ (13)

where $I_{a}$ and $I_{b}$ are given by

$\displaystyle\\!\\!\\!\\!\\!\\!\\!I_{a}$ $\displaystyle=$
$\displaystyle\frac{\lambda\mathrm{e}^{-\Delta
d}\left(\mathrm{erf}\left(d+1\right)\right)}{2\Delta}+\frac{\lambda\left(1-\mathrm{erf}\left(d\sqrt{\Delta+1}\right)\right)}{2\Delta\sqrt{\Delta+1}},$
(14a) $\displaystyle\\!\\!\\!\\!\\!\\!\\!I_{b}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{\pi\sigma_{Y}^{2}}{2}}\rho\left(1-\mathrm{erf}\left(\sqrt{2\sigma_{Y}a}d\right)\right),$
(14b)

where $\lambda=2\sigma_{Y}^{2}\rho\sqrt{\pi a}$, $\rho$ is given in (11a), and
$d=(x-\mu_{Y})/(2\sigma_{Y}^{2}\sqrt{a})$. From (13), we approximate the CDF
of $\gamma^{*}=\bar{\gamma}R^{2}$ as

$\displaystyle F_{\gamma^{*}}(y)$ $\displaystyle=$
$\displaystyle\mathrm{Pr}\left(\gamma^{*}\leq y\right)\approx
F_{\tilde{R}}\left(\sqrt{y/\bar{\gamma}}\right).$ (15)

Figure 2: PDF and CDF of SNR ($\gamma^{*}$) for $\bar{\gamma}=-10$dB. The
combinations of $M$ and $N$ for Case 1 to Case 4 are set to $\\{M=64,N=32\\}$,
$\\{M=64,N=64\\}$, $\\{M=144,N=64\\}$, and $\\{M=64,N=128\\}$.

Remark 1: We plot the approximated PDF and CDF of $\gamma^{*}$ by using the
analysis in (12) and (15), respectively, in Fig. 2. Monte-Carlo simulations
are also plotted in the same figure for various $M$ and $N$ to verify the
accuracy of our approximations. From Fig. 2, we observe that our analytical
approximations for the PDF (12) and CDF (15) of $\gamma^{*}$ are accurate even
for moderately large values for $M$ and $N$.

$\displaystyle\mathcal{R}_{lb}=\mathrm{log}_{2}\left(1+\frac{\bar{\gamma}\left(\xi_{u}+\sigma_{Y}^{2}+2\mu_{u}\mu_{Y}+\mu_{u}^{2}+\mu_{Y}^{2}\right)^{3}}{\sum\nolimits_{n\in
C_{4}}\binom{4}{n}\left(2\xi_{u}\right)^{n/2}\Gamma\left(n/2+1\right)\frac{\psi}{2\sqrt{\pi}}\sum\nolimits_{i\in
C_{n}}\binom{n}{i}\left({2\sigma_{Y}^{2}}\right)^{(n-i)/2}\mu_{Y}^{i}I\left(n-i,\frac{-\mu_{Y}}{2\sigma_{Y}^{2}}\right)}\right)$
(17)
$\displaystyle\hat{\mathcal{R}}_{ub}=\mathrm{log}_{2}\left(1+\bar{\gamma}\left(\xi_{u}+{(\mu_{Y}\sin(\tau))}/{\tau}\left[2\mu_{u}+{(\mu_{Y}\sin(\tau))}/{\tau}\right]+{4\sigma_{Y}^{2}}/({16-\pi^{2}})\left[4-{\pi^{2}\sin(\tau)^{2}}/{(4\tau^{2})}\right]\right)\right)$
(19)

##  IV Performances Analysis

###  IV-A Outage probability

An outage event occurs when the optimal received SNR (7) falls below a
threshold SNR ($\gamma_{th}$). To this end, we define the the outage
probability of the proposed system model as
$P_{out}=P_{r}\left(\gamma\leq\gamma_{th}\right)$. From (15), we can compute
the a tight approximation for the outage probability as $P_{out}\approx
F_{\gamma^{*}}(\gamma_{th})$.

###  IV-B Average achievable rate

The average achievable rate of the proposed system can be defined as
$\mathcal{R}=\mathbb{E}\\!\left[{\mathrm{log}_{2}\left(1+\gamma^{*}\right)}\right]$.
The exact derivation of this expectation in $\mathcal{R}$ appears
mathematically intractable. Thus, we resort to tight upper/lower bounds for
$\mathcal{R}$ as $\mathcal{R}_{lb}\lesssim\mathcal{R}\lesssim\mathcal{R}_{ub}$
by invoking the Jensen’s inequality [18]. Next, we derive $\mathcal{R}_{ub}$
as (see Appendix C)

$\displaystyle\mathcal{R}_{ub}=\mathrm{log}_{2}\left(1+\bar{\gamma}\left(\xi_{u}+\sigma_{Y}^{2}+2\mu_{u}\mu_{Y}+\mu_{u}^{2}+\mu_{Y}^{2}\right)\right).$
(16)

We derive $\mathcal{R}_{lb}$ as given in (17) at the top of the next page.

##  V Impact of discrete phase-shift adjustments

Due to the hardware limitation, the adoption of continuous phase-shift
adjustments for passive reflective elements at the IRS is practically
challenging. Thus, we investigate the feasibility of adopting discrete phase-
shifts for the proposed set-up via phase-shift quantization. It is assumed
that a limited number of discrete phase-shifts is available to select at the
$n$th reflector such that $\hat{\theta}_{n}^{*}=\pi\varsigma/2^{B-1}$, where
$B$ denotes the number of quantization bits, $\varsigma=\underset{q\in\\{0,\pm
1,\cdots,\pm 2^{B-1}\\}}{\mathrm{argmin}}|{\theta}_{n}^{*}-\pi q/2^{B-1}|$,
and $\theta_{n}^{*}$ is the optimal phase-shift in Section II-B. Then, we can
define the error of the continuous and quantized phase-shifts as
$\varepsilon_{n}={\theta}_{n}^{*}-\hat{\theta}_{n}^{*}$. For a large number of
quantization levels, $\varepsilon_{n}$ can shown to be uniformly distributed
as $\varepsilon_{n}\sim\mathcal{U}\left[-\tau,\tau\right)$ with
$\tau=\pi/2^{B}$ [19]. The signal and error $\varepsilon_{n}$ becomes
uncorrelated for a high number of quantization levels [19]. Thus, the optimal
SNR in (7) can be rewritten with discrete phase-shift as

$\displaystyle\\!\\!\\!\\!\\!\\!\hat{\gamma}^{*}\\!=\\!\bar{\gamma}\left|\lambda_{u}\\!+\\!\sum_{n\in{\mathcal{N}}}\beta_{n}\lambda_{{g}_{n}}\lambda_{{h}_{n}}\mathrm{e}^{j\varepsilon_{n}}\right|^{2}\\!=\\!\bar{\gamma}\left((\lambda_{u}\\!+\\!Y_{R})^{2}\\!+\\!Y_{I}^{2}\right)\\!,$
(18)

where
$Y_{R}=\sum_{n\in{\mathcal{N}}}\beta_{n}\lambda_{{g}_{n}}\lambda_{{h}_{n}}\cos(\varepsilon_{n})$
and
$Y_{I}=\sum_{n\in{\mathcal{N}}}\beta_{n}\lambda_{{g}_{n}}\lambda_{{h}_{n}}\sin(\varepsilon_{n})$.
By following steps similar to those in Appendix C, an upper bound for the
achievable rate with phase-shift quantization errors
$(\hat{\mathcal{R}}_{ub})$ can be derived by using (18) as shown in (19).

##  VI Numerical Results

The system parameters for our simulations are given below:
$\zeta_{v}=\left(d_{0}/d_{v}\right)^{\kappa}\times 10^{\varphi_{v}/10}$ is
used to model large-scale fading, where $v\in\\{u_{m},h_{mn},g_{n}\\}$ for
$m\in\mathcal{M}$ and $n\in\mathcal{N}$. The transmission distance between
nodes is denoted by $d_{v}$, $d_{0}=1$ m is a reference distance, the path-
loss exponent is $\kappa=2.8$, and log-normal shadow fading is captured by
$10^{\varphi_{v}/10}$ with $\varphi_{v}\sim(0,8)$ [20]. In our system
topology, the IRS and $D$ are in positioned at fixed locations and $250\,$m
apart, while the APs are uniformly distributed over an area of $1000\times
1000$ $\mathrm{m}^{2}$. The amplitudes of reflection coefficients are set to
$\beta_{n}=0.9$ for $n\in\mathcal{N}$, which is a typical assumption for IRSs
[10, 16].

Figure 3: The outage probability for different $M$ and $N$ and $\gamma_{th}=0$
dB. The combinations of $M$ and $N$ for Case-1 to Case-6 are set to
$\\{M=36,N=16\\}$, $\\{M=36,N=32\\}$, $\\{M=16,N=64\\}$, $\\{M=36,N=64\\}$,
$\\{M=64,N=64\\}$, and $\\{M=36,N=128\\}$.

In Fig. 3, we plot the outage probability as a function of the average
transmit SNR ($\bar{\gamma}$) for different combinations of distributed APs
$(M)$ and reflective elements $(N)$ at the IRS. For comparison purposes, we
also plot the outage probability for the APs-to-$D$ direct transmission
(without using an IRS) for $M=64$ in the same figure. We use our closed-form
derivation in (15) to plot the analytical outage probability approximations,
and we plot the exact counterparts through Monte-Carlo simulation. The latter
is used to verify the accuracy/tightness of our outage probability
approximations. According to Fig. 3, the tightness of our outage analysis
improves with as $M$ or/and $N$ increase. The reason for this is that large
$M$ or/and $N$ improves the accuracy of CLT. Moreover, the outage probability
can be reduced by either increasing $M$ or/and $N$. For example, at an average
SNR of $-5\,$dB, the outage probability can be reduced by $99.9$% by doubling
$N$ from $16$ (Case-1) to $32$ (Case-2) while keeping $M=36$. Moreover, by
increasing $M,N$ from $\\{M=36,N=32\\}$ in Case-2 to $\\{M=64,N=64\\}$ in
Case-4, the average SNR required to achieve an outage probability of $10^{-3}$
can be reduced by $155.6\%$ dB. From Fig. 3, we observe that the proposed IRS-
aided cell-free set-up outperforms the APs-to-$D$ direct transmission. For
instance, the set-up without IRS needs an average transmit SNR of $18$ dB to
reach an outage probability of $10^{-2}$, which is about $177.6\%$ increase
over the transmit SNR requirement for the Case-5 with the IRS-aided set-up for
the same number of APs $(M=64)$. Thus, the co-existence of IRSs within a cell-
free set-up can be beneficial in reducing the system outage probability.

Figure 4: The average achievable rate for $N\in\\{16,32,64,128,256\\}$ and
$M=64$.

In Fig. 4, we study the average achievable rate of the proposed system as a
function of the average transmit SNR ($\bar{\gamma}$) for
$N\in\\{16,32,64,128,256\\}$. We also compare the achievable rates of APs-
to-$D$ direct transmission and the IRS-aided transmission. The upper and lower
bounds for the achievable rates are plotted by using our analysis in (16) and
(17), respectively. We again validate the accuracy of our analysis through
Monte-Carlo simulations of the exact achievable rate. The tightness of our
upper/lower rate bounds is clearly depicted in enlarged portion of Fig. 4. We
observe that the rate gains can be achieved by increasing the number of
reflective elements in the IRS. Fig. 4 also illustrates that an IRS can be
embedded within a cell-free set-up to boost the achievable gains. For
instance, an IRS with $N=16$ provides a rate gain of about $180$ % compared to
the APs-to-$D$ transmission without an IRS at an average transmit SNR of $0$
dB.

Figure 5: The impact of discrete phase-shifts with phase-shift quantization on
the average achievable rate for different $M$ and $N$. The combinations of $M$
and $N$ for Case-1 to Case-4 are set to $\\{M=36,N=32\\}$, $\\{M=64,N=32\\}$,
$\\{M=36,N=64\\}$, and $\\{M=64,N=64\\}$.

In Fig. 5, we investigate the impact of discrete phase-shifts and the number
of quantization bits ($B$) by plotting the percentage rate ratio
$(\mathcal{R}_{ub}^{per})$ against the average transmit SNR for different
combinations of $M$ and $N$. The phase-shift quantization errors are uniformly
distributed: $\mathcal{U}\left[-\pi/2^{B},\pi/2^{B}\right)$. The percentage
rate ratio is defined as follows:
$\mathcal{R}_{ub}^{per}=\hat{\mathcal{R}}_{ub}/\mathcal{R}_{ub}\times 100\%$,
where $\hat{\mathcal{R}}_{ub}$ and $\mathcal{R}_{ub}$ are the upper bounds of
the average achievable rate with and without phase-shift quantization errors
given in (19) and (16), respectively. Monte-Carlo simulation curves are also
generated to validate our analysis. Fig. 5 shows that the impact of phase-
shift quantization errors vanishes when a higher $B$ is used. For instance, we
can recover more than $98\%$ of the average rate when $4$ bit quantization is
used at the IRS compared to the system with continuous phase-shift
adjustments. As per Fig. 5, $\mathcal{R}_{ub}^{per}$ improves in the high SNR
regime. For example, by varying $B$ as 1, 2, and 4 bits, the average rate can
be recovered more than $90\%$, $98\%$, and almost $100\%$, respectively, at a
transmit SNR of $20\,$dB. Fig. 5 shows that a higher number of $M,N$ is also
beneficial for recovering the achievable rate in the moderate-to-large
transmit SNR regime.

##  VII Conclusion

In this paper, the feasibility of adopting an IRS embedded within a cell-free
set-up has been explored. The optimal received SNR through multiple
distributed APs with an IRS-aided channel has been statistically characterized
by deriving the tight PDF and CDF approximations. This probabilistic SNR
analysis has been used to derive tight approximations/bounds for the outage
probability and the average achievable rate in closed-form. The impairments of
discrete phase-shifts with equalization errors have been explored. The
accuracy of our performance analysis of the proposed system set-up has been
verified by providing Monte-Carlo simulations. We observe from our numerical
results that IRS-aided cell-free system set-ups may be used to reduce the
outage probability and boost the achievable rates of next-generation wireless
systems.

##  Appendix A The derivation of PDF of $\tilde{R}$ in (10)

By using the fact that $\lambda_{u}$ and $\tilde{Y}$ are independent random
variables, we derive the PDF of $\tilde{R}$ as

$\displaystyle f_{\tilde{R}}(x)$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}f_{u}(u)f_{\tilde{Y}}(x-u)du$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!=2a\rho\mathrm{e}^{-\frac{(x-\mu_{Y})^{2}}{2\sigma_{Y}^{2}}}\int_{0}^{\infty}u\mathrm{e}^{-au^{2}+bu}du$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!=2a\rho\mathrm{e}^{-\frac{(x-\mu_{Y})^{2}}{2\sigma_{Y}^{2}}}\mathrm{e}^{\frac{b^{2}}{4a}}\int_{0}^{\infty}u\mathrm{e}^{-a\left(u-\frac{b}{2a}\right)^{2}}du$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\stackrel{{\scriptstyle(a)}}{{=}}2a\rho\mathrm{e}^{-\frac{(x-\mu_{Y})^{2}}{2\sigma_{Y}^{2}}}\mathrm{e}^{\frac{b^{2}}{4a}}\\!\\!\left(\underbrace{\int_{-b/2a}^{\infty}\\!\\!\\!t\mathrm{e}^{-at^{2}}dt}_{I_{1}}+\frac{b}{2a}\underbrace{\int_{-b/2a}^{\infty}\\!\\!\\!\mathrm{e}^{-at^{2}}dt}_{I_{2}}\right),$

where $b=(x-\mu_{Y})/\sigma_{Y}^{2}$. The step $(a)$ is obtained by letting
$t=u-b/2a$. Then, we can evaluate $I_{1}$ in (A) as

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!I_{1}$ $\displaystyle=$
$\displaystyle\int_{-b/2a}^{\infty}\\!\\!\\!t\mathrm{e}^{-at^{2}}dt\stackrel{{\scriptstyle(b)}}{{=}}\left[-\mathrm{e}^{-at^{2}}/2a\right]_{-b/2a}^{\infty}=\mathrm{e}^{-b^{2}/2a},$
(20)

where the step $(b)$ is computed by using [17, Eqn. 2.33.12]. Next, we
evaluate $I_{2}$ as

$\displaystyle I_{2}$ $\displaystyle=$
$\displaystyle\int_{-b/2a}^{\infty}\mathrm{e}^{-at^{2}}dt\stackrel{{\scriptstyle(c)}}{{=}}\left[\frac{\sqrt{\pi}\mathrm{erf}\left(\sqrt{a}t\right)}{2\sqrt{a}}\right]_{-b/2a}^{\infty}$
(21) $\displaystyle=$
$\displaystyle\frac{\sqrt{\pi}}{2\sqrt{a}}\left(1-\mathrm{erf}\left(\frac{-b}{2\sqrt{a}}\right)\right),$

where the step $(c)$ is due to [17, Eqn. 2.33.16]. We substitute (20) and (21)
into (A) to obtain the PDF of $\tilde{R}$ in (10).

##  Appendix B The derivation of CDF of $\tilde{R}$ in (13)

We substitute (10) into (13) to derive $I_{a}$ as

$\displaystyle I_{a}\\!$ $\displaystyle=$
$\displaystyle\\!\sqrt{\pi}\rho\\!\\!\\!\int_{x}^{\infty}\\!\\!\\!\left(\\!\frac{u-\mu_{Y}}{2\sigma_{Y}^{2}\sqrt{a}}\\!\right)\\!\mathrm{e}^{\\!-\Delta\left(\\!\frac{u-\mu_{Y}}{2\sigma_{Y}^{2}\sqrt{a}}\\!\right)^{\\!2}}\\!\\!\left(\\!\mathrm{erf}\left(\frac{u-\mu_{Y}}{2\sigma_{Y}^{2}\sqrt{a}}\right)\\!+\\!1\\!\right)du$
(22) $\displaystyle\stackrel{{\scriptstyle(d)}}{{=}}$
$\displaystyle\lambda\int_{d}^{\infty}t\mathrm{exp}\left(-\Delta
t^{2}\right)\left(\mathrm{erf}\left(t\right)+1\right)dt$
$\displaystyle\stackrel{{\scriptstyle(e)}}{{=}}$
$\displaystyle\lambda\left[\frac{-\mathrm{e}^{-\Delta
t^{2}}(\mathrm{erf}\left(t\right)+1)}{2\Delta}\right]_{d}^{\infty}+\lambda\int_{d}^{\infty}\frac{\mathrm{e}^{-t^{2}(\Delta+1)}}{2\Delta}dt$
$\displaystyle\stackrel{{\scriptstyle(f)}}{{=}}$
$\displaystyle\frac{\lambda\mathrm{e}^{-\Delta
d}\left(\mathrm{erf}\left(d+1\right)\right)}{2\Delta}+\frac{\lambda\left(1-\mathrm{erf}\left(d\sqrt{\Delta+1}\right)\right)}{2\Delta\sqrt{\Delta+1}},$

where $\lambda=2\sigma_{Y}^{2}\rho\sqrt{\pi a}$ and
$d=(x-\mu_{Y})/(2\sigma_{Y}^{2}\sqrt{a})$. The step $(d)$ is obtained by
through $t=(u-\mu_{Y})/2\sigma_{Y}^{2}\sqrt{a}$. The step $(e)$ is written by
invoking part-by-part integration, while the step $(f)$ is due to [17, Eqn.
2.33.16]. Next, we compute $I_{b}$ as

$\displaystyle I_{b}$ $\displaystyle=$
$\displaystyle\rho\int_{x}^{\infty}\mathrm{e}^{-\left(\frac{u-\mu_{Y}}{2\sigma_{Y}^{2}}\right)^{2}}du\stackrel{{\scriptstyle(g)}}{{=}}\sqrt{2\sigma_{Y}^{2}}\rho\int_{\sqrt{2\sigma_{Y}a}d}^{\infty}\mathrm{e}^{-t^{2}}dt$
(23) $\displaystyle\stackrel{{\scriptstyle(h)}}{{=}}$
$\displaystyle\sqrt{\frac{\pi\sigma_{Y}^{2}}{2}}\rho\left(1-\mathrm{erf}\left(\sqrt{2\sigma_{Y}a}d\right)\right),$

where the step $(g)$ is due to a changing of dummy variable as
$t=(u-\mu_{Y})/(2\sigma_{Y}^{2})$, and the step $(h)$ is resulted due to [17,
Eqn. 2.33.16].

##  Appendix C The derivation of $\mathcal{R}_{lb}$ and $\mathcal{R}_{ub}$ in
(17) and (16)

First, by invoking Jensen’s inequality, $\mathcal{R}_{lb}$ and
$\mathcal{R}_{ub}$ can be defined as

$\displaystyle\mathcal{R}_{lb}$ $\displaystyle=$
$\displaystyle\mathrm{log}_{2}\left(1+\left(\mathbb{E}\\!\left[{1/\tilde{\gamma}^{*}}\right]\right)^{-1}\right),$
(24a) $\displaystyle\mathcal{R}_{ub}$ $\displaystyle=$
$\displaystyle\mathrm{log}_{2}\left(1+\mathbb{E}\\!\left[{\tilde{\gamma}^{*}}\right]\right).$
(24b)

Then, we evaluate the expectation term in (24b) as

$\displaystyle\mathbb{E}\\!\left[{\tilde{\gamma}^{*}}\right]$ $\displaystyle=$
$\displaystyle\mathbb{E}\\!\left[{\bar{\gamma}\tilde{R}^{2}}\right]=\bar{\gamma}\mathbb{E}\\!\left[{(\lambda_{u}+\tilde{Y})^{2}}\right]$
(25) $\displaystyle=$ $\displaystyle\bar{\gamma}\sum_{n\in
C_{2}}\\!\\!\binom{2}{n}\mathbb{E}\\!\left[{\lambda_{u}^{(2-n)}}\right]\mathbb{E}\\!\left[{\tilde{Y}^{n}}\right]$
$\displaystyle=$
$\displaystyle\bar{\gamma}\left(\xi_{u}+\mu_{u}^{2}+\sigma_{Y}^{2}+\mu_{Y}^{2}+2\mu_{u}\mu_{Y}\right),$

where $\mu_{u}=\sqrt{\pi\xi_{u}/2}$. Moreover, $\mu_{Y}$ and $\sigma_{Y}^{2}$
are given in (9a) and (9b), receptively. By substituting (25) into (24b),
$\mathcal{R}_{ub}$ can be computed as (16). Next, we can write the expectation
term in (24a) as

$\displaystyle\mathbb{E}\\!\left[{1/\tilde{\gamma}^{*}}\right]={1}/{\mathbb{E}\\!\left[{\tilde{\gamma}^{*}}\right]}+{\mathbb{V}\mathrm{ar}\\!\left[{\tilde{\gamma}^{*}}\right]}/{\left(\mathbb{E}\\!\left[{\tilde{\gamma}^{*}}\right]\right)^{3}},$
(26)

where $\mathbb{E}\\!\left[{\tilde{\gamma}^{*}}\right]$ is defined in (25) and
$\mathbb{V}\mathrm{ar}\\!\left[{\tilde{\gamma}^{*}}\right]=\bar{\gamma}^{2}\mathbb{E}\\!\left[{\tilde{R}^{4}}\right]-\left(\mathbb{E}\\!\left[{\tilde{\gamma}^{*}}\right]\right)^{2}$.
Then, we can compute $\mathbb{E}\\!\left[{\tilde{R}^{4}}\right]$ as follows:

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\mathbb{E}\\!\left[{\tilde{R}^{4}}\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}\\!\left[{(\lambda_{u}+\tilde{Y})^{4}}\right]\\!=\sum_{n\in
C_{4}}\\!\binom{4}{n}\mathbb{E}\\!\left[{\lambda_{u}^{(4-n)}}\right]\mathbb{E}\\!\left[{\tilde{Y}^{n}}\right],$
(27)

where the $n$th moment of $\lambda_{u}^{n}$ is denoted by
$\mathbb{E}\\!\left[{\lambda_{u}^{n}}\right]$. We compute
$\mathbb{E}\\!\left[{\lambda_{u}^{n}}\right]$ as

$\displaystyle\mathbb{E}\\!\left[{\lambda_{u}^{n}}\right]$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}x^{n}f_{u}(x)dx=\int_{0}^{\infty}\frac{x^{n+1}}{\xi_{u}}\mathrm{exp}\left(-\frac{x^{2}}{2\xi_{u}}\right)dx$
(28) $\displaystyle\stackrel{{\scriptstyle(i)}}{{=}}$
$\displaystyle\left(2\xi_{u}\right)^{n/2}\Gamma\left(n/2+1\right),$

where the step $(m)$ is evaluated from [17, Eqn. 2.33.10] and
$\Gamma(t)=\int_{0}^{\infty}x^{t}\mathrm{e}^{-x}dx$ is the Gamma function [17,
Eqn. 8.310.1]. Then, we evaluate $\mathbb{E}\\!\left[{\tilde{Y}^{n}}\right]$
for $n\in C_{4}^{\prime}$ as

$\displaystyle\mathbb{E}\\!\left[{\tilde{Y}^{n}}\right]$ $\displaystyle=$
$\displaystyle\frac{\psi}{\sqrt{2\pi\sigma_{Y}^{2}}}\int_{0}^{\infty}y^{n}\mathrm{e}^{-\frac{(y-\mu_{Y})^{2}}{2\sigma_{Y}^{2}}}dy$
$\displaystyle\stackrel{{\scriptstyle(j)}}{{=}}$
$\displaystyle\frac{\psi}{\sqrt{\pi}}\int_{{-\mu_{Y}}/{\sqrt{2\sigma_{Y}^{2}}}}^{\infty}\left(\sqrt{2\sigma_{Y}^{2}}t+\mu_{Y}\right)^{n}\mathrm{e}^{-t^{2}}dt$
$\displaystyle\stackrel{{\scriptstyle(k)}}{{=}}$
$\displaystyle\\!\frac{\psi}{2\sqrt{\pi}}\sum\limits_{i\in
C_{n}}\\!\\!\binom{n}{i}\left({2\sigma_{Y}^{2}}\right)^{\frac{n-i}{2}}\mu_{Y}^{i}I\\!\left(n\\!-\\!i,\frac{-\mu_{Y}}{2\sigma_{Y}^{2}}\right),$

where the step $(j)$ is due to a changing of the dummy variable, the step
$(k)$ is obtained by expanding
$\left(\\!\sqrt{2\sigma_{Y}^{2}}t\\!+\\!\mu_{Y}\\!\right)^{\\!\\!n}\\!$ based
on $n\\!$ value. Moreover, $\\!I(\\!\cdot,\cdot\\!)\\!$ is given as

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!I\\!\left(m,t\right)\\!=\\!\begin{cases}(-1)^{m}\gamma\left(\frac{m+1}{2},t^{2}\right)+\Gamma\left(\frac{m+1}{2}\right),&\text{for}\,\,t\leq
0,\\\ \Gamma\left(\frac{m+1}{2},t^{2}\right),&\text{otherwise},\end{cases}$
(30)

where $\gamma(\lambda,x)=\int_{0}^{x}\mathrm{e}^{-t}t^{\lambda-1}dt$ is the
lower incomplete Gamma function [17, Eqn. 8.350.1]. Finally,
$\mathcal{R}_{lb}$ is derived as (17).

## References

  * [1] H. Q. Ngo _et al._ , “Cell-Free Massive MIMO: Uniformly Great Service for Everyone,” in _IEEE 16th Int. Workshop on Signal Process. Adv. in Wireless Commun. (SPAWC)_ , June 2015, pp. 201–205.
  * [2] ——, “Cell-Free Massive MIMO versus Small Cells,” _IEEE Trans. Wireless Commun._ , vol. 16, no. 3, pp. 1834–1850, Mar. 2017.
  * [3] M. D. Renzo _et al._ , “Smart Radio Environments Empowered by Reconfigurable AI Meta-Surfaces: An idea whose time has come,” _EURASIP J. Wireless Commun. Net._ , May 2019.
  * [4] C. Liaskos _et al._ , “A New Wireless Communication Paradigm through Software-Controlled Metasurfaces,” _IEEE Commun. Mag._ , vol. 56, no. 9, pp. 162–169, 2018.
  * [5] D. L. Galappaththige, D. Kudathanthirige, and G. Amarasuriya Aruma Baduge, “Performance Analysis of Distributed Intelligent Reflective Surface Aided Communications,” in _IEEE Global Commun. Conf. (GLOBECOM)_ , May 2020, pp. 1–6, (submitted).
  * [6] H. Q. Ngo _et al._ , “On the Total Energy Efficiency of Cell-Free Massive MIMO,” _IEEE Trans. Green Commun. Netw._ , vol. 2, no. 1, pp. 25–39, 2018\.
  * [7] E. Nayebi _et al._ , “Precoding and Power Optimization in Cell-Free Massive MIMO Systems,” _IEEE Trans. Wireless Commun._ , vol. 16, no. 7, pp. 4445–4459, July 2017.
  * [8] D. L. Galappaththige and G. Amarasuriya, “Cell-Free Massive MIMO with Underlay Spectrum-Sharing,” in _IEEE Int. Conf. on Commun. (ICC)_ , 2019, pp. 1–7.
  * [9] Ö. Özdogan, E. Björnson, and E. G. Larsson, “Intelligent Reflecting Surfaces: Physics, Propagation, and Pathloss Modeling,” _IEEE Wireless Commun. Lett._ , vol. 9, no. 5, pp. 581–585, 2020.
  * [10] Q. Wu and R. Zhang, “Intelligent Reflecting Surface Enhanced Wireless Network via Joint Active and Passive Beamforming,” _IEEE Trans. Wireless Commun._ , pp. 1–1, 2019.
  * [11] Y. Han _et al._ , “Large Intelligent Surface-Assisted Wireless Communication Exploiting Statistical CSI,” _IEEE Trans. Veh. Technol._ , vol. 68, no. 8, pp. 8238–8242, Aug 2019.
  * [12] J. Chen _et al._ , “Intelligent Reflecting Surface: A Programmable Wireless Environment for Physical Layer Security,” _IEEE Access_ , vol. 7, pp. 82 599–82 612, 2019.
  * [13] S. Abeywickrama, R. Zhang, and C. Yuen, “Intelligent Reflecting Surface: Practical Phase Shift Model and Beamforming Optimization,” in _IEEE Int. Conf. on Commun. (ICC)_ , 2020, pp. 1–6.
  * [14] Z. Ding and H. Vincent Poor, “A Simple Design of IRS-NOMA Transmission,” _IEEE Commun. Lett._ , vol. 24, no. 5, pp. 1119–1123, 2020\.
  * [15] A. Papoulis and S. U. Pillai, _Probability, Random Variables, and Stochastic Processes_ , 4th ed. McGraw Hill, 2002.
  * [16] Q. Wu and R. Zhang, “Towards Smart and Reconfigurable Environment: Intelligent Reflecting Surface Aided Wireless Network,” _IEEE Commun. Mag._ , vol. 58, no. 1, pp. 106–112, 2020.
  * [17] I. Gradshteyn and I. Ryzhik, _Table of Integrals, Series, and Products_ , 7th ed. Academic Press, 2007\.
  * [18] Q. Zhang _et al._ , “Power Scaling of Uplink Massive MIMO Systems with Arbitrary-Rank Channel Means,” _IEEE J. Sel. Areas Signal Process._ , vol. 8, no. 5, pp. 966–981, Oct. 2014.
  * [19] S. Haykin and M. Moher, _Communication Systems_ , 5th ed. Wiley India Pvt. Limited, 2009.
  * [20] T. L. Marzetta _et al._ , _Fundamentals of Massive MIMO_. Cambridge University Press, Cambridge, UK, 2016\.


# Probing criticality with deep learning in relativistic heavy-ion collisions

Yige Huang  Key Laboratory of Quark and Lepton Physics (MOE) & Institute of
Particle Physics,Central China Normal University, Wuhan 430079, China  Long-
Gang Pang<EMAIL_ADDRESS>Key Laboratory of Quark and Lepton Physics (MOE)
& Institute of Particle Physics,Central China Normal University, Wuhan 430079,
China  Xiaofeng Luo<EMAIL_ADDRESS>Key Laboratory of Quark and Lepton
Physics (MOE) & Institute of Particle Physics,Central China Normal University,
Wuhan 430079, China  Xin-Nian Wang<EMAIL_ADDRESS>Key Laboratory of Quark and
Lepton Physics (MOE) & Institute of Particle Physics,Central China Normal
University, Wuhan 430079, China  Nuclear Science Division, Lawrence Berkeley
National Laboratory, Berkeley, CA 94720, USA

###### Abstract

Systems with different interactions could develop the same critical behaviour
due to the underlying symmetry and universality. Using this principle of
universality, we can embed critical correlations modeled on the 3D Ising model
into the simulated data of heavy-ion collisions, hiding weak signals of a few
inter-particle correlations within a large particle cloud. Employing a point
cloud network with dynamical edge convolution, we are able to identify events
with critical fluctuations through supervised learning, and pick out a large
fraction of signal particles used for decision-making in each single event.

##  I Introduction

Quantum Chromodynamics (QCD) is the fundamental theory of the strong
interaction. Exploring the phase structure of strongly interacting QCD matter
is one of the main goals of heavy-ion collision experiment Fukushima and
Hatsuda (2011); Bzdak  _et al._ (2020); Luo and Xu (2017). Lattice QCD Aoki
_et al._ (2009); Ding  _et al._ (2019, 2015) predicts a smooth crossover
transition from normal hadronic phase to Quark-Gluon Plasma (QGP) around
temperature $T_{c}$=156 MeV at vanishing baryon chemical potential ($\mu_{B}$
= 0 MeV). At finite baryon density region, QCD-based models calculations Shi
_et al._ (2014); Gao and Liu (2016); Fischer (2019); Fu  _et al._ (2020)
indicate that there is a possible QCD critical point (CP), which is the end
point of the first-order phase transition boundary between the hadronic matter
and QGP.

Searching for the CP is one of the most important goals in beam energy scan
(BES) program at the Relativistic Heavy-ion Collider (RHIC) Fukushima and
Hatsuda (2011); Bzdak  _et al._ (2020); Luo and Xu (2017). Many theoretical
and experimental efforts have been made to locate the CP Stephanov (2004,
2006); Luo and Xu (2017). One avenue is to classify the smooth crossover and
first order phase transition using the information from the final state
particle spectra and collective flow Hofmann  _et al._ (1976); Stoecker and
Greiner (1986); Brachmann  _et al._ (2000a, b); Csernai and Rohrich (1999);
Ivanov  _et al._ (2002); Rischke  _et al._ (1995); Stoecker (2005); Csernai
_et al._ (2005); Nara  _et al._ (2017, 2018a, 2018b); Paech  _et al._ (2003).
This method looks for the consequences of the softening of the equation of
state since the pressure gradients are much smaller in a medium with a first
order phase transition than a smooth crossover transition, which leads to
slower fluid acceleration and smaller transverse momenta of final state
particles. Another avenue is to search for the enhanced fluctuations when the
system goes through the critical point. These includes, for example,
fluctuations of conserved charges Stephanov (2009, 2011); Aggarwal  _et al._
(2010); Adamczyk  _et al._ (2014a, b, 2018); Adam  _et al._ (2021); Abdallah
_et al._ (2021), hydrodynamic fluctuations Nahrgang  _et al._ (2011); Herold
_et al._ (2013); Plumberg and Kapusta (2017), fluctuations caused by spinodal
instabilities Li and Ko (2016); Scavenius  _et al._ (2001); Palhares and Fraga
(2010); Herold  _et al._ (2014); Li and Ko (2017); Chomaz  _et al._ (2004);
Randrup (2004); Sasaki  _et al._ (2007); Steinheimer and Randrup (2012, 2013);
Steinheimer  _et al._ (2014) and enhanced light nuclei yield ratio due to
baryon density fluctuations Sun  _et al._ (2018); Yu  _et al._ (2020); Sun
_et al._ (2021); Zhao  _et al._ (2021).

Many critical phenomena in systems with different interactions can develop the
same critical behaviour with a universality that is dictated by the symmetry
of the systems and can be described by same critical exponents Wilson and
Kogut (1974). Lee and Yang proved that the Ising model in a magnetic field and
a lattice gas are mathematically equivalent Lee and Yang (1952). Employing
this universality, one can therefore map the QCD equation of state to that
given by a 3-dimensional Ising model with the same universality class Lee and
Yang (1952); Stephanov (2004); Pradeep and Stephanov (2019); Karthein  _et
al._ (2021); Teaney (2021); Bluhm  _et al._ (2020) to study the QCD phase
diagram. The divergence of the correlation length near the critical point will
lead to the critical opalescence and scaling invariant, which means that the
systems are self-similar when the resolution changes. One thus expects that
particles from the freeze-out hyper-surface close to the critical point have
multi-particle fractal structure in the momentum space Bialas and Peschanski
(1988); Satz (1989); Hwa (1990); Antoniou  _et al._ (2001); Wu  _et al._
(2020). Experimentally, intermittency analysis has been proposed to probe the
self-similarity and density fluctuations in heavy-ion collisions. Though a
non-trivial intermittency phenomenon is observed recently by the NA61/SHINE
experiment at CERN SPS Anticic  _et al._ (2015); Davis (2020); Davis  _et al._
(2019) in Ar+Sc collisions at 150 AGeV, the magnitude of background
fluctuations is big and the power law scaling is not fully established. No
intermittency signal is observed in C+C, Pb+Pb and Be+Be collisions with
similar collision energies. Critical Monte Carlo simulations suggest a maximum
critical proton fraction smaller than $0.3$% in Be+Be collision, indicating
that traditional intermittency analysis may fail in looking for the weak
signal of self-similarity, if the fraction of CMC particless is small compared
with uncorrelated background . It is interesting to explore whether the state-
of-the-art deep learning can help to identify the weak intermittency signal
from each event of heavy ion collisions.

Recently deep learning has been used to study the QCD equation of states by
classifying phase transition types, using convolution neural network Pang  _et
al._ (2018); Pang (2021); Du  _et al._ (2020); Kvasiuk  _et al._ (2020) and
point cloud network Steinheimer  _et al._ (2019); Kuttan  _et al._ (2020). In
heavy ion collisions at low energies, auto-encoder with a single latent
variable is also used to study the order parameter of the nuclear liquid-gas
phase transition Wang  _et al._ (2020). In these studies, deep learning is
powerful in mapping momentum or charge distributions of particles to the type
of QCD phase transitions. In this study, we will train a dynamical edge
convolution network plus a point cloud network to identify weak intermittency
signals of critical fluctuations, from exotic uncorrelated background
particles. Employing Critical Monte Carlo (CMC) Antoniou  _et al._ (2001); Wu
_et al._ (2020), we encode the self-similarity in the inter-particle distances
in momentum space. Further, we assume that only a small fraction of particles
have intermittency which does not change the single particle distribution.

This paper is organized as follows. In Sec.II, we present the JAM transport
model which is used to generate data on multiple particle production in heavy
ion collisions. The CMC is used to generate intermittency signals of critical
fluctuations and the deep neural network is used for both classification and
tagging. In Sec. III, the prediction accuracy is compared for point cloud
network and dynamical edge convolution neural network. We also show the
performance of signal-particle tagging. In Sec. IV, we discuss and summarize
the findings and the implications of the present work.

##  II Method

Probing critical fluctuations in heavy-ion collisions is a typical inverse
problem. The information of criticality should be transmitted through the
dynamical evolution of the dense medium in heavy-ion collisions and get
encoded in the final state hadrons that are recorded by detectors. In the
forward process, relativistic hydrodynamics as well as hadronic transport
model are widely used to generate single particle distribution and multi-
hadron correlations. In the present study, we use a hadronic transport model
JAM Nara  _et al._ (2000); Nara (2019) to generate background events without
critical fluctuations. On the other hand, to introduce critical fluctuations,
the so called Critical Monte-Carlo (CMC) model Antoniou  _et al._ (2001); Wu
_et al._ (2020) is applied to generate a series of correlated particle
momentum, which will be used to replace the momentum of particles in JAM
events.

In the inverse process, a point cloud network and a dynamical edge convolution
network are trained to identify critical fluctuations from large amount of
uncorrelated background particles. The traditional intermittency analysis is
also carried out to probe the encoded critical signals in the JAM events and
validate the effectiveness of the deep learning method.

###  II.1 The JAM and Critical Monte-Carlo model

JAM model is a hadronic transport model to simulate heavy-ion collisions Sorge
(1995, 1997); Bass  _et al._ (1998); Bleicher  _et al._ (1999); Kahana  _et
al._ (1996); Li and Ko (1998); Lin  _et al._ (2005); Nara  _et al._ (2000);
Nara (2019); Weil  _et al._ (2016). It simulates the complicated process from
initial stage nuclear collisions to multiple particle production and final
state hadronic interactions. Independent binary collisions among hadrons
including produced ones are modeled using the vacuum hadron-hadron scattering
cross section. In the present study, the mean field mode of JAM model is used
to generate background events without including the critical fluctuations.

To simulate events involving critical fluctuations, Critical Monte-Carlo (CMC)
model Antoniou  _et al._ (2001, 2006); Wu  _et al._ (2020) is used to generate
a series of correlated particle momentum according to a power law function:

$f(\Delta p)=A\Delta p^{-\alpha}$ (1)

where $\Delta p$ is the distance of two CMC particles along an axis in
momentum space. $\nu=1/6$ is an index related to the universality class of
Ising model, and we let $\alpha=1+\nu$. $a$ and $b$ are the minimum and
maximum of $\Delta p$, and in out study, we set $a=2\times
10^{-7}\mathrm{GeV/c}$ and $b=2\mathrm{GeV/c}$. $A=(\nu
a^{\nu}b^{\nu})/(b^{\nu}-a^{\nu})$, is the normalization coefficient which is
independent of $\Delta p$. In this study, we only consider 2D momentum space
($p_{y},p_{y}$). The Levy flight random walk algorithm proposes the next step
with strides respecting the distribution $f(\Delta p)=A\Delta p^{-\alpha}$ for
$\Delta p_{x}$ and $\Delta p_{y}$ independently, and in this way, two sequence
of $p_{x}$ and $p_{y}$ of CMC particles are generated whose adjacent
differences $\Delta p$ obey the power law distribution. The self-similarity or
intermittency is thus encoded in these CMC particles, which is related to the
observed large local density fluctuations associated with the critical point.

For such a probability density function $f(\Delta p)=A\Delta p^{-1-\nu}$
within a range of (a, b), it is possible to derive its cumulative distribution
function:

$F(\Delta p)=\frac{b^{\nu}(\Delta p^{\nu}-a^{\nu})}{\Delta
p^{\nu}(b^{\nu}-a^{\nu})}$ (2)

where $F(\Delta p)$ is the cumulative distribution function of random variable
$\Delta p$, $F(\Delta p)=\int_{a}^{b}{f(\Delta p)\mathrm{d}\Delta p}$. And one
can then calculate the inverse function of $F(\Delta p)$:

$\Delta p(F)=(\frac{a^{\nu}b^{\nu}}{b^{\nu}-b^{\nu}F+a^{\nu}F})^{1/\nu}$ (3)

By randomly picking up a $F$ respecting to uniform distribution between 0 and
1, and using Eq. 3, one can obtain a $\Delta p$.

###  II.2 Data set preparation

We generate about $2.2\times 10^{5}$ events of Au+Au central collisions at
$\sqrt{s_{\mathrm{NN}}}$ = 27 GeV with impact parameters $b<3\ \mathrm{fm}$.
Each event consists of hundreds of charged particles including pion, kaon and
proton. The transverse momentum $p_{x}$ and $p_{y}$ are considered as two
features of each particle. Therefore, each event stores one particle cloud in
2-dimensional momentum space. $2\times 10^{5}$ events are used to form the
training set, while the number of events for validation and test are $1\times
10^{3}$ and $2\times 10^{4}$, respectively. For each JAM event, a
corresponding CP event is created that encodes the critical fluctuation
signals from CMC model. As a result, $4.4\times 10^{5}$ events in total are
used in our study. To avoid data pollution, event with critical fluctuations
and its corresponding JAM event are always put in the same data category. In
this case, if one JAM event is in the training data, the event with critical
fluctuations associated with that JAM event is also put in the training data.
We will refer to these events with critical fluctuations as CP events and
these particles encoded with the critical fluctuations as CMC particles. Since
the CMC model only generates the momentum correlation pattern and does not
include the information of specific particle species, we don’t distinguish
between the types of particles when performing the replacement of particle in
a JAM event.

For a given JAM event, we use replacing rate $\eta=N_{CMC}/N_{JAM}$ to
describe the multiplicity ratio of CMC events to JAM events, the number of CMC
particles introduced into its corresponding CP event can reflect how strongly
the critical signal is encoded. In our study, two kinds of CP events with
$\eta=5\%$ and $\eta=10\%$, respectively, are prepared. The detailed replacing
procedures are listed below:

  1. 1.

Randomly select a particle in the chosen JAM event, use its $(p_{x},p_{y})$ as
the starting momentum for generating the CMC event.

  2. 2.

Fill a histogram $H$ of the transverse momentum distribution from the
generated CMC event. Remark the maximum magnitude of this histogram as $M$.

  3. 3.

Loop over the particles in the JAM events. For each particle, find its
corresponding $p_{T}$ bin in $H$, record the content of $H$ in the $p_{T}$ bin
as $f$.

  4. 4.

Get a random number $y$ in range from $0$ to $M$ respecting to uniform
distribution. If $y\leq f$, randomly select a CMC particle in the $p_{T}$ bin
and replace this JAM particle with it; and if $y>f$, give up this JAM particle
and go back to step 3 to find next JAM particle.

  5. 5.

Repeat step 3 to 4 until all the CMC particles are used or all the JAM
particles are looped.

By applying such algorithm, it is possible to keep the $p_{T}$ spectra of the
substituted JAM particles close to that of the introduced CMC particles, hence
the $p_{T}$ spectra of the JAM event and the corresponding CP event are quit
similar. Even if there has a fluctuation of $p_{T}$ distribution, the overall
$p_{T}$ spectrum will not be greatly affected due to the small fraction of CMC
particles (5% or 10%) in the CP event. Considering the momentum resolution of
experimental detector, we introduced a uncertainty for momentum of each
particle in JAM event with a smearing as $\delta p_{i}\approx\pm 0.05p_{i}$,
where $i=x,y$. The smearing operation will be done after the JAM and CP events
are generated.

###  II.3 Intermittency analysis

Local density fluctuations near the QCD critical point can be probed by
intermittency analysis of scaled factorial moments Wu  _et al._ (2020) in
relativistic heavy-ion collisions. The scaled factorial moments (SFM)Wu  _et
al._ (2020) are defined as follows,

$F_{q}(M)=\frac{\langle\frac{1}{M^{D}}\sum^{M^{D}}_{i=1}{n_{i}(n_{i}-1)\cdot\cdot\cdot(n_{i}-q+1)}\rangle}{\langle\frac{1}{M^{D}}\sum^{M^{D}}_{i=1}n_{i}\rangle^{q}}$
(4)

where $M$ is the number of grids in momentum space with equal size, $D$ is the
dimension, $i$ is the number of particles in the $i$th momentum-grid, and $q$
is the order of the SFM method.

When $M$ is large, the power law dependence of SFM on the number of
partitioned bins implies a self-similar correlations in the studied
systemBialas and Peschanski (1986, 1988).

$F_{q}(M)\approx(M^{D})^{\phi_{q}}$ (5)

The intermittency index $\phi_{q}$ can characterize the strength of
intermittency behavior and is related to the anomalous fractal dimension of
the systemDe Wolf  _et al._ (1996). And there are studies show that using
intermittency measurement together with the estimated freeze-out parameters
can estimate the possible critical region of the QCD CEPAntoniou and Diakonos
(2019).

Figure 1: The second order scaled factorial moments analysis for uncorrelated
JAM events and events with critical fluctuations. The upper-panel shows the
absolute values of SFM for JAM events and events with 5% and 10% CMC
particles. To avoid the overlap of markers, results of critical events are
slightly shifted horizontally for a clearer visualization. The lower-panel
shows the ratios between critical and normal JAM events. No significant
differences are observed for the absolute SFM values and their ratios.

In the present study, the second order SFM ($q=2$) in two dimensional space
($D=2$) are studied for $M=$ 2, 4, 8, 16, 32, 50. As we take the experimental
detectors into consideration, in SFM calculation, we only take no more than 50
grids for each dimension in a range of plus-minus 2.5 $\mathrm{GeV/c}$ to keep
$p_{T}$ resolution to be like experimental options and at about 0.1 GeV/c.

As shown in Figure. 1, the intermittency analysis using the SFM method Anticic
_et al._ (2015); Davis (2020); Davis  _et al._ (2019); Wu  _et al._ (2020) can
not differentiate CP events with 5% and 10% CMC particles that carry critical
fluctuations from uncorrelated JAM events.

###  II.4 Dynamical edge convolution neural network

Figure 2: Dynamical edge convolution neural network with point cloud module
for both classification and tagging. The edge convolution block looks for k
nearest neighbors of each particle to obtain a latent representation of that
very particle, with short or long range correlations encoded deeply in. The
representation of each particle are used in two tasks. One is the
classification task to identify critical fluctuations from uncorrelated
background events. The other is the tagging task to label correlated particles
used for decision making.

A graph-based dynamical edge convolution neural network is trained for our
multi-task learning. The input to the neural network are the particle cloud of
each event, which consists of a list of particles with their information on
$(p_{x},p_{y})$. The output of the neural network corresponds to two tasks.
The first task is the binary classification which requires true labels of each
single event for supervised learning, with CP indicating events with critical
fluctuations and JAM indicating events without. The second task is the
particle tagging which requires true labels of each single particle, with 0 or
1 to indicate whether the particle is generated using Critical Monte Carlo
model.

Shown in Figure. 2 is the architecture of our neural network. There are two
kNN plus dynamical edge convolution blocks connecting to the input layer. In
the first block, kNN is used to find the k-nearest neighbors of each particle
in $(p_{x},p_{y})$ space. A fully connected network is used to learn edge
features $\phi(\vec{p}_{i},\vec{p}_{j})$ between the $i$’th particle and its
$j$’th neighbor. This module is shared by all its neighbors of particle $i$ to
produced edge features and that explains the name ”edge convolution”. The
information of particle $i$ together with its edge features are feed to the
second block. Edge convolution layer would not only make use of the features
of input neuron itself, but also take the relevance between the clustered
units near that neuron into consideration, thus it can effectively capture the
correlation information between particles.

The second kNN find the k-nearest neighbors of each particle in feature space.
It is thus possible to correlate particles that are far away in momentum
space. The neighbors of each particle change dynamically when the distances
are computed in feature space, that is why the method is called ”dynamical
edge convolution”.

The features of each particle together with its ”local” information are
flattened and feed to a fully connected neural network to get a high
dimensional latent variable for each particle. The latent variable provides a
high dimensional representation of each particle. The above neural network is
also shared by all particles and is called 1D convolution neural network
(CNN). Finally, the latent variables of each particle are used for two
different tasks. The module of ”Classification” task is shown in the lower
right corner. A global max pooling gets the maximum values of each feature
among all particles. This symmetric permutation operation learns the global
feature of each particle cloud and is used to determine whether it is a CP or
JAM event. The module of ”Tagging” task is shown on the right of Figure. 2. A
1D CNN with one output neuron is used to tag each particle in the particle
cloud. This module provides interpretation on whether the correlated particles
are used to identify events with critical fluctuations. We have labeled
correlated CMC particles as ”signal” and uncorrelated JAM particles as
”noise”. Binary cross entropy is used to compute the differences between the
tagging output and the true labels of each particle. The loss values of
tagging module is added to the total loss with a weighting factor $10^{-3}$
such that the network focus more on ”classification” task.

For comparison, we also train a point-cloud network without the kNN and
dynamical edge convolution blocks shown in Figure. 2. The $(p_{x},p_{y})$ of
each particle is directly feed to 1D CNN with 256, 128 and 64 channels
respectively for classification. Global average pooling layer is used in this
simple point-cloud network as it performs better here. Without kNN and
dynamical edge convolution, the network can not capture much local information
for intermittency identification.

##  III Results and discussion

###  III.1 Classification accuracy

Shown in the Figure. 3 are the training (solid lines) and validation (dashed
lines) accuracy as a function of training epochs. Both training and validation
accuracy increase as the model is trained longer with more epochs. The
validation accuracy reaches a maximum of 99.3%, which means that deep learning
is able to classify each single event with high accuracy, for uncorrelated JAM
events and events mixed with 90% uncorrelated JAM particles and 10% CMC
particles ($\eta=10\%$). For a smaller replacing rate ($\eta=5\%$), both
validation and training accuracy decrease as compared with ($\eta=10\%$),
whose maximum value is about $93.3\%$. Note: the smeared 5% and 10% both got
93.3% acc. for validation set, while the 10% one got higher score for test
set. The validation accuracy is slightly higher than training accuracy caused
by the dropout and batch normalization layers used in the network. These two
kinds of layers are known to be able to increase the generalization of the
network by introducing noise during training.

Shown in Table. 1 are the testing accuracy of four different configurations.
Using the dynamical edge convolution plus point cloud network we constructed
in this study, the testing accuracy are $97.7\%$ for $10\%$ replacing rate and
$92.8\%$ for $5\%$ replacing rate, which are not quite far away from the
validation accuracy. Removing the dynamical edge convolution block, we have
tested the performance of the point cloud network with varying numbers of
layers and neurons per layer to get the best testing accuracy. The testing
accuracy decreases to $84.8\%$ for $10\%$ replacing rate and $83.4\%$ for
$5\%$ replacing rate.

Another test set is prepared to make sure that the network make their decision
based on multi-particle correlation in the CMC particles. In this test set, 5%
or 10% particles of a JAM event are replaced by same amount of particles
sampled randomly from many other events, one particle from each event to
eliminate the two particle correlation in the replaced particles. If our
network trained to identify CMC particles is fooled to classify these mixed
events as CMC events, it means that the network learns the missing correlation
in the replaced particles as compared with original JAM particles. In
practice, our trained network treat these mixed events as JAM events, which is
a proof that the network make their predictions using signals of CMC
particles.

Figure 3: The training and validation accuracy as a function of epochs. The
training accuracy is in solid lines, for replacing rate $5\%$ (blue) and
$10\%$(red). The validation accuracy is in dashed lines for replacing rate
$5\%$ and $10\%$.

Testing accuracy  $\eta$ Edge-Conv Point-Cloud Net $5\%$ 92.8% 83.4% $10\%$
97.7% 84.8%

Table 1: The testing accuracy for dynamical edge convolution network and a
simple point cloud network.

###  III.2 Interpretability: tagging

To figure out how does the network make its decision in identifying critical
fluctuations from the background, we have added a tagging layer to the neural
network. To quantify the tagging performance, we introduce two metrics as
follows,

$r_{\rm c}=\frac{N_{C}}{N_{C}+N_{M}},\quad\;r_{\rm
t}=\frac{N_{C}}{N_{C}+N_{W}}$ (6)

where $r_{\rm c}$ is the catching rate defined as the ratio between the number
of correctly tagged particles $N_{C}$ and total number of signal particles
$N_{C}+N_{M}$, where $N_{M}$ is the number of signal particles missed by the
tagging module. $r_{\rm t}$ is the tagging rate defined as the ratio between
the number of correctly tagged particles $N_{C}$ and the total number of
tagged particles $N_{C}+N_{W}$, where $N_{W}$ is the number of wrongly tagged
uncorrelated particles.

The average catching rates $r_{\rm c}=73.6\%$ for $\eta=5$% and $r_{\rm
c}=75.9$% for $\eta=10$% indicate that the network may use about $3/4$ of the
correlated particles to make its decision. On the other hand, the tagging rate
$r_{t}=94.5$% for $\eta=5$% and $r_{t}=95.4$% for $\eta=10$% are much higher
than catching rate $r_{c}$. This result tells us that the tagging module can
label CMC particles quite precisely.

Since both edge convolution and the following 1D convolution layers of tagging
module perform the same transformation for each particle, we can reversely
track the tensor of labeled particles in the hidden feature space in the
forward propagation process of neural network. For each input CP event, by
checking the feature space after passing edge convolution layer, for a total
of $N$ CMC particles well tagged, we find the $k$ nearest particles in the
feature space corresponding to the feature vector of each particle, and count
the number $M$ of CMC particles that were also well tagged. The proportion of
those well tagged CMC particles from kNN to the total number of these kNN
particles can then be calculated as $\frac{M}{k\times N}=94\%$. This result
indicates that, the feature space transformation guided by edge convolution
can aggregate CMC particles into a cluster in the new feature space, and then
the tagging module can label them through the subsequent 1D convolution
layers.

Figure 4: The upper subplots show the comparison of JAM event and its
corresponding CP event, in which the grey dots are the unchanged JAM
particles, and the red ones are the critical particles introduced by CMC
events. The lower subplots are labeled results of tagging network, and the red
dots refer to particles which were tagged correctly, while the blue ones are
JAM particles labeled as CMC ones, while the grey dots are unlabeled
particles. The graphs on the left show an example of $\eta=5$%, while the ones
on the right show an example of $\eta=10$%. Although the CMC clusters in the
two examples shown are all distributed on the right side of phase space, the
location of CMC particles are not restricted indeed and they can be on any
corner of the plot.

Figure 4 demonstrates the output of the tagging module. In the upper subplots,
grey dots represent unchanged JAM particles and red dots represent all the CMC
particles in two testing events. The corresponding tagging output for these
two events are shown in the two lower subplots, where the red dots represent
CMC particles correctly tagged by the network while the blue ones are JAM
particles but incorrectly tagged as CMC particles. In average, $3/4$ of CMC
particles are recognized by the tagging module. And as discussed before, the
incorrectly tagged particles are much fewer than correctly tagged CMC
particles. The two figures in the left are for $5\%$ replacing rate while the
ones on the right are for $10\%$ replacing rate.

Figure 5 shows the SFM calculation of $\eta=5\%$ CP events and the SFM of
tagged particles of them, the former ones event have no increment with the
increase of $M^{2}$ while the tagged ones present slight power law. This
result reflects that the tagging module can somehow extract the encoded
intermittency information.

Figure 5: The ’Mixed’ labeled red diamond markers represent the SFM results of
all particles from $\eta$=5% CP events, while the ’Tagged’ labeled blue square
markers stand for the SFM of tagged part of those events. As $M^{2}$ increase,
the red diamonds have a flat performance, and the blue squares show a
increment.

##  IV Summary and outlook

In summary, we have constructed a dynamical edge convolution plus point cloud
network to identify the weak intermittency signal from the experimental data
of heavy-ion collisions. We have demonstrated that such a state-of-the-art
deep learning network enables us to achieve a testing accuracy 92.8% if only
5% of JAM particles in each event are replaced by correlated CMC particles.
The performance increases to $97.7$% if the replacing rate of correlated
particles increases to 10%. Removing the dynamical edge convolution block will
decrease the performance by a large margin. Using tagging module, we further
demonstrate that the network can use around $3/4$ of correlated particles to
make their decision. At the same time, only about 5% of uncorrelated
background particles are incorrectly tagged as CMC particles.

We observe that the network can identify self-similarity or scaling invariant
from uncorrelated background. This is important for experimental data analysis
since only one indication of intermittency is observed in Ar + Sc collisions
whereas several other systems with similar collision energies fail. Different
from previous theoretical studies, we preserve the single particle
distribution while introducing a small fraction of particles with multi
particle fractal structure. This is more realistic but also difficult for the
traditional intermittency analysis. Based on our study, deep learning shows
strong pattern recognition ability in identifying weak intermittency signals
associated with critical phenomena. The method developed in this study can be
applied to probe the critical fluctuations in heavy-ion collisions and can
also be used to explore the criticality of other systems.

## Acknowledgement

We thank Jin Wu for helpful discussions on the critical monte carlo model.
This work is supported by the National Key Research and Development Program of
China (Grant No. 2020YFE0202002 and 2018YFE0205201), the National Natural
Science Foundation of China under Grant Nos. 12122505, 11935007, 11221504,
11890711, 11861131009 and 12075098, and by the Director, Office of Energy
Research, Office of High Energy and Nuclear Physics, Division of Nuclear
Physics, of the U.S. Department of Energy (DOE) under grant No. DE-
AC02-05CH11231, by the U.S. National Science Foundation under No. OAC- 2004571
within the X-SCAPE Collaboration. Computations are performed at Nuclear
Science Computer Center at CCNU (NSC3). LG Pang and YG Huang also acknowledge
the support provided by Huawei Technologies Co., Ltd.

## References

  * Fukushima and Hatsuda (2011) Kenji Fukushima and Tetsuo Hatsuda, “The phase diagram of dense QCD,” Rept. Prog. Phys. 74, 014001 (2011), arXiv:1005.4814 [hep-ph] .
  * Bzdak  _et al._ (2020) Adam Bzdak, Shinichi Esumi, Volker Koch, Jinfeng Liao, Mikhail Stephanov, and Nu Xu, “Mapping the Phases of Quantum Chromodynamics with Beam Energy Scan,” Phys. Rept. 853, 1–87 (2020), arXiv:1906.00936 [nucl-th] .
  * Luo and Xu (2017) Xiaofeng Luo and Nu Xu, “Search for the QCD Critical Point with Fluctuations of Conserved Quantities in Relativistic Heavy-Ion Collisions at RHIC : An Overview,” Nucl. Sci. Tech. 28, 112 (2017), arXiv:1701.02105 [nucl-ex] .
  * Aoki  _et al._ (2009) Y. Aoki, Szabolcs Borsanyi, Stephan Durr, Zoltan Fodor, Sandor D. Katz, Stefan Krieg, and Kalman K. Szabo, “The QCD transition temperature: results with physical masses in the continuum limit II.” JHEP 06, 088 (2009), arXiv:0903.4155 [hep-lat] .
  * Ding  _et al._ (2019) H. T. Ding _et al._ (HotQCD), “Chiral Phase Transition Temperature in ( 2+1 )-Flavor QCD,” Phys. Rev. Lett. 123, 062002 (2019), arXiv:1903.04801 [hep-lat] .
  * Ding  _et al._ (2015) Heng-Tong Ding, Frithjof Karsch, and Swagato Mukherjee, “Thermodynamics of strong-interaction matter from Lattice QCD,” Int. J. Mod. Phys. E 24, 1530007 (2015), arXiv:1504.05274 [hep-lat] .
  * Shi  _et al._ (2014) Chao Shi, Yong-Long Wang, Yu Jiang, Zhu-Fang Cui, and Hong-Shi Zong, “Locate QCD Critical End Point in a Continuum Model Study,” JHEP 07, 014 (2014), arXiv:1403.3797 [hep-ph] .
  * Gao and Liu (2016) Fei Gao and Yu-xin Liu, “QCD phase transitions via a refined truncation of Dyson-Schwinger equations,” Phys. Rev. D 94, 076009 (2016), arXiv:1607.01675 [hep-ph] .
  * Fischer (2019) Christian S. Fischer, “QCD at finite temperature and chemical potential from Dyson–Schwinger equations,” Prog. Part. Nucl. Phys. 105, 1–60 (2019), arXiv:1810.12938 [hep-ph] .
  * Fu  _et al._ (2020) Wei-jie Fu, Jan M. Pawlowski, and Fabian Rennecke, “QCD phase structure at finite temperature and density,” Phys. Rev. D 101, 054032 (2020), arXiv:1909.02991 [hep-ph] .
  * Stephanov (2004) Mikhail A. Stephanov, “QCD phase diagram and the critical point,” Prog. Theor. Phys. Suppl. 153, 139–156 (2004), arXiv:hep-ph/0402115 .
  * Stephanov (2006) M. A. Stephanov, “QCD phase diagram: An Overview,” PoS LAT2006, 024 (2006), arXiv:hep-lat/0701002 .
  * Hofmann  _et al._ (1976) J. Hofmann, Horst Stoecker, Ulrich W. Heinz, W. Scheid, and W. Greiner, “Possibility of Detecting Density Isomers in High Density Nuclear MACH Shock Waves,” Phys. Rev. Lett. 36, 88–91 (1976).
  * Stoecker and Greiner (1986) Horst Stoecker and W. Greiner, “High-Energy Heavy Ion Collisions: Probing the Equation of State of Highly Excited Hadronic Matter,” Phys. Rept. 137, 277–392 (1986).
  * Brachmann  _et al._ (2000a) J. Brachmann, S. Soff, A. Dumitru, Horst Stoecker, J. A. Maruhn, W. Greiner, L. V. Bravina, and D. H. Rischke, “Antiflow of nucleons at the softest point of the EoS,” Phys. Rev. C 61, 024909 (2000a), arXiv:nucl-th/9908010 .
  * Brachmann  _et al._ (2000b) J. Brachmann, A. Dumitru, Horst Stoecker, and W. Greiner, “The Directed flow maximum near c(s) = 0,” Eur. Phys. J. A 8, 549–552 (2000b), arXiv:nucl-th/9912014 .
  * Csernai and Rohrich (1999) L. P. Csernai and D. Rohrich, “Third flow component as QGP signal,” Phys. Lett. B 458, 454 (1999), arXiv:nucl-th/9908034 .
  * Ivanov  _et al._ (2002) Yu. B. Ivanov, E. G. Nikonov, W. Noerenberg, A. A. Shanenko, and V. D. Toneev, “Directed flow of baryons in heavy ion collisions,” Acta Phys. Hung. A 15, 117–130 (2002), arXiv:nucl-th/0011004 .
  * Rischke  _et al._ (1995) Dirk H. Rischke, Yaris Pursun, Joachim A. Maruhn, Horst Stoecker, and Walter Greiner, “The Phase transition to the quark - gluon plasma and its effects on hydrodynamic flow,” Acta Phys. Hung. A 1, 309–322 (1995), arXiv:nucl-th/9505014 .
  * Stoecker (2005) Horst Stoecker, “Collective flow signals the quark gluon plasma,” Nucl. Phys. A 750, 121–147 (2005), arXiv:nucl-th/0406018 .
  * Csernai  _et al._ (2005) L. P. Csernai, A. Anderlik, Cs. Anderlik, V. K. Magas, E. Molnar, A. Nyiri, D. Rohrich, and K. Tamosiunas, “The 3rd flow component as a QGP signal,” Acta Phys. Hung. A 22, 181–186 (2005), arXiv:hep-ph/0405277 .
  * Nara  _et al._ (2017) Yasushi Nara, Harri Niemi, Jan Steinheimer, and Horst Stöcker, “Equation of state dependence of directed flow in a microscopic transport model,” Phys. Lett. B 769, 543–548 (2017), arXiv:1611.08023 [nucl-th] .
  * Nara  _et al._ (2018a) Yasushi Nara, Harri Niemi, Akira Ohnishi, Jan Steinheimer, Xiaofeng Luo, and Horst Stöcker, “Enhancement of elliptic flow can signal a first order phase transition in high energy heavy ion collisions,” Eur. Phys. J. A 54, 18 (2018a), arXiv:1708.05617 [nucl-th] .
  * Nara  _et al._ (2018b) Yasushi Nara, Jan Steinheimer, and Horst Stoecker, “The enhancement of v4 in nuclear collisions at the highest densities signals a first-order phase transition,” Eur. Phys. J. A 54, 188 (2018b), arXiv:1809.04237 [nucl-th] .
  * Paech  _et al._ (2003) K. Paech, Horst Stoecker, and A. Dumitru, “Hydrodynamics near a chiral critical point,” Phys. Rev. C 68, 044907 (2003), arXiv:nucl-th/0302013 .
  * Stephanov (2009) M. A. Stephanov, “Non-Gaussian fluctuations near the QCD critical point,” Phys. Rev. Lett. 102, 032301 (2009), arXiv:0809.3450 [hep-ph] .
  * Stephanov (2011) M. A. Stephanov, “On the sign of kurtosis near the QCD critical point,” Phys. Rev. Lett. 107, 052301 (2011), arXiv:1104.1627 [hep-ph] .
  * Aggarwal  _et al._ (2010) M. M. Aggarwal _et al._ (STAR), “Higher Moments of Net-proton Multiplicity Distributions at RHIC,” Phys. Rev. Lett. 105, 022302 (2010), arXiv:1004.4959 [nucl-ex] .
  * Adamczyk  _et al._ (2014a) L. Adamczyk _et al._ (STAR), “Energy Dependence of Moments of Net-proton Multiplicity Distributions at RHIC,” Phys. Rev. Lett. 112, 032302 (2014a), arXiv:1309.5681 [nucl-ex] .
  * Adamczyk  _et al._ (2014b) L. Adamczyk _et al._ (STAR), “Beam-energy dependence of charge separation along the magnetic field in Au+Au collisions at RHIC,” Phys. Rev. Lett. 113, 052302 (2014b), arXiv:1404.1433 [nucl-ex] .
  * Adamczyk  _et al._ (2018) L. Adamczyk _et al._ (STAR), “Collision Energy Dependence of Moments of Net-Kaon Multiplicity Distributions at RHIC,” Phys. Lett. B 785, 551–560 (2018), arXiv:1709.00773 [nucl-ex] .
  * Adam  _et al._ (2021) J. Adam _et al._ (STAR), “Nonmonotonic Energy Dependence of Net-Proton Number Fluctuations,” Phys. Rev. Lett. 126, 092301 (2021), arXiv:2001.02852 [nucl-ex] .
  * Abdallah  _et al._ (2021) Mohamed Abdallah _et al._ (STAR), “Cumulants and correlation functions of net-proton, proton, and antiproton multiplicity distributions in Au+Au collisions at energies available at the BNL Relativistic Heavy Ion Collider,” Phys. Rev. C 104, 024902 (2021), arXiv:2101.12413 [nucl-ex] .
  * Nahrgang  _et al._ (2011) Marlene Nahrgang, Stefan Leupold, Christoph Herold, and Marcus Bleicher, “Nonequilibrium chiral fluid dynamics including dissipation and noise,” Phys. Rev. C 84, 024912 (2011), arXiv:1105.0622 [nucl-th] .
  * Herold  _et al._ (2013) Christoph Herold, Marlene Nahrgang, Igor Mishustin, and Marcus Bleicher, “Chiral fluid dynamics with explicit propagation of the Polyakov loop,” Phys. Rev. C 87, 014907 (2013), arXiv:1301.1214 [nucl-th] .
  * Plumberg and Kapusta (2017) Christopher Plumberg and Joseph I. Kapusta, “Hydrodynamic fluctuations near a critical endpoint and Hanbury-Brown–Twiss interferometry,” Phys. Rev. C 95, 044910 (2017), arXiv:1702.01368 [nucl-th] .
  * Li and Ko (2016) Feng Li and Che Ming Ko, “Spinodal instabilities of baryon-rich quark-gluon plasma in the Polyakov–Nambu–Jona-Lasinio model,” Phys. Rev. C 93, 035205 (2016), arXiv:1601.00026 [nucl-th] .
  * Scavenius  _et al._ (2001) O. Scavenius, A. Dumitru, E. S. Fraga, J. T. Lenaghan, and A. D. Jackson, “First order chiral phase transition in high-energy collisions: Can nucleation prevent spinodal decomposition?” Phys. Rev. D 63, 116003 (2001), arXiv:hep-ph/0009171 .
  * Palhares and Fraga (2010) Leticia F. Palhares and Eduardo S. Fraga, “Droplets in the cold and dense linear sigma model with quarks,” Phys. Rev. D 82, 125018 (2010), arXiv:1006.2357 [hep-ph] .
  * Herold  _et al._ (2014) Christoph Herold, Marlene Nahrgang, Igor Mishustin, and Marcus Bleicher, “Formation of droplets with high baryon density at the QCD phase transition in expanding matter,” Nucl. Phys. A 925, 14–24 (2014), arXiv:1304.5372 [nucl-th] .
  * Li and Ko (2017) Feng Li and Che Ming Ko, “Spinodal instabilities of baryon-rich quark matter in heavy ion collisions,” Phys. Rev. C 95, 055203 (2017), arXiv:1606.05012 [nucl-th] .
  * Chomaz  _et al._ (2004) Philipe Chomaz, Maria Colonna, and Jorgen Randrup, “Nuclear spinodal fragmentation,” Phys. Rept. 389, 263–440 (2004).
  * Randrup (2004) Jorgen Randrup, “Spinodal decomposition during the hadronization stage at RHIC?” Phys. Rev. Lett. 92, 122301 (2004), arXiv:hep-ph/0308271 .
  * Sasaki  _et al._ (2007) C. Sasaki, B. Friman, and K. Redlich, “Density fluctuations in the presence of spinodal instabilities,” Phys. Rev. Lett. 99, 232301 (2007), arXiv:hep-ph/0702254 .
  * Steinheimer and Randrup (2012) Jan Steinheimer and Jorgen Randrup, “Spinodal amplification of density fluctuations in fluid-dynamical simulations of relativistic nuclear collisions,” Phys. Rev. Lett. 109, 212301 (2012), arXiv:1209.2462 [nucl-th] .
  * Steinheimer and Randrup (2013) Jan Steinheimer and Jorgen Randrup, “Spinodal density enhancements in simulations of relativistic nuclear collisions,” Phys. Rev. C 87, 054903 (2013), arXiv:1302.2956 [nucl-th] .
  * Steinheimer  _et al._ (2014) Jan Steinheimer, Jørgen Randrup, and Volker Koch, “Non-equilibrium phase transition in relativistic nuclear collisions: Importance of the equation of state,” Phys. Rev. C 89, 034901 (2014), arXiv:1311.0999 [nucl-th] .
  * Sun  _et al._ (2018) Kai-Jia Sun, Lie-Wen Chen, Che Ming Ko, Jie Pu, and Zhangbu Xu, “Light nuclei production as a probe of the QCD phase diagram,” Phys. Lett. B 781, 499–504 (2018), arXiv:1801.09382 [nucl-th] .
  * Yu  _et al._ (2020) Ning Yu, Dingwei Zhang, and Xiaofeng Luo, “Search for QCD critical point by transverse velocity dependence of anti-deuteron to deuteron ratio,” Chin. Phys. C 44, 014002 (2020), arXiv:1812.04291 [nucl-th] .
  * Sun  _et al._ (2021) Kai-Jia Sun, Che Ming Ko, Feng Li, Jun Xu, and Lie-Wen Chen, “Enhanced yield ratio of light nuclei in heavy ion collisions with a first-order chiral phase transition,” Eur. Phys. J. A 57, 313 (2021), arXiv:2006.08929 [nucl-th] .
  * Zhao  _et al._ (2021) Wenbin Zhao, Kai-jia Sun, Che Ming Ko, and Xiaofeng Luo, “Multiplicity scaling of light nuclei production in relativistic heavy-ion collisions,” Phys. Lett. B 820, 136571 (2021), arXiv:2105.14204 [nucl-th] .
  * Wilson and Kogut (1974) K. G. Wilson and John B. Kogut, “The Renormalization group and the epsilon expansion,” Phys. Rept. 12, 75–199 (1974).
  * Lee and Yang (1952) T. D. Lee and Chen-Ning Yang, “Statistical theory of equations of state and phase transitions. 2. Lattice gas and Ising model,” Phys. Rev. 87, 410–419 (1952).
  * Pradeep and Stephanov (2019) Maneesha Sushama Pradeep and Mikhail Stephanov, “Universality of the critical point mapping between Ising model and QCD at small quark mass,” Phys. Rev. D 100, 056003 (2019), arXiv:1905.13247 [hep-ph] .
  * Karthein  _et al._ (2021) J. M. Karthein, D. Mroczek, A. R. Nava Acuna, J. Noronha-Hostler, P. Parotto, D. R. P. Price, and C. Ratti, “Strangeness-neutral equation of state for QCD with a critical point,” Eur. Phys. J. Plus 136, 621 (2021), arXiv:2103.08146 [hep-ph] .
  * Teaney (2021) Derek Teaney, “Dynamics of Critical Fluctuations in Nucleus-Nucleus Collisions,” Nucl. Phys. A 1005, 121750 (2021).
  * Bluhm  _et al._ (2020) Marcus Bluhm _et al._ , “Dynamics of critical fluctuations: Theory – phenomenology – heavy-ion collisions,” Nucl. Phys. A 1003, 122016 (2020), arXiv:2001.08831 [nucl-th] .
  * Bialas and Peschanski (1988) A. Bialas and Robert B. Peschanski, “Intermittency in Multiparticle Production at High-Energy,” Nucl. Phys. B 308, 857–867 (1988).
  * Satz (1989) Helmut Satz, “Intermittency and Critical Behavior,” Nucl. Phys. B 326, 613–618 (1989).
  * Hwa (1990) Rudolph C. Hwa, “Fractal Measures in Multiparticle Production,” Phys. Rev. D 41, 1456 (1990).
  * Antoniou  _et al._ (2001) N. G. Antoniou, Y. F. Contoyiannis, F. K. Diakonos, A. I. Karanikas, and C. N. Ktorides, “Pion production from a critical QCD phase,” Nucl. Phys. A 693, 799–824 (2001), arXiv:hep-ph/0012164 .
  * Wu  _et al._ (2020) Jin Wu, Yufu Lin, Yuanfang Wu, and Zhiming Li, “Probing QCD critical fluctuations from intermittency analysis in relativistic heavy-ion collisions,” Phys. Lett. B 801, 135186 (2020), arXiv:1901.11193 [nucl-th] .
  * Anticic  _et al._ (2015) T. Anticic _et al._ (NA49), “Critical fluctuations of the proton density in A+A collisions at 158$A$ GeV,” Eur. Phys. J. C 75, 587 (2015), arXiv:1208.5292 [nucl-ex] .
  * Davis (2020) Nikolaos Davis (NA61/SHINE), “Searching for the critical point of strongly interacting matter in nucleus-nucleus collisions at CERN SPS,” PoS EPS-HEP2019, 305 (2020).
  * Davis  _et al._ (2019) Nikolaos Davis, Nikolaos Antoniou, and Fotios K. Diakonos (Na61/Shine), “Recent results from proton intermittency analysis in nucleus-nucleus collisions from NA61/SHINE at CERN SPS,” PoS CORFU2018, 154 (2019).
  * Pang  _et al._ (2018) Long-Gang Pang, Kai Zhou, Nan Su, Hannah Petersen, Horst Stöcker, and Xin-Nian Wang, “An equation-of-state-meter of quantum chromodynamics transition from deep learning,” Nature Commun. 9, 210 (2018), arXiv:1612.04262 [hep-ph] .
  * Pang (2021) Long-Gang Pang, “Machine learning for high energy heavy ion collisions,” Nucl. Phys. A 1005, 121972 (2021).
  * Du  _et al._ (2020) Yi-Lun Du, Kai Zhou, Jan Steinheimer, Long-Gang Pang, Anton Motornenko, Hong-Shi Zong, Xin-Nian Wang, and Horst Stöcker, “Identifying the nature of the QCD transition in relativistic collision of heavy nuclei with deep learning,” Eur. Phys. J. C 80, 516 (2020), arXiv:1910.11530 [hep-ph] .
  * Kvasiuk  _et al._ (2020) Yu. Kvasiuk, E. Zabrodin, L. Bravina, I. Didur, and M. Frolov, “Classification of Equation of State in Relativistic Heavy-Ion Collisions Using Deep Learning,” JHEP 07, 133 (2020), arXiv:2004.14409 [nucl-th] .
  * Steinheimer  _et al._ (2019) Jan Steinheimer, Longgang Pang, Kai Zhou, Volker Koch, Jørgen Randrup, and Horst Stoecker, “A machine learning study to identify spinodal clumping in high energy nuclear collisions,” JHEP 12, 122 (2019), arXiv:1906.06562 [nucl-th] .
  * Kuttan  _et al._ (2020) Manjunath Omana Kuttan, Kai Zhou, Jan Steinheimer, Andreas Redelbach, and Horst Stoecker, “An equation-of-state-meter for CBM using PointNet,” JHEP 21, 184 (2020), arXiv:2107.05590 [hep-ph] .
  * Wang  _et al._ (2020) Rui Wang, Yu-Gang Ma, R. Wada, Lie-Wen Chen, Wan-Bing He, Huan-Ling Liu, and Kai-Jia Sun, “Nuclear liquid-gas phase transition with machine learning,” Phys. Rev. Res. 2, 043202 (2020), arXiv:2010.15043 [nucl-th] .
  * Nara  _et al._ (2000) Y. Nara, N. Otuka, A. Ohnishi, K. Niita, and S. Chiba, “Study of relativistic nuclear collisions at AGS energies from p + Be to Au + Au with hadronic cascade model,” Phys. Rev. C 61, 024901 (2000), arXiv:nucl-th/9904059 .
  * Nara (2019) Yasushi Nara, “JAM: an event generator for high energy nuclear collisions,” EPJ Web Conf. 208, 11004 (2019).
  * Sorge (1995) H. Sorge, “Flavor production in Pb (160-A/GeV) on Pb collisions: Effect of color ropes and hadronic rescattering,” Phys. Rev. C 52, 3291–3314 (1995), arXiv:nucl-th/9509007 .
  * Sorge (1997) H. Sorge, “Soft transverse expansion in pb(158 agev) on pb collisions: preequilibrium motion or first order phase transition?” Physics Letters B 402, 251–256 (1997).
  * Bass  _et al._ (1998) S.A. Bass, M. Belkacem, M. Bleicher, M. Brandstetter, L. Bravina, C. Ernst, L. Gerland, M. Hofmann, S. Hofmann, J. Konopka, G. Mao, L. Neise, S. Soff, C. Spieles, H. Weber, L.A. Winckelmann, H. Stöcker, W. Greiner, Ch. Hartnack, J. Aichelin, and N. Amelin, “Microscopic models for ultrarelativistic heavy ion collisions,” Progress in Particle and Nuclear Physics 41, 255–369 (1998).
  * Bleicher  _et al._ (1999) M Bleicher, E Zabrodin, C Spieles, S A Bass, C Ernst, S Soff, L Bravina, M Belkacem, H Weber, H Stöcker, and W Greiner, “Relativistic hadron-hadron collisions in the ultra-relativistic quantum molecular dynamics model,” Journal of Physics G: Nuclear and Particle Physics 25, 1859–1896 (1999).
  * Kahana  _et al._ (1996) S. H. Kahana, D. E. Kahana, Y. Pang, and T. J. Schlagel, “Modeling relativistic heavy ion collisions at the AGS,” Ann. Rev. Nucl. Part. Sci. 46, 31–70 (1996).
  * Li and Ko (1998) Bao-An Li and C.M. Ko, “Excitation functions of stopping power and flow in relativistic heavy-ion collisions,” Nuclear Physics A 630, 556–562 (1998), nucleus-Nucleus Collisions.
  * Lin  _et al._ (2005) Zi-Wei Lin, Che Ming Ko, Bao-An Li, Bin Zhang, and Subrata Pal, “A Multi-phase transport model for relativistic heavy ion collisions,” Phys. Rev. C 72, 064901 (2005), arXiv:nucl-th/0411110 .
  * Weil  _et al._ (2016) J. Weil _et al._ , “Particle production and equilibrium properties within a new hadron transport approach for heavy-ion collisions,” Phys. Rev. C 94, 054905 (2016), arXiv:1606.06642 [nucl-th] .
  * Antoniou  _et al._ (2006) N. G. Antoniou, F. K. Diakonos, A. S. Kapoyannis, and K. S. Kousouris, “Critical opalescence in baryonic QCD matter,” Phys. Rev. Lett. 97, 032002 (2006), arXiv:hep-ph/0602051 .
  * Bialas and Peschanski (1986) A. Bialas and Robert B. Peschanski, “Moments of Rapidity Distributions as a Measure of Short Range Fluctuations in High-Energy Collisions,” Nucl. Phys. B 273, 703–718 (1986).
  * De Wolf  _et al._ (1996) E. A. De Wolf, I. M. Dremin, and W. Kittel, “Scaling laws for density correlations and fluctuations in multiparticle dynamics,” Phys. Rept. 270, 1–141 (1996), arXiv:hep-ph/9508325 .
  * Antoniou and Diakonos (2019) Nikolaos G. Antoniou and Fotios K. Diakonos, “Ising-QCD phenomenology close to the critical point,” J. Phys. G 46, 035101 (2019), arXiv:1802.05857 [hep-ph] .