# Methodology for Analyzing Proton Multiplicity Fluctuations with Azimuthal
Partitions in Heavy-Ion Collisions

Dylan Neff<EMAIL_ADDRESS>Department of Physics and Astronomy, University
of California, Los Angeles, California 90095, USA  IRFU, CEA, Université
Paris-Saclay, 91191 Gif-sur-Yvette, France  Zhongling Ji  Department of
Physics and Astronomy, University of California, Los Angeles, California
90095, USA  Physics Department, Syracuse University, Syracuse, NY 13244, USA
Roli Esha  Department of Physics and Astronomy, Stony Brook University, Stony
Brook, NY 11790 USA  Gang Wang  Department of Physics and Astronomy,
University of California, Los Angeles, California 90095, USA  Huan Zhong Huang
Department of Physics and Astronomy, University of California, Los Angeles,
California 90095, USA  Key Laboratory of Nuclear Physics and Ion-beam
Application (MOE), and Institute of Modern Physics, Fudan University,
Shanghai-200433, People’s Republic of China

###### Abstract

A primary objective in high-energy heavy-ion collisions is to investigate the
phase transition between confined and deconfined color matter. Complementary
to the cumulants of conserved charges integrated over the full azimuth, we
introduce a novel experimental approach to explore particle fluctuations in
azimuthal partitions, which are potentially sensitive to the first-order phase
transition in heavy-ion collisions. We evaluate proton multiplicity ($N_{w}$)
fluctuations in azimuthal partitions of width $w$ to quantitatively estimate
the clustering tendency among these protons. The $\Delta\sigma^{2}$ observable
is defined as the normalized difference between the variance of the $N_{w}$
distribution and the binomial baseline. We demonstrate the feasibility and
characteristics of this observable through simulations using the AMPT and
MUSIC+FIST models. We also use a Gaussian correlation model to illustrate that
the dependence of $\Delta\sigma^{2}$ on $w$ can be parameterized to accurately
extract the strength and the range of the input interaction among protons.

##  I Introduction

Cumulants, particularly higher-order ones, of net-baryon distributions are
expected to be sensitive to the order of the phase transition between the
deconfined parton system and normal nuclear matter in high-energy heavy-ion
collisions [1, 2, 3]. In experiments, neutrons are typically undetectable, and
the net-proton number is measured as a proxy for the net-baryon number.
Considerable effort has been devoted to the measurements of proton and net-
proton cumulants, as demonstrated, for instance, in Refs. [4, 5]. Yet, the
accurate measurement of cumulants in experiments and their comparison with
theoretical predictions have both posed formidable difficulties. Higher-order
cumulants are not only more sensitive to the pertinent physics, but also more
susceptible to detector anomalies, where a few problematic events can have an
outsized impact [6]. Moreover, the statistical uncertainty in cumulant
measurements rises with order, making higher-order analyses more challenging.
On the other hand, data interpretation must consider other effects, such as
the finite volume and lifetime of the collision system [7], initial volume
fluctuations [8], baryon number transport to midrapidities, and distinctions
between net-proton and net-baryon cumulants [9].

In this paper, we propose an alternative approach to investigate proton
fluctuations and to explore its possible sensitivity to the nature of the
phase transition. During a hypothetical first-order phase transition, the
surface tension between regions of different stable phases is anticipated to
produce clumps of distinct baryon densities in the cooling medium [10]. These
clumps in coordinate space could manifest in momentum space through the radial
expansion of the system and, upon hadronization, may result in correlations
among final-state particles. Accordingly, we can investigate the potential
increase in the clustering of final-state protons as a signal of the first-
order phase transition, though the size of this signal may be significantly
reduced when mapping from coordinate space to momentum space.

Figure 1: Example event with six total protons, including three within a
$120^{\circ}$-partition (shaded area). Red arrows represent the directions of
proton momenta.

To examine proton clustering, we analyze the event-by-event distribution of
proton multiplicity ($N_{w}$) within an azimuthal partition of width $w$ (in
degrees) randomly placed on each event. Figure 1 illustrates a mock event with
$N_{120}=3$, where the total proton multiplicity is $N=6$. Figure 2 provides
an example of the $N_{120}$ distribution for events with $N=6$ in the 0–5%
centrality range of Au+Au collisions at $\sqrt{s_{NN}}=39$ GeV from a
multiphase transport (AMPT) model [11]. If protons within an event are
completely uncorrelated, the probability of a proton falling into a
$w$-partition is $P_{0}=\frac{w}{2\pi}$, and the $N_{w}$ values for events
with $N$ protons follow a binomial distribution with $N$ trials and the
probability of success $P_{0}$. To quantify proton clustering, we define the
$\Delta\sigma^{2}$ observable,

$\Delta\sigma^{2}\equiv\frac{\sigma^{2}(N_{w})-\sigma^{2}_{\text{binomial}}}{N(N-1)},$
(1)

where $\sigma^{2}(N_{w})$ is the variance of the measured $N_{w}$
distribution, and the baseline $\sigma^{2}_{\text{binomial}}=NP_{0}(1-P_{0})$
is the binomial variance expected from uncorrelated protons. The normalization
factor $N(N-1)$ will be explained in Sec. II.1.

Figure 2: Example $N_{120}$ distribution for 500 events with $N=6$ in the 0–5%
centrality range of Au+Au collisions at $\sqrt{s_{NN}}=39$ GeV from AMPT
simulations.

When an attractive interaction forms proton clusters in azimuth, corresponding
voids will also emerge. Therefore, the randomly placed partitions will sample
both excess clusters and excess voids, resulting in a $N_{w}$ distribution
broader than the binomial reference and hence a positive $\Delta\sigma^{2}$.
Conversely, a repulsive interaction tends to distribute protons more
uniformly, leading to fewer clusters or voids, and producing a negative
$\Delta\sigma^{2}$. Figure 3 depicts the $N_{w}$ distributions with
significant attraction (a) and repulsion (b) from a Gaussian correlation
simulation, displaying narrower and wider widths, respectively, than the
binomial distribution (red lines).

The novel $\Delta\sigma^{2}$ observable possesses several strengths that we
aim to leverage. The normalization scheme renders it insensitive to
homogeneous detector inefficiencies, while an event-mixing technique mitigates
contributions from azimuthal-dependent detector inefficiencies and other
detector effects. Furthermore, multiple sampling over the same event (event
resampling) optimally utilizes the statistics of the data. More importantly,
beyond technical aspects, azimuthal partitions unveil a finer structure of the
collision event than an integrated conserved charge.

In Sec. II, we outline the comprehensive methodology for measuring proton
$\Delta\sigma^{2}$. We proceed to study the dependence of this observable on
$N$ and $w$ in simulations using the AMPT and MUSIC+FIST models [12, 13, 14,
15, 16, 17] in Sec. III. In Sec. IV, we develop a Gaussian correlation model
to examine the observable’s sensitivity to the strength and the range of
proton correlations. A summary and outlook are provided in Sec. V. Additional
details and tests on momentum conservation are included in Appendix A.

Figure 3:  Simulations of the $N_{w}$ distributions with significant
attraction (a) and repulsion (b) from a Gaussian correlation model discussed
in Sec. IV. Red lines represent the binomial distributions expected from
uncorrelated protons.

##  II Methodology

In the subsequent analyses, we select the kinematic region for protons
similarly to the RHIC Beam Energy Scan phase I (BES-I) data, with $|y|<0.5$
and $0.4<p_{T}<2$ GeV/$c$. Events are divided into separate classes based on
$N$ and collision centrality, and a $N_{w}$ distribution is constructed for
each class. In this section, we will address the technical details related to
normalization, event mixing, correction for elliptic flow, and event
resampling.

###  II.1  $\Delta\sigma^{2}$ Observable

We begin with $N$ protons, each independently distributed in azimuth, with the
probability density function (PDF) $p(\varphi)$. The integrated probability
for each proton falling in the azimuthal window $[\psi,\psi+w]$ is
$P(\psi)=\int_{\psi}^{\psi+w}p(\varphi)d\varphi$. Here, $\psi$ represents the
orientation of the PDF, which we take to be evenly distributed. The variance
of the $N_{w}$ distribution is

$\begin{split}\sigma^{2}(N_{w})=&E_{\psi}\left[E_{N_{w}|\psi}[N_{w}^{2}]\right]-(E_{\psi}\left[E_{N_{w}|\psi}[N_{w}]\right])^{2}\\\
=&E_{\psi}\left[NP(\psi)\left(1-P(\psi)\right)+\left(NP(\psi)\right)^{2}\right]-(NP_{0})^{2}\\\
=&NP_{0}(1-P_{0})+N(N-1)\left(E_{\psi}[P^{2}(\psi)]-P_{0}^{2}\right)\\\
=&\sigma_{\text{binomial}}^{2}+N(N-1)\left[\int_{0}^{2\pi}\frac{P^{2}(\psi)}{2\pi}\,d\psi-
P_{0}^{2}\right],\end{split}$ (2)

where $\sigma^{2}_{\text{binomial}}\equiv NP_{0}(1-P_{0})$ is the binomial
variance for a uniform $p(\varphi)$ with $P_{0}\equiv\frac{w}{2\pi}$, and
$E_{X|Y}\left[f(X)\right]\equiv\sum_{X}P(X|Y)f(X)$ denotes the conditional
expectation. Rearranging Eq. (2), we can express the observable
$\Delta\sigma^{2}$ from Eq. (1) as the normalized deviation from the binomial
baseline:

$\displaystyle\Delta\sigma^{2}$ $\displaystyle\equiv$
$\displaystyle\frac{\sigma^{2}(N_{w})-\sigma_{\text{binomial}}^{2}}{N(N-1)}$
(3) $\displaystyle=$
$\displaystyle\int_{0}^{2\pi}\frac{P^{2}(\psi)}{2\pi}\,d\psi-P_{0}^{2}.$

While the derivation is based on independent protons that only produce a non-
negative $\Delta\sigma^{2}$ [18], we show in Sec. IV that this observable is
sensitive to proton interactions that could yield both positive and negative
correlations.

###  II.2 Event Mixing

Mixed events are commonly employed in heavy-ion collision analyses, as they
account for detector effects present in all events while eliminating any
intra-event correlation between particles. In particular, mixed events reflect
azimuthal-dependent inefficiencies, whose contribution to $\Delta\sigma^{2}$
must be subtracted from same events,

$\Delta\sigma^{2}_{\text{mix-
corrected}}=\Delta\sigma^{2}_{\text{same}}-\Delta\sigma^{2}_{\text{mixed}}.$
(4)

While not vital to the model studies in this paper, which do not include
detector effects, the event-mixing procedure is crucial in real-data analyses.
In such analyses, a mixed event with $N$ protons can be generated by sampling
each proton from a distinct event. The sampling is performed $N$ times from
events with collision location and collision centrality similar to the current
event under study.

###  II.3 Elliptic Flow Correction

In heavy-ion collisions, the geometric anisotropy of the initial overlap
region can be transformed through hydrodynamic expansion [19] into the
anisotropy of particle emission in momentum space. These collective motions
mimic the clustering behavior in the $\Delta\sigma^{2}$ observable and,
consequently, contribute as background to the physics related to the phase
transition. In a specific kinematic region, the collective motions are usually
quantified by decomposing the $\varphi$ distribution of produced particles in
each collision with a Fourier series [20]:

$\frac{1}{N}\frac{dN}{d\varphi}=\frac{1}{2\pi}\big{(}1+\sum_{n=1}^{\infty}2v_{n}\cos
n\Delta\varphi\big{)},$ (5)

where $\Delta\varphi$ is the azimuthal angle of a particle relative to the
reaction plane (spanned by impact parameter and beam momenta). The
coefficients $v_{n}\equiv\langle\cos n\Delta\varphi\rangle$ are experimentally
obtainable by averaging over particles of interest and over events.
Conventionally, we refer to $v_{1}$ as directed flow, $v_{2}$ as elliptic
flow, and $v_{3}$ as triangular flow.

Applying the $p(\varphi)$ distribution in Eq. (5), we obtain the $v_{n}$
contribution to $\Delta\sigma^{2}$ in Eq. (3),

$\Delta\sigma^{2}_{v_{n}}=\frac{2v_{n}^{2}}{\pi^{2}n^{2}}\sin^{2}\frac{nw}{2}.$
(6)

Directed flow is typically small in magnitude at midrapidities and exhibits an
odd function of rapidity, allowing us to neglect its contribution to a
symmetric rapidity region. The factor of $n^{2}$ in the denominator of Eq. (6)
rapidly diminishes the contribution from higher-order harmonics, irrespective
of their $v_{n}$ magnitudes. Therefore, we will only focus on the correction
for elliptic flow, and subtract the corresponding contribution according to
Eq. (6). Note that the correction depends on $w$ and that the contribution
vanishes for partitions with $w=180^{\circ}$.

###  II.4 Event Resampling

Instead of randomly placing a single azimuthal partition on each event, we can
place multiple partitions on a single event and append these multiplicities to
the $N_{w}$ distribution. This procedure, referred to as “event resampling”,
optimally utilizes the statistics in each event at the cost of generating
distributions composed of entries that are not fully independent. The
resampling is grounded on the assumption that, after correcting for flow and
detector effects, the relevant physics has an equal probability in all
azimuthal directions, regardless of the scale of $w$. In other words, there is
no preferred azimuthal direction for a potential phase transition signal.
Event resampling converges to the true moments of a distribution faster than
taking a single sample (partition), and according to simulations in Ref. [18],
a Poisson block bootstrap procedure offers accurate estimates of the
statistical uncertainties on the moments. We will, therefore, default to
sampling 72 partitions per event in the following analyses.

##  III Application to AMPT and MUSIC+FIST Models

In this section, we highlight certain features of the $\Delta\sigma^{2}$
observable using simulated data from the AMPT and MUSIC+FIST models. As
$\Delta\sigma^{2}$ involves two independent variables, $N$ and $w$, we delve
into the dependence of $\Delta\sigma^{2}$ on each of them.

###  III.1 Model Description

The AMPT model [11] is an event generator that simulates the microscopic
interactions in relativistic heavy-ion collisions. The initial conditions are
managed by the HIJING model [21, 22], with the creation of hadrons, strings,
and minijet partons from initial scattering. We employ the string melting
version of this model [23], in which the initial-state hadrons and minijets
are converted into partons, whose interactions are described by Zhang’s parton
cascade (ZPC) [24]. At freeze-out, hadronization is governed by a quark
coalescence model, and a relativistic transport (ART) model [25] is then
exploited to process the final-state hadronic scattering. The string melting
version of AMPT reasonably well reproduces particle spectra and elliptic flow
in heavy-ion collisions at RHIC and the LHC [26]. A recent update of AMPT with
an improved quark coalescence algorithm [27] is used in this analysis.

MUSIC+FIST is a combination of the MUSIC hydrodynamic model [12, 13, 14, 15,
16] and the FIST sampler [17]. The collision system takes the initial
conditions from the Monte Carlo Glauber model [28] and undergoes the
hydrodynamic evolution in MUSIC. The medium expands until reaching a
“switching” energy density, $\varepsilon_{\text{SW}}=0.26$ GeV/fm3. The FIST
sampler then generates final-state particles based on the produced freeze-out
hypersurface. This model pre-determines the number of each particle species to
sample in each event, ensuring global conservation of conserved quantum
numbers, but at the cost of not conserving momentum. An excluded-volume effect
has been incorporated into the original FIST model [29], ensuring that upon
sampling the hypersurface, no two baryons are sampled within a volume of
radius $r=1$ fm. This radius is an adjustable parameter.

Figure 4: (a) AMPT simulations of $\sigma^{2}(N_{120})$ vs $N$ in the 0–5%
most central Au+Au collisions at $\sqrt{s_{NN}}=$ 39 GeV using 72 partition
samples per event, for both same events (blue circles) and mixed events (green
squares). The red line corresponds to the binomial variance,
$NP_{0}(1-P_{0})$. (b) $\Delta\sigma^{2}$ vs $N$, calculated from
$\sigma^{2}(N_{120})$ in panel (a) by subtracting the binomial variance and
then dividing by $N(N-1)$. Figure 5: $\Delta\sigma^{2}(N_{120})$ vs $N$ in
0–5% Au+Au collisions at different beam energies from AMPT (red circles),
default MUSIC+FIST (blue squares), and MUSIC+FIST with an excluded-volume
effect (purple triangles). The horizontal dashed lines represent the weighted
averages of each data set, with the bands about the lines representing
statistical uncertainties. Figure 6: $\langle\Delta\sigma^{2}\rangle$,
obtained by averaging the $\Delta\sigma^{2}$ values from Fig. 5 over $N$, as a
function of beam energy, from AMPT (red circles), default MUSIC+FIST (blue
squares), and MUSIC+FIST with an excluded-volume effect (purple triangles).
Figure 7: AMPT simulations of $\langle\Delta\sigma^{2}\rangle$ as a function
of the average reference multiplicity for each centrality class in Au+Au
collisions at various beam energies. Results for each of the six energies are
plotted in a different color and connected by straight lines.

###  III.2 $N$ Dependence

Figure 4(a) shows the AMPT simulations of $\sigma^{2}(N_{120})$ vs $N$ in the
0–5% most central Au+Au collisions at $\sqrt{s_{NN}}=$ 39 GeV using 72
partition samples per event, for both same events (blue circles) and mixed
events (green squares). The expected binomial variance (red line)
$NP_{0}(1-P_{0})$ is closely followed by the mixed-event variance, slightly
above the same-event one. Figure 4(b) presents the corresponding
$\Delta\sigma^{2}$ values as a function of $N$.
$\Delta\sigma^{2}_{\text{mixed}}$ is consistent with zero, as AMPT lacks
azimuthal inhomogeneities across events, unlike real data with detector
effects. On the other hand, $\Delta\sigma^{2}_{\rm same}$ is significantly
negative, implying an effectively repulsive interaction among protons in the
AMPT events. In the subsequent analyses, we skip the correction with mixed
events, and only correct $\Delta\sigma^{2}_{\rm same}$ for elliptic flow, as
described in Section II.3, to obtain the final observable $\Delta\sigma^{2}$.

Figure 8: $\langle\Delta\sigma^{2}\rangle$ as a function of $w$ in 0–5% Au+Au
collisions at various beam energies from AMPT (red), default MUSIC+FIST
(blue), and MUSIC+FIST with an excluded-volume effect (purple). The points are
fit to a quadratic function described in Eq. (8).

Figure 5 shows the flow-corrected $\Delta\sigma^{2}(N_{120})$ vs $N$ in 0–5%
Au+Au collisions at various beam energies from AMPT (red circles), default
MUSIC+FIST (blue squares), and MUSIC+FIST with an excluded-volume effect
(purple triangles). At all energies, the AMPT results are significantly
negative, indicating repulsive interactions. The $\Delta\sigma^{2}$ values
from both versions of MUSIC+FIST are much closer to zero, implying a weaker
correlation between protons than in AMPT. In each data set, $\Delta\sigma^{2}$
remains roughly constant with respect to $N$, supporting the normalization
scheme in Eq. (3), which is defined such that $\Delta\sigma^{2}$ should be
independent of $N$. This constancy justifies averaging $\Delta\sigma^{2}$ over
$N$, as represented by the horizontal lines.

Figure 6 shows $\langle\Delta\sigma^{2}\rangle$, obtained by averaging the
$\Delta\sigma^{2}$ values from Fig. 5 over $N$, as a function of beam energy.
In the default MUSIC+FIST model, $\langle\Delta\sigma^{2}\rangle$ is minimally
but significantly negative, with no evident energy dependence. The excluded-
volume version of MUSIC+FIST exhibits markedly more negative outcomes compared
with the default version. Since the excluded-volume effect is expected to
produce an effectively repulsive interaction, the more negative
$\langle\Delta\sigma^{2}\rangle$ values corroborate that our observable is
indeed sensitive to repulsion. The AMPT results seem to reveal much stronger
repulsive interactions than MUSIC+FIST, and the repulsion strength intensifies
at lower collision energies.

To explore the disparity between MUSIC+FIST and AMPT, as well as the beam-
energy dependence of the latter, we depict in Fig. 7
$\langle\Delta\sigma^{2}\rangle$ as a function of the average reference
multiplicity in each centrality class for the investigated beam energies.
Reference multiplicity is defined by the number of charged hadrons within the
pseudorapidity range of $|\eta|<0.5$ in an event. The AMPT points at all
energies appear to follow a universal curve, with the potential exception of
the 7.7 GeV data being slightly less negative than other energies. This
decreasing trend toward lower multiplicities explains the dependence on
center-of-mass energy in Fig. 6, because for the fixed 0–5% centrality,
collisions at lower energies produce fewer final-state particles. We attribute
this pronounced multiplicity dependence primarily to global momentum
conservation. In a peripheral collision, where only a few particles are
produced, each particle bears greater responsibility for conserving momentum.
This tendency results in a more uniform distribution of particles in azimuth
and, consequently, a negative $\Delta\sigma^{2}$. Each additional particle
introduces an extra degree of freedom, allowing for a more flexible means to
satisfy global momentum conservation, and the negative correlation is diluted.
This effect approaches zero as more particles are added to the event. As
mentioned in Sec. III.1, the FIST sampler does not conserve momentum.
Therefore, the MUSIC+FIST results do not manifest the same trend as AMPT in
Fig. 6.

This hypothesis of momentum conservation is substantiated by simulations
presented in Appendix A. In these simulations, particles are generated with
random momenta which are slightly rotated until the total momentum of the
event approaches zero. The resulting $\langle\Delta\sigma^{2}\rangle$ as a
function of $M$, the total number of particles in the event, mirrors the trend
observed in AMPT.

###  III.3 $w$ Dependence

The preceding discussion has been confined to a partition width of
$w=120^{\circ}$. Intuitively, varying the partition width may offer
sensitivity to the dynamical scale associated with the measured correlation.
We will investigate this sensitivity in Sec. IV with a Gaussian correlation
model. In this subsection, we first illustrate the dependence of
$\langle\Delta\sigma^{2}\rangle$ on $w$ and then discuss its parametrization.

Figure 8 shows $\langle\Delta\sigma^{2}\rangle$ as a function of $w$ in 0–5%
Au+Au collisions at various beam energies from AMPT (red circles), default
MUSIC+FIST (blue squares), and MUSIC+FIST with an excluded-volume effect
(purple triangles). The apparent symmetry about $180^{\circ}$ arises from the
mathematical nature of Eq. (3). Substituting $w$ with $(2\pi-w)$ is equivalent
to replacing $P_{0}$ with $(1-P_{0})$ and, in the integral, replacing
$P(\psi)$ with $[1-P(\psi)]$. Under this operation, Eq. (3) remains unchanged.
The proof entails the following identity,

$\displaystyle\int_{0}^{2\pi}P(\psi)d\psi$ $\displaystyle=$
$\displaystyle\int_{0}^{2\pi}d\psi\int_{\psi}^{\psi+w}p(\varphi)d\varphi$ (7)
$\displaystyle=$
$\displaystyle\int_{0}^{w}d\varphi\int_{0}^{2\pi}p(\varphi+\psi)d\psi$
$\displaystyle=$ $\displaystyle\int_{0}^{w}d\varphi=w.$

Figure 9: The general scale $b$ extracted from the $w$ dependence of
$\langle\Delta\sigma^{2}\rangle$ in Fig. 8 using the quadratic fit in Eq. (8)
as a function of reference multiplicity from AMPT. Figure 10: The inverse
curvature $z$ extracted from the $w$ dependence of
$\langle\Delta\sigma^{2}\rangle$ in Fig. 8 using the quadratic fit in Eq. (8)
as a function of reference multiplicity from AMPT.

We parameterize the $w$ dependence of $\langle\Delta\sigma^{2}\rangle$ with a
quadratic function,

$\langle\Delta\sigma^{2}\rangle=b\left[1-\left(\frac{w-180}{z}\right)^{2}\right],$
(8)

where $b$ and $z$ quantify the general scale and the inverse curvature,
respectively, supposedly reflecting the strength and the range of the
correlation between protons. We extract these two parameters from the AMPT
simulations in Fig. 8, and present $b$ and $z$ as a function of reference
multiplicity in Figs. 9 and 10, respectively. The universal trend in Fig. 9
looks akin to that in Fig. 7, with a larger scale, since $b$ corresponds to
the maximum $\langle\Delta\sigma^{2}\rangle$ magnitude occurring at
$w=180^{\circ}$. Representing the general scale with $b$ removes the
arbitrariness of the choice in $w$. In Fig. 10, there is no apparent
dependence of $z$ on either beam energy or reference multiplicity. We
anticipate applying this methodology to real data and comparing the outcomes
with model simulations to deepen our understanding of the pertinent physics.

##  IV Analysis with Gaussian Correlation Model

In this section, we examine the sensitivity of the $\Delta\sigma^{2}$
observable to proton interactions. We devise a simple toy model incorporating
both positive and negative Gaussian-type correlations among protons. We will
demonstrate that our proposed method can effectively extract the input
correlation.

###  IV.1 Model Description

In this model, we sample protons one by one, each ($\varphi_{i}$) from a
specific azimuthal probability distribution determined by the positions of all
previously generated protons. The first proton follows a uniform distribution
in the interval $[0,2\pi)$. The PDF for the second proton is a combination of
a uniform distribution and a Gaussian centered at $\varphi_{1}$,
$P(\varphi_{2})\propto
1+Ae^{-\frac{1}{2}\left(\frac{\varphi_{2}-\varphi_{1}}{\sigma}\right)^{2}}$,
with appropriate normalization. The third proton has
$P(\varphi_{3})\propto\left[1+Ae^{-\frac{1}{2}\left(\frac{\varphi_{3}-\varphi_{1}}{\sigma}\right)^{2}}\right]\left[1+Ae^{-\frac{1}{2}\left(\frac{\varphi_{3}-\varphi_{2}}{\sigma}\right)^{2}}\right]$.
This process continues, and the PDF for the $n$-th proton receives a
contribution from all previous protons:

$P(\varphi_{n>2})\propto\prod_{i=1}^{n-1}{\left[1+Ae^{-\frac{1}{2}\left(\frac{\varphi_{n}-\varphi_{i}}{\sigma}\right)^{2}}\right]}.$
(9)

Since Gaussian distributions extend from $-\infty$ to $\infty$, $\varphi_{i}$
could fall outside the $[0,2\pi)$ range. In that case, we fold it back into
$[0,2\pi)$.

The two parameters in this Gaussian correlation model, $A$ and $\sigma$,
represent the strength and the range, respectively, of the correlation among
protons. Parameter $A$ can take either positive or negative values,
corresponding to positive (clustering) or negative (repulsion) correlations.
The range of $A$, $[-1,+\infty)$, is asymmetric about zero. $A$ being $-1$
prevents any proton from being placed exactly on top of a pre-existing proton.
Conversely, $A$ must approach $+\infty$ to ensure each proton is placed
directly atop the first one. This asymmetry implies that simulation results
with positive and negative values of $A$ should not be compared on the same
footing, except for small magnitudes of $A$.

###  IV.2 $\Delta\sigma^{2}$ Analysis

Utilizing the Gaussian correlation model, we create proton events with $N$
ranging from 2 to 80, encompassing the maximum multiplicity observed in the
AMPT data. Subsequently, we analyze these events using the same approach as
conducted in Sec. III. We systematically explore different combinations of $A$
and $\sigma$ parameters to understand the observable’s behavior across the
model’s parameter space. In this simulation, there is no need for correction
regarding either mixed events or elliptic flow. We use 72 partition samples
per event.

Figure 11 shows the simulations of $\Delta\sigma^{2}(N_{120})$ as a function
of $N$ from the Gaussian correlation model, with the input parameters $A=\pm
0.002$, $\pm 0.006$, and $\pm 0.01$, and $\sigma=1$. All attractive sets with
positive $A$ produce positive $\Delta\sigma^{2}$ values, whereas all repulsive
sets with negative $A$ render negative $\Delta\sigma^{2}$ values. In addition,
larger $|A|$ values lead to larger magnitudes of $\Delta\sigma^{2}$. Thus,
$\Delta\sigma^{2}(N_{120})$ appears to be correlated with the strength of the
proton correlation. Similar to Fig. 5, $\Delta\sigma^{2}$ remains constant
against $N$. We will use $\langle\Delta\sigma^{2}\rangle$, averaged over $N$,
in the following analyses.

Figure 11: Simulations of $\Delta\sigma^{2}(N_{120})$ as a function of $N$
from the Gaussian correlation model. The input parameters are $A=\pm 0.002$,
$\pm 0.006$, and $\pm 0.01$, and $\sigma=1$. Figure 12:
$\langle\Delta\sigma^{2}\rangle$, obtained by averaging over $N$, vs the width
of the azimuthal partition, $w$, from the Gaussian correlation model. The
input parameters are $A=\pm 0.002$ and $\pm 0.008$, and $\sigma=0.4$ and 0.8.

Figure 12 shows the simulations of $\langle\Delta\sigma^{2}\rangle$ as a
function of $w$ from the Gaussian correlation model with input parameters
$A=\pm 0.002$ and $\pm 0.008$, and $\sigma=0.4$ and 0.8. With the same $w$,
the magnitude of $\langle\Delta\sigma^{2}\rangle$ increases with $|A|$,
regardless of whether the proton interaction is attractive or repulsive. These
trends can be reasonably well described by the quadratic function in Eq. (8),
where the sign of the obtained $b$ parameter matches the sign of the input
$A$, and a larger $|A|$ corresponds to a larger $|b|$. The fit curves with the
same $\sigma$ input seem to intersect the $x$-axis around the same location in
Fig. 12, and therefore, they have similar $z$ values.

Figure 13: The inverse curvature $z$ extracted from the $w$ dependence of
$\langle\Delta\sigma^{2}\rangle$ using the quadratic fit in Eq. (8) as a
function of input parameter $\sigma$ from the Gaussian correlation model. The
input parameter $A$ is chosen to be $\pm 0.002$, $\pm 0.006$, and $\pm 0.01$.
The dashed curve represents a simple function, $z=12/\sigma+131$, describing
the relation between $z$ and $\sigma$. Some markers have been slightly shifted
along the $x$-axis to improve clarity.

Figure 13 shows the inverse curvature $z$ extracted from the $w$ dependence of
$\langle\Delta\sigma^{2}\rangle$ using the quadratic fit in Eq. (8) as a
function of input parameter $\sigma$ from the Gaussian correlation model, with
input parameter $A=\pm 0.002$, $\pm 0.006$, and $\pm 0.01$. Note that the
units of $z$ and $\sigma$ are degree and radian, respectively. The extracted
$z$ seems to be solely dependent on the $\sigma$ value, independent of input
parameter $A$. As illustrated by the dashed curve, $z$ follows a simple
function of $\sigma$, $z=12/\sigma+131$. This supports the idea that $z$ holds
a meaningful physical interpretation, reflecting the correlation range.
According to the realistic AMPT simulations in Fig. 10, the $z$ values
typically fall at or below $145^{\circ}$, corresponding to $\sigma\geq 0.8$ in
Fig. 13.

Figure 14 presents the general scale $b$ extracted from the $w$ dependence of
$\langle\Delta\sigma^{2}\rangle$ using the quadratic fit in Eq. (8) as a
function of input parameter $|A|$ from the Gaussian correlation model, with
$\sigma$ ranging from 0.4 to 1.2. The sign of $b$ corresponds to positive or
negative $A$, indicating attractive or repulsive proton interactions. In cases
where $\sigma\geq 0.8$, $b$ follows a linear function of $A$ passing through
$(0,0)$, regardless of the magnitude of $\sigma$. Hence, in realistic
scenarios, $b$ serves as an indicator of the correlation strength in this
model. When $\sigma\leq 0.8$, the $b$ values deviate from the universal trend,
but the linear relation between $b$ and $A$ still holds. If real data yield
$z$ values corresponding to small $\sigma$, the measured $b$ can still be
employed to determine the correlation strength, and its sign indicates whether
the interaction is attractive or repulsive.

Figure 14: The general scale $b$ extracted from the $w$ dependence of
$\langle\Delta\sigma^{2}\rangle$ using the quadratic fit in Eq. (8) as a
function of input parameter $|A|$ from the Gaussian correlation model. The
positive (negative) $b$ values correspond to attractive (repulsive)
interactions with positive (negative) $A$. The input parameter $\sigma$ varies
from 0.4 to 1.2. The dashed lines represent linear functions passing through
the point of origin.

One concern regarding the Gaussian correlation model is the sequential
placement of protons. While this procedure mimics the interaction between
protons, the protons may not be treated equally. As each proton is generated
according to the positions of all previous ones, protons placed earlier appear
to be influenced by fewer other protons. However, this effect is expected to
be negligible for two reasons. Firstly, the two angles in the Gaussian
distribution in Eq. (9) are mathematically symmetric, indicating that the
implemented interaction between two protons is mutual. Secondly, the
$\Delta\sigma^{2}$ observable randomly samples and resamples all protons
without any specific sequence or bias.

##  V Conclusion

To explore the nature of the phase transition in heavy-ion collisions, we
advocate a new approach to investigate proton multiplicity fluctuations in
azimuthal partitions. We introduce a $\Delta\sigma^{2}$ observable to quantify
the azimuthal clustering of final-state protons. The robustness of the
measurement against detector inefficiency is ensured through the normalization
scheme and the mixed-event correction. Additionally, the correction for
elliptic flow eliminates the contribution of collective motion to the
clustering signal. Leveraging the event resampling technique allows for
optimal utilization of available statistics.

We have applied the proposed methodology to simulations using the AMPT and the
MUSIC+FIST models for Au+Au collisions at center-of-mass energy ranging from
7.7 to 62.4 GeV. The resulting $\Delta\sigma^{2}$ values remain independent of
the total proton multiplicity ($N$), validating the effectiveness of the
normalization scheme. The discrepancy in $\Delta\sigma^{2}$ between the two
versions of MUSIC+FIST highlights the observable’s sensitivity to the
excluded-volume effect. The variation in AMPT results with collision energy
and event multiplicity demonstrate general features due to global momentum
conservation.

With both models, we have also investigated the dependence of
$\langle\Delta\sigma^{2}\rangle$, averaged over $N$, on the azimuthal
partition width $w$. We characterize the $w$ dependence with a quadratic
function, utilizing parameters $b$ and $z$ to quantify the general scale and
the inverse curvature, respectively. With a Gaussian correlation model, we
illustrate that the extracted $b$ and $z$ values reflect the strength and the
range, respectively, of the input interaction between protons.

The application of the proposed approach to the RHIC BES data will provide new
insight into the nature of the quark-hadron phase transition in the QCD phase
diagram. In particular, measurements of the innovative observable
$\Delta\sigma^{2}$ will allow us to probe the clustering phenomenon among
protons. While addressing the overall scale of $\Delta\sigma^{2}$ necessitates
further evaluation of the momentum conservation effect, the inverse curvature
shows little dependence on the scale magnitude and reflects mostly the
correlation width in the Gaussian correlation model. The correlation width in
azimuth could be linked to the underlying physics of the phase transition,
making its beam-energy dependence of great experimental interest.

###### Acknowledgements.

The authors thank Zhiwan Xu, Yunshan Cheng, Thomas Marshall, Aditya Dash, and
Xiatong Wu for many fruitful discussions. Z. J., G. W., and H. H. are
supported by the U.S. Department of Energy under Grant No. DE-FG02-88ER40424
and by the National Natural Science Foundation of China under Contract
No.1835002. We acknoledge the support from the Office of Nuclear Physics in
the Office of Science of the Department of Energy and the Commissariat ‘a
l’Energie Atomique (France).

##  Appendix A Impact of Momentum Conservation

The conservation of momentum in each event is expected to induce a negative
correlation between particles. In an extreme case where only two particles are
produced back to back, their momenta are perfectly anti-correlated. If more
particles emerge in the final state, each becomes less crucial for conserving
the total momentum. Consequently, as the number of particles approaches
infinity, the anti-correlation between them vanishes. While we comprehend the
behavior in extreme cases, the precise functional form of this effect remains
elusive.

We have employed a simple model to examine the impact of momentum conservation
on the $\Delta\sigma^{2}$ observable. In each event, we generate $M$
particles, with each momentum component of each particle randomly sampled from
a uniform distribution within $(0,2)$ GeV/$c$. We calculate the total momentum
of the event as the sum of individual momenta. Each particle is then slightly
rotated away from this total momentum, the rotation amount being proportional
to the particle’s momentum magnitude and the event’s total momentum magnitude.
After all particles undergo these slight rotations, we recalculate the total
momentum, and the process is repeated. As the rotations diminish with each
iteration, the total momentum converges to zero. We then randomly label 40% of
these particles as protons, and further apply the kinematic selection with
$|y|<0.5$ and $p<2$ GeV/$c$.

Figure 15: Proton $\langle\Delta\sigma^{2}\rangle$, averaged over $N$, as a
function of total particle number, $M$, from a simple model that conserves
total momentum. The fit function suggests that a simple form of $\propto 1/M$
can well describe the trend.

Figure 15 displays proton $\langle\Delta\sigma^{2}\rangle$, averaged over $N$,
as a function of $M$ from the model simulation. As anticipated, the effect of
momentum conservation is evident in notably negative
$\langle\Delta\sigma^{2}\rangle$ values, which become more negative for events
with lower multiplicities. We fit this trend with a power law function plus a
constant, and the best fit suggests a $1/M$ dependence without an offset. This
trend mirrors observations in the AMPT simulations. The straightforward
functional relationship indicates that the impact of momentum conservation on
the $\langle\Delta\sigma^{2}\rangle$ observable may be analytically
manageable, implying the potential for a correction procedure.

While the simulation results are enticing, we note that the momentum
conservation model presented here is oversimplified. Several factors, such as
particle spectra, the proton fraction, and the particle momentum distribution
may impact the $\langle\Delta\sigma^{2}\rangle$ values. Additionally, the
model lacks many physics aspects involved in heavy-ion collisions, and it is
beyond the scope of this study to construct a realistic model for such
collisions.

## References

  * Stephanov  _et al._ [1999] M. Stephanov, K. Rajagopal, and E. Shuryak, Phys. Rev. D 60, 114028 (1999).
  * Stephanov [2011] M. A. Stephanov, Phys. Rev. Lett. 107, 052301 (2011).
  * Asakawa  _et al._ [2009] M. Asakawa, S. Ejiri, and M. Kitazawa, Phys. Rev. Lett. 103, 262301 (2009).
  * Adam  _et al._ [2021] J. Adam, L. Adamczyk, J. R. Adams, _et al._ (STAR Collaboration), Phys. Rev. Lett. 126, 092301 (2021).
  * Aboona  _et al._ [2023] B. E. Aboona, J. Adam, L. Adamczyk, _et al._ (STAR Collaboration), Phys. Rev. Lett. 130, 082301 (2023).
  * Chen  _et al._ [2024] L. Chen, Y.-Y. Zhao, Y. Cheng, _et al._ , Phys. Rev. C 109, 034911 (2024).
  * Tang  _et al._ [2023] S. Tang, S. Wu, and H. Song, Phys. Rev. C 108, 034901 (2023).
  * Gorenstein and Gaździcki [2011] M. I. Gorenstein and M. Gaździcki, Phys. Rev. C 84, 014904 (2011).
  * Kitazawa and Asakawa [2012] M. Kitazawa and M. Asakawa, Phys. Rev. C 85, 021901 (2012).
  * Steinheimer and Randrup [2013] J. Steinheimer and J. Randrup, Phys. Rev. C 87, 054903 (2013).
  * Lin  _et al._ [2005] Z.-W. Lin, C. M. Ko, B.-A. Li, _et al._ , Phys. Rev. C 72, 064901 (2005).
  * Schenke  _et al._ [2010] B. Schenke, S. Jeon, and C. Gale, Phys. Rev. C 82, 014903 (2010).
  * Schenke  _et al._ [2011] B. Schenke, S. Jeon, and C. Gale, Phys. Rev. Lett. 106, 042301 (2011).
  * Paquet  _et al._ [2016] J.-F. m. c. Paquet, C. Shen, G. S. Denicol, _et al._ , Phys. Rev. C 93, 044906 (2016).
  * Denicol  _et al._ [2018] G. S. Denicol, C. Gale, S. Jeon, _et al._ , Phys. Rev. C 98, 034916 (2018).
  * Huovinen and Petersen [2012] P. Huovinen and H. Petersen, The European Physical Journal A 48, 171 (2012).
  * Vovchenko  _et al._ [2022] V. Vovchenko, V. Koch, and C. Shen, Phys. Rev. C 105, 014904 (2022).
  * Neff [2023] D. Neff,  _Exploring the QCD Phase Diagram at STAR: Measurement of Proton Fluctuations in Azimuthal Partitions_ , Ph.D. thesis, UCLA (2023).
  * Heinz and Snellings [2013] U. Heinz and R. Snellings, Annual Review of Nuclear and Particle Science 63, 123 (2013).
  * Poskanzer and Voloshin [1998] A. M. Poskanzer and S. A. Voloshin, Phys. Rev. C 58, 1671 (1998).
  * Wang and Gyulassy [1991] X.-N. Wang and M. Gyulassy, Phys. Rev. D 44, 3501 (1991).
  * Gyulassy and Wang [1994] M. Gyulassy and X.-N. Wang, Computer Physics Communications 83, 307 (1994).
  * Lin and Ko [2002] Z.-w. Lin and C. M. Ko, Phys. Rev. C 65, 034904 (2002).
  * Zhang [1998] B. Zhang, Computer Physics Communications 109, 193 (1998).
  * Li and Ko [1995] B.-A. Li and C. M. Ko, Phys. Rev. C 52, 2037 (1995).
  * Lin [2014] Z.-W. Lin, Phys. Rev. C 90, 014904 (2014).
  * He and Lin [2017] Y. He and Z.-W. Lin, Phys. Rev. C 96, 014910 (2017).
  * Miller  _et al._ [2007] M. L. Miller, K. Reygers, S. J. Sanders, and P. Steinberg, Annual Review of Nuclear and Particle Science 57, 205 (2007).
  * Vovchenko [2022] V. Vovchenko, Phys. Rev. C 106, 064906 (2022).


# Neural Network Quantisation
for Faster Homomorphic Encryption ††thanks: This project has received funding
from the European Research Council (ERC) under the European Union’s Horizon
2020 research and innovation programme (grant agreement nr. 101020005) and
from the Defense Advanced Research Projects Agency (DARPA) under Contract No.
HR0011-21-C-0034. Michiel Van Beirendonck is funded by FWO as Strategic Basic
(SB) PhD fellow (project number 1SD5621N). Jan-Pieter D’Anvers is funded by
FWO (Research Foundation – Flanders) as junior post-doctoral fellow (contract
number 133185).

Wouter Legiest, Furkan Turan, Michiel Van Beirendonck, Jan-Pieter D’Anvers and
Ingrid Verbauwhede  COSIC, KU Leuven, Belgium
<EMAIL_ADDRESS>

###### Abstract

Homomorphic encryption (HE) enables calculating on encrypted data, which makes
it possible to perform privacy-preserving neural network inference. One
disadvantage of this technique is that it is several orders of magnitudes
slower than calculation on unencrypted data. Neural networks are commonly
trained using floating-point, while most homomorphic encryption libraries
calculate on integers, thus requiring a quantisation of the neural network. A
straightforward approach would be to quantise to large integer sizes (e.g.
$32\text{\,}\mathrm{bit}$) to avoid large quantisation errors. In this work,
we reduce the integer sizes of the networks, using quantisation-aware
training, to allow more efficient computations. For the targeted MNIST
architecture proposed by Badawi et al. [1], we reduce the integer sizes by 33%
without significant loss of accuracy, while for the CIFAR architecture, we can
reduce the integer sizes by 43%. Implementing the resulting networks under the
BFV homomorphic encryption scheme using SEAL, we could reduce the execution
time of an MNIST neural network by 80% and by 40% for a CIFAR neural network.

###### Index Terms:

convolutional neural networks, quantisation, privacy-preserving machine
learning, fully homomorphic encryption

##  I Introduction

Homomorphic encryption (HE) allows performing calculations on encrypted data.
This technique enables applications where data is processed in untrusted
environments (e.g. a cloud environment) while ensuring that this environment
does not learn anything about the data itself. As such, it is a promising
technique to make privacy-preserving machine learning possible.

A downside of HE is that it significantly increases the size of encrypted
data. As a result, encrypted operations are typically several orders of
magnitude slower than their unencrypted counterparts. This work tries to
accelerate neural network inference under homomorphic encryption by using
quantisation techniques to reduce the data size and, thus, the computational
cost.

Neural network frameworks generally use a floating-point representation to
represent network parameters and intermediate variables. However, HE systems
like BFV [2] encode only integers, requiring an additional conversion step to
convert the floating-point neural network parameters to the integer HE
variables. While it is possible to design neural networks that work solely
with integer representations, previous works have only studied such networks
in a non-HE related context [3, 4, 5].

In addition, this conversion is an essential step before porting it to
hardware. For instance, a plaintext $32\text{\,}\mathrm{bit}$ floating-point
addition is $30\times$ more energy-consuming 111Energy consumption using a
$45\text{\,}\mathrm{nm}$ CMOS technology. than an $8\text{\,}\mathrm{bit}$
integer equivalent [6]. By using the conversion, we can select smaller HE
parameters that lead to limited resource use and better management of corner
cases. Therefore, making the behaviour of the system faster and more
predictable in general.

As calculations in these non-HE integer-only networks are performed, the sizes
of the integer variables increase. The intermediate values are commonly scaled
down to a smaller number after each layer to keep these integer-only networks
manageable. This means the most significant bits are held after each
operation, while the least significant are discarded.

Unfortunately, these reduction operations are based on division or shift
operations, which natively are not supported in HE schemes, so downscaling
cannot easily be performed. Therefore, in neural network HE inferences, the
intermediate values will grow throughout the inference and the final
calculations will need to operate on very large integers. For instance, when
all of the weights of a neural network are converted to
$32\text{\,}\mathrm{bit}$, a 10-layer CIFAR network will produce integers with
bit-sizes up to $614\text{\,}\mathrm{bit}$. The maximum bit-length of these
output integers will be denoted as the _‘final integer width’_ (FIW), and we
will show that this value significantly affects the overall computation cost.

Gilad-Bachrach et al. [7] implemented the first artificial feedforward neural
network under homomorphic encryption using the HE scheme YASHE [8]. Note that
an attack proposed by Albrecht et al. [9] reduced the security level of this
scheme and is therefore considered broken in practice. Gilad-Bachrach et al.
[7] proposed a specialised, HE-focussed _CryptoNets_ architecture for the
MNIST dataset [10].

One of the downsides of the CPU implementations of CryptoNets is the high
latency of $250\text{\,}\sec$ for an MNIST image. It was improved by Brutzkus
et al. [11] with the Low-Latency CryptoNets (LoLa) architecture. Using the BFV
scheme, optimisations in the underlying HE library SEAL and a different
approach to representing the ciphertext data, a latency of $0.29\text{\,}\sec$
was reached, an improvement of $93\times$ relative to CryptoNets. In addition,
Brutzkus et al. [11] proposes variants of the LoLa network for processing the
CIFAR-10 dataset [12]. They report an accuracy of 74.1% and a latency of
$730\text{\,}\sec$.

Badawi et al. [1] implemented the BFV scheme on GPUs. They propose two
architectures, one smaller for MNIST and one more extensive for CIFAR.
Accordingly, their CIFAR network boasts an accuracy of 77.55% and a latency of
$304.43\text{\,}\sec$.

In this work, we improve upon the state-of-the-art HE neural networks by
considering advanced neural network quantisation techniques. We first
investigate post-training quantisation, a method typically used in the state-
of-the-art, and show that there is a limit to how many intermediate variables
can be scaled down without significantly affecting accuracy. We then show that
quantisation-aware training can indeed be used to substantially scale down
these intermediate variables without a similar accuracy penalty. In the end,
we reduced the final integer width with 33% for MNIST and 43% for CIFAR,
allowing a speedup with factors 80% and 40%, respectively, over typical 8-bit
post-training quantisation networks as used in the state-of-the-art.

##  II Preliminaries

###  II-A Homomorphic encryption

Homomorphic encryption enables performing arithmetic operations on encrypted
data. Take the following example: consider an asymmetric encryption system
with two integer numbers $x$ and $y$. They can be encrypted by the using
encryption key $\mathsf{pk}$ to $c_{x}=\mathsf{Enc}(\mathsf{pk},x)$ and
$c_{y}=\mathsf{Enc}(\mathsf{pk},y)$. These two ciphertexts are sent to a
server that cannot be trusted. The server can perform an operation $\diamond$
on both ciphertext $c_{xy}=c_{x}\diamond c_{y}$, which is the equivalent of
doing an addition on the plaintexts. The result of this operation is then sent
back to the user, who can decrypt this message. Using the decryption key
$\mathsf{sk}$, the resulting plaintext message $z$ is obtained by
$z=\mathsf{Dec}(\mathsf{sk},c_{xy})$. The message $z$ will have a value of
$z=x+y$. Altogether, the server does not obtain any information about the
integers $x$ and $y$ while possessing the unencrypted data.

A limitation of this form of encryption is that it only allows certain
operations, i.e. addition or multiplication of two ciphertexts. Execution of
non-linear functions is normally performed using a polynomial approximation
that uses only addition, subtraction, and multiplication. Moreover, a division
in the HE schemes CKKS and BFV is theoretically possible, but it is costly and
thus avoided in practice [13].

The biggest problem with the lack of a division operation is that variables
will grow during computation. For example, when multiplying two
$8\text{\,}\mathrm{bit}$ integers, the result becomes roughly
$16\text{\,}\mathrm{bit}$. In unencrypted neural network implementations, this
variable can be divided with a power of two to get back to
$8\text{\,}\mathrm{bit}$, making it more manageable for the next layer.
However, such an operation is not possible in encrypted neural network
inference. This leads to large intermediate and output integers. The maximum
bit-length of these output integers will be denoted as the ‘final integer
width’ (FIW).

The HE scheme must be instantiated with more extensive parameters to
accommodate these larger variables, which comes at a significant cost. Once a
specific variable size is reached, additional techniques are required to allow
large representations. More specifically, to ensure a correct representation
in the plaintext space during inference, a residue numeral system (RNS),
following the Chinese remainder theorem, is used to divide the large numbers
into several smaller numbers. This leads to several smaller HE instances that
could be run in parallel. Since each instance consumes computing resources,
decreasing the variable sizes can significantly reduce the number of RNS
instances and, thus, the computational cost.

###  II-B Neural network

A neural network is a machine learning technique consisting of a network of
small interconnected computation units called neurons. These neurons can be
adapted, which enables the network to ‘learn’ a specific, human-like task such
as classifications of images. A neuron will take a number of inputs, perform a
weighted sum over these inputs, and output a function of the result of this
sum. Neurons are grouped to form layers, and different behaviour can be
obtained depending on their configuration.

Since a typical division or non-linear function cannot be executed trivially
under FHE, we use a slightly adapted version of the classical neural network
layers. A dense or convolutional layer is representable under FHE. However,
the activation function is approximated by a $f(x)=x^{2}$ square function,
resulting in a _Square layer_. Moreover, the _scaled average pooling_ layer is
replaced by an equivalent where the inputs are summed, but the division is
omitted.

###  II-C Architectures

This work uses the two architectures developed by Badawi et al. [1] for
homomorphic inference. These networks are used as test cases to research the
effect of quantisation on homomorphic encryption inference. Both architectures
omit the last (Sigmoid) activation function since it only maps the output to
the unit interval. For a detailed description of the network, we refer the
reader to the paper of Badawi et al. [1].

The first architecture used in this paper focuses on the MNIST dataset [10].
It is based on the HCNN [1] and consists of two convolutional, two square
activation and one dense layer. The authors stated an accuracy of 99% for this
architecture.

The second architecture is designed to classify the more complex CIFAR-10
dataset [12]. The 10-layered network especially uses the scaled average
pooling and square layer. The originally proposed HCNN architecture was
slightly modified in our implementation by not using padding, as this only
results in reduced accuracy: our floating-point model model obtains an
accuracy of 73.28%, while the original HCNN reports 77.8%.

##  III Post-Training Quantisation (PTQ)

Usually, floating-point numbers with single or double precision are used to
represent the weights and biases of a network. However, it is possible to
convert these numbers to $8\text{\,}\mathrm{bit}$ integers without a notable
reduction in accuracy [14]. A further reduction in the representation might
have a more detrimental effect on the neural network accuracy. Converting an
existing (floating-point) neural network into a quantised (integer) version is
called post-training quantisation (PTQ).

During PTQ quantisation, a real value $r\in[\alpha,\beta]$ is converted to a
$b$-bit integer $q$. The process is determined by two factors: the zero-point
$Z$ and the scale factor $S$, using the following formula:

$\displaystyle q$ $\displaystyle=\left\lfloor\frac{r}{S}+Z\right\rceil.$ (1)

Dequantisation can be done through the formula $r=S(q-Z)$, where the quantised
value is converted back to its original scale.

The scale factor $S$ determines the quantisation step size. The zero-point $Z$
is the quantised value $q$ corresponding to the real value $r=0$ and positions
the range of representable numbers optimally.

When $Z\neq 0$, we say the quantisation is asymmetric or affine. This
quantisation explicitly uses the zero point, often set at
$Z=-\alpha\cdot(2^{b}-1)/(\beta-\alpha)$. A second option is a symmetric
quantisation, which reduces the overhead of dealing with the zero point by
setting it to zero. Commonly, the values are mapped to a signed symmetric
interval $[\alpha,\beta]=[-2^{b-1},2^{b-1}-1]$, although an unsigned interval
is also possible. Symmetric quantisation is a more limited but easier-to-
handle quantisation technique.

To evaluate the effect of quantisation, we first determined the distribution
of the parameters of both MNIST and CIFAR networks, plotted in Figure 1. Since
both networks possess a symmetric distribution, the mean of the values is
zero, and thus symmetric quantisation is the best candidate to convert the
signed real numbers for both networks. To determine the ideal scale factor for
the weights, three candidates are tested:

The first scale factor $S=1/(2^{b-1}-1)$ only considers the bit width. No
account is taken of the size or distribution of the real numbers. The second
scaling factor $S=\max(|\boldsymbol{W}|)/(2^{b-1}-1)$ considers the largest
absolute value that can be represented in the quantised interval, and extreme
values that are outside the quantisation range are quantised to the edges of
the quantisation interval.

To understand the influence of these different scale factors on the network,
we build a Python framework that evaluates the effect of post-training
quantisation on the accuracy and FIW. The framework takes a neural network,
converts each of the weights to an integer representation and then executes a
neural network inference. We process each of the $10\;000$ images in the test
set for these experiments to determine the accuracy and FIW. The maximum of
all individual final integer widths and corresponding accuracies are reported
in Table I.

Furthermore, we also reduce the sizes of the input coefficients. This way, we
can obtain an even lower final integer width. In all the experiments, the
MNIST data is scaled down from its typical $8\text{\,}\mathrm{bit}$ to
$2\text{\,}\mathrm{bit}$. However, since CIFAR images are more complex, the
same reduction could lead to unacceptable accuracies. Therefore we chose not
to reduce the CIFAR dataset.

The results in Table I show that for both networks, we can quantise until
$8\text{\,}\mathrm{bit}$ without an accuracy drop. When we want to use lower
quantisation, the accuracy starts to drop. One of the reasons is that in these
cases, many of the weights are quantised to zero, which causes much of
information to ‘disappear’ and results in a diminished FIW.

Figure 1: Overview of the weight distribution of the MNIST and CIFAR architecture.  TABLE I: Results of the quantised model using Brevitas and post-training quantisation with different scale factors for the architectures for both MNIST and CIFAR. |  | PTQ – $S=\frac{\max(|\boldsymbol{W}|)}{2^{b-1}-1}$ | PTQ – $S=\frac{\beta-\alpha}{2^{b-1}-1}$ | PTQ – $S=\frac{1}{2^{b-1}-1}$ | QAT
---|---|---|---|---|---
Dataset | Quantisation | Acc [%] | FIW [bit] | Acc [%] | FIW [bit] | Acc [%] | FIW [bit] | Acc [%] | FIW [bit]
| $32\text{\,}\mathrm{bit}$ | 98.43 | 237 | 98.43 | 231 | 98.43 | 233 | 98.41 | 238
MNIST | $8\text{\,}\mathrm{bit}$ | 98.41 | 69 | 98.44 | 63 | 98.44 | 64 | 98.3 | 70
$3\text{\,}\mathrm{bit}$ | 94.39 | 32 | 44.03 | 26 | 79.04 | 27 | 98.3 | 38
| $2\text{\,}\mathrm{bit}$ | 14.4 | 20 | 11.35 | 3 | 11.53 | 6 | 98.46 | 29
| $32\text{\,}\mathrm{bit}$ | 73.09 | 583 | 73.09 | 571 | 73.09 | 570 | 73.28 | 614
CIFAR | $8\text{\,}\mathrm{bit}$ | 73.0 | 202 | 73.18 | 187 | 73.09 | 186 | 73.04 | 205
| $4\text{\,}\mathrm{bit}$ | 52.24 | 135 | 18.67 | 123 | 9.96 | 128 | 72.49 | 153

##  IV Quantisation-aware training (QAT)

In the previous section, we showed that neural networks could be quantised to
$8\text{\,}\mathrm{bit}$ integers, but the accuracy reduces for a more drastic
quantisation. To reduce the bit width of the network further, we can make the
network aware of the quantisation during its training. Before training is
started, the quantisation technique and parameters are chosen and introduced
into the training graphs as ‘fake quantisation’ nodes, which simulate the low-
precision behaviour of the quantisation. These nodes quantise a real input
using Equation 1 and immediately perform a dequantisation immediately
afterwards, thus injecting an error that the quantisation would cause.
Depending on the quantisation used, this method can result in a network with
approximately the same accuracy as a full-precision network while using low-
precision parameters.

Using Brevitas [15], we trained the same networks as in the post-training
quantisation experiments of the previous section. Brevitas is a library to
develop and train quantisation-aware hardware-ready networks. We used its
‘weights-only quantisation process’, in which a quantisation error is
exclusively injected in the weights.

The results can be seen on the right in Table I. For lower bit widths, the
accuracy of the QAT is significantly better than the PTQ case, remaining
approximately the same as the full precision network. One of the reasons is
that quantisation-aware training will prevent parameter sparsity and ensure
that each parameter is used correctly.

TABLE II: Results of the quantisation-aware training and homomorphic inference. |  Network |  |  Quant |  Acc. [%] |  FIW [$\text{\,}\mathrm{bit}$] |  Seq. time [min.] |  No. of inst.
---|---|---|---|---|---|---|---
|  CryptoNets |  [7] |  5-$10\text{\,}\mathrm{bit}$ |  99 |  80 | - |  2
MNIST |  HCNN |  [1] |  $4\text{\,}\mathrm{bit}$ |  99 |  43 | - |  1
|  Our Work |  |  $32\text{\,}\mathrm{bit}$ |  98.6 |  238 |  98.45 |  7
|  Our Work |  |  $8\text{\,}\mathrm{bit}$ |  98.51 |  70 |  41.63 |  3
|  Our Work |  |  $4\text{\,}\mathrm{bit}$ |  98.65 |  45 |  12.2 |  1
|  Our Work |  |  $2\text{\,}\mathrm{bit}$ |  98.46 |  29 |  8.78 |  1
|  LoLa |  [11] |  8-$9\text{\,}\mathrm{bit}$ |  74.1 |  93 | - |  4
CIFAR |  HCNN |  [1] |  $8\text{\,}\mathrm{bit}$ |  77.55 |  218 | - |  10
|  Our Work |  |  $32\text{\,}\mathrm{bit}$ |  73.28 |  614 |  12801 |  30
|  Our Work |  |  $8\text{\,}\mathrm{bit}$ |  73.04 |  205 |  4267 |  10
|  Our Work |  |  $4\text{\,}\mathrm{bit}$ |  72.49 |  153 |  3413 |  8
|  Our work |  |  $2\text{\,}\mathrm{bit}$ |  69.14 |  124 |  2560 |  6

Table II compares our QAT network to earlier works. Notable is that the
CryptoNets and HCNN implementations use PTQ techniques, but no QAT techniques.
Using QAT, we can quantise the weights to as low as $2\text{\,}\mathrm{bit}$,
giving the network a much lower FIW with minimal to no drop in accuracy.
Compared to a full precision network, i.e. quantising the parameters to
$32\text{\,}\mathrm{bit}$ integers, the FIW is reduced with a factor of 8.2
for the MNIST network and a factor of 5 for the CIFAR network. Compared to the
numbers presented by HCNN, our smallest networks have a 33% and 43% smaller
final integer width for MNIST and CIFAR, respectively, while boasting similar
accuracy.

##  V Evaluation

In this section, we evaluate our newly developed quantised neural networks by
implementing them using the Pyfhel [16] library, which is a software package
that provides python-bindings for Microsoft’s SEAL library [17]. An encrypted
inference is executed using the integer-based BFV scheme. All of our tests are
run using Python 3.9.13, Brevitas 0.7.1, Pyfhel 3.3.1 (using SEAL 3.7), on an
Intel Xeon Silver 4208 CPU.

One of the most compelling optimisations in certain HE schemes was the
introduction of batching or packing, as described by Smart and Vercauteren
[18]. It provides a way to pack multiple plaintext messages into a single
ciphertext as if it were a vector of plaintexts. In our implementation, we use
batching to pack each input channel into a single ciphertext. A single
ciphertext is used for a (black-and-white) MNIST image, and three ciphertexts
are needed for an (RGB) CIFAR image.

Due to batching, we cannot implement a dot-based matrix-vector multiplication
since we need access to the individual elements of a ciphertext. Therefore
rotation-based versions of each neural network layer are implemented based on
previous works. Dathathri et al. [19] propose an algorithm to calculate a
single convolution kernel on a subset of the input data. We adapted the
algorithm further to enable us to simultaneously apply an input kernel to a
complete channel. As proposed by Juvekar et al. [20], a rotation-based
algorithm is used to execute a matrix-vector multiplication for the dense
layer. This algorithm will perform the multiplication using only vector
addition, multiplications and rotations.

When converting to an almost binary size ($2\text{\,}\mathrm{bit}$), extra
sparsity is introduced, which we use to reduce the latency. Before encoding
the weight during the encrypted inference, we check for a zero vector. If
there is one, all the associated operations using this vector can be omitted.
This results in a speedup of 28% between a $4\text{\,}\mathrm{bit}$ and
$2\text{\,}\mathrm{bit}$ network.

###  V-A Homomorphic encryption parameter selection

To determine and select suitable HE parameters, we first analyse the final
integer width that determines whether we need multiple instances. The SEAL
library limits the maximum size of the plaintext modulus to
$60\text{\,}\mathrm{bit}$ for performance reasons. Due to the outcomes of our
QAT experiments, we need to represent larger plaintext spaces and use a
residue numeral system (RNS). An overview of the used HE parameters is given
in Table III.

TABLE III: Used HE parameters Network | Quantisation | N | $\log\;q$ |  Plaintext modulus
---|---|---|---|---
MNIST | $8\text{\,}\mathrm{bit}$ | $2^{14}$ | 389 |  35184371138561,
$4\text{\,}\mathrm{bit}$ | $2^{14}$ | 389 |  35184371138561
$2\text{\,}\mathrm{bit}$ | $2^{14}$ | 389 |  1073643521
CIFAR | $8\text{\,}\mathrm{bit}$ | $2^{15}$ | 825 |  Same as $4\text{\,}\mathrm{bit}$ \+
8257537, 6946817
$4\text{\,}\mathrm{bit}$ | $2^{15}$ | 825 |  Same as $2\text{\,}\mathrm{bit}$ \+
5767169, 6750209
$2\text{\,}\mathrm{bit}$ | $2^{15}$ | 825 |  1376257, 1769473, 2424833, 2752513, 3604481, 3735553

###  V-B Results

We report the sequential times for various quantisation on the right in Table
II. To account for the number of instances, the ‘sequential time’ is given,
corresponding to the total time when each instance is executed sequentially
and reflecting the use of computing resources. For the CIFAR network, the work
of Badawi et al. [1] uses ten instances, each possessing a plaintext size of
around $21\text{\,}\mathrm{bit}$. Using the same sizes, our smallest network
($2\text{\,}\mathrm{bit}$) only requires six instances.

The MNIST architecture’s smallest network is 80% faster than the best
$8\text{\,}\mathrm{bit}$ PTQ network. This is due to the smaller FIW and
because it uses the additional sparsity of the weights. As for the CIFAR
architecture, we obtain a 40% speedup compared to the $8\text{\,}\mathrm{bit}$
PTQ network, which is equal to the quantisation used by HCNN.

##  VI Conclusion

The absence of a division operation in some fully homomorphic encryption
schemes implies that variables keep growing during computations. In this work,
we tested two main quantisation techniques to reduce the size of the internal
variables, which in turn affects computation cost. We first looked at the
limitations of post-training quantisation and showed that there is a lower
limit to the quantisation (in our case $8\text{\,}\mathrm{bit}$) before the
accuracy significantly drops. To further reduce the variable sizes, we
developed a quantisation-aware training framework. We reduced the final
integer width with 33% for MNIST and 43% for CIFAR, compared to the state-of-
the-art HCNN architecture. In our experiments, the quantisation aware
training, allowing for reducing the network weights up to
$2\text{\,}\mathrm{bit}$, boasts an 80% and 40% speedup for the MNIST and
CIFAR network, respectively, over typical $8\text{\,}\mathrm{bit}$ weights
obtained with post-training quantisation.

## References

  * Badawi et al. [2021] A. A. Badawi, C. Jin, J. Lin, C. F. Mun, S. J. Jie, B. H. M. Tan, X. Nan, K. M. M. Aung, and V. R. Chandrasekhar, “Towards the AlexNet moment for homomorphic encryption: HCNN, the first homomorphic CNN on encrypted data with GPUs,” _IEEE Trans. Emerg. Top. Comput._ , vol. 9, no. 3, pp. 1330–1343, 2021.
  * Fan and Vercauteren [2012] J. Fan and F. Vercauteren, “Somewhat practical fully homomorphic encryption,” _IACR Cryptol. ePrint Arch._ , p. 144, 2012.
  * Jacob et al. [2018] B. Jacob, S. Kligys, B. Chen, M. Zhu, M. Tang, A. G. Howard, H. Adam, and D. Kalenichenko, “Quantization and training of neural networks for efficient integer-arithmetic-only inference,” in _CVPR_. Computer Vision Foundation / IEEE Computer Society, 2018, pp. 2704–2713.
  * Moons et al. [2016] B. Moons, B. D. Brabandere, L. V. Gool, and M. Verhelst, “Energy-efficient ConvNets through approximate computing,” in _WACV_. IEEE Computer Society, 2016, pp. 1–8.
  * Vandersteegen et al. [2021] M. Vandersteegen, K. Van Beeck, and T. Goedemé, “Integer-only cnns with 4 bit weights and bit-shift quantization scales at full-precision accuracy,” _Electronics_ , vol. 10, no. 22, 2021. [Online]. Available: https://www.mdpi.com/2079-9292/10/22/2823
  * Horowitz [2014] M. Horowitz, “1.1 computing’s energy problem (and what we can do about it),” in _ISSCC_. IEEE, 2014, pp. 10–14.
  * Gilad-Bachrach et al. [2016] R. Gilad-Bachrach, N. Dowlin, K. Laine, K. E. Lauter, M. Naehrig, and J. Wernsing, “CryptoNets: Applying neural networks to encrypted data with high throughput and accuracy,” in _ICML_ , ser. JMLR Workshop and Conference Proceedings, vol. 48. JMLR.org, 2016, pp. 201–210.
  * Bos et al. [2013] J. W. Bos, K. E. Lauter, J. Loftus, and M. Naehrig, “Improved security for a ring-based fully homomorphic encryption scheme,” in _IMACC_ , ser. Lecture Notes in Computer Science, vol. 8308. Springer, 2013, pp. 45–64.
  * Albrecht et al. [2016] M. R. Albrecht, S. Bai, and L. Ducas, “A subfield lattice attack on overstretched NTRU assumptions - cryptanalysis of some FHE and graded encoding schemes,” in _CRYPTO (1)_ , ser. Lecture Notes in Computer Science, vol. 9814. Springer, 2016, pp. 153–178.
  * LeCun et al. [1998] Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner, “Gradient-based learning applied to document recognition,” _Proc. IEEE_ , vol. 86, no. 11, pp. 2278–2324, 1998.
  * Brutzkus et al. [2019] A. Brutzkus, R. Gilad-Bachrach, and O. Elisha, “Low latency privacy preserving inference,” in _ICML_ , ser. Proceedings of Machine Learning Research, vol. 97. PMLR, 2019, pp. 812–821.
  * Krizhevsky and Hinton [2009] A. Krizhevsky and G. Hinton, “Learning multiple layers of features from tiny images,” University of Toronto, Toronto, Ontario, Tech. Rep. 0, 2009.
  * Dowlin et al. [2017] N. Dowlin, R. Gilad-Bachrach, K. Laine, K. E. Lauter, M. Naehrig, and J. Wernsing, “Manual for using homomorphic encryption for bioinformatics,” _Proc. IEEE_ , vol. 105, no. 3, pp. 552–567, 2017.
  * Han et al. [2016] S. Han, H. Mao, and W. J. Dally, “Deep compression: Compressing deep neural network with pruning, trained quantization and huffman coding,” 2016.
  * Pappalardo [2021] A. Pappalardo, “Xilinx/brevitas,” 2021. [Online]. Available: https://doi.org/10.5281/zenodo.3333552
  * Ibarrondo and Viand [2021] A. Ibarrondo and A. Viand, “Pyfhel: Python for homomorphic encryption libraries,” in _WAHC ’21: Proceedings of the 9th on Workshop on Encrypted Computing & Applied Homomorphic Cryptography, Virtual Event, Korea, 15 November 2021_. WAHC@ACM, 2021, pp. 11–16. [Online]. Available: https://doi.org/10.1145/3474366.3486923
  * [17] “Microsoft SEAL (release 4.0),” https://github.com/Microsoft/SEAL, Mar. 2022, microsoft Research, Redmond, WA.
  * Smart and Vercauteren [2011] N. P. Smart and F. Vercauteren, “Fully homomorphic SIMD operations,” _IACR Cryptol. ePrint Arch._ , p. 133, 2011.
  * Dathathri et al. [2019] R. Dathathri, O. Saarikivi, H. Chen, K. Laine, K. E. Lauter, S. Maleki, M. Musuvathi, and T. Mytkowicz, “CHET: an optimizing compiler for fully-homomorphic neural-network inferencing,” in _PLDI_. ACM, 2019, pp. 142–156.
  * Juvekar et al. [2018] C. Juvekar, V. Vaikuntanathan, and A. P. Chandrakasan, “GAZELLE: A low latency framework for secure neural network inference,” in _USENIX Security Symposium_. USENIX Association, 2018, pp. 1651–1669.