# FlowPrecision: Advancing FPGA-Based Real-Time Fluid Flow Estimation with
Linear Quantization

Tianheng Ling, Julian Hoever, Chao Qian, Gregor Schiele  Intelligent Embedded
Systems Lab, University of Duisburg-Essen, 47057 Duisburg, Germany
{tianheng.ling, chao.qian<EMAIL_ADDRESS><EMAIL_ADDRESS>due.de

###### Abstract

In industrial and environmental monitoring, achieving real-time and precise
fluid flow measurement remains a critical challenge. This study applies linear
quantization in FPGA-based soft sensors for fluid flow estimation,
significantly enhancing Neural Network model precision by overcoming the
limitations of traditional fixed-point quantization. Our approach achieves up
to a 10.10% reduction in Mean Squared Error and a notable 9.39% improvement in
inference speed through targeted hardware optimizations. Validated across
multiple data sets, our findings demonstrate that the optimized FPGA-based
quantized models can provide efficient, accurate real-time inference, offering
a viable alternative to cloud-based processing in pervasive autonomous
systems.

###### Index Terms:

Flow Estimation, Soft Sensors, In-situ AI, Quantization, FPGA-based Inference,
Embedded Systems

##  I Introduction

Pervasive in-situ AI systems have immense potential in various applications,
such as industrial and environmental monitoring. One critical task in these
settings is accurately measuring fluid flow, which is challenging,
particularly in harsh environments such as drilling pipes or sewers. Direct-
contact physical sensors often degrade quickly and require frequent
maintenance, while non-contact sensors may lack precision or may be cost-
prohibitive. An emerging solution to this challenge is to use _soft sensors_ ,
which utilize data from cheap, non-contact sensors to estimate flow rate
accurately.

Recent advancements in flow estimation have utilized physical formulas,
statistical models, and deep learning algorithms. Among them, Tomperi et al.
[2] employed multivariate linear regression with cost-efficient distance
sensors for wastewater flow estimation. Chhantyal et al. [8] integrated data
fusion with Artificial Neural Networks (ANNs) for ultrasonic level sensor-
based flow estimation. Noori et al. [4] combined shallow ANNs with radar
sensors, and Liang et al. [7] developed models for complex flow scenarios
using ultrasonic sensors. However, these approaches primarily focus on cloud-
based soft sensors, susceptible to timing and network dependency issues,
potentially affecting service reliability.

In response to the need for autonomy and reliability in soft sensor
deployment, our study shifts the focus to deploying ANN-based soft sensor
algorithms directly on Internet of Things (IoT) devices. This local analysis
approach mitigates network dependencies but introduces new challenges due to
the limited computational resources available on IoT devices. These
constraints necessitate the optimization of ANN models for efficient
execution, a process we address through model quantization. Specifically, we
propose quantizing ANN models for deployment on embedded Field Programmable
Gate Arrays (FPGAs) using linear quantization with adaptive scalers. This
method aligns quantization parameters with the actual data distributions of
each tensor, thereby enhancing model precision and balancing energy
efficiency. Our key contributions include:

  * •

We apply linear quantization to Multi-Layer Perceptron (MLP) models for fluid
flow estimation, reducing test loss by up to 10.10%.

  * •

We conduct an ablation study to validate our approach against precision loss
issues in fixed-point quantization.

  * •

We test our method across three data sets, demonstrating its broad
applicability and robustness.

  * •

We develop a specialized hardware accelerator that enhances model precision
while maintaining energy efficiency and acceptable resource utilization
overhead.

The reminder content is organized as follows: Section II reviews pertinent
literature, setting the stage for our research. Section III delves into the
technical aspects of fixed-point and linear quantization. Section IV details
our approach and implementation. Section V describes the experimental setup,
and Section VI presents and discusses our findings. Finally, Section VII
concludes the paper and outlines future research.

##  II Related Work

Accurate flow estimation using soft sensors requires high-bit quantization to
ensure precision, mainly due to its regression-oriented nature. Ling et al.
[6] advanced in this area by quantizing an MLP model for FPGA platforms,
employing a uniform power-of-2 scaling factor in their 8-bit fixed-point
quantization approach. While their approach showcased FPGA deployments’ speed
and energy efficiency, it faltered in precision when compared to linear
quantized models on Microcontroller Units (MCUs). The primary limitation of
fixed-point quantization stems from its assumption of uniform tensor
distribution, leading to notable precision loss amidst real-world data
variability. In our pursuit to enhance model precision for soft sensor
applications on FPGAs, we also examined the work of Jain et al. [11]. Their
introduction of per-tensor scaling considerably ameliorated compatibility
across diverse hardware platforms and more effectively adapted to the
variability inherent in tensor distributions. Beyond refining the fixed-point
quantization, linear quantization is a formidable alternative. This technique
applies a zero point and scaling factor to convert floating-point tensors into
integer tensors linearly. Gaining traction in prominent deep learning
frameworks like TensorFlow and PyTorch, linear quantization is epitomized by
the Brevitas library [10] in PyTorch, showcasing enhanced model precision
through its application. tgbhu8hu86tg

Our work builds upon these works, aiming to refine model quantization for FPGA
deployments to enhance the precision of soft sensors in flow estimation. We
employ adaptive scaling, dynamically adjusting quantization parameters in
alignment with the actual data distributions within models. This approach not
only overcomes the constraints of traditional power-of-2 scaling but also aims
to elevate numerical precision. To counterbalance the computational overhead
inherent in adaptive scaling, we incorporated targeted FPGA optimizations,
ensuring energy efficiency and performance efficacy.

##  III Fundamentals

Quantization in numerical representation involves converting continuous real
numbers from the domain $\mathbb{R}$ into discrete values within the quantized
domain $\mathbb{Q}$. Consider a tensor $X$, composed of real-valued elements
$x$, and its quantized counterpart $X_{q}$. As depicted in Equation 1, the
scale factor $S$, a floating-point parameter, defines the relationship between
$X$ and $X_{q}$, while the zero point $Z$ is a fixed-point parameter
representing zero in $X$. The rounding operation approximates $x$ to the
nearest integer, and the clamp function ensures $x_{q}$ remains within the
bounds of a $b$-bit signed integer $(-2^{b-1},2^{b-1}-1)$. De-quantization,
detailed in Equation 2, reverses this process, converting $X_{q}$ back into an
approximate real-valued tensor $X^{\prime}$, using the same scale factor $S$
and zero point $Z$.

$\displaystyle x\mapsto x_{q}$
$\displaystyle=\text{clamp}\left(\text{round}(\frac{x}{S})+Z,-2^{b-1},2^{b-1}-1\right)$
(1) $\displaystyle x_{q}\mapsto x^{\prime}$ $\displaystyle=S\cdot(x_{q}-Z)$
(2)

###  III-A Fixed-point Quantization

Fixed-point quantization, commonly used in FPGA deployments, is known for
simplicity and deterministic nature. It assigns a fixed bit-width to each
number, divided into integer and fractional parts. This format is represented
as $(a,b)$, where $a$ and $b$ denote bits for fractional values and total bit
count, respectively. The scale factor $S$ is preset to $2^{-a}$, with no zero
point $Z$. Despite its utility, fixed-point quantization has limitations. It
offers a narrower range than floating-point numbers, risking overflow and
precision loss. It may allocate bits inefficiently, covering values outside
the needed range and reducing precision. Additionally, it presumes a symmetric
data distribution around zero, a condition rarely met in real-world data,
leading to quantization errors. These drawbacks highlight the necessity for a
more adaptable quantization method to meet varying precision demands in
practical applications.

###  III-B Linear Quantization

In addressing the limitations of fixed-point quantization, our approach
utilizes linear quantization with adaptive parameters, as per [9]. This method
dynamically computes the scaling factors $S$ and zero points $Z$ during
training, customizing them to each tensor’s distribution. Equations 3 and 4
illustrate how these parameters are computed based on the tensor’s observed
minimum $\alpha$ and maximum $\beta$ values. This adaptive approach ensures a
more accurate numerical representation, especially for data distributions that
are not symmetric around zero, thereby enhancing model performance by
minimizing quantization errors and optimizing bit precision utilization.

$\displaystyle S$ $\displaystyle=\frac{\alpha-\beta}{2^{b}-1}$ (3)
$\displaystyle Z$
$\displaystyle=\text{clamp}\left(\text{round}((2^{b-1}-1)-\frac{\alpha}{S}),-2^{b-1},2^{b-1}-1\right)$
(4)

##  IV Our Approach

Our approach utilizes an MLP model, similar to those in [4, 6], comprising a
hidden layer with a ReLU activation function and an output layer. Adapted for
FPGA deployment, our method facilitates more flexible linear quantization
within the PyTorch framework. Consistent with the methodologies outlined in
[6], we vary the number of neurons in the hidden layer (10, 30, 60, and 120)
to explore the effects of model complexity on performance. Next, we will
describe our approach from both software and hardware implementation
perspectives.

###  IV-A Customized Software Implementation

####  IV-A1 QAT with Linear Quantization

Quantization-Aware Training (QAT) is a favored quantization strategy that
enhances ANNs’ resilience to quantization errors. This technique simulates the
quantization process during training, incorporating both quantization and de-
quantization steps. The non-differentiable nature of rounding operations in
backpropagation is addressed using the Straight-Through Estimator [3]. We
transition from the theoretical underpinnings of QAT to its practical
application in our model, identifying the primary quantization objects, as
detailed in Table LABEL:tab:quantization_parameters. These include the inputs
$X$, weights $W^{1}$, biases $B^{1}$, and activations $A_{1}$ for the hidden
layer, and inputs $A_{2}$, weights $W_{2}$, biases $B_{2}$, and final outputs
$Y$ for the output layer. Unlike fixed-point quantization with static scale
factors, our approach dynamically calculates the quantization parameters after
each training batch, aligning them with the actual data distribution to
minimize quantization errors. The specific quantization parameters for each
object are concisely summarized in Table LABEL:tab:quantization_parameters.
For example, the quantization scale factor for the hidden layer’s inputs is
$S_{X}$, with a corresponding zero point of $Z_{X}$. Notably, all biases have
no zero points for computational simplicity. Although our implementation
supports mixed-precision quantization, all quantization objects in this study
are uniformly quantized at 8 bits for model consistency and simplified
analysis.

TABLE I: Quantization Description

Layers |  | Quantization
---
Objects
| Quantization
---
Parameters
Hidden Layer | $X$ |  $S_{X}$, $Z_{X}$
$W^{1}$ |  $S_{W^{1}}$, $Z_{W^{1}}$
$B^{1}$ | $S_{B^{1}}$
$A^{1}$ |  $S_{A^{1}}$, $Z_{A^{1}}$
ReLU | $A^{1}$ |  $S_{A^{1}}$, $Z_{A^{1}}$
$A^{2}$ |  $S_{A^{1}}$, $Z_{A^{1}}$
Output Layer | $A^{2}$ |  $S_{A^{1}}$, $Z_{A^{1}}$
$W^{2}$ |  $S_{W^{2}}$, $Z_{W^{2}}$
$B^{2}$ | $S_{B^{2}}$
$Y$ |  $S_{Y}$, $Z_{Y}$

In our approach, we adapt the ReLU function’s implementation. During QAT, we
employ PyTorch’s standard torch.nn.functional.relu function. For subsequent
integer-only inference, the quantization parameters of the ReLU function’s
inputs are inherited from the outputs of the preceding layer, i.e., scale
$S_{A^{1}}$ and zero point $Z_{A^{1}}$. We use the same quantization
parameters for the outputs of the ReLU function to ensure uniformity and
coherence in quantization across inputs and outputs. Similarly, for the inputs
of the output layer, the inputs $A_{2}$ inherit their quantization parameters
directly from the hidden layer’s outputs $A_{1}$. This approach maintains a
consistent and coherent quantization method throughout the model, facilitating
the integrity of the quantization process across varying layers.

####  IV-A2 Enhanced Integer-only Model Inference

For deployment on resource-constrained FPGAs, which do not support floating-
point operations, we have implemented an integer-only model inference
approach. This ensures that all computations and data transfers between layers
are conducted using integers, preserving the model’s efficiency and precision.

##### Integer-only Linear Layers

We employ the integer-only linear layer calculation method proposed by Benoit
et al. [5]. This method is applied starting from the hidden layer of our
model. The computation leverages quantization parameters determined during the
QAT phase, which are detailed in Equation 5. To streamline the calculation, we
approximate the bias term as $S_{B^{1}}B^{1}_{q}\approx
S_{X}S_{W^{1}}B^{1*}_{q}$, leading to the formulation presented in Equation 6.
Notably, the term $\frac{S_{X}S_{W_{1}}}{S_{A^{1}}}$ in this equation is the
sole floating-point component. To maintain our commitment to integer-only
computations, we employ bit-shift operations for approximating this term, as
elaborated in Equation 8. Specifically, we first convert the floating-point
term into a positive integer $M_{0}$, followed by a right bit-shift operation
by $n$ positions to closely approximate the original floating-point value.
This integer-only computation approach is also consistently applied to the
output layer, as indicated in Equation 9, ensuring a uniform integer-only
processing methodology throughout our model.

$\footnotesize S_{A^{1}}(A^{1}_{q}-Z_{A^{1}})\approx
S_{X}(X_{q}-Z_{X})S_{W^{1}}(W^{1}_{q}-Z_{W^{1}})+S_{B^{1}}B^{1}_{q}$ (5)
$\footnotesize S_{A^{1}}(A^{1}_{q}-Z_{A^{1}})\approx
S_{X}(X_{q}-Z_{X})S_{W^{1}}(W^{1}_{q}-Z_{W^{1}})+S_{X}S_{W^{1}}B^{1*}_{q}$ (6)
$\footnotesize
A^{1}_{q}\approx\frac{S_{X}S_{W^{1}}}{S_{A^{1}}}((X_{q}-Z_{X})(W^{1}_{q}-Z_{W^{1}})+B^{1*}_{q})+Z_{A^{1}}$
(7) $\footnotesize\frac{S_{X}S_{W^{1}}}{S_{A^{1}}}\approx 2^{-n}M_{0}$ (8)
$\footnotesize
Y_{q}\approx\frac{S_{A^{1}}S_{W^{2}}}{S_{Y}}((A^{2}_{q}-Z_{A^{1}})(W^{2}_{q}-Z_{W^{2}})+B^{2*}_{q})+Z_{Y}$
(9)

##### Integer-only ReLU

In our implementation, we align with TensorFlow Lite’s approach for the
integer-only ReLU function. Specifically, we use $Z_{A^{1}}$ as a threshold,
as formulated in Equation 10. This threshold is designed to align with the
zero value in the floating-point space, ensuring that the behavior of our
integer-only ReLU function closely mirrors that of its floating-point
counterpart. This method is critical in maintaining the precision of the model
during quantization.

$\footnotesize A^{2}_{q}\approx\max(Z_{A^{1}},A^{1}_{q})$ (10)

###  IV-B Optimized Model Inference on FPGAs

To effectively integrate linear quantization within an MLP accelerator,
ensuring efficiency on par with fixed-point quantization, we have implemented
several key optimizations.

####  IV-B1 Linear Layer Optimization

Our approach to optimizing the linear layer architecture draws inspiration
from the designs presented in [1]. We have adapted their VHDL templates to
meet our specific needs for integer-only inference.

##### Configurable Parameters

The VHDL template in our study incorporates a range of configurable
parameters: $M_{0}$, $n$, $Z_{X}$, $Z_{W}$, and $Z_{Y}$. These parameters are
crucial for achieving integer-only inference and align closely with those in
the software quantization framework, ensuring seamless integration between
hardware and software components. Focusing on the output layer as a case in
point, the parameters $M_{0}$ and $n$ are instrumental for bit-shifting
operations that hone in on the precise approximation of the scale factor
$\frac{S_{A^{1}}S_{W^{2}}}{S_{Y}}$. Here, $Z_{X}$, $Z_{W}$, and $Z_{Y}$
meticulously align with $Z_{A^{1}}$, $Z_{W^{2}}$, and $Z_{Y}$ respectively,
thus maintaining the integrity of our quantization approach throughout the
computational process.

##### Pipelined Matrix Multiplication

To exploit the parallelization capabilities of FPGAs, we refine our approach
for linear quantization and integer-only inference. Our optimizations include
a) subtracting zero points before multiplication and accumulation (MAC)
operations and b) redesigning the MAC unit with a pipelined architecture to
accommodate complex scaling needs. These algorithmic modifications, detailed
in Algorithm 1, ensure efficient parallel task execution, reducing processing
delays and enabling higher clock frequencies. Specifically, steps 7 to 9,
involving data fetching and zero-point subtraction, are intentionally
separated from the main MAC operations. Furthermore, the bit-shifting for
approximating $M$ is strategically placed after the main computation loop,
ensuring efficient scaling application just before populating the output
buffer $Y$. Significantly, steps 3 to 5 and 12 to 14 are designed to be
executed efficiently within a single 100MHz clock cycle, even though they are
in distinct stages.

Algorithm 1 MAC Algorithm in the Linear Layer

0: $x$ is an $K$-element vector, $W$ is an $J\times K$ matrix B is an
$J$-element vector

1: Initialization: $sum\leftarrow 0$, $j\leftarrow 0$

2: repeat

3: $k\leftarrow 0$

4: Load: $W[j][0]$, $x[0]$, $B[j]$

5: $sum\leftarrow sum+B[j]$

6: repeat

7: Load $W[j][k+1]$, $x[k+1]$

8: $w\leftarrow W[j][k]-Z_{W}$, $x\leftarrow x[k]-Z_{X}$

9: $sum\leftarrow sum+w\cdot x$

10: $k\leftarrow k+1$

11: until $k=K$

12: $y\leftarrow(sum\cdot M_{0})>>n$

13: Store: $Y[j]\leftarrow y+Z_{Y}$

14: $j\leftarrow j+1$

15: until $j=J$

15: $Y$

####  IV-B2 ReLU Optimization

The ReLU function optimization is meticulously calibrated to handle inputs on
an element-wise basis, ensuring a straightforward and delay-free operation.
Including the configurable parameter $Z_{A}$ sets a threshold for the input
tensor, enabling the comparator to adjust the output as needed. If the input
is less than $Z_{A}$, the output is set to this threshold value. Otherwise, it
retains the input value. This efficient logic ensures minimal operational
delay, allowing immediate output updates following input changes.

####  IV-B3 Network Component Integration

By integrating the linear layer and ReLU optimizations, we construct a
comprehensive network component to realize the MLP model. This component
sequentially interconnects layers with activation functions, reflecting the
data flow established in the software implementation. The templates for
assembling the MLP models are thoroughly validated through ongoing
simulations. Our hardware implementation is verified to produce outputs
consistent with those from the software implementation, as demonstrated by
regular tests conducted on the test data set.

##  V Experiments

###  V-A Data Sets and Task

In this study, we utilized three distinct data sets. The specifics of each
data set are summarized in Table II. DS1, an open data set referenced from
[8], simulates drilling conditions using an advanced flow loop system. This
system includes a mud tank and a Venturi construction, with high-fidelity data
captured at a 10kHz sampling rate from three level sensors. The data from
these sensors serves as input data, while a precise Coriolis mass flow meter
provides the target data, as depicted in Figures 1. To further strengthen our
model evaluation, we incorporated two additional data sets, DS2 and DS3,
generously provided by Viumdal, a co-author of the paper [8]. These data sets
are characterized by upward and downward flow trends, introducing additional
complexities and enriching the comprehensiveness of our model’s assessment.
All three data sets, DS1, DS2, and DS3, are divided using the same proportion:
75% of the data is allocated for training, and the remaining 25% is equally
split between validation and testing. All data sets were normalized to a range
between 0 and 1. To enhance the robustness of our model, we have employed
seven-fold cross-validation across all data sets.

TABLE II: Flow Estimation Data Sets Data Sets | Description
---|---
DS1 | 1800 samples with upward trend only
DS2 | 4439 samples with upward & downward trends
DS3 | 4985 samples with upward & downward trends
Figure 1: Data Set DS1: Measurements With Upward Trend Figure 2: Distribution
of Test MSEs Across Various Hidden Neurons

###  V-B Experimental Settings

Our model training utilized the Adam optimizer with hyperparameters set as:
$\beta_{1}=0.9$, $\beta_{2}=0.98$, and $\epsilon=10^{-9}$. We started with an
initial learning rate of 0.001, employing a scheduler with a step size of 3
and a decay factor $\gamma$ of 0.5 to adjust it. The Mean Squared Error (MSE)
was chosen as the loss function for model training and evaluation. We
conducted 100 experiments, each running for 500 epochs and with an early
stopping mechanism to avoid overfitting. Notably, the test MSE reported in
this study is calculated using de-normalized outputs and target data.

TABLE III: Models with Various Configurations Models | Description
---|---
M-Float | A 32-Bit Floating-Point Model
M-Fixed | An (6,8) Fixed-Point Quantized Model Using QAT
M-Linear | An 8-Bit Integer Quantized Model Using QAT

To elucidate the distinctions among the models with varying configurations, we
present Table III serving as a concise summary. The M-Float model, adhering to
the standard 32-bit floating-point configuration, closely mirrors the
precision levels reported in [4]. This model serves as a reference point,
helping us gauge the extent of precision loss attributable to quantization.
The M-Fixed model utilizes a conventional fixed-point (6,8) configuration
during QAT, indicating an 8-bit total width with 6 bits allocated for the
fractional part. This model is referenced as the baseline from related work
[6] for contextual purposes. Furthermore, the M-Linear model, employing linear
quantization during QAT, is converted into an 8-bit integer format,
illustrating our approach to quantization.

###  V-C Target Hardware and Model Deployment

We deployed the MLP model on a Spartan-7 FPGA for comparison, following the
process described in [6]. The process began with converting the trained models
within PyTorch framework into VHDL templates. These templates were then
tailored to embed the model’s weights and biases into ROM instances and
integrate other quantization parameters, forming the basis for the RTL
designs. Subsequently, these designs were synthesized into bitfiles using the
Vivado software suite, enabling deployment on the FPGA. We assessed the
model’s timing and power consumption using Vivado’s analytical tools, which
allow for a comparison with the estimates provided in [6], even without access
to their customized hardware.

##  VI Results and Evaluation

###  VI-A Test Loss Assessment Across Quantization Methods

In our first experiment, we conducted a comparative analysis of linear and
fixed-point quantization methods, assessing their impact on models with
varying complexities. Additionally, the M-Float models served as references.
This analysis encompassed three model types, each tested with varying numbers
of hidden neurons (10, 30, 60, and 120) in the hidden layer, resulting in 12
distinct model variations. Each model underwent 7-fold cross-validation on the
DS1 data set, with each fold involving 100 independent training runs. The
distribution of test MSEs is illustrated in Figure 2. Our findings revealed
that the M-Float models, featuring floating-point parameters and activations,
consistently achieved the lowest test MSEs, surpassing the quantized models.
Models with higher counts of hidden neurons tended to yield lower test MSEs.
Among the quantized models, the M-Linear models generally exhibited lower
average and minimum test MSEs compared to the M-Fixed models. Moreover, the
M-Linear models demonstrated greater stability, as indicated by their smaller
variance in test MSEs, unlike the M-Fixed models. These results, consistent
across various folds and model complexities, underscore the effectiveness and
reliability of linear quantization for improving model performance.

###  VI-B Ablation Study of Linear Quantization

Following the same experimental setup, we conducted an ablation study on data
set DS1 at Fold 0 to evaluate the impact of linear quantization on model
performance. In this study, each linear layer was independently quantized
using either fixed-point (F) or linear (L) quantization methods, resulting in
16 distinct model configurations. The performance of these models is
summarized in Table IV. When analyzing models with 10 hidden neurons, as shown
in the first four rows of Table IV, we observed a consistent increase in MSE
when switching from linear to fixed-point quantization. This finding, in line
with our initial experiment’s results, underscores the critical role of linear
quantization in reducing MSE. This pattern holds true for models with a higher
number of hidden neurons, further affirming the effectiveness of linear
quantization in enhancing model precision across various model sizes.

TABLE IV: Ablation Study of Linear Quantization | Hidden
---
Neurons
| Quantization Configuration
---
(Hidden Layer / Output Layer)
MSE
10 | L / L | 73.75
L / F | 82.00
F / L | 83.37
F / F | 85.86
30 | L / L | 66.10
L / F | 74.83
F / L | 72.25
F / F | 78.36
60 | L / L | 59.60
L / F | 66.36
F / L | 67.93
F / F | 71.82
120 | L / L | 62.87
L / F | 65.38
F / L | 67.36
F / F | 68.14

###  VI-C Stability Assessment on Linear Quantization

In our further experiments using data sets DS2 and DS3, we adhered to the same
experimental setup to assess the generalizability of linear quantization. The
summarized outcomes in Table V consistently highlight the superiority of
linear quantization over fixed-point quantization, echoing our earlier
observations with data set DS1. Particularly noteworthy in data set DS2 was
the up to 9.70% reduction in test MSE achieved by linear quantization in
models with 10 hidden neurons, and a notable decrease of 4.42% even in the
more complex model configuration with 120 hidden neurons. The trends observed
in data set DS3 mirrored these improvements, evident across all neuron
configurations, with the largest model configuration exhibiting a 9.40%
reduction in MSE. These consistent results across varied data sets firmly
establish the robustness and efficacy of linear quantization, significantly
enhancing the performance of quantized models.

TABLE V: Model Performance Across Data Sets Data Sets |  Hidden Neurons | MSE
---|---|---
M-Float | M-Fixed | M-Linear
DS2 | 10 | 75.38 | 82.72 | 74.70($\downarrow$9.70%)
30 | 70.74 | 77.30 | 75.55($\downarrow$2.26%)
60 | 67.73 | 74.32 | 72.34($\downarrow$2.66%)
120 | 65.51 | 73.09 | 69.86 ($\downarrow$4.42%)
DS3 | 10 | 59.49 | 64.41 | 60.03 ($\downarrow$6.80%)
30 | 55.53 | 61.14 | 57.63 ($\downarrow$5.74%)
60 | 51.17 | 61.68 | 58.22($\downarrow$5.61%)
120 | 51.60 | 59.24 | 53.67($\downarrow$9.40%)

###  VI-D Performance on FPGA Deployment

Our final set of experiments assessed the performance of models deployed on
FPGAs, with a focus on inference time, power consumption, and energy
efficiency. As detailed in Table VI, the M-Linear model consistently
demonstrated a modest reduction in inference time compared to the M-Fixed
model, across various hidden neuron configurations on DS1. Notably, this
increased processing speed was accompanied by an elevated power consumption,
often exceeding a 10% increase. This trade-off between inference speed and
power consumption is attributed to the pipelined parallelization in the linear
layer implementation. The resource utilization, outlined in Table VII, shows
that the M-Linear model requires more Look-Up Tables (LUTs) and Block RAM
(BRAM), especially as the neuron count increases. The higher use of Digital
Signal Processing slices (DSPs) is necessary for the complex scaling logic,
particularly evident at higher model complexities, as detailed in Algorithm 1,
specifically at step 12. The uptick in resource use correlates with the
model’s complexity, indicating that linear quantization enhances inference
speed through more intensive use of FPGA resources. Despite the increased
power consumption, the energy cost per inference of the M-Linear model is only
up to 5.71% higher than the M-Fixed model. Notably, this energy consumption is
still quite lower, by 56.8$\times$, compared to model inference on MCUs, as
reported in [6].

TABLE VI: Model Performance on XC7S15 FPGA | Models
---
| Hidden
---
Neurons
| Inference
---
Time ($\mu$s)
| Power
---
(mW)
| Energy
---
($\mu$J/inference)
M-Fixed | 10 | 1.04 | 28 | 0.03
30 | 3.04 | 29 | 0.09
60 | 6.04 | 29 | 0.18
120 | 12.04 | 29 | 0.35
M-Linear | 10 | 1.01 ($\downarrow$2.88%) | 31 ($\uparrow$10.71%) | 0.03 (0.0%)
30 | 2.81 ($\downarrow$7.57%) | 32 ($\uparrow$10.34%) | 0.09 (0.0%)
60 | 5.51 ($\downarrow$8.78%) | 33 ($\uparrow$13.79%) | 0.18 (0.0%)
120 | 10.91 ($\downarrow$9.39%) | 34 ($\uparrow$17.24%) | 0.37 ($\uparrow$5.71%)
TABLE VII: Resource Utilization on XC7S15 FPGA Models | Hidden Neurons |  LUT(10400) |  BRAM(10) |  DSP(20)
---|---|---|---|---
M-Fixed | 10 | 355 | 0 | 0
30 | 373 | 0 | 0
60 | 387 | 0.5 | 0
120 | 375 | 1 | 0
M-Linear | 10 | 467 | 0 | 2
30 | 463 | 0.5 | 2
60 | 479 | 1 | 2
120 | 518 | 1.5 | 2

##  VII Conclusion and Future Work

This study has demonstrated the significant advantages of linear quantization
in MLP models for FPGA deployment. Our experimental findings have shown that
adjusting quantization parameters substantially reduces MSE and improves model
stability across various neuron configurations and data sets. While linear
quantized models demand increased power consumption and resource utilization,
notable gains in inference speed justify these trade-offs, presenting a
compelling option for applications where rapid computational responses are
essential. Looking forward, we will concentrate on further refining linear
quantization techniques. Our primary objective is to limit the precision loss
in quantized models to a maximum of 5% compared to their full-precision
counterparts. In addition, we plan to explore quantization at lower bit depths
to deploy more complex MLP models for flow estimation on FPGA platforms. This
endeavor aims to enhance both the efficiency and effectiveness of these models
in real-time flow estimation.

Acknowledgments. The authors gratefully acknowledge the financial support
provided by the Federal Ministry of Economic Affairs and Climate Protection of
Germany for the RIWWER project (01MD22007C).

## References

  * [1] C. Qian, T. Ling, and G. Schiele, “Elasticai: Creating and deploying energy-efficient deep learning accelerator for pervasive computing,” in _2023 IEEE International Conference on Pervasive Computing and Communications Workshops and other Affiliated Events (PerCom Workshops)_. IEEE, 2023, pp. 297–299.
  * [2] J. Tomperi, P. M. Rossi, and M. Ruusunen, “Estimation of wastewater flowrate in a gravitational sewer line based on a low-cost distance sensor,” _Water Practice & Technology_, p. wpt2022171, 2022.
  * [3] P. Yin, J. Lyu, S. Zhang, S. Osher, Y. Qi, and J. Xin, “Understanding straight-through estimator in training activation quantized neural nets,” _arXiv preprint arXiv:1903.05662_ , 2019.
  * [4] N. Noori, T. Waag, H. Viumdal, R. Sharma, M. Jondahl, and A. Jinasena, “Non-newtonian fluid flow measurement in open venturi channel using shallow neural network time series and non-contact level measurement radar sensors,” in _SPE Norway Subsurface Conference_. OnePetro, 2020.
  * [5] B. Jacob, S. Kligys, B. Chen, M. Zhu, M. Tang, A. Howard, H. Adam, and D. Kalenichenko, “Quantization and training of neural networks for efficient integer-arithmetic-only inference,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2018.
  * [6] T. Ling, C. Qian, and G. Schiele, “On-device soft sensors: Real-time fluid flow estimation from level sensor data,” 2023.
  * [7] H. Liang, C. Song, Z. Li, and H. Yang, “Volume flow rate calculation model of non-full pipe multiphase flow based on ultrasonic sensors,” _Physics of Fluids_ , vol. 35, no. 3, 2023.
  * [8] K. Chhantyal, H. Viumdal, S. Mylvaganam, and G. Elseth, “Ultrasonic level sensors for flowmetering of non-newtonian fluids in open venturi channels: Using data fusion based on artificial neural network and support vector machines,” in _2016 IEEE Sensors Applications Symposium (SAS)_. IEEE, 2016, pp. 1–6.
  * [9] R. Krishnamoorthi, “Quantizing deep convolutional networks for efficient inference: A whitepaper,” _arXiv preprint arXiv:1806.08342_ , 2018.
  * [10] A. Pappalardo. (2023) Xilinx/brevitas. [Online]. Available: https://doi.org/10.5281/zenodo.3333552
  * [11] S. Jain, A. Gural, M. Wu, and C. Dick, “Trained quantization thresholds for accurate and efficient fixed-point inference of deep neural networks,” _Proceedings of Machine Learning and Systems_ , pp. 112–128, 2020.


# Slicing Pomerons in ultraperipheral collisions using forward neutrons from
nuclear breakup

M. Alvioli  Consiglio Nazionale delle Ricerche, Istituto di Ricerca per la
Protezione Idrogeologica, via Madonna Alta 126, I-06128, Perugia, Italy
Istituto Nazionale di Fisica Nucleare, Sezione di Perugia, via Pascoli 23c,
I-06123, Perugia, Italy  V. Guzey  University of Jyvaskyla, Department of
Physics, P.O. Box 35, FI-40014 University of Jyvaskyla, Finland and Helsinki
Institute of Physics, P.O. Box 64, FI-00014 University of Helsinki, Finland
M. Strikman  Pennsylvania State University, University Park, PA, 16802, USA

###### Abstract

We argue that measurements of forward neutrons from nuclear breakup in
inclusive high energy photon-nucleus ($\gamma A$) scattering provide a novel
complementary way to study small-$x$ dynamics of QCD in heavy-ion
ultraperipheral collisions (UPCs). Using the leading twist approximation to
nuclear shadowing, we calculate the distribution over the number of
evaporation neutrons produced in $\gamma Pb$ collisions at the LHC. We
demonstrate that it allows one to determine the distribution over the number
of wounded nucleons (inelastic collisions), which constrains the mechanism of
nuclear shadowing of nuclear parton distributions.

Heavy-ion scattering, ultraperipheral collisions, nuclear shadowing

##  I Introduction and motivation

Understanding of the QCD dynamics of hard high energy interactions and the
structure of nuclei and nucleons is one of the main directions of theoretical
and experimental studies at the Large Hadron Collider (LHC) and the
Relativistic Heavy Ion Collider (RHIC). Of particular interest is the limit of
very small momentum fractions $x$, when the linear Dokshitzer-Gribov-Lipatov-
Altarelli-Parisi (DGLAP) approximation is expected to break down
Dokshitzer:1978hw ; Gribov:1983ivg  and a regime close to the black disk limit
(BDL) Frankfurt:2001nt  may set in. Its observation is one of the prime
objectives of the planned Electron-Ion Collider (EIC) at Brookhaven National
Laboratory Accardi:2012qut ; AbdulKhalek:2021gbh , which will ultimately reach
$x$ $\sim$ 10-3 for momentum transfers of a few GeV. At the same time, it was
pointed out some time ago that ultraperipheral collisions (UPCs) of two ions
at the LHC, where a photon emitted by one of the nuclei interacts with the
other nucleus, allow one to probe down to $x\sim 10^{-5}-10^{-4}$, depending
on a particular reaction channel and the detector geometry Strikman:2005yv ;
Baltz:2007kq .

Electron-nucleus collisions at the EIC and UPCs of heavy ions at the LHC
present two options for studying small-$x$ dynamics, which are largely
complementary. At the LHC practically all data are collected for one heavy
nucleus and one cannot directly access the dependence of cross sections on the
virtuality of the probe. At the same time, one can reach very small $x$, which
is provided by a wide rapidity coverage of the LHC detectors and the large
invariant photon-nucleus collision energies exceeding by a factor of 100 the
design energies at the EIC.

Over the last decade the data taken in the LHC and RHIC kinematics discovered
a significant nuclear suppression of coherent $J/\psi$ photoproduction in Pb-
Pb and Au-Au UPCs compared to the impulse approximation prediction
ALICE:2013wjo ; ALICE:2012yye ; CMS:2016itn ; ALICE:2021gpt ; ALICE:2019tqa ;
LHCb:2022ahs ; LHCb:2021bfl ; CMS:2023snh ; ALICE:2023jgu ; STAR:2023gpk .
When interpreted in the leading twist approximation (LTA) Frankfurt:2011cs ,
it amounts to strong gluon nuclear shadowing Guzey:2013xba ; Guzey:2013qza :

$R_{A}^{g}(x,Q^{2})=\frac{g_{A}(x,Q^{2})}{Ag_{N}(x,Q^{2})}<1\,,$ (1)

where $g_{A}(x,Q^{2})$ and $g_{N}(x,Q^{2})$ are the nucleus and nucleon gluon
densities, respectively. Typical numbers reported by the LHC experiments
ALICE:2013wjo ; ALICE:2012yye ; CMS:2016itn ; ALICE:2021gpt ; ALICE:2019tqa ;
LHCb:2022ahs ; LHCb:2021bfl ; CMS:2023snh ; ALICE:2023jgu  correspond to

$\displaystyle R^{g}_{\rm Pb}(x=10^{-3},Q^{2}_{\rm eff}=3\,{\rm GeV}^{2})$
$\displaystyle\approx$ $\displaystyle 0.6\,,$ $\displaystyle R^{g}_{\rm
Pb}(x=10^{-4},Q^{2}_{\rm eff}=3\,{\rm GeV}^{2})$ $\displaystyle\approx$
$\displaystyle 0.5\,,$ (2)

with a similar suppression extending down to $x$ $\sim$ 10-5. Here $Q_{\rm
eff}$ is the effective resolution scale determined by the charm quark mass.
These values of $R^{g}_{\rm Pb}$ agree very well with the LTA predictions for
nuclear shadowing made more than 10 years ago Frankfurt:2011cs . Note that
this interpretation of the $J/\psi$ UPC data is complicated at the next-to-
leading order (NLO) of the perturbative expansion in powers of $\log Q^{2}$
(perturbative QCD) due to large cancellations between the leading-order (LO)
and NLO gluon terms, which leaves a numerically important quark contribution
Eskola:2022vpi ; Eskola:2022vaf . A way to stabilize the perturbation series
and restore the gluon dominance in this process on the proton target was
suggested in Jones:2015nna ; Jones:2016ldq .

Other hard UPC processes considered in the literature are inclusive
Strikman:2005yv ; Guzey:2018dlm  and diffractive Guzey:2016tek  dijet
photoproduction, which are still to produce final results ATLAS:2017kwa ;
ATLAS:2022cbd , timelike Compton scattering Pire:2008ea ; Schafer:2010ud ;
Peccini:2021rbt , and heavy quark photoproduction Klein:2002wm ;
Goncalves:2009ey ; Goncalves:2017zdx . Overall these findings confirm the
conclusion of Baltz:2007kq  that UPCs provide a very effective tool to access
the small-$x$ dynamics of the strong interactions and the nuclear structure in
hard, semi-hard and soft regimes.

So far the experimental studies of UPCs have mainly been focusing on coherent
and incoherent production of light and heavy vector mesons. In this paper, we
would like to outline several possible directions of future UPC studies, which
were not discussed in the review Baltz:2007kq  due to the lack of experimental
confirmations of large nuclear shadowing. We explore for the first time
possibilities of testing the small-$x$ nuclear shadowing dynamics by measuring
the rates of forward neutron production from nuclear breakup in the zero
degree calorimeters (ZDCs) at the LHC. Our numerical studies demonstrate that
the number of produced neutrons is correlated with the number of wounded
nucleons (inelastic photon-nucleon interactions), which presents a
complementary way to study the mechanism of nuclear shadowing of nuclear
parton distributions. In particular, this allows one to effectively probe the
centrality (impact parameter) dependence of nuclear shadowing by
experimentally controlling the number of unitarity cuts of diffractive
exchanges (Pomerons), which build up the leading twist nuclear shadowing.

This paper is organized as follows. In Sec. II we briefly summarize
expectations based on applications of the Abramovski-Gribov-Kancheli (AGK)
theorem to photon-nucleus scattering, model the distribution over the number
of wounded nucleons and estimate its average value in the current UPC
kinematics. Section III presents our predictions for the distributions over
the number of emitted forward neutrons from nuclear breakup in inelastic
photon-nucleus scattering and its connection to parameters of the leading
twist approximation of nuclear shadowing. Our conclusions and outlook are
given Sec. IV.

##  II Abramovski-Gribov-Kancheli (AGK) cutting rules, nuclear shadowing and
the number of wounded nucleons in $\gamma A$ scattering

It was demonstrated by Abramovski, Gribov and Kancheli in 1973
Abramovsky:1973fm  that different unitary cuts of the diagrams corresponding
to multi-Pomeron (color singlet) exchanges result in different multiplicities
of produced particles in the central rapidity region and that the absorptive
part of the amplitude can be expressed in terms of a small number of cut
diagrams, which are related by combinatorial factors; these are the so-called
Abramovski-Gribov-Kancheli (AGK) cutting rules or the AGK cancellation. For
the interpretation of the AGK rules in QCD and other effective field theories,
see Refs. Treleani:1994at ; Jalilian-Marian:2004vhw ; Gelis:2006yv ;
Kovner:2006wr ; Nikolaev:2006mx .

The application of the AGK cutting rules to photon-nucleus scattering allows
one to express the nuclear shadowing correction to the total nuclear cross
section $\sigma_{\rm tot}^{\gamma A}$ in terms of the diffractive cross
section on individual nucleons Frankfurt:1991nx ; Frankfurt:1998ym . Indeed,
the nuclear amplitude can be constructed as a series corresponding to the
interactions with $k=1,2,\dots,A$ nucleons of the target, which are mediated
by Pomeron exchanges. The interaction with one nucleon gives the impulse
approximation (IA),

$\sigma_{\rm tot}^{\gamma A(\rm IA)}\,=\,A\,\sigma_{\rm tot}^{\gamma N}\,,$
(3)

where $\sigma_{\rm tot}^{\gamma N}$ is the total photon-nucleon cross section.
For the interaction with two nucleons, the nuclear cross section is given by a
sum of cut two-Reggeon amplitudes corresponding to both Reggeons uncut
($\sigma_{0}$), one cut Reggeon ($\sigma_{1}$), and two cut Reggeons
($\sigma_{2}$),

$\sigma_{\rm tot}^{\gamma A(k=2)}=\sigma_{0}+\sigma_{1}+\sigma_{2}\,.$ (4)

Figure 1 shows the diagrams corresponding to the three terms in Eq. (4). The
AGK rules relate these cross sections as follows

$\displaystyle\sigma_{1}$ $\displaystyle=$ $\displaystyle-4\,\sigma_{0}\,,$
$\displaystyle\sigma_{2}$ $\displaystyle=$ $\displaystyle 2\,\sigma_{0}\,,$
(5)

where $\sigma_{1}$ is the absorptive (screening) correction to the IA cross
section, and $\sigma_{2}$ corresponds to the inelastic interaction with both
nucleons. In Eq. (5), we neglected the small real part of the amplitude for
simplicity.

Substituting Eq. (5) in Eq. (4), one finds for the nuclear shadowing
correction $\Delta\sigma_{\rm tot}^{\gamma A}$ in the approximation, where one
keeps the interaction with only one and two nucleons,

$\Delta\sigma_{\rm tot}^{\gamma A}=\sigma_{\rm tot}^{\gamma A}-\sigma_{\rm
tot}^{\gamma A(\rm IA)}=\sigma_{\rm tot}^{\gamma A(k=2)}=-\sigma_{0}\,.$ (6)

This relation constitutes the foundation of the Gribov model of nuclear
shadowing for the total cross section of pion-deuteron scattering, which
relates the nuclear shadowing correction to the pion-nucleon diffractive cross
section Gribov:1968jf . Note that the AGK rules were derived later and were
not explicitly used in Gribov:1968jf . For the generalization of this approach
for the calculation of total cross sections of photon scattering on heavier
nuclei, see Gribov:1968gs ; Piller:1999wx ; Adeluyi:2006xy ; Adeluyi:2007wk .

Figure 1: Unitary cuts of the two-Reggeon forward amplitude, which gives rise
to the nuclear shadowing correction due to the interaction with $k=2$ nucleons
of the nuclear target: $(a)$ the diffractive cut with both Reggeons uncut and
the cross section $\sigma_{0}$, $(b)$ one cut Reggeon corresponding to a
single multiplicity and the cross section $\sigma_{1}$, and $(c)$ two cut
Reggeons corresponding to a double multiplicity and the cross section
$\sigma_{2}$.

Equation (6) is general and does not require an expansion over twists or the
partonic interpretation. However, in lepton-proton deep inelastic scattering
(DIS), $\sigma_{0}$ is a leading twist cross section, which by the virtue of
the QCD collinear factorization theorem for hard diffraction Collins:1997sr
can be interpreted in terms of diffractive parton distribution functions
(PDFs). As a result, it allows one to formulate the leading twist
approximation for small-$x$ nuclear PDFs of individual flavors (valence and
sea quarks and gluons) Frankfurt:2011cs ; Frankfurt:1998ym ; Frankfurt:2003zd
.

Another important application of the AGK cutting rules involves the total
hadron-nucleus inelastic cross section defined as the difference between the
total and total elastic (coherent plus incoherent) cross sections, which is
obtained using the unitary form of the Glauber theory Bertocchi:1976bq .
Generalizing this result to photon-nucleus scattering by taking into account
that the photon takes part in the strong interactions by means of its hadronic
fluctuations, the photon-nucleus total inelastic cross section can be
presented in the following form Alvioli:2016gfo

$\sigma_{\rm inel}^{\gamma A}=\sum_{\nu=1}^{A}\sigma_{\nu}\,,$ (7)

where

$\sigma_{\nu}=\frac{A!}{(A-\nu)!\nu!}\int d^{2}\vec{b}\int d\sigma
P_{\gamma}(\sigma)(\sigma_{\rm inel}T_{A}(\vec{b}))^{\nu}(1-\sigma_{\rm
inel}T_{A}(\vec{b}))^{A-\nu}\,.$ (8)

In Eq. (8), $\vec{b}$ is the impact parameter (transverse coordinate) of the
interacting nucleon, $T_{A}(\vec{b})=\int dz\rho_{A}(\vec{b},z)$, where
$\rho_{A}(\vec{b},z)$ is the nuclear density normalized to unity, and
$\sigma_{\rm inel}=0.85\,\sigma$ is the inelastic cross section for the
interaction of a hadronic fluctuation of the photon with a target nucleon,
which is based on the estimate that the $\rho$ meson-nucleon elastic cross
section constitutes approximately 15% of the total one. Unlike the sign-
alternating Glauber series for the nuclear total cross section Franco:1965wi ,
the cross sections $\sigma_{\nu}$ in Eqs. (7) and (8) are positive and
represent the physical process, where $\nu$ nucleons undergo inelastic
scattering, while the remaining $A-\nu$ nucleons provide absorption. In the
literature, one uses the term “wounded nucleons” Bialas:1976ed  and the
notation $\nu=N_{\rm coll}$.

The distribution $P_{\gamma}(\sigma)$ gives the probability density for
hadronic fluctuations of the real photon to interact with nucleons with the
cross section $\sigma$ Alvioli:2016gfo ; Frankfurt:2022jns . While the shape
of $P_{\gamma}(\sigma)$ cannot be calculated from the first principles, one
can reliably model it using constraints on its first moments and the
small-$\sigma$ and large-$\sigma$ limits. Indeed, the total photon-proton
cross section $\sigma_{\gamma p}(W)$ and the cross section of photon
diffractive dissociation on the proton $d\sigma_{\gamma p\to Xp}(W,t=0)/dt$
constrain the first two moments of $P_{\gamma}(\sigma)$ as follows

$\displaystyle\sigma_{\gamma p}(W)$ $\displaystyle=$ $\displaystyle\int
d\sigma P_{\gamma}(\sigma)\sigma\,,$ $\displaystyle\frac{d\sigma_{\gamma
p}(W,t=0)}{dt}$ $\displaystyle=$ $\displaystyle\frac{1}{16\pi}\int d\sigma
P_{\gamma}(\sigma)\sigma^{2}\,,$ (9)

where $W$ is the invariant photon-nucleon center-of-mass-energy. Further, in
the small-$\sigma$ limit, one can express $P_{\gamma}(\sigma)$ in terms of the
quark-antiquark component of the photon light-cone wave function and the color
dipole cross section, which results in $P_{\gamma}(\sigma)\propto 1/\sigma$.
In the opposite limit of large $\sigma$, the photon behaves as a superposition
of the $\rho$, $\omega$ and $\phi$ vector mesons in the spirit of the vector
meson dominance model and, hence, $P_{\gamma}(\sigma)$ can be modeled using
hadronic (cross section) fluctuations in $\rho$ mesons, which in turn are
related to those for pions. Finally, the small-$\sigma$ and large-$\sigma$
regimes can be smoothly interpolated. Note that this matching is achieved
best, when the light quark masses $m_{q}$ are taken to be those of the
constituent quarks, $m_{q}\sim 300$ MeV. For details, see Alvioli:2016gfo ;
Frankfurt:2022jns .

The left panel of Fig. 2 presents the resulting distribution
$P_{\gamma}(\sigma)$ as a function of $\sigma$ at $W=100$ GeV. Its shape and
normalization are constrained by the procedure outlined above, the
parametrization of $\sigma_{\gamma p}(W)$ ParticleDataGroup:2014cgo , and the
data on $d\sigma_{\gamma p\to Xp}(W,t=0)/dt$ Chapin:1985mf . Since the $W$
dependence of $P_{\gamma}(\sigma)$ is weak, the presented distribution is
applicable to a wide range of energies probed in heavy-ion UPCs at the LHC.
Note that the distribution $P_{\gamma}(\sigma)$ parametrizes the so-called
resolved photon contribution to photon-induced scattering and does contain the
direct photon contribution.

Using Eqs. (7) and (8), one can readily define the $P(\nu)$ probability
distribution for the number of wounded nucleons $\nu$ in inelastic photon-
nucleus scattering as follows Alvioli:2016gfo ,

$P(\nu)=\frac{\sigma_{\nu}}{\sum_{\nu=1}^{A}\sigma_{\nu}}\,,$ (10)

where $\sigma_{\nu}$ is given by Eq. (8). To calculate it, we use the Monte
Carlo generator for nucleon configurations in complex nuclei Alvioli:2009ab ,
which also includes nucleon-nucleon correlations in the nucleus wave function
Alvioli:2008rw ; Alvioli:2007zz , and the Gribov-Glauber model for photon-
nucleus scattering, where the hadronic structure of the photon is described by
the distribution $P_{\gamma}(\sigma)$ Alvioli:2016gfo . The resulting
distribution $P(\nu)$ as a function of $\nu$ for lead (Pb) is shown by the
curve labeled “Color Fluctuations” in the right panel of Fig. 2.

Figure 2: (a) The $P_{\gamma}(\sigma)$ probability density for hadronic
fluctuations of the real photon to interact with nucleons with the cross
section $\sigma$ at $W=100$ GeV. (b) The $P(\nu)$ distribution over the number
of wounded nucleons (inelastic interactions) $\nu$ in inelastic photon-nucleus
(Pb) scattering. The three curves correspond to the three models for
$P_{\gamma}(\sigma)$ parametrizing hadronic fluctuations of the real photon,
see text for details. The insert emphasizes the region of large $\nu$. The
figures are adopted from Alvioli:2016gfo .

Our modeling of small-$\sigma$ hadronic fluctuations of the photon is based on
the quark-antiquark component of the photon wave function and corresponds to a
weak nuclear shadowing in disagreement with the observed strong gluon nuclear
shadowing, see Eqs. (1) and (2). To take the latter into account, we model the
nuclear suppression of the dipoles with $\sigma\leq\sigma_{0}=20$ mb by the
factor of $R^{g}_{A}$, which leads to the modified expression for the
distribution $\tilde{P}_{\gamma}(\sigma)$,

$\tilde{P}_{\gamma}(\sigma)=\left[R^{g}_{A}(x,Q^{2}_{\rm
eff})\theta(\sigma_{0}-\sigma)+\theta(\sigma-\sigma_{0})\right]P_{\gamma}(\sigma)\,,$
(11)

where $x=Q^{2}_{\rm eff}/W^{2}$ and $P_{\gamma}(\sigma)$ is shown in Fig. 2
(left panel). The distribution $P(\nu)$ corresponding to $\sigma_{\nu}$
calculated using $\tilde{P}_{\gamma}(\sigma)$ is given by the curve
“Generalized CF” in the right panel of Fig. 2.

Finally, to test the importance of cross section (color) fluctuations in the
real photon, we calculated $P(\nu)$ neglecting these fluctuations and using

$P_{\gamma}(\sigma)=\delta(\sigma-25\,{\rm mb})$ (12)

in Eq. (8). The result is given by the curve “Glauber” in Fig. 2. One can see
from this figure that color fluctuations significantly increase the
distribution $P(\nu)$ at small and especially large $\sigma$; the latter is
emphasized in the insert.

In the total inelastic photon-nucleus cross section, the AGK cancellations
manifest themselves as the observation that the average number of wounded
nucleons $\langle N_{\rm coll}\rangle$ is inversely proportional to the
nuclear shadowing factor. Generalizing the result of Bertocchi:1976bq  for
hadron-nucleus scattering to the case of photon-induced scattering, one
obtains

$\langle N_{\rm
coll}\rangle\equiv\sum_{\nu=1}^{A}P(\nu)\nu=\frac{\sum_{\nu=1}^{A}\nu\sigma_{\nu}}{\sum_{\nu=1}^{A}\sigma_{\nu}}=\frac{A\sigma_{\rm
inel}^{\gamma N}}{\sigma_{\rm inel}^{\gamma A}}\,,$ (13)

where $\sigma_{\rm inel}^{\gamma N}$ is the photon-nucleon inelastic cross
section. Considering a particular hard process in inelastic photon-nucleus
scattering that probes the nuclear gluon distribution, e.g., inclusive
charmonium (bottomonium) production $\gamma+A\to J/\psi(\Upsilon)+X$ or
inclusive heavy-quark dijet production $\gamma+A\to Q{\bar{Q}}+X$, one obtains
using Eq. (13):

$\langle N_{\rm coll}\rangle\approx\frac{1}{R^{g}_{\rm Pb}(x,Q^{2})}\lesssim
2\,.$ (14)

In this estimate we used the following considerations. In general, the effect
of nuclear shadowing for the nuclear inelastic cross section is somewhat
larger than that for the total cross section. However, the theoretical
uncertainties of the leading twist approximation Frankfurt:2011cs  largely
mask the differences and make the shadowing effects approximately equal
(within uncertainties) in the two cases. Finally, in the last step, we used
the results of Eq. (2).

Figure 3 shows predictions of the leading twist approximation for the average
number of wounded nucleons $\langle N_{\rm coll}\rangle=1/R^{g}_{\rm
Pb}(x,Q^{2})$ in inelastic photon-nucleus (Pb) scattering as a function of $x$
at $Q^{2}=3$, 20, and 1000 GeV2. These values of $Q^{2}$ correspond to
photoproduction $J/\psi$, $\Upsilon$, and high-$p_{T}$ dijets, respectively.
One can see from the figure that in the discussed kinematics, the average
number of wounded nucleons is modest and the series in Eq. (13) converges
rather rapidly. In particular, we have checked that it is saturated by first
six terms with a 5% precision. Note, however, that the convergence slows down
in the limit of small $x$.

Figure 3: The LTA predictions for the average number of wounded nucleons
$\langle N_{\rm coll}\rangle$ in inelastic photon-nucleus (Pb) scattering as a
function of $x$ at $Q^{2}=3$, 20, and 1000 GeV2.

Measurements of $\langle N_{\rm coll}\rangle$ present a new method to study
nuclear shadowing in inelastic photon-nucleus scattering. Unlike the
observables used so far, the constraint of Eq. (14) indicates that one can
perform a “Pomeron surgery” of nuclear shadowing by cutting a small number of
Pomeron exchanges controlling the number of inelastic interactions with target
nucleons. As a result, it gives an opportunity for an experimental
determination of a small number of parameters quantifying nuclear shadowing,
which leads to a systematic improvement of its theoretical description. We
elaborate on it in the following section.

##  III The distributions over the number of wounded nucleon and forward
neutrons from nuclear breakup

While the average number of inelastic interactions $\langle N_{\rm
coll}\rangle$ encodes information on the energy and scale dependence of
nuclear shadowing, its dependence on the impact parameter of the collision is
averaged over. To obtain a more microscopic description of nuclear shadowing
and restore its impact parameter dependence, one needs to determine not only
$\langle N_{\rm coll}\rangle$, but also the entire distribution over the
number of wounded nucleons. This can be done using experimental data on the
neutron emission resulting from nucleus fragmentation in a given UPC process,
e.g., in inclusive quarkonium or dijet photoproduction in heavy-ion UPCs with
an additional condition of $Xn$ neutrons in the zero degree calorimeter (ZDC)
on the nuclear target side ATLAS:2017kwa ; ATLAS:2022cbd .

Very little is known about the dynamics of neutron emission in high energy
scattering off heavy nuclei. The ALICE collaboration measured the distribution
over the number of collisions $\langle N_{\rm coll}\rangle$ in proton-nucleus
scattering as determined by the energy release ($E_{T}$) at central rapidities
and the neutron multiplicity as a function of $E_{T}$. It was observed that
$\langle N_{\rm coll}(E_{T})\rangle$ is linearly proportional to number of
evaporation neutrons $\langle M_{n}(E_{T})\rangle$ for the same $E_{T}$ bins
at least up to $\langle N_{\rm coll}\rangle\sim 10$, for details, see
ALICE:2021poe . Note that in our case, $\langle N_{\rm coll}\rangle$ is much
lower, see Eq. (14).

Another important observation made by E665 experiment at Fermilab is that in
muon-nucleus deep inelastic scattering (DIS) in coincidence with detection of
slow neutrons, $\mu^{-}+A\to n+X$, the average neutron multiplicity $\langle
M_{n}\rangle$ for the lead target is E665:1995utr

$\langle M_{n}\rangle\approx 5\,.$ (15)

This result has been understood in the framework of cascade models of nuclear
DIS Strikman:1998cc ; Larionov:2018igy , where soft neutrons are produced
either directly in DIS on a bound nucleon or through statistical decay (de-
excitation) of the excited residual nucleus, leading to neutron
evaporation111The geometry of the heating is very different from the case of
$AA$ collisions, where in each collision a large postion of nucleons in each
nucleus interacts and de-excitation of the spectators only occurs close to the
interaction surface Alvioli:2010yk .. The latter mechanism depends strongly on
the hadron formation time: to describe the energy spectrum of emitted
neutrons, one has to assume that only nucleons and pions with momenta $\leq 1$
GeV/c could be involved in final state interactions. A similar conclusion was
reached by M. Baker (private communication, 2022) using the BeAGLE Monte Carlo
generator Chang:2022hkt .

This suggests the following space-time picture of forward neutron production
in high energy photon-nucleus scattering. Well before the target, the incoming
photon fluctuates into long-lived hadronic components, which pass through the
nucleus and interact inelastically with several nucleons. This leads to the
creation of holes in the nucleus (particle-hole excitations in terminology of
a nuclear shell model), which de-excite and cool the nucleus by evaporating
neutrons. It also produces a number of soft particles with the momenta less
than 1 GeV/c, which in turn generate more neutrons.

The nucleon fragmentation weakly depends on the incident energy due to Feynman
scaling and, hence, it is not significantly affected by the energy
conservation constraint, which is important in the case, when one uses hadron
multiplicities at central rapidities ATLAS:2015hkr . In this case, the energy
transferred to the rest of the nucleus, which heats the residual nuclear
system, is proportional to $\langle N_{\rm coll}\rangle$. Since the Fermilab
data E665:1995utr  corresponds to the average momentum fraction $\langle
x\rangle=0.015$, where the nuclear shadowing effect is small, one finds that
$\langle N_{\rm coll}\rangle$ $\approx$ 1, see Eq. (13). Thus, every inelastic
photon-nucleon interaction results on average in 5 forward neutrons.

Figure 4: (a) The probability distribution (frequency) for Pb to produce $N$
forward neutrons assuming their Poisson distribution with $\nu=\langle N_{\rm
coll}\rangle$=1, 2, and 3 independent emissions (wounded nucleons) and the
$\langle M_{n}\rangle=5$ average multiplicity for a single inelastic
scattering. (b) The probability distribution of forward neutron emission
$P_{\rm comb}(N)$ (17) in the model combining the distribution over the number
of wounded nucleons, which is calculated using three models for hadronic
fluctuations in the photon (see Fig. 2), with the Poisson distribution of
produced neutrons.

To test this hypothesis, we perform two numerical studies. First, we consider
a simple model, which assumes that the probability density of neutron emission
is given by the Poisson distribution and that each hole created in the target
nucleus generates independently on average $\langle M_{n}\rangle$ neutrons.
Therefore, the neutron probability distribution for $\nu=\langle N_{\rm
coll}\rangle$ wounded nucleons is

$P_{\rm Poisson}(N;\lambda=\nu\langle M_{n}\rangle)=\frac{(\nu\langle
M_{n}\rangle)^{N}e^{-\nu\langle M_{n}\rangle}}{N!}\,,$ (16)

where $N$ is the number of produced neutrons (neutron multiplicity). The
resulting probability density (frequency) for Pb as a function of $N$ is shown
in the left panel of Fig. 4. In this estimate, we used $\langle
M_{n}\rangle=5$, see Eq. (15) and Ref. E665:1995utr , and $\nu=\langle N_{\rm
coll}\rangle$ = 1, 2 or 3 independent neutron emissions (wounded nucleons).
One can see from the figure that since the distributions for different $\nu$
are peaked at different values of $N$ and do not significantly overlap, the
measurement of the forward neutron multiplicity can be used to reliably
separate contributions of different numbers of wounded nucleons $\langle
N_{\rm coll}\rangle$. Our analysis also shows that this separation becomes
even cleaner with an increase of $\langle M_{n}\rangle$ (not shown here).

In the second study, we combine the distribution over the number of wounded
nucleons $P(\nu)$ that we discussed in Sec. II, see Eq. (10) and Fig. 2 (right
panel), with the Poisson distribution of produced neutrons, see Eq. (16) and
Fig. 4 (left panel). The resulting probability distribution of forward
neutrons is given by the following expression,

$P_{\rm comb}(N)=\sum_{\nu=1}^{A}P(\nu)P_{\rm Poisson}(N;\nu\langle
M_{n}\rangle)\,.$ (17)

The right panel of Fig. 4 presents $P_{\rm comb}(N)$ as a function of forward
neutrons $N$ for $\langle M_{n}\rangle=5$. The three curves correspond to the
three models for the hadronic (color) fluctuations of the real photon, see the
right panel of Fig. 2. One can see from the figure that cross section
fluctuations in the real photon noticeably affect the shape of the neutron
distribution: its maximum is shifted toward smaller $N$ compared to the
“Glauber” result and is also somewhat suppressed by the leading twist
shadowing in the “Generalized CF” case.

Our numerical studies suggest that one can examine details of the theoretical
description of small-$x$ nuclear shadowing in QCD using the distribution of
forward neutrons from nuclear breakup emitted in a given hard UPC process,
which is directly correlated with the number of inelastic interactions
(wounded nucleons). By studying the dispersion of this distribution, one can
single out the individual contributions of $\nu$ = 1, 2, 3, $\dots$ wounded
nucleons. In particular, using Eq. (8) in the approximation that the series in
Eq. (13) is saturated by first few terms, which corresponds to the limit of
small-to-modest nuclear shadowing, one obtains Alvioli:2016gfo ;
Alvioli:2014sba

$\sigma_{2}(x,Q^{2})\equiv\frac{\langle\sigma_{\rm
inel}^{2}\rangle}{\langle\sigma_{\rm inel}\rangle}=\frac{\langle N_{\rm
coll}-1\rangle}{\frac{(A-1)}{2}\int d^{2}\vec{b}\,T_{A}^{2}(\vec{b})}\,,$ (18)

and

$\sigma_{3}(x,Q^{2})\equiv\frac{\langle\sigma_{\rm
inel}^{3}\rangle}{\langle\sigma_{\rm
inel}^{2}\rangle}=\frac{\langle\sigma_{\rm
inel}^{3}\rangle}{\langle\sigma_{\rm inel}\rangle}\frac{\langle\sigma_{\rm
inel}^{2}\rangle}{\langle\sigma_{\rm inel}^{2}\rangle}=\frac{(\langle N_{\rm
coll}-2\rangle)(\langle N_{\rm coll}-1\rangle)}{\langle N_{\rm
coll}-1\rangle}\frac{\frac{(A-1)}{2}\int
d^{2}\vec{b}\,T_{A}^{2}(\vec{b})}{\frac{(A-2)(A-1)}{2}\int
d^{2}\vec{b}\,T_{A}^{3}(\vec{b})}\,,$ (19)

where $\langle\sigma_{\rm inel}^{n}\rangle=\int d\sigma
P_{\gamma}(\sigma)\sigma_{\rm inel}^{n}$. The cross sections
$\sigma_{2}(x,Q^{2})$ and $\sigma_{3}(x,Q^{2})$ are essential ingredients of
the LTA approach Frankfurt:2011cs . While $\sigma_{2}(x,Q^{2})$ is determined
using the HERA data on inclusive diffraction in lepton-proton DIS,
$\sigma_{3}(x,Q^{2})$ is model-dependent and leads to theoretical
uncertainties of this approach. Thus, an independent determination
$\sigma_{2}(x,Q^{2})$ and $\sigma_{3}(x,Q^{2})$ using UPCs with forward
neutrons will supply new constraints on these quantities.

Note that to reach a high accuracy in such an analysis, one needs to calibrate
the theoretical description against the kinematics, where only one target
nucleon is struck, e.g., $\gamma+A\to{\rm 2\ jets}+X$ or quasi-elastic
$J/\psi$ production for $x_{A}$ $\geq$ 0.01, where the effect of nuclear
shadowing is small, and $\langle N_{\rm coll}\rangle$ $\approx$ 1.

##  IV Conclusions and outlook

In this paper, we advertise measurements of forward neutrons from nuclear
breakup in inclusive high energy photon-nucleus scattering in heavy-ion UPCs,
e.g., charmonium (bottomonium) production $\gamma+A\to J/\psi(\Upsilon)+X$ or
heavy-quark dijet production $\gamma+A\to Q{\bar{Q}}+X$, as a novel way to
study the QCD dynamics at small $x$. The key quantity is the number of
inelastic photon-nucleon interactions (the number of wounded nucleons): its
average value $\langle N_{\rm coll}\rangle$ is proportional to inverse of the
gluon nuclear shadowing and its first moments are sensitive to details of the
leading twist mechanism of nuclear shadowing. Our numerical analysis has
demonstrated that the number of forward neutrons from nuclear breakup detected
in the ZDC on the nuclear target side is rather unambiguously proportional to
the number of wounded nucleons, which provides a practical opportunity for
novel studies of nuclear shadowing.

On top of providing stringent tests of the dynamics of leading twist shadowing
of gluon PDFs, it would be possible to explore effects related to proximity to
the black disk limit of the strong interaction. For example, one can study
fragmentation of leading hadrons in $\gamma A$ scattering and look for
suppression of their multiplicity as a function of Feynman $x_{F}$ and $W$ as
well as for broadening of their transverse momentum distribution
Frankfurt:2001nt . These effects should be more pronounced for central
collisions characterized by an enhanced activity in the ZDC. It should be
possible to construct from the data an analog of the $R_{CP}$ ratio, which
would probe the density dependence of fragmentation. It would also be useful
to construct similar quantities for low $p_{T}$ charm production.

Another interesting application is for multiparton interactions in proton-
nucleus ($pA$) scattering. It was argued in Alvioli:2019kcy  that the single
and double scattering can be separated using their dependence on the impact
parameter: the former is proportional to $A$, while the latter $\propto
A^{4/3}$. However, since both hard interactions are typically detected in a
limited range of rapidities $|y|\leq 3-4$, centrality is difficult to
determine from the transverse energy $E_{T}$ signal because multiparton
interactions also contribute to $E_{T}$. The use of forward neutrons in ZDCs
would alleviate this problem.

One should point out that the neutrons detected in ZDCs can be a promising
complementary way to determine centrality of various photon-nucleus and
proton-nucleus inelastic collisions expanding the use of ZDCs beyond their
current use in vector meson diffractive production and for determining of
centrality of the heavy-ion collisions. The main advantage of using forward
neutrons rather than the transverse energy $E_{T}$ for the determination of
centrality is a much larger distance in rapidity between the rapidity of the
hard process and that of the process used for determination of the centrality.

The methods presented in this paper can be readily generalized to the case of
virtual photons and allow one to predict the distribution over the number of
forward neutrons in inelastic photon-nucleus scattering at the EIC.

###### Acknowledgements.

The research of V.G. was funded by the Academy of Finland project 330448, the
Center of Excellence in Quark Matter of the Academy of Finland (projects
346325 and 346326), and the European Research Council project
ERC-2018-ADG-835105 YoctoLHC. The research of M.S. was supported by the US
Department of Energy Office of Science, Office of Nuclear Physics under Award
No. DE- FG02-93ER40771.

## References

  * (1) Y. L. Dokshitzer, D. Diakonov and S. I. Troian, Phys. Rept. 58, 269-395 (1980)
  * (2) L. V. Gribov, E. M. Levin and M. G. Ryskin, Phys. Rept. 100, 1-150 (1983)
  * (3) L. Frankfurt, V. Guzey, M. McDermott and M. Strikman, Phys. Rev. Lett. 87, 192301 (2001) [arXiv:hep-ph/0104154 [hep-ph]].
  * (4) A. Accardi, J. L. Albacete, M. Anselmino, N. Armesto, E. C. Aschenauer, A. Bacchetta, D. Boer, W. K. Brooks, T. Burton and N. B. Chang, et al. Eur. Phys. J. A 52, no.9, 268 (2016) [arXiv:1212.1701 [nucl-ex]].
  * (5) R. Abdul Khalek, A. Accardi, J. Adam, D. Adamiak, W. Akers, M. Albaladejo, A. Al-bataineh, M. G. Alexeev, F. Ameli and P. Antonioli, et al. Nucl. Phys. A 1026, 122447 (2022) [arXiv:2103.05419 [physics.ins-det]].
  * (6) M. Strikman, R. Vogt and S. N. White, Phys. Rev. Lett. 96, 082001 (2006) [arXiv:hep-ph/0508296 [hep-ph]].
  * (7) A. J. Baltz, G. Baur, D. d’Enterria, L. Frankfurt, F. Gelis, V. Guzey, K. Hencken, Y. Kharlov, M. Klasen and S. R. Klein, et al. Phys. Rept. 458, 1-171 (2008)
  * (8) E. Abbas et al. [ALICE], Eur. Phys. J. C 73, no.11, 2617 (2013) [arXiv:1305.1467 [nucl-ex]].
  * (9) B. Abelev et al. [ALICE], Phys. Lett. B 718, 1273-1283 (2013) [arXiv:1209.3715 [nucl-ex]].
  * (10) V. Khachatryan et al. [CMS], Phys. Lett. B 772, 489-511 (2017) [arXiv:1605.06966 [nucl-ex]].
  * (11) S. Acharya et al. [ALICE], Eur. Phys. J. C 81, no.8, 712 (2021) [arXiv:2101.04577 [nucl-ex]].
  * (12) S. Acharya et al. [ALICE], Phys. Lett. B 798, 134926 (2019) [arXiv:1904.06272 [nucl-ex]].
  * (13) R. Aaij et al. [LHCb], JHEP 07, 117 (2022) [arXiv:2107.03223 [hep-ex]].
  * (14) R. Aaij et al. [LHCb], JHEP 06, 146 (2023) [arXiv:2206.08221 [hep-ex]].
  * (15) A. Tumasyan et al. [CMS], Phys. Rev. Lett. 131, no.26, 262301 (2023) [arXiv:2303.16984 [nucl-ex]].
  * (16) S. Acharya et al. [ALICE], JHEP 10, 119 (2023) [arXiv:2305.19060 [nucl-ex]].
  * (17) [STAR], [arXiv:2311.13632 [nucl-ex]].
  * (18) L. Frankfurt, V. Guzey and M. Strikman, Phys. Rept. 512, 255-393 (2012) [arXiv:1106.2091 [hep-ph]].
  * (19) V. Guzey, E. Kryshen, M. Strikman and M. Zhalov, Phys. Lett. B 726, 290-295 (2013) [arXiv:1305.1724 [hep-ph]].
  * (20) V. Guzey and M. Zhalov, JHEP 10, 207 (2013) [arXiv:1307.4526 [hep-ph]].
  * (21) K. J. Eskola, C. A. Flett, V. Guzey, T. Löytäinen and H. Paukkunen, Phys. Rev. C 106, no.3, 035202 (2022) [arXiv:2203.11613 [hep-ph]].
  * (22) K. J. Eskola, C. A. Flett, V. Guzey, T. Löytäinen and H. Paukkunen, Phys. Rev. C 107, no.4, 044912 (2023) [arXiv:2210.16048 [hep-ph]].
  * (23) S. P. Jones, A. D. Martin, M. G. Ryskin and T. Teubner, J. Phys. G 43, no.3, 035002 (2016) [arXiv:1507.06942 [hep-ph]].
  * (24) S. P. Jones, A. D. Martin, M. G. Ryskin and T. Teubner, Eur. Phys. J. C 76, no.11, 633 (2016) [arXiv:1610.02272 [hep-ph]].
  * (25) V. Guzey and M. Klasen, Phys. Rev. C 99, no.6, 065202 (2019) [arXiv:1811.10236 [hep-ph]].
  * (26) V. Guzey and M. Klasen, JHEP 04, 158 (2016) [arXiv:1603.06055 [hep-ph]].
  * (27) [ATLAS], “Photo-nuclear dijet production in ultra-peripheral Pb+Pb collisions,” ATLAS-CONF-2017-011.
  * (28) [ATLAS], “Photo-nuclear jet production in ultra-peripheral Pb+Pb collisions at $\sqrt{s}_{\text{NN}}=5.02$ TeV with the ATLAS detector,” ATLAS-CONF-2022-021.
  * (29) B. Pire, L. Szymanowski and J. Wagner, Phys. Rev. D 79, 014010 (2009) [arXiv:0811.0321 [hep-ph]].
  * (30) W. Schafer, G. Slipek and A. Szczurek, Phys. Lett. B 688, 185-191 (2010) [arXiv:1003.0610 [hep-ph]].
  * (31) G. M. Peccini, L. S. Moriggi and M. V. T. Machado, Phys. Rev. D 103, no.5, 054009 (2021) [arXiv:2101.08338 [hep-ph]].
  * (32) S. R. Klein, J. Nystrand and R. Vogt, Phys. Rev. C 66, 044906 (2002) [arXiv:hep-ph/0206220 [hep-ph]].
  * (33) V. P. Goncalves, M. V. T. Machado and A. R. Meneses, Phys. Rev. D 80, 034021 (2009) [arXiv:0905.2067 [hep-ph]].
  * (34) V. P. Gonçalves, G. Sampaio dos Santos and C. R. Sena, Nucl. Phys. A 976, 33-45 (2018) [arXiv:1711.04497 [hep-ph]].
  * (35) V. A. Abramovsky, V. N. Gribov and O. V. Kancheli, Yad. Fiz. 18, 595-616 (1973), Sov. J. Nucl. Phys. 18, 308-317 (1974).
  * (36) D. Treleani, Int. J. Mod. Phys. A 11, 613-654 (1996)
  * (37) J. Jalilian-Marian and Y. V. Kovchegov, Phys. Rev. D 70, 114017 (2004) [erratum: Phys. Rev. D 71, 079901 (2005)] [arXiv:hep-ph/0405266 [hep-ph]].
  * (38) F. Gelis and R. Venugopalan, Nucl. Phys. A 776, 135-171 (2006) [arXiv:hep-ph/0601209 [hep-ph]].
  * (39) A. Kovner and M. Lublinsky, JHEP 11, 083 (2006) [arXiv:hep-ph/0609227 [hep-ph]].
  * (40) N. N. Nikolaev and W. Schafer, Phys. Rev. D 74, 074021 (2006) [arXiv:hep-ph/0607307 [hep-ph]].
  * (41) L. L. Frankfurt and M. I. Strikman, Phys. Lett. B 382, 6-12 (1996)
  * (42) L. Frankfurt and M. Strikman, Eur. Phys. J. A 5, 293-306 (1999) [arXiv:hep-ph/9812322 [hep-ph]].
  * (43) V. N. Gribov, Sov. Phys. JETP 29, 483-487 (1969) [Zh. Eksp. Teor. Fiz. 56, 892-901 (1969)]
  * (44) V. N. Gribov, Sov. Phys. JETP 30, 709-717 (1970) [Zh. Eksp. Teor. Fiz. 57, 1306-1323 (1969)]
  * (45) G. Piller and W. Weise, Phys. Rept. 330, 1-94 (2000) [arXiv:hep-ph/9908230 [hep-ph]].
  * (46) A. Adeluyi and G. Fai, Phys. Rev. C 74, 054904 (2006) [arXiv:hep-ph/0610214 [hep-ph]].
  * (47) A. Adeluyi and T. Nguyen, Phys. Rev. C 75, 054911 (2007)
  * (48) J. C. Collins, Phys. Rev. D 57, 3051-3056 (1998) [erratum: Phys. Rev. D 61, 019902 (2000)] [arXiv:hep-ph/9709499 [hep-ph]].
  * (49) L. Frankfurt, V. Guzey and M. Strikman, Phys. Rev. D 71, 054001 (2005) [arXiv:hep-ph/0303022 [hep-ph]].
  * (50) L. Bertocchi and D. Treleani, J. Phys. G 3, 147 (1977)
  * (51) M. Alvioli, L. Frankfurt, V. Guzey, M. Strikman and M. Zhalov, Phys. Lett. B 767, 450-457 (2017) [arXiv:1605.06606 [hep-ph]].
  * (52) V. Franco and R. J. Glauber, Phys. Rev. 142, 1195-1214 (1966)
  * (53) A. Bialas, M. Bleszynski and W. Czyz, Nucl. Phys. B 111, 461-476 (1976)
  * (54) L. Frankfurt, V. Guzey, A. Stasto and M. Strikman, Rept. Prog. Phys. 85, no.12, 126301 (2022) [arXiv:2203.12289 [hep-ph]].
  * (55) K. A. Olive et al. [Particle Data Group], Chin. Phys. C 38, 090001 (2014)
  * (56) T. J. Chapin, R. L. Cool, K. A. Goulianos, K. A. Jenkins, J. P. Silverman, G. R. Snow, H. Sticker, S. N. White and Y. H. Chou, Phys. Rev. D 31, 17-30 (1985)
  * (57) M. Alvioli, H. J. Drescher and M. Strikman, Phys. Lett. B 680, 225-230 (2009) [arXiv:0905.2670 [nucl-th]].
  * (58) M. Alvioli, C. Ciofi degli Atti, I. Marchino, V. Palli and H. Morita, Phys. Rev. C 78, 031601 (2008) [arXiv:0807.0873 [nucl-th]].
  * (59) M. Alvioli, C. Ciofi degli Atti and H. Morita, Phys. Rev. Lett. 100, 162503 (2008)
  * (60) S. Acharya et al. [ALICE], JHEP 08, 086 (2022) [arXiv:2107.10757 [nucl-ex]].
  * (61) M. R. Adams et al. [E665], Phys. Rev. Lett. 74, 5198-5201 (1995) [erratum: Phys. Rev. Lett. 80, 2020-2021 (1998)]
  * (62) M. Strikman, M. G. Tverskoii and M. B. Zhalov, Phys. Lett. B 459, 37-42 (1999) [arXiv:nucl-th/9806099 [nucl-th]].
  * (63) A. B. Larionov and M. Strikman, Phys. Rev. C 101, no.1, 014617 (2020) [arXiv:1812.08231 [hep-ph]].
  * (64) M. Alvioli and M. Strikman, Phys. Rev. C 83, 044905 (2011) [arXiv:1008.2328 [nucl-th]].
  * (65) W. Chang, E. C. Aschenauer, M. D. Baker, A. Jentsch, J. H. Lee, Z. Tu, Z. Yin and L. Zheng, Phys. Rev. D 106, no.1, 012007 (2022) [arXiv:2204.11998 [physics.comp-ph]].
  * (66) G. Aad et al. [ATLAS], Eur. Phys. J. C 76, no.4, 199 (2016) [arXiv:1508.00848 [hep-ex]].
  * (67) M. Alvioli, L. Frankfurt, V. Guzey and M. Strikman, Phys. Rev. C 90, 034914 (2014) [arXiv:1402.2868 [hep-ph]].
  * (68) M. Alvioli, M. Azarkin, B. Blok and M. Strikman, Eur. Phys. J. C 79, no.6, 482 (2019) [arXiv:1901.11266 [hep-ph]].