# Deblurring for Nuclei: 3D Characteristics of Heavy-Ion Collisions

Pawel Danielewicz<EMAIL_ADDRESS>Facility for Rare Isotope Beams and
Department of Physics and Astronomy, Michigan State University, East Lansing,
Michigan 48824, USA  Mizuki Kurata-Nishimura  RIKEN Nishina Center, Hirosawa
2-1, Wako, Saitama 351-0198, Japan

###### Abstract

Observables from nuclear and high-energy experiments can be degraded by
detector performance and/or methodology in extracting the observables, such as
of the final-state characteristics of heavy-ion collisions in relation to a
coarsely estimated reaction-plane direction. We propose the use of deblurring
methods, such as in optics, to correct for observable degradation. Our main
focus is the restoration of triple-differential particle distributions in
heavy-ion collisions. We demonstrate that these could be extracted from
collision measurements following the Richardson-Lucy deblurring method from
optics. We illustrate basic features of the restoration methodology in a
schematic model assuming either ideal or more realistic particle detection.
The inferred three-dimensional (3D) distributions for collisions may easier to
interpret in terms of collision dynamics and sought properties of bulk matter
than the currently employed Fourier coefficients, that combine information
from different azimuthal angles relative to the reaction plane.

reaction plane; Richardson-Lucy algorithm; resolution

##  I Introduction

Observables from nuclear and high-energy experiments can be degraded by
detector performance and/or methodology in extracting the observables, such as
tied to limited control over the initial state of a reaction. With the latter
respect, the three-dimensional (3D) particle distributions relative to the
reaction plane in heavy-ion collisions are of much interest, but the reaction
plane direction can be only coarsely assessed, so effectively only blurred
distributions could be directly determined. The characteristics at a fixed
impact parameter are of interest for the heavy-ion collisions, but again the
impact parameter can be assessed only coarsely. Direct reactions of fast rare
isotopes are of much interest, but fast rare isotopes are produced in flight
with their beam consequently spread out in space, angle and energy,
potentially erasing details in observables. There may be different strategies
to deal with such situations, such as to regress in the observables,
particularly to those that are less degraded by experimental limitations or
easier to correct. Else, one can resort to theory and filter the theoretical
predictions through the features of the experimental setup and methodology.
The problem there is that nuclear theory is usually approximate at one level
or another. With observables being degraded and being compared against
uncertain reference, important insights might get lost. As an alternative to
the practice so far, we propose the use of deblurring techniques to cope with
any degraded observables. The deblurring has a long record in optics, being
employed there to correct images for flaws of optical apparatus. With the
deblurring as an option, one may have a wider range of observables to choose
from and may let those observables to stand on their own, prior to any theory
comparisons, inspiring the latter rather than being resigned to them. Of our
primary interest is the inference of 3D characteristics of the final states of
heavy-ion collisions. However, different strategy details may be relevant for
other contexts in nuclear physics and potentially in high-energy physics.

In heavy-ion collisions, any directional aspects of observables have been
primarily tied to the incident beam direction, as in particular evidenced in
the use of rapidity along the beam and/or transverse momentum magnitude in
presenting data. However, the symmetry around the beam direction is broken in
the initial state by the relative displacement of the nuclear centers, by
impact parameter. The transverse displacement and the beam axis define the
reaction plane for the collision. The symmetry breaking has important
consequences for the collisions, such as in setting the preferred direction of
transverse momentum transfer in the system and in the emergence of participant
and spectator products in more energetic collisions. At lower energies, heavy
projectile and target residues emerge from a collision largely aligned with
the reaction plane and other characteristics of the final state tend to be
related to these. At higher energies, much could be learned on reaction of
strongly interacting matter to compression, from product emission in
correlation to the reaction plane. However, at the same time, estimating the
reaction plane becomes a challenge at those higher energies since the system
largely vaporizes. Our considerations here will largely refer to the latter
situation. At still higher energies, Lorentz dilation freezes fluctuations
present the initial state, that can then compete with the breaking of the
symmetry around the beam axis by the reaction plane. Even there our proposed
methodology may have applications, though it might not provide the most
interesting information one may seek in the context.

Importantly, the way in which the symmetry breaking in the initial state gets
reflected in transverse anisotropies in particle emission, that may represent
transverse collective motion and shadowing, can provide significant
information about the colliding system, including temporal progress of the
collision and bulk properties of the matter such as equation of state and
viscosity. For that reason, strategies not only to identify the reaction plane
direction, but also quantify the azimuthal anisotropies have been intensely
pursued over decades, starting with telling whether the transverse collective
motion is there at all [1, 2], through estimation of the reaction plane
direction and its uncertainty [3, 4], to detailed quantification of
anisotropies in terms of Cartesian moments [3, 5] and eventually Fourier
coefficients [6, *demoulins_collective_1989, 8, 9, 10, 11]. Importance of
correcting for detector inefficiencies and circumventing correlations
unrelated to spatial anisotropies, that could distort inferences, was
recognized [5, 10, 11].

We will be connecting to those prior developments here at different level,
whether directly, when addressing heavy-ion collisions, or in inspiration,
behind the considerations in a simplified model. Of specific importance is
that correlation of particle emission pattern with the reaction plane can be
exploited for estimating the direction of the reaction plane [3]. For low
event statistics, transverse Cartesian moments of the distributions can be
evaluated. Because of the uncertainty in the reaction plane determination,
values of the moments, Cartesian or Fourier, need to be renormalized, on a
moment by moment basis, and this may be done in a self-consistent manner [3,
6, *demoulins_collective_1989, 8, 9]. However, for high event statistics, the
whole distribution associated with the reaction plane could be addressed with
the Fourier coefficients [12]. Taking another look, if such a distribution
were evaluated directly using estimated directions for the reaction plane in
collision events, the distribution would come out blurred, because of the
estimated direction of the plane fluctuates around the true direction from one
event to another. This brings up analogy to the situation in optics where a
photo is taken with a camera that shakes. For such situations and other, where
distortions of image occur due to understood distortions of light intensity,
deblurring techniques have been developed in the optical contexts [13]. We
propose to use such techniques in the context of central heavy-ion collisions
to arrive at three-dimensional (3D) distributions tied to the reaction plane.
The interpretations of the latter distributions may be easier than of Fourier
coefficients, especially those high, given that these coefficients tend to
quantify information incidentally lumped together, through the reference to
the same value of transverse momentum magnitude. We moreover envision the use
of such techniques to other situations in nuclear or high-energy physics where
observables are distorted due to limitations of methodology or detector
performance.

We start out by reviewing the Richardson-Lucy (RL) deblurring algorithm [14,
15], that is popular in optics, based on Bayesian considerations, and well
suited, in our opinion, for nuclear applications including heavy-ion
collisions. We then turn to a schematic one-dimensional (1D) model, where a
system moving at unknown velocity emits particles according to a single-
particle distribution that is forward-backward symmetric in the c.m. frame. A
number of particles is measured and the goal is to determine the velocity
distribution of particles relative to the center of mass, even though the
center of mass is not known on the case by case basis. The center of mass
velocity may be estimated from the emitted particles, but it straggles
relative to the true velocity. With this, the distribution relative to the
center of mass evaluated from the emitted particles alone is smeared out
compared to the true distribution. We demonstrate that the original
distribution may be restored from the simulated observations by combining the
central-limit theorem with the RL deblurring algorithm. Next we repeat the
model assuming that particles are measured with a detector for which the
detection window is comparable with the spread of the particle distribution.
In this case, the central-limit theorem cannot be easily employed and the
restoration problem becomes significantly nonlinear. We demonstrate that in
that more involved situation, the original velocity distribution may be still
restored from the measured distribution following a self-consistent RL
algorithm. We thereafter turn to simulated distributions for heavy-ion
collisions. We assume a typical situation recognized as a good representation
of the final-state emission for semicentral collisions at few hundred
MeV/nucl, with a local equilibrium combined with sideward, radial and elliptic
flows. With a moderate number of light charged particles registered, with mass
$A\leq 4$, it is possible to deblur their distributions associated with the
reaction plane, even when these distributions vary by an order of magnitude or
more in the transverse directions. We complement the outline of the strategy
for data analysis with results from transport theory, illustrating what kind
of information could be accessed in the deblurred distributions, that might
not be easily seen in the azimuthal moments for the distributions.

##  II Deblurring

The blurring problem may be stated in the form of the equation:

$n({\zeta})=\int\text{d}{\xi}\,P({\zeta}|{\xi})\,{\mathcal{N}}({\xi})\,.$ (1)

Here, ${\mathcal{N}}({\xi})$ is the distribution in coordinates ${\xi}$ that
faithfully characterizes the measured system and $P({\zeta}|{\xi})$ is the
probability density that the system is measured around $\zeta$ when it is
actually at $\bf\xi$. Finally, $n({\bf\zeta})$ is the distribution attributed
to the system in effect of the direct measurements. The goal of deblurring is
to determine ${\mathcal{N}}({\xi})$ from the measured $n({\zeta})$, given
knowledge about $P$. While the problem appears stated as linear, in practical
situations $P$ may have a dependence on ${\mathcal{N}}$, too.

In the reaction-plane problem, the blurring is due to the fact that the
azimuthal angle relative to the estimated reaction plane is not equal to the
angle relative to the true reaction plane. The deviation may vary from one
collision event to another and even particle to a particle, within the
strategy of estimating the reaction plane orientation [3, 6]. In detector
efficiency problems, $P$ might account for particle being misidentified or
missed or whole event disregarded consequently under a trigger. Depending on
the situation, the probability density integrates to 1 or not,
$\int_{Z}\text{d}\zeta\,P({\zeta}|{\xi})=P_{Z}(\xi)\leq 1$. For generality
further on, we will consider the possibility of weighing the domain $Z$ with a
function $W(\zeta)$, with the net weighted probability being then
$P_{WZ}(\xi)=\int_{Z}\text{d}\zeta\,W(\zeta)\,P({\zeta}|{\xi})$. The function
$W$ may be then used to effectively shrink $Z$ such as in disregarding
inferior quality events.

For inferring ${\mathcal{N}}$ directly from measured $n$, it might be tempting
to invert (1) directly, but that strategy is likely to amplify short-
wavelength noise always present in $n$. In optical contexts alternative
methods have been developed of which the Richardson-Lucy (RL) deconvolution
algorithm [14, 15], based on Bayesian considerations that we next lay out, is
particularly popular.

Let $Q({\xi}|{\zeta})$ denote complementary probability density to $P$, that
the system is at $\xi$ when it is measured at $\zeta$. Then the probability
that the system is within $\text{d}\xi$ while measured within $\text{d}\zeta$
can be expressed in two different ways:

$Q({\xi}|{\zeta})\,n({\zeta})\,\text{d}\zeta\,\text{d}\xi=P({\zeta}|{\xi})\,{\mathcal{N}}({\xi})\,\text{d}\xi\,\text{d}\zeta\,.$
(2)

This yields

$Q({\xi}|{\zeta})=\frac{P({\zeta}|{\xi})\,{\mathcal{N}}({\xi})}{\int\text{d}{\xi^{\prime}}\,P({\zeta}|{\xi^{\prime}})\,{\mathcal{N}}({\xi^{\prime}})}\,,$
(3)

and

${\mathcal{N}}({\xi})=\frac{\int\text{d}{\zeta}\,Q({\xi}|{\zeta})\,W(\zeta)\,n({\zeta})}{P_{WZ}(\xi)}\,.$
(4)

The RL method solves the last two equations iteratively:

$\displaystyle Q^{(r)}({\xi}|{\zeta})$ $\displaystyle=$
$\displaystyle\frac{P({\zeta}|{\xi})\,{\mathcal{N}}^{(r)}({\xi})}{\int\text{d}\xi^{\prime}\,P({\zeta}|{\xi^{\prime}})\,{\mathcal{N}}^{(r)}({\xi^{\prime}})}\,,$
(5) $\displaystyle{\mathcal{N}}^{(r+1)}({\xi})$ $\displaystyle=$
$\displaystyle\frac{\int\text{d}{\zeta}\,Q^{(r)}({\xi}|{\zeta})\,W(\zeta)\,n({\zeta})}{P_{WZ}(\xi)}\,,$
(6)

where $r$ is iteration step index. Any implicit dependence of $P$ on
${\mathcal{N}}$ may be handled in the equations iteratively. One notable
feature of the equations is that distributions that start as nonnegative stay
such during the iterations. Combination of the equations above yields

${\mathcal{N}}^{(r+1)}({\xi})=A^{(r)}(\xi)\cdot{\mathcal{N}}^{(r)}({\xi})\,,$
(7)

where

$A^{(r)}(\xi)=\frac{\int\text{d}\zeta\,\frac{n(\zeta)}{n^{(r)}(\zeta)}\,W(\zeta)\,P(\zeta|\xi)}{\int\text{d}\zeta^{\prime}\,W(\zeta^{\prime})\,P(\zeta^{\prime}|\xi)}\,,$
(8)

and

$n^{(r)}(\zeta)=\int\text{d}{\xi}\,P({\zeta}|{\xi})\,{\mathcal{N}}^{(r)}({\xi})\,.$
(9)

If the variable space is discretized, then the equations become

${\mathcal{N}}_{i}^{(r+1)}=A_{i}^{(r)}\cdot{\mathcal{N}}_{i}^{(r)}\,,$ (10)
$A_{i}^{(r)}=\sum_{j}\frac{n_{j}}{n_{j}^{(r)}}\,T_{ji}\,,$ (11)

with

$n_{j}^{(r)}=\sum_{i}P_{ji}\,{\mathcal{N}}_{i}^{(r)}\,,$ (12)

and

$T_{ji}=W_{j}\,P_{ji}\Big{/}\sum_{j^{\prime}}W_{j^{\prime}}\,P_{j^{\prime}i}\,.$
(13)

To provide more discussion of the conventions and mathematics here, in the
standard RL method, the variation of emphasis over the $Z$ domain is not
employed, i.e., $W\equiv 1$. We underscore the possibility of a different
emphasis across $Z$ for the sake of flexibility when the certainty in
$n(\zeta)$ varies across the domain. The variation in $W$ might also be used
to test the robustness of conclusions. Further on in this work we employ
$W_{j}=1$, except for rare cases of $n_{j}^{(r)}=0$ when we put $W_{j}=0$.
Within the equations behind the RL method, Eq. (1) represents forward mapping,
while (4) represents backward. Depending on the context, the matrix $P_{ji}$
acting on a vector from the left in (12) may be called a blurring matrix or
forward transfer matrix. The sum of elements in a column is bounded by 1,
$\sum_{j}P_{ji}\leq 1$, with the sum being exactly 1 when the probability is
preserved. The matrix $T_{ji}$ acting on a vector from the right in (11) may
be called a backward transfer matrix. The sum of elements in a column of that
matrix is exactly 1, $\sum_{j}T_{ji}=1$.

When $n$ and ${\mathcal{N}}$ extend over the same domain, it is common to use
$n$ to start the iterations for ${\mathcal{N}}$. It is common for the
iterations to progress quickly first, but then gradually stall, with seesaw
instabilities developing over many iterations. Remedies include an
acceleration [16] by using a power $\nu>1$ for the amplification factor $A$,
terminating the iterations after a moderate number of steps, and/or using a
regularization factor [16, 17]:

$I^{(r)}=\frac{1}{1-\lambda{\bf
H}\cdot{\bf\nabla}\,({\bf\nabla}{\mathcal{N}}^{(r)}/|{\bf\nabla}{\mathcal{N}}^{(r)}|)}\,.$
(14)

Here, ${\bf H}$ is the discretization vector within the domain of
${\mathcal{N}}$ and $\lambda$ is a small factor that prevents build-up of any
seesaw pattern in effect of the iterations. I.e., Eq. (10) can be replaced
with

${\mathcal{N}}_{i}^{(r+1)}({\xi})=\big{[}A_{i}^{(r)}(\xi)\big{]}^{\nu}\cdot
I_{i}^{(r)}\cdot{\mathcal{N}}_{i}^{(r)}({\xi})\,,$ (15)

in an accelerated regularized LR scheme. Excessive acceleration powers,
$\nu>2$, have been show to lead to an instability on its own in the
iterations.

Obviously other methods of inferring $n$ from ${\mathcal{N}}$ are possible,
such as in decomposing them in orthogonal functions, including harmonics [8,
12]. The danger in the latter case is of combining functions with varying sign
where noise is handled independently. In the case of functions that change by
orders of magnitude, the danger is of arriving at a poorer representation for
${\mathcal{N}}$ in the regions of low values than even provided by $n$.

##  III 1D Model

We next present a schematic and thus potentially transparent 1D model. It may
at first seem disconnected from the issue of reaction-plane deblurring. In
reality, it is closely related.

We consider a projectile moving at an unknown velocity $V$ that varies from an
event to an event. The projectile de-excites emitting $N=10$ identical
particles sampled from a uniform velocity distribution from -1 to +1 relative
to $V$, in some arbitrary velocity units. The particle velocities are measured
and the goal is to determine, from event statistics, the distribution of
particles in the projectile frame. For uncorrelated emission, the rms of the
sought distribution is obviously easy to estimate from the average difference
of velocities between any two particles squared:

$\langle(v_{1}^{\text{lab}}-v_{2}^{\text{lab}})^{2}\rangle=\langle(v_{1}-v_{2})^{2}\rangle=2\langle
v^{2}\rangle\,.$ (16)

Here, $v$ are velocities relative to the projectile and $v^{\text{lab}}$ are
in the laboratory system: $v^{\text{lab}}=v+V$.

In assessing velocity of a particle relative to the projectile, the velocity
of the projectile can be estimated with the average velocity of $N-1=9$
remaining particles, in a similar manner as when reaction plane is estimated
from remaining particle in an event [3]. The so assessed velocity $V^{\prime}$
fluctuates though around the true velocity $V$. The probability distribution
from simulated 10,000 events is shown as a solid curve in Fig. 1. The central-
limit theorem states that for a large particle number in one event such a
distribution approaches a Gaussian

$\frac{\text{d}P}{\text{d}V^{\prime}}\simeq\sqrt{\frac{(N-1)}{2\pi\langle
v^{2}\rangle}}\,\text{exp}\bigg{[}-\frac{(N-1)(V^{\prime}-V)^{2}}{2\langle
v^{2}\rangle}\bigg{]}\,.$ (17)

We compare that Gaussian distribution to the one from simulations and it is
apparent that the two distributions are nearly indistinguishable in the region
where they are significant.

Figure 1: Probability density for system velocity $V^{\prime}$, estimated with
the average velocity for $N-1=9$ particles sampled from a symmetric uniform
velocity distribution centered around the true velocity $V$. The solid (red)
line shows the distribution from a direct simulation, calculated with $H=0.1$
binning, and the dashed (blue) line - the distribution estimated with the
central-limit theorem.

Because of the fluctuations of the reference $V^{\prime}$, around the true
velocity $V$, the inferred velocities $v^{\prime}$ relative to the projectile
are inaccurate, $v^{\prime}$ = $v+V-V^{\prime}$, and the inferred distribution
in velocity gets smeared out:

$\frac{\text{d}N}{\text{d}v^{\prime}}(v^{\prime})=\int\text{d}V^{\prime}\,\frac{\text{d}P}{\text{d}V^{\prime}}(V^{\prime}-V)\,\frac{\text{d}N}{\text{d}v}(v)\,.$
(18)

In Fig. 2, we show the original distribution of particles relative to the
projectile in the simulations, together with that from simulated measurements
and the smearing is evident. In collecting the statistics, we use $H=0.1$ bin
size. The third distribution represented in Fig. 2 is one from deblurring
using the regularized and accelerated RL algorithm. In 1D, the regularization
factor is

$I_{i}^{(r)}=\begin{cases}\frac{1}{1-\lambda}\,,&\mbox{if
}\,\,{\mathcal{N}}_{i}^{(r)}<{\mathcal{N}}_{i-1,i+1}^{(r)}\,,\\\
\frac{1}{1+\lambda}\,,&\mbox{if
}\,\,{\mathcal{N}}_{i}^{(r)}>{\mathcal{N}}_{i-1,i+1}^{(r)}\,,\\\
1\,,&\mbox{otherwise}\,.\end{cases}$ (19)

We employ $\lambda=0.01$ and accelerating power $\nu=1.99$ and start with
${\mathcal{N}}^{(1)}=n$. Moreover, we carry out the iterations in (15) by
using the central-limit blurring function (17) with (16). In principle, it
could be possible to update the blurring functions during iterations to make
it consistent with the inferred d$N$/d$v$, for finite $N$, if desired accuracy
of restoration called for that. At the level of the statistics we employ, the
iterations that employ our more basic procedure stabilize at a distribution
that is largely consistent with the original.

Figure 2: Distribution of particles in velocity within the projectile frame
for 1D model: original represented by the solid (red) line, inferred from
measurements represented by the dashed (blue) line and restored through
deblurring and represented by the dash-dot (green) line. In the simulation,
10,000 of events with $N=10$ particles each were generated. A binning of
$H=0.1$ was used for the distributions. The restoration was done assuming the
central-limit blurring function.

##  IV Model with Detection Inefficiencies

Here, we expand the 1D model of the previous section assuming detection
inefficiencies. This continuous a connection to the reaction-plane deblurring,
though maybe the particular challenges are overemphasized compared to a
typical case of the reaction-plane determination. The overemphasis tests the
robustness of the restoration procedures which can be useful for other
applications in nuclear physics.

Again in the model we assume that a projectile moves at unknown velocity $V$
that varies from an event to event. The projectile de-excites emitting $N$
identical particles sampled from a uniform velocity distribution from -1 to +1
relative to $V$, in the arbitrary velocity units. The particle velocities are
measured and the goal is to determine, from event statistics, the distribution
of particles in the projectile frame. Now, however, we assume that the
particles are detected using a detector with efficiency $E(v^{\text{lab}})$
varying across the detection window
$|v^{\text{lab}}|<v^{\text{lab}}_{\text{max}}=2$:

$E(v^{\text{lab}})=\begin{cases}\frac{1}{2}+\frac{1}{2}\cos{\Big{(}\pi\,v^{\text{lab}}/v^{\text{lab}}_{\text{max}}\Big{)}}\,,&\mbox{if}\,\,|v^{\text{lab}}|<v^{\text{lab}}_{\text{max}}\,,\\\
0\,,&\mbox{otherwise}\,.\\\ \end{cases}$ (20)

The efficiency is illustrated in Fig. 3 with the superimposed distribution of
emitted particles for an exemplary projectile velocity $V=0.5$. Given that the
number $n$ of detected particles is generally reduced compared to the number
$N$ of emitted, $n\leq N$, we assume now $N=30$ particles and demand that a
minimum $n_{\text{min}}=10$, $n\geq n_{\text{min}}$ is detected for the events
to be analysed with the goal of finding $\text{d}N/\text{d}v$ in the
projectile frame.

Figure 3: Detection efficiency $E$ (solid line) as a function of particle
laboratory velocity $v^{\text{lab}}$, in the 1D model with inefficiencies,
with superimposed particle distribution $dN/dv^{\text{lab}}$ (dotted line) for
a projectile moving at velocity $V=0.5$ (short vertical dashed line).

In the sampling of the events, we distribute the projectile velocity $V$
uniformly in the laboratory frame over the range
$|V|<V_{\text{max}}=v_{\text{max}}+v^{\text{lab}}_{\text{max}}+H/2$. Here,
$H=0.1$ continues to be our bin size employed in discretizing distributions.
The precise value of $V_{\text{max}}$ is ultimately of a minor importance, as
the probability of accepting an event for analysis becomes very small when
$|V|$ gets close to $v^{\text{lab}}_{\text{max}}$ terminating the detector
acceptance. Varying $V_{\text{max}}$ in some range above then largely amounts
to adding or subtracting events that do not get analysed. Average distribution
of particles relative to the projectile velocity, for the original events and
for the simulated measurements, is illustrated in Fig. 4. Note that the
average in Fig. 4 is over the simulated events, that may or may not be
analysed. Any extra events that do not get analysed, due to
$n<n_{\text{min}}$, contribute a drop in the overall normalization of
$\text{d}n/\text{d}v^{\prime}$ distribution found in analysing the events,
without altering the shape of that distribution. That normalization eventually
drops out in the restoration.

Figure 4: Distribution of particles in velocity within the projectile frame
for 1D model with inefficiencies: original represented by the solid (red)
line, inferred from measurements represented by the dashed (blue) line and
restored through deblurring and represented by the dash-dot (green) line. The
simulated events originally contained $N=30$ particles. In the measurement the
particles were assumed to be accepted with the efficiency of Fig. 3 and the
events were processed if the number of accepted particles exceeded a minimum:
$n\geq n_{\text{min}}=10$. A binning of $H=0.1$ was used for the distributions
and in the restoration matrices.

In determining a contribution to $\text{d}n/\text{d}v^{\prime}$ from each
analysed event, we progress as before and for any particular particle we
estimate the projectile velocity $V^{\prime}$ as the average velocity of the
remaining $n-1$ particles in an event. In spite of the starting particle
number being significantly higher than before and only the low cut-off for
observed particle number $n_{\text{min}}$ coinciding with the prior starting
number, the shape of the measured distribution in Fig. 4 is now much more
degraded compared to the original, than in the previous model consideration of
Fig. 2. The impact of the significant, rapidly changing inefficiency on the
measurement, that can be expressed with conditional probabilities, and on
restoration, is discussed below.

On the technical side, the probability densities for determining the system
velocity to be $V^{\prime}$, while the original is $V$, and for the particle
velocity to be $v^{\prime}$ in the system frame, when the original is $v$, do
not depend anymore on just difference of density arguments, as in previous
example, but on each argument separately, and these densities are not
identical to each other anymore. In the measurement, when the projectile
velocity $V$ pushes against a region of reduced detection efficiency, one can
expect the velocity $V^{\prime}$, determined from decay products, to be biased
towards the region of higher efficiency. This is indeed seen in terms of
conditional probability density
$\frac{\text{d}P}{\text{d}V^{\prime}}(V^{\prime}|V)$, the examples of which
are shown in Fig. 5. Also expected drop in the overall probability of an event
getting accepted for the analysis can be seen as $V$ gets closer to the
particle detection limits. As far as the width of the distributions for a
given $V$ is concerned, the quickly changing inefficiency has a dual effect,
of increasing the width due to a reduced accepted particle number and
decreasing due to shrinking of the accepted particle distribution. In the end,
the width is only modestly changing with $V$ and it is actually shrunk
compared to the previous case with 100% detection but fewer particles at the
start.

Figure 5: Conditional probability density for determining velocity
$V^{\prime}$, from decay products in the 1D model with inefficiencies, at
three indicated velocities $V$ of the projectile. Those original velocities
$V$ are marked with the short vertical dashed lines by the abscissa. Following
the strategy of using the remainder of the system as a reference for any
single particle, in generating the illustrated density, $N-1=29$ particles can
be sampled for an incomplete event and the set analysed if at least
$n_{\text{min}}-1=9$ particles pass the detector acceptance.

The conditional probability for the velocity relative to projectile frame,
$\frac{\text{d}P}{\text{d}v^{\prime}}(v^{\prime}|v)$, can be represented as a
convolution of the efficiency function $E$ with the conditional probability
for the projectile velocity,
$\frac{\text{d}P}{\text{d}V^{\prime}}(V^{\prime}|V)$:

$\frac{\text{d}P}{\text{d}v^{\prime}}(v^{\prime}|v)=\frac{1}{2V_{\text{max}}}\int_{-V_{\text{max}}}^{V_{\text{max}}}\text{d}V\,E(V+v)\,\frac{\text{d}P}{\text{d}V^{\prime}}(V+v-v^{\prime}|V)\,.$
(21)

The conditional probability density at exemplary values of the original
velocity $v$, from the model simulations, is illustrated in Fig. 6. To
underscore, that density is established no matter what the value of
$\text{d}N/\text{d}v$ at $v$ is, even when that value is zero. Though then
that last distribution might seem to disappear from the ongoing consideration,
it is obviously there in $\frac{\text{d}P}{\text{d}v^{\prime}}(v^{\prime}|v)$,
as the estimation of the projectile velocity $V^{\prime}$ depends on the
overall $\text{d}N/\text{d}v$, within
$\frac{\text{d}P}{\text{d}V^{\prime}}(V^{\prime}|V)$ in (21). In the model in
the previous section, that inherent dependence could be largely described in
terms of the easily estimated dispersion for the distribution, but here this
dependence is generally more complex. The measured distribution of particles
in relative velocity, $\text{d}n/\text{d}v^{\prime}$, is related to the
original distribution $\text{d}N/\text{d}v$ with, cf. Eq. (1),

$\frac{\text{d}n}{\text{d}v^{\prime}}=\int\text{d}v\,\frac{\text{d}P}{\text{d}v^{\prime}}(v^{\prime}|v)\,\frac{\text{d}N}{\text{d}v}\,.$
(22)

Figure 6: Conditional probability density for finding velocity $v^{\prime}$
relative to the projectile, for a particle measured in an event in the 1D
model with inefficiencies, at three indicated true velocities $v$ relative to
the projectile. Those original velocities $v$ are marked with the short
vertical dashed lines by the abscissa.

In comparing the conditional probability density for relative velocity
$\frac{\text{d}P}{\text{d}v^{\prime}}(v^{\prime}|v)$, of Fig. 6 to the density
for projectile velocity $\frac{\text{d}P}{\text{d}V^{\prime}}(V^{\prime}|V)$,
of Fig. 5, we can see that the density for the relative velocity is much
wider, exhibits nearly as much bias as the other for finite velocity values,
and it more clearly exhibits skewness. The width and bias in the conditional
$\frac{\text{d}P}{\text{d}v^{\prime}}(v^{\prime}|v)$ have their origin in the
widths and biases of contributing conditional
$\frac{\text{d}P}{\text{d}V^{\prime}}(V^{\prime}|V)$ in Eq. (21). The strong
increase in the widths for
$\frac{\text{d}P}{\text{d}v^{\prime}}(v^{\prime}|v)$ can be attributed to the
range of $\frac{\text{d}P}{\text{d}V^{\prime}}(V^{\prime}|V)$ contributions,
with different $V$ and changing bias. Regarding the bias, it may be observed
in Eq. (21) and otherwise, that $v^{\prime}-v=-(V-V^{\prime})$, so the bias
contributed to $\frac{\text{d}P}{\text{d}v^{\prime}}(v^{\prime}|v)$ by events
is opposite to the bias contributed to
$\frac{\text{d}P}{\text{d}V^{\prime}}(V^{\prime}|V)$. The fact that the bias
is primarily negative in the forward direction for
$\frac{\text{d}P}{\text{d}v^{\prime}}(v^{\prime}|v)$ in Fig. 6, means that the
forward relative velocity will be primarily measured for events where the
projectile velocity is in the negative direction, and the other way around.
This is consistent with the singled out velocity and reference particle set
fitting in an optimal manner within the detector acceptance window, one
pushing the front and the other the rear of that window. The skewness can be
tied to the fact that large widths become comparable to the pace at which
efficiency is changing. Both the probability for analysing an instance of a
given relative velocity $v$,
$P(v)=\int\text{d}v^{\prime}\,\frac{\text{d}P}{\text{d}v^{\prime}}(v^{\prime}|v)$,
and net velocity $V$,
$P(V)=\int\text{d}V^{\prime}\,\frac{\text{d}P}{\text{d}V^{\prime}}(V^{\prime}|V)$,
fall with increase in magnitude of the respective velocity, see Figs. 5 and 6.
Overall, the large width of the conditional probability density
$\frac{\text{d}P}{\text{d}v^{\prime}}(v^{\prime}|v)$, larger for central $v$
than in the model of the preceding section, the bias, not there before in the
model, and the fall of the analysis probability with $|v|$, again not there
before, all contribute to much stronger deterioration of the shape of
distribution in relative velocity in the measurements compared to the original
shape, cf. Figs. 2 and 4.

For restoration of $\text{d}N/\text{d}v$, Eq. (22) is discretized. The
difficulty now is that the conditional probability density at iteration $r$,
$\frac{\text{d}P^{(r)}}{\text{d}v^{\prime}}(v^{\prime}|v)$, depends in a
nontrivial manner on ${\text{d}N^{(r)}}/{\text{d}v}$. We start the iterations
of Sec. II with ${\text{d}N^{(1)}}/{\text{d}v}=\langle
E\rangle^{-1}\,{\text{d}n}/{\text{d}v^{\prime}}$, at $v^{\prime}=v$, and we
compute $\frac{\text{d}P^{(1)}}{\text{d}v^{\prime}}(v^{\prime}|v)$ through a
simulation. Here, $\langle E\rangle$ is the typical detector acceptance that
we conservatively take as average efficiency over the detection window, i.e.,
$\langle E\rangle=0.5$. Simulations for determining
$\frac{\text{d}P^{(r)}}{\text{d}v^{\prime}}(v^{\prime}|v)$ are numerically
costly compared to the matrix multiplications in the RL algorithm and this is
likely to be the case in similar situations of restoration. We tame increase
in the cost by updating
$\frac{\text{d}P^{(r)}}{\text{d}v^{\prime}}(v^{\prime}|v)$ only once in a
while during the iterations, infrequently enough for the RL iterations on
their own to converge in-between. As is illustrated in Fig. 4, despite of the
strong detector inefficiencies, strongly distorted measured distribution, and
greater complexity in the restoration, the procedure can restore the original
single-particle distribution sufficiently well for typical practical
applications.

##  V Reaction-Plane Deblurring

The problem of the blurring of particle distributions and other
characteristics of heavy-ion collisions lies in attributing azimuthal angle
$\Phi^{\prime}$ to the reaction plane that is actually directed at angle
$\Phi$. With this, the angles relative to the reaction plane attributed to the
particles are inaccurate, $\phi^{\prime}=\phi+\Phi-\Phi^{\prime}$, and
measured distributions are going to be smeared out:

$\frac{\text{d}^{3}\,N}{p_{\perp}\,\text{d}p_{\perp}\,\text{d}y\,\text{d}\phi^{\prime}}=\int\text{d}\Phi^{\prime}\,\frac{\text{d}P}{\text{d}\Phi^{\prime}}(\Phi^{\prime}-\Phi)\,\frac{\text{d}^{3}\,N}{p_{\perp}\,\text{d}p_{\perp}\,\text{d}y\,\text{d}\phi}(\phi)\,.$
(23)

Owing to the asymmetries in emission tied with the reaction plane, the
direction of the reaction plane is estimated with a combination of particle
momenta, most commonly sum of weighted transverse momenta, ${\bf
q}_{\nu}=\omega_{\nu}\,{\bf p}_{\perp\nu}$. When the direction of the reaction
plane is used as a reference for a certain particle $\mu$, that particular
particle is omitted from the reference, to avoid a self-correlation. I.e., the
reaction plane direction is estimated with the direction of the vector [3, 18]

${\bf Q}_{\mu^{\prime}}=\sum_{\nu\neq\mu}\omega_{\nu}\,{\bf
p}_{\perp\nu}=\sum_{\nu\neq\mu}{\bf q}_{\nu}\,.$ (24)

The direction is obviously the same as for the average vector ${\bf q}$ of the
remaining particles

$\overline{\bf q_{\mu^{\prime}}}=\frac{1}{N-1}\sum_{\nu\neq\mu}{\bf
q}_{\nu}\,,$ (25)

where $N$ is particle multiplicity. When ${\bf q}_{\nu}$ are largely
uncorrelated, except through the reaction-plane geometry, the central-limit
theorem implies the distribution of the average vector $\overline{\bf q}$ in
relation to the reaction plane approaches Gaussian for large $N$:

$\frac{\text{d}^{2}P}{\text{d}^{2}\overline{q}}\simeq\frac{N-1}{2\pi\,\sigma_{x}\,\sigma_{y}}\,\text{exp}\bigg{\\{}-\frac{N-1}{2}\bigg{[}\frac{\big{(}\overline{q}_{x}-\langle
q_{x}\rangle\big{)}^{2}}{\sigma_{x}^{2}}+\frac{\overline{q}_{y}^{2}}{\sigma_{y}^{2}}\bigg{]}\bigg{\\}}\,.$
(26)

Here, index $x$ is for a component along the reaction-plane direction and $y$
– along the perpendicular direction. The dispersions $\sigma$ and average
$\langle q_{x}\rangle$ refer to single-particle emission in relation to the
reaction plane,

$\begin{split}\sigma_{x}^{2}&=\langle(q_{x}-\langle
q_{x}\rangle)^{2}\rangle=\langle q_{x}^{2}\rangle-\langle
q_{x}\rangle^{2}\,,\\\ \sigma_{y}^{2}&=\langle q_{y}^{2}\rangle\,,\end{split}$
(27)

Under the condition that the correlations tied to the reaction plane dominate
[5], the parameters for the probability distribution can be determined from
constructs out of the ${\bf q}$ contributions averaged over particle sets and
events:

$\begin{split}\langle{\bf q}_{1}\cdot{\bf q}_{2}\rangle&=\langle
q_{x}\rangle^{2}\,,\\\ \langle 2({\bf q}_{1}\cdot{\bf q}_{2})({\bf
q}_{2}\cdot{\bf q}_{3})-({\bf q}_{1}\cdot{\bf
q}_{3})\,q_{2}^{2}\rangle&=\langle q_{x}\rangle^{2}\,(\langle
q_{x}^{2}\rangle-\langle q_{y}^{2}\rangle)\,.\end{split}$ (28)

Here, the indices 1, 2 and 3 refer to any three particles contributing to
${\bf Q}$ in an event.

Now, we can express $\overline{\bf q}$ in terms of azimuthal deviation from
the true direction of the reaction plane, $\Delta\Phi=\Phi^{\prime}-\Phi$

$\overline{q}_{x}=\overline{q}\cos{\Delta\Phi}\,,\hskip
20.00003pt\overline{q}_{y}=\overline{q}\sin{\Delta\Phi}\,,$ (29)

and write the probability density (26) as

$\frac{\text{d}^{2}P}{\text{d}^{2}\overline{q}}=\frac{\text{d}^{2}P}{\overline{q}\,\text{d}\overline{q}\,\text{d}\Phi^{\prime}}\simeq{\mathcal{A}}\,\text{exp}\big{[}-{\mathcal{B}}(\Delta\Phi)\,\overline{q}^{2}+2{\mathcal{C}}(\Delta\Phi)\,\overline{q}\big{]}\,.$
(30)

Integration over magnitude of $\overline{q}$ yields the central-limit
distribution of the estimated reaction plane relative to the true plane:

$\frac{\text{d}P}{\text{d}\Phi^{\prime}}=\int\text{d}\overline{q}\,\overline{q}\,\frac{\text{d}^{2}P}{\text{d}^{2}\overline{q}}\simeq\alpha(\Delta\Phi)\,\Big{\\{}1+\sqrt{\pi}\,\beta(\Delta\Phi)\,\exp{\big{(}\beta^{2}(\Delta\Phi)\big{)}}\big{[}1+\text{erf}\big{(}\beta(\Delta\Phi)\big{)}\big{]}\Big{\\}}\,.$
(31)

Here, the terms $\beta$ and $\alpha$, both dependent on the deviation
$\Delta\Phi$ from the true reaction plane, are

$\beta(\Delta\Phi)=\frac{{\mathcal{C}}(\Delta\Phi)}{\sqrt{{\mathcal{B}}(\Delta\Phi)}}=\sqrt{\frac{N-1}{2\big{(}\sigma_{y}^{2}\cos^{2}{\Delta\Phi}+\sigma_{x}^{2}\sin^{2}{\Delta\Phi}\big{)}}}\,\frac{\sigma_{y}}{\sigma_{x}}\,\langle
q_{x}\rangle\cos{\Delta\Phi}\approx\sqrt{\frac{N-1}{\langle
q^{2}\rangle}}\,\langle q_{x}\rangle\,\cos{\Delta\Phi}\,.$ (32)

and

$\alpha(\Delta\Phi)=\frac{\mathcal{A}}{2{\mathcal{B}}(\Delta\Phi)}=\frac{\sigma_{x}\,\sigma_{y}\,\exp{\big{[}-\frac{(N-1)\langle
q_{x}\rangle^{2}}{2\sigma_{x}^{2}}\big{]}}}{2\pi\,\big{(}\sigma_{y}^{2}\,\cos^{2}{\Delta\Phi}+\sigma_{x}^{2}\,\sin^{2}{\Delta\Phi}\big{)}}\approx\frac{1}{2\pi}\,\exp{\bigg{[}-\frac{(N-1)\langle
q_{x}\rangle^{2}}{\langle q^{2}\rangle}\bigg{]}}\,,$ (33)

see also Ref. [8]. The final approximations in Eqs. (32) and (33) pertain to
the case of weak azimuthal anisotropies in the distributions relative to the
reaction plane. These underscore that the main variation with angle in (31) is
tied to the factor of $\cos{\Delta\Phi}$ in (32).

An example of the distribution (31) is shown in Fig. 8. We next illustrate the
envisioned procedure of reaction plane deblurring in the case of the final
state modelled using a local-equilibrium model that incorporates common
concepts for intemediate-energy heavy-ion collisions.

###  V.1 Deblurring Example: Local-Equilibrium Model

We illustrate deblurring for a 3D distribution associated with the reaction
plane, deducing
$\text{d}^{3}\,N/p_{\perp}\,\text{d}p_{\perp}\,\text{d}y\,\text{d}\phi$ from
$\text{d}^{3}\,N/p_{\perp}\,\text{d}p_{\perp}\,\text{d}y\,\text{d}\phi^{\prime}$,
c.f. Eq. (23), by simulating the final state of a collision in a local
equilibrium model. As a reference beam energy for the collision, we take
$300\,\text{MeV/nucl}$ and consider $A\leq 4$ charged products in the final
state.

Our goal is to mimic a state that is consistent at a very coarse level with
the observations [18, 19]. We have no intention to dive into any controversies
around the freeze-out in collisions, whether the colliding systems literarily
go through an equilibrium when expanding into vacuum or whether they
separately freeze out chemically and kinematically. Specifically, we assume
that the system freeze-out can be described in terms of a local equilibrium,
at a uniform freeze-out density $\rho_{f}$ and temperature $T_{f}$, combined
with collective motion for the freeze-out location. Within the latter, the
longitudinal expansion dominates, but also present are a radial transverse
expansion, with a weak elliptical modulation, and finally a weak sideward
flow, typical for mid-central collisions. We take the freeze-out density at
$\rho_{f}=\rho_{0}/6$ and, for simplicity, a neutron-proton symmetric system.
Assuming that the local kinetic energy for nucleons is a fraction of the
kinetic energy available per nucleon, we arrive at the freeze-out temperature
$T_{f}\simeq 29\,\text{MeV}$. To fix the attention, we assume that the average
number of protons included in the measurements, whether as free or in
clusters, is $\langle Z_{t}\rangle=50$. Then, $\rho_{f}$ and $T_{f}$ give us
charged particle multiplicities, $\langle N_{p}\rangle=19.0$, $\langle
N_{d}\rangle=17.8$, $\langle N_{t}\rangle=\langle N_{h}\rangle=3.7$, and
$\langle N_{\alpha}\rangle=1.0$. The $A\geq 3$ multiplicities are low, for a
given proton multiplicity, compared to the data from the specific beam energy
region [19]. However, our main goal is to explore the situation with different
species used in reaction plane determination, not a literal reproduction of
relative yield data. We sample multiplicities for the species from Poisson
distributions for the average specie multiplicities and we combine a momentum
sampled from a local equilibrium distribution for the specific species with a
collective velocity boost that itself combines the aforementioned persistence
of longitudinal motion, radial transverse expansion and elliptic and sideward
flows, all of the latter common between the species. In Fig. 7, we show
exemplary characteristics of protons and deuterons in the generated events.
Second-order coefficients of azimuthal anisotropy have low values in this
simulation, in fact for protons quite comparable to statistical errors even
for a large number of events, and are not shown.

Figure 7: Characteristics of protons and deuterons within the local
equilibrium model. The lines represent differential distributions in
normalized rapidity, $\text{d}N_{X}/\text{d}y_{R}$, and symbols represent
first-order azimuthal-asymmetry coefficients,
$v_{1X}=\langle\cos{(\phi-\Phi)}\rangle_{X}$, all vs $y_{R}$. Here, the
rapidity is in the center of mass and scaled with the rapidity of the beam,
$y_{R}=y/y_{\text{Beam}}$.

In each of the events in the sample, when considering particle $\mu$, we
estimate the reaction plane direction with the remaining particles in the
event. We use

$\omega_{\nu}=\begin{cases}\text{sgn}\,y_{R}\,,&\text{if $|y_{R}|>\delta$},\\\
0\,,&\text{if $|y_{R}|<\delta$},\end{cases}$ (34)

in Eq. (24), with $\delta=0.17$ [3]. The distribution of estimated plane
directions, relative to the true direction, is shown in Fig. 8 for this choice
of weights. More optimal weights could be chosen, yielding a more narrow
distribution [20, 18]. Our goal here, though, is of presenting deblurring
opportunities and not of making of making everything optimal - actual
experimental analysis will likely face other tensions. In addition to the
distribution found in the simulations, we show in Fig. 8 the central-limit
result (31), with widths from Eqs. (27) and (28), and $N$ replaced by the
average $\langle N\rangle$ for the simulation. Remarkably, the curves cannot
be distinguished by eye, which bodes well for using central-limit results in
practice, in lieu of any complete self-consistent simulation. Notably, testing
of the proximity to the central-limit can be carried out directly in the
experiment, by randomly dividing events into subevents of nearly equal
multiplicity [3] and comparing the relative distribution of estimated
directions of the reaction plane to that from Eq. (31) with $(N-1)$ replaced
by $\langle N\rangle/4$.

Figure 8: Probability density for the azimuthal angle $\Phi^{\prime}$ of the
reaction plane in the local equilibrium model, when estimated with the
direction of the vector $\overline{\bf q}$ out of $N-1$ particles in an event,
relative to the true angle $\Phi$. Here, $N$ is the total particle number in
an event. The solid (red) line shows the density from a direct simulation in
the local equilibrium model and the dashed (blue) line shows the density
estimated with the central-limit Eq. (31). In using the last equation,
$\langle q_{x}\rangle$, $\sigma_{x}^{2}$ and $\sigma_{y}^{2}$ in $\alpha$ and
$\beta$ were obtained from averages over the event simulations, following the
strategy laid out in the text. In addition, the multiplicity $N$, actually
varying from an event to event, was replaced with the single average value
$\langle N\rangle$ over the simulated events.

We next illustrate the blurring that occurs when attempting to measure the
triple differential distribution. To the blurred distribution we apply
deblurring. Specifically, we examine the spectrum of deuterons at a moderately
forward rapidity, $y_{R}=0.5$. We choose deuterons rather than protons in the
illustration, because the push of the latter in the reaction plane is
relatively meager in the parameterization of the flows we chose. Figure 9
shows three versions of the transverse momentum spectra in the reaction plane
($p_{y}=0$): as recorded for the event sample using the known direction of the
reaction plane (diamonds), as recorded (crosses) when using the direction
estimated with a construct from particle momenta, Eqs. (24) and (25), and as
obtained from deblurring the blurred spectrum (circles). The blurring function
for the latter could be determined self-consistently, through additional
simulations during the deblurring iterations, but given Fig. 8 we instead use
the central-limit formula (31). The discretization step in the RL algorithm is
the same for the blurring function and the distribution and we use
$H=15\degree$. The values from restoration settle after about 5 RL iterations
with $\nu=1.99$ and $\lambda=0.005$. The deblurring is carried out for each
transverse momentum bin separately and the inset in Fig. 9 shows the case of a
high central bin momentum, $p_{\perp}=1.17\,\text{GeV}/c$, with statistical
fluctuations evident.

Figure 9: Triple differential deuteron distribution at $y_{R}=0.5$, in the
local equilibrium model. The main figure shows the distribution in the
reaction plane, as a function of the momentum component in the reaction plane.
The inset shows the distribution at $p_{\perp}=1.17\,\text{GeV}/c$, as a
function of azimuthal angle about the reaction plane. In each case, three
versions of the distribution are displayed: one determined directly using the
known direction of the reaction plane in the modelled events (diamonds), one
obtained by estimating the direction of the plane with the vector ${\bf Q}$
(crosses), out of momenta of remaining particles in the simulated events, Eq.
(24), and one arrived by applying the deblurring to the measured triple-
differential distribution, with the RL algorithm applied to Eq. (23) and the
blurring function from the central limit theorem, Eq. (31).

It is apparent in Fig. 9 that the restored spectrum well coincides with the
original, as the respective representing symbols largely overlap. It is
apparent, for the high momenta there, that the restoration can reproduce
spectra varying by an order of magnitude over the azimuthal angle. Much of the
variation over angle is already there in the simulated measured spectrum, see
the inset. The restoration just improves details in that variation and does so
mostly on a local scale in the azimuthal angle. Have we chosen different
weights for $\bf Q$ in Eq. (34), the blurring function in Fig. 8 could have
been more narrow, and there would be less needed improvement left for the
deblurring.

For practitioners of the analysis of heavy-ion collisions, it should be
apparent that the considerations in the 1D model earlier are just veiled
versions of the reaction-plane analysis, stripped of the dimension that plays
no direct role in the restoration and stripped of the periodicity in the other
dimension. With this, the steps that must be taken to treat detector
inefficiencies consistently with the rest of the reaction-plane analysis
should be apparent from the preceding section. Particle misidentifications
would give rise to matrix elements in the transfer matrix coupling not only
different particle identities, but also rapidities. Interestingly, it might
even be possible to handle a situation where the particles could only be
identified on a statistical level [5].

###  V.2 Why 3D Characteristics?

It may be asked why engage in deblurring, if procedures exist for determining
azimuthal Fourier coefficients, as a functions of transverse momentum [9, 10,
11, 12]. The lowest nonvanishing coefficients usually have a straightfoward
physical meaning and may be easier to determine for low event statistics than
refined distributions. However, a multitude of coefficients of different
orders, changing with rapidity and transverse momentum, combining information
from different azimuthal directions, can be a challenge at operational level.
To illustrate how the view of reactions could be expanded by examining
distributions that include azimuthal angle relative to the reaction plane, and
potentially 3D distributions, we reach to the pBUU transport model [21]. The
latter and other transport models have been extensively used to simulate a
variety of heavy-ion collisions in the beam energy region from few tens of
MeV/nucl to few GeV/nucl, describing different observables at quantitative and
semi-quantitative level [22]. With this predictions of such models might be
considered pertinent, but as azimuthal and more broadly 3D spectra have not
been accessible experimentally, they were normally not considered in the
theory either.

Multiple transport simulations have been carried out in the context of
$270\,\text{MeV/nucl}$ Sn + Sn experiments performed recently at RIKEN [23].
In Fig. 10, we show exemplary results from such simulations, proton and
neutron spectra from 132Sn + 124Sn collisions at $b=3.3\,\text{fm}$, within
the reaction plane at $y_{R}=0.5$. With projectile and target of different
mass, the rapidity is taken here in the nucleon-nucleon center of mass and is
normalized to the beam. There are similarities and differences between the
spectra in the conventional local equilibrium picture, such as in Fig. 9, and
from transport, such as in Fig. 10, when the spectra are examined in the
reaction plane. Both types of spectra have maxima shifted in the positive
direction of the reaction plane. The one in the conventional picture is
largely symmetric about the fairly flat maximum and close to parabolic in the
logarithmic scale. The maxima in the spectra from transport are sharp, the
spectra are largely piecewise exponential in momentum, have knees and there is
visible asymmetry between the sides extending into the positive and negative
sides of the reaction plane. The slope for the positive side, right after the
maximum, is sharper than for the negative side. This could be due to a larger
fraction of spectator matter towards the positive side of the reaction plane
in space, for particles moving forward, than towards the negative. The knees
in the spectra are at different transverse momenta on the two sides of the
reaction plane. A sharp maximum could be produced by Coulomb interactions, but
these cannot explain the maximum for neutrons, that seems to be even sharper
than for protons.

Figure 10: Triple differential distributions in the reaction plane at
$y_{R}=0.5$, for neutrons (circles) and protons (diamonds), from the pBUU
simulations of 132Sn + 124Sn collisions at $270\,\text{MeV/nucl}$ and
$b=3.3\,\text{fm}$. The rapidity $y_{R}$ is here in the nucleon-nucleon center
of mass and normalized to the beam.

When the transport spectra are averaged over azimuthal angle, the maxima in
the spectra move to $p_{\perp}=0$ and soften. The knees soften too. Overall,
the two types of spectra, from transport and the conventional picture, become
qualitatively much closer to each other with the azimuthal angle-averaging
than without. When azimuthal asymmetries in the spectra are explored in terms
of Fourier coefficients, usually just one or two lowest ones in the particular
energy regime, details such as in Fig. 10 are beyond resolution.

Differences in the spectra for the sides of the system with different
participant-spectator composition on different sides of the reaction plane may
be accessible only through the 3D examination of those spectra. Shift of a
maximum in the spectrum within the reaction plane may offer the only
opportunity to examine its shape, as a maximum can be difficult to assess when
it coincides with beam direction. Comparison of the knees in relation to the
maximum, between experiment and theory may help to clarify their origin and
clarify the level of understanding of the collisions within theory, such as of
the elementary collisions taking place on shell.

Turning to the second example, heavy-ion collisions compress nuclear matter to
densities above normal at conditions approaching thermal equilibrium. With
this they represent opportunity to learn about the nuclear equation of state
(EOS). The current interest is in the component of EOS, so-called symmetry
energy that describes energy change with the change in relative neutron-proton
asymmetry, at different net nucleonic densities. Energies of more massive
nuclei, and other data, constrain the symmetry energy at subnormal densities
and especially at the density representing an average for the nuclei,
$\rho\simeq 2\rho_{0}/3$, at a value of about $25.4\,\text{MeV}$ [24]. From
there on, the region in density opens up where the symmetry energy is poorly
constrained and this includes the pace of symmetry-energy variation around
$\rho_{0}$. Yet when different parametrizations of the symmetry energy are
explored in transport simulations, that pass near the asserted value, it is
very hard to find sensitivity to the symmetry energy at $\rho\gtrsim\rho_{0}$,
in nearly any predicted observable. One reason are low asymmetries for the
nuclear systems as a whole, additionally depleted in the center of the matter
due to migration of the asymmetry to the surface. There is some advantage in
going to heavier systems, as surface to volume ratio drops and Coulomb
interactions help to maintain significant interior asymmetry, both in the
initial state and in reaction dynamics [25]. Access to the 3D information
should help too. Indeed, if one examines nucleonic spectra in direction
perpendicular to the reaction plane, in $250\,\text{MeV/nucl}$ 208Pb + 208Pb
collisions at $b=4\,\text{fm}$, at $y_{R}=0$, see Fig. 11, one can observe
definite changes in the ratio of neutron to proton yields, when the symmetry
energy evolves from moderately soft to stiff, with its slope parameter at
$\rho_{0}$ changing from $L=38.7\,\text{MeV}$ to $L=105.5\,\text{MeV}$. Both
parametrizations of the symmetry pass by the consensus value at $\rho\sim
2\rho_{0}/3$.

Figure 11: Ratio of neutron yield to proton yield in $250\,\text{MeV/nucl}$
208Pb + 208Pb collisions at $b=4\,\text{fm}$, at midrapidity, in the direction
perpendicular to the reaction plane, as a function of transverse momentum. The
two symbol sets represent the results for two different parametrizations of
symmetry energy: diamonds - stiff with $L=105.5\,\text{MeV}$ and circles -
moderately soft with $L=38.7\,\text{MeV}$. The two parametrizations are chosen
so that they yield the symmetry energy at $\rho\sim 2\rho_{0}/3$ consistent
with the conclusions from binding energies of heavy nuclei and other
determinations pertaining to subnormal densities.

The fact that a stiff symmetry energy raises the neutron-proton ratio at high
transverse momentum agrees with intuitive expectations. However, the fact that
it lowers the ratio at low momenta is somewhat surprising. Likely the
unexpected features of both Figs. 10 or 11 are tied to evolution of particle
emission with time. The ratio being higher than the overall for the system of
1.54 in Fig. 11, for most transverse momenta, irrespective of the EOS, may be
due to Coulomb interactions that push protons out not only in the transverse
directions, but also along the beam axis.

In practical comparisons of transport theory to data it is common to make sure
that many more rudimentary aspects of the data are understood on equal footing
with those more sophisticated. I.e., other aspects of measured and calculated
distributions would need to be simultaneously checked and more calculations
would need to be carried out, than here, to draw credible physics conclusions.

##  VI Conclusions

We have proposed to apply deblurring, adopted from optics, to nuclear
observables subject to degradation, whether on account of the inference method
or detector performance, such as 3D characteristics relative to the reaction
plane in heavy-ion collisions. In some instances, such as that of the reaction
plane determination or the schematic 1D model, and sufficiently many particles
detected, the blurring function may be determined, relatively faithfully,
following the central limit theorem. In a more general situation, the blurring
function may need to be determined through simulations, possibly carried out
self-consistently with the deblurring, as for the 1D model with inefficiency
that cannot be linearized, or through independent measurements, as in
assessing detector performance. In heavy-ion collisions, the deblurring, such
as with the Richardson-Lucy method, can work locally in the azimuthal angle
only improving the 3D characteristics of collision events that are already
resolved to some extent with the estimated reaction plane. This is a potential
improvement over the complete Fourier decomposition of 3D characteristics in
the azimuth. Importantly, estimation of the reaction plane can be made with
particles measured in a different detector than the one used to detect the
particle or particles for which the 3D characteristics are desired.

In analysis of actual heavy-ion data, other correlations, than those
associated with the reaction plane, will be present in the final states, in
particular tied to the total momentum conservation, interactions at low
relative velocity and sequential decays. Strategies to deal with these have
been developed in other contexts [5, 26, 12] and first the base needs to be
established ignoring these correlations, which has been attempted here.

Finally, we hope that the deblurring, revealing 3D distributions, could help
extract novel physics information from the heavy-ion collisions, that we tried
to illustrate. Beyond single-particle distributions tied to the reaction
plane, we hope that the strategy could be used for low-velocity correlations
[27].

###### Acknowledgements.

The authors benefited from discussions with Scott Pratt and with members of
the S$\pi$RIT Collaboration. This work was supported by the U.S. Department of
Energy Office of Science under Grant DE-SC0019209.

## References

  * Gyulassy  _et al._ [1982] M. Gyulassy, K. A. Frankel, and H. Stöcker, Physics Letters B 110, 185 (1982).
  * Danielewicz and Gyulassy [1983] P. Danielewicz and M. Gyulassy, Physics Letters B 129, 283 (1983).
  * Danielewicz and Odyniec [1985] P. Danielewicz and G. Odyniec, Physics Letters B 157, 146 (1985), LBL Report 18600 (1984).
  * Ollitrault [1992] J.-Y. Ollitrault, Physical Review D 46, 229 (1992).
  * Danielewicz  _et al._ [1988] P. Danielewicz, H. Ströbele, G. Odyniec, D. Bangert, R. Bock, R. Brockmann, J. W. Harris, H. G. Pugh, W. Rauch, R. E. Renfordt, A. Sandoval, D. Schall, L. S. Schroeder, and R. Stock, Physical Review C 38, 120 (1988).
  * Demoulins  _et al._ [1990] M. Demoulins, D. L’Hôte, J. P. Alard, J. Augerat, R. Babinet, N. Bastid, F. Brochard, C. Cavata, N. De Marco, P. Dupieux, H. Fanet, Z. Fodor, L. Fraysse, P. Gorodetzky, J. Gosset, T. Hayashino, M. C. Lemaire, A. Le Merby, B. Lucas, J. Marroncle, G. Montarou, M. J. Parizet, J. Poitou, C. Racca, W. Schimmerling, Y. Terrien, and O. Valette, Physics Letters B 241, 476 (1990).
  * Demoulins [1989] M. Demoulins,  _Collective Flow Studies in Central Collisions Between Nuclei at Several Hundreds of MeV per Nucleon_ , Ph.D., L’Universite de Paris-Sud, Orsay (1989), Report CEA-N-2628.
  * Voloshin and Zhang [1996] S. Voloshin and Y. Zhang, Zeitschrift für Physik C Particles and Fields 70, 665 (1996).
  * Poskanzer and Voloshin [1998] A. M. Poskanzer and S. A. Voloshin, Physical Review C 58, 1671 (1998).
  * Borghini  _et al._ [2001] N. Borghini, P. M. Dinh, and J.-Y. Ollitrault, Physical Review C 64, 054901 (2001).
  * Bilandzic  _et al._ [2011] A. Bilandzic, R. Snellings, and S. Voloshin, Physical Review C 83, 044913 (2011).
  * HADES Collaboration  _et al._ [2020] HADES Collaboration, J. Adamczewski-Musch, O. Arnold, C. Behnke, A. Belounnas, A. Belyaev, J. Berger-Chen, A. Blanco, C. Blume, M. Böhmer, P. Bordalo, S. Chernenko, L. Chlad, I. Ciepal, C. Deveaux, J. Dreyer, E. Epple, L. Fabbietti, O. Fateev, P. Filip, P. Fonte, C. Franco, J. Friese, I. Fröhlich, T. Galatyuk, J. Garzón, R. Gernhäuser, O. Golosov, M. Golubeva, R. Greifenhagen, F. Guber, M. Gumberidze, S. Harabasz, T. Heinz, T. Hennino, S. Hlavac, C. Höhne, R. Holzmann, A. Ierusalimov, A. Ivashkin, B. Kämpfer, T. Karavicheva, B. Kardan, I. Koenig, W. Koenig, M. Kohls, B. Kolb, G. Korcyl, G. Kornakov, F. Kornas, R. Kotte, A. Kugler, T. Kunz, A. Kurepin, A. Kurilkin, P. Kurilkin, V. Ladygin, R. Lalik, K. Lapidus, A. Lebedev, L. Lopes, M. Lorenz, T. Mahmoud, L. Maier, A. Malige, M. Mamaev, A. Mangiarotti, J. Markert, T. Matulewicz, S. Maurus, V. Metag, J. Michel, D. Mihaylov, S. Morozov, C. Müntz, R. Münzer, L. Naumann, K. Nowakowski, Y. Parpottas, V. Pechenov, O. Pechenova, O. Petukhov, K. Piasecki, J. Pietraszko, W. Przygoda, K. Pysz, S. Ramos, B. Ramstein, N. Rathod, A. Reshetin, P. Rodriguez-Ramos, P. Rosier, A. Rost, A. Rustamov, A. Sadovsky, P. Salabura, T. Scheib, H. Schuldes, E. Schwab, F. Scozzi, F. Seck, P. Sellheim, I. Selyuzhenkov, J. Siebenson, L. Silva, U. Singh, J. Smyrski, Y. Sobolev, S. Spataro, S. Spies, H. Ströbele, J. Stroth, C. Sturm, O. Svoboda, M. Szala, P. Tlusty, M. Traxler, H. Tsertos, E. Usenko, V. Wagner, C. Wendisch, M. Wiebusch, J. Wirth, D. Wójcik, Y. Zanevsky, and P. Zumbruch, Physical Review Letters 125, 262301 (2020).
  * Vankawala  _et al._ [2015] F. Vankawala, A. Ganatra, and A. Patel, International Journal of Computer Applications 116, 15 (2015).
  * Richardson [1972] W. H. Richardson, Journal of the Optical Society of America 62, 55 (1972).
  * Lucy [1974] L. B. Lucy, The Astronomical Journal 79, 745 (1974).
  * Remmele and Hesser [2009] S. Remmele and J. Hesser, in  _Bildverarbeitung für die Medizin 2009_, Informatik aktuell, edited by H.-P. Meinzer, T. M. Deserno, H. Handels, and T. Tolxdorff (Springer, Berlin, Heidelberg, 2009) pp. 400–404.
  * Dey  _et al._ [2006] N. Dey, L. Blanc‐Feraud, C. Zimmer, P. Roux, Z. Kam, J.-C. Olivo‐Marin, and J. Zerubia, Microscopy Research and Technique 69, 260 (2006).
  * Andronic  _et al._ [2006] A. Andronic, J. Łukasik, W. Reisdorf, and W. Trautmann, The European Physical Journal A - Hadrons and Nuclei 30, 31 (2006).
  * Reisdorf  _et al._ [2010] W. Reisdorf, A. Andronic, R. Averbeck, M. L. Benabderrahmane, O. N. Hartmann, N. Herrmann, K. D. Hildenbrand, T. I. Kang, Y. J. Kim, M. Kis̆, P. Koczoń, T. Kress, Y. Leifels, M. Merschmeyer, K. Piasecki, A. Schüttauf, M. Stockmeier, V. Barret, Z. Basrak, N. Bastid, R. C̆aplar, P. Crochet, P. Dupieux, M. Dz̆elalija, Z. Fodor, P. Gasik, Y. Grishkin, B. Hong, J. Kecskemeti, M. Kirejczyk, M. Korolija, R. Kotte, A. Lebedev, X. Lopez, T. Matulewicz, W. Neubert, M. Petrovici, F. Rami, M. S. Ryu, Z. Seres, B. Sikora, K. S. Sim, V. Simion, K. Siwek-Wilczyńska, V. Smolyankin, G. Stoicea, Z. Tymiński, K. Wiśniewski, D. Wohlfarth, Z. G. Xiao, H. S. Xu, I. Yushmanov, and A. Zhilin, Nuclear Physics A 848, 366 (2010).
  * Danielewicz [1995] P. Danielewicz, Physical Review C 51, 716 (1995).
  * Danielewicz [2000] P. Danielewicz, Nuclear Physics A 673, 375 (2000).
  * Wolter  _et al._ [2021] H. Wolter, M. Colonna, D. Cozma, P. Danielewicz, C. M. Ko, R. Kumar, A. Ono, M. B. Tsang, J. Xu, Y.-X. Zhang, Z.-Q. Feng, T. Gaitanos, A. Le Fevre, N. Ikeno, Y. Kim, S. Mallik, P. Napolitani, D. Oliinychenko, T. Ogawa, M. Papa, J. Su, R. Wang, Y.-J. Wang, J. Weil, F.-S. Zhang, G.-Q. Zhang, Z. Zhang, J. Aichelin, L.-W. Chen, H. Elfner, C. Hartnack, S. Hashimoto, S. Jeon, K. Kim, M. Kim, B.-A. Li, C.-H. Lee, Q.-F. Li, Z.-X. Li, U. Mosel, Y. Nara, K. Niita, A. Ohnishi, T. Sato, T. Song, A. Sorensen, and N. Wang (2021), submitted to Progress in Particle and Nuclear Physics.
  * Jhang  _et al._ [2021] G. Jhang, J. Estee, J. Barney, G. Cerizza, M. Kaneko, J. W. Lee, W. G. Lynch, T. Isobe, M. Kurata-Nishimura, T. Murakami, C. Y. Tsang, M. B. Tsang, R. Wang, D. S. Ahn, L. Atar, T. Aumann, H. Baba, K. Boretzky, J. Brzychczyk, N. Chiga, N. Fukuda, I. Gasparic, B. Hong, A. Horvat, K. Ieki, N. Inabe, Y. J. Kim, T. Kobayashi, Y. Kondo, P. Lasko, H. S. Lee, Y. Leifels, J. Łukasik, J. Manfredi, A. B. McIntosh, P. Morfouace, T. Nakamura, N. Nakatsuka, S. Nishimura, R. Olsen, H. Otsu, P. Pawłowski, K. Pelczar, D. Rossi, H. Sakurai, C. Santamaria, H. Sato, H. Scheit, R. Shane, Y. Shimizu, H. Simon, A. Snoch, A. Sochocka, Z. Sosin, T. Sumikama, H. Suzuki, D. Suzuki, H. Takeda, S. Tangwancharoen, H. Toernqvist, Y. Togano, Z. G. Xiao, S. J. Yennello, J. Yurkon, Y. Zhang, M. Colonna, D. Cozma, P. Danielewicz, H. Elfner, N. Ikeno, C. M. Ko, J. Mohs, D. Oliinychenko, A. Ono, J. Su, Y. J. Wang, H. Wolter, J. Xu, Y.-X. Zhang, and Z. Zhang, Physics Letters B 813, 136016 (2021).
  * Lynch and Tsang [2021] W. G. Lynch and M. B. Tsang, arXiv:2106.10119 [astro-ph, physics:nucl-ex, physics:nucl-th] (2021).
  * Stone  _et al._ [2021] J. R. Stone, P. Danielewicz, and Y. Iwata (2021), submitted to Physics Letters B.
  * Borghini  _et al._ [2002] N. Borghini, P. M. Dinh, J.-Y. Ollitrault, A. M. Poskanzer, and S. A. Voloshin, Physical Review C 66, 014901 (2002).
  * Lisa  _et al._ [2000] M. A. Lisa, N. N. Ajitanand, J. M. Alexander, M. Anderson, D. Best, F. P. Brady, T. Case, W. Caskey, D. Cebra, J. L. Chance, P. Chung, B. Cole, K. Crowe, A. C. Das, J. E. Draper, M. L. Gilkes, S. Gushue, M. Heffner, A. S. Hirsch, E. L. Hjort, L. Huo, M. Justice, M. Kaplan, D. Keane, J. C. Kintner, J. Klay, D. Krofcheck, R. A. Lacey, J. Lauret, H. Liu, Y. M. Liu, R. McGrath, Z. Milosevich, G. Odyniec, D. L. Olson, S. Y. Panitkin, C. Pinkenburg, N. T. Porile, G. Rai, H. G. Ritter, J. L. Romero, R. Scharenberg, L. Schroeder, B. Srivastava, N. T. B. Stone, T. J. M. Symons, R. Wells, J. Whitfield, T. Wienold, R. Witt, L. Wood, and W. N. Zhang, Physics Letters B 496, 1 (2000).