\section{Introduction}
In recent years, the dynamics of heavy ion collisions at intermediate
energy has been extensively investigated within the framework of transport theories, such as the Nordheim approach, in which the Vlasov equation for the
one-body phase space density, $f({\bf r}, {\bf p},t)$, is supplemented with a Pauli-blocked
Boltzmann collision term \cite{Bertsch_rep,Aldo_rep}.
The basic ingredients that enter the resulting transport
equation, often called  Boltzmann-Uehling-Uhlenbeck (BUU) equation,
are
the self-consistent mean-field potential and the two-body scattering cross sections.
These transport models hence describe the time evolution of the reduced
one-body density in phase-space and,
consequently, they are suited for
the description of one-body observables, such as inclusive particle spectra
in nuclear collisions, average collective flows and excitations.
However, they cannot provide a reliable description of fluctuation phenomena,
such as multi-fragmentation processes, i.e. the break-up of excited nuclear
systems into many pieces. In fact,
neither fluctuations of one-body observables nor many-body correlations can be
addressed with this class of mean-field models.
Hence suitable extensions, including fluctuations of the one-body density, have to
be considered.

An intense theoretical work
on fluctuations in nuclear dynamics has started in the past years,
also stimulated by the availability of
large amounts of experimental data on fragment formation in intermediate
energy heavy ion collisions and the possibility to observe a liquid-gas
phase transition \cite{Frankland,Moretto,Dagostino,Chomaz,EPJA,rep,Ono}.
In order to introduce fluctuations in transport theories,
a number of different avenues have been taken, that can
be essentially reconducted to two different classes of models.
One is the class of molecular dynamics models \cite{Ono,md,feld,ONOab,Pap,FFMD} while the other
kind is represented by stochastic mean-field approaches
\cite{Ayik,Randrup,rep}.

In molecular dynamics models the many-body state is represented
by a simple product wave function, with or without antisymmetrization.
The single particle wave functions are assumed to have a fixed Gaussian shape.
In this way, though nucleon wave functions are supposed to be
independent (mean-field approximation), the use of localized wave packets
induces many-body correlations both in mean-field propagation
and hard two body scattering (collision integral), which
is treated stochastically.
Hence this way to introduce many-body correlations
and produce a trajectory branching is essentially based on the use of empirical gaussian
wave packets.
If wave functions were allowed to assume any shape, the method would
become identical to standard mean-field descriptions.
While the wave function localization appears appropriate
 to describe final fragmentation channels,
where each single particle wave function should be localized within a
fragment,
the use of fixed shape localized wave packets in the full dynamics
could affect the correct description of
one-body effects, such as spinodal instabilities and zero sound
propagation \cite{FFMD,Akira_comparison}.

On the other side, in the so-called stochastic mean-field approaches,
the stochastic extension of the transport treatment for the one-particle
density is obtained by introducing a stochastic term representing the
fluctuating part of the collision integral \cite{Ayik,Randrup}, in close analogy
with the Langevin equation for a brownian motion.
This can be derived as the next-order correction, in the equation describing
the time evolution of $f$, with respect to the standard
average collision integral, leading to the Boltzmann-Langevin (BL) equation.
Thus, the system is still described solely in terms of the reduced one-body density
$f$, but this function experiences a stochastic time evolution in response
to the random effect of the fluctuating collision term.
In this way density fluctuations are introduced, that are amplified when
instabilities or bifurcations occur in the dynamics.
This procedure is suitable also for addressing multifragmentation phenomena,
since fragments can be associated with the regions where the spacial density
becomes larger, which finally can be reconstructed by sampling the
one-body distribution function.

A specific method for solving the Boltzmann-Langevin equation by direct
numerical simulation was introduced in Refs.\cite{fluct_Chomaz,fluct_burgio}.
In this numerical implementation, the one-particle density
$f({\bf r}, {\bf p})$ is
represented on a lattice of grid points in phase space and the collision integral is treated
by considering all possible transitions between phase space cells, adding a noise term whose
features are related to the average rate of transitions between two specified initial cells and
two final cells.
The numerical implementation of this method has only been possible in two dimensions (2D) because it
requires too large computer resources in 3D.

Hence several approximate solutions of the BL equations have been formulated,
mostly based on the projection of the BL noise only on a given dynamical variable
(such as the local quadrupole tensor of the momentum distribution)
or on  ${\bf r}$  space only \cite{Suraud,Salvo}.
More tractable fluctuating terms, such as a stochastic force added to the mean-field potential,
have also been proposed and extensively applied to
multifragmentation studies \cite{BOB,EPJA,rep}.

However, the implementation of the full structure in phase space of the original BL term
can still be considered as an important goal to reach.
In fact, this would allow to treat a more general class of phenomena, where
the correct description of fluctuations and correlations in  ${\bf p}$ space
is essential (such as particle production and fragment velocity correlations for instance).
Moreover, also in the multi-fragmentation mechanism, that is dominated
by spacial density fluctuations, a more accurate representation of the
full phase space dynamics, including fluctuations, would allow to improve the description of the fragment
kinematical properties and correlations.

We also stress the general interest of this effort. Indeed transport phenomena occur in many physical
systems, for which a more precise description of the time evolution of the one-body distribution
function, including the effect of many-body correlations,
would be important.

A first attempt to introduce a fluctuating collision term in a 3D transport approach
was made by Bauer et al. \cite{fluct_bauer}.
This method can be implemented relatively easily into standard transport codes that adopt  the
scattering of pseudo-particles (or test particles) as a
method of solution of the collision integral
and it consists essentially in forcing similar two-body collisions to occur
for neighboring test particles, defined according to a given distance in phase space,
so that effectively two nucleons are involved in each particular
collision event.
The distance should reproduce the phase-space shape of the nucleon wave function.
In this way the random nature of the two-body scattering, that in the standard
codes applies only to test particles and is washed-out when using a huge number of
them, is transferred to entire nucleons.
However, in the procedure proposed in Ref.\cite{fluct_bauer} the Pauli blocking is checked only for the collision of the
two original test particles and not for the entire swarm affected,
leading to some unpleasant features in the description of fermionic systems.
Indeed
the Pauli-blocking violation introduces important inaccuracies in
the fluctuations of  the one-body density.

In the present manuscript we present a new method to reconstruct the phase space nucleon wave function in mean-field approaches,
in such a way that the Pauli-blocking is checked for the entire cloud of moved particles.
This will improve the description of the
fluctuation variance,
approaching the one
expected for fermionic systems.
We will also pay special attention to the definition of the phase-space metric that would
optimize the value of this variance.

The paper is organized as follows:
We will first recall the main ingredients of the BL theory, in order to connect
the formalism with the numerical implementation adopted (Section 2).
Then we will discuss in more detail
the methods that have been proposed so far to solve the
Boltzmann-Langevin equation (Section 3).
The new procedure that we follow to build fluctuations is presented in
Section 4. Several results demonstrating and analyzing in detail the method
are discussed from Section 5 to 8.
Conclusions and perspectives are drawn in Section 9.

\section{The Boltzmann-Langevin equation}

Within the semi-classical framework,
the stochastic transport equation of motion
for the one-body distribution function $f$
can be expressed in the following form,
\begin{equation}\label{BLE}
\dot{f}\ \equiv \ {\partial\over\partial t}f\ -\ \{h[f],f\}\
=\ K[f]\ =\ \bar{K}[f]\ +\ \delta K[f]\ ,
\end{equation}
where the left side describes the collisionless propagation
of the individual particles in their common self-consistent one-body field,
while the right side expresses the effect of the residual binary collisions.

The collision term $K[f]$
has a stochastic character.
For example, the distance a particle travels in the medium
before colliding is stochastic,
as is the resulting scattering angle.
At the Boltzmann level,
Eq.(\ref{BLE}) includes only
the {\em average} part of the collision term, $\bar{K}[f]$.
Usually, in nuclear systems,
the quantum statistics is taken into account
by adding suitable Fermi blocking
in the single-particle final states, leading to the Uehling-Uhlenbeck
collision term \cite{Nordheim1928,UUPR43}.
The resulting nuclear Boltzmann equation
exists in many implementations that differ with respect to
both the physical input
(such as the types of constituents included,
the form of their effective Hamiltonian,
and their differential interaction cross sections)
and the numerical methods employed
(whether of test particle or lattice type)
and various names have been employed in the literature, including
BUU (for Boltzmann-Uehling-Uhlenbeck),
VUU (for Vlasov-Uehling-Uhlenbeck),
Landau-Vlasov and BNV (for Boltzmann-Nordheim-Vlasov) \cite{Bertsch_rep,Aldo_rep,Bao_an,seb,Bao_IBUU}.
The  Boltzmann-Langevin treatment
includes also the {\em stochastic} part of the collision term,
$\delta K[f]
\cite{Ayik,Randrup}.

In the simple physical
scenario where the residual interaction can be considered as binary
collisions that are well localized in space and time, the average part is
given by
\cite{Nordheim1928,UUPR43}:
\begin{equation}
\bar{K}({\bf r},{\bf p}_1)\ =\ {g}\sum_{234}W(12;34)\left[ \bar{f}%
_{1}\bar{f}_{2}f_{3}f_{4}-f_{1}f_{2}\bar{f}_{3}\bar{f}_{4}\right] \ ,
\end{equation}
where $f_{i}$ is a short-hand notation for $f({\bf r},{\bf p}_i,t)$
and $\bar{f}\equiv 1-f$ is the associated Fermi blocking factor.
$g$ is the degeneracy factor.

The basic transition rate is simply related to the differential cross section
for the corresponding two-body scattering process,
\begin{equation}
W(12;34)\ =\ v_{12}\left({d\sigma \over d\Omega }\right)_{12\to34}
\delta({\bf p}_1+{\bf p}_2-{\bf p}_3-{\bf p}_4)\ ,
\end{equation}
being $v_{12}\equiv |\bold{v}_1-\bold{v}_2|$ the relative velocity,
and it thus has corresponding symmetry properties,
\begin{equation}
W(12;34)\ =\ W(21;34)\ =\ W(34;12)\ .  \label{Wsym}
\end{equation}

Since it arises from the same elementary two-body processes, the stochastic
part of the collision term is fully determined by the basic transition rate
as well, as a manifestation of the fluctuation-dissipation theorem.
With the collisions assumed to be local in space and time,
the correlation function for the fluctuating part of the collision term
 has the following form,
\begin{equation}\label{dKdK}
\prec \delta K({\bf r},{\bf p_1},t)\ \delta K({\bf r}',{\bf p_1}',t')\succ =
C({\bf p_1},{\bf p_1}',{\bf r},t)\ \delta ({\bf r}-{\bf r}')\ \delta (t-t')\ ,
\end{equation}
where $\prec\cdot\succ$ denotes the average with respect to
the ensemble of possible trajectories resulting from
the current one-body density $f$.
Furthermore, for elastic scattering,
the correlation kernel is given by \cite{Ayik}:
\begin{equation}
C({\bf p_1},{\bf p_1}',{\bf r},t)=\delta
_{11'}\sum_{234}W(12;34)F(12;34)
\nonumber
\end{equation}
\begin{equation}
+ \sum_{34}\left[ W(11';34)F(11';34)-2W(13;1'4)F(13;1'4)\right] \ ,
\end{equation}

with the short-hand notations $\delta_{11'}\equiv h^D\delta({\bf p}_1-{\bf p}_1')$
and $F(12;34)\equiv f_{1}f_{2}\bar{f}_{3}\bar{f}_{4}+%
\bar{f}_{1}\bar{f}_{2}f_{3}f_{4}$
($D$ is the dimension of the space considered).
The symmetry properties (\ref{Wsym}) of the transition rate
ensure that the following sum rules hold,
\begin{eqnarray}
\sum_1 C({\bf p}_1,{\bf p}_2,{\bf r},t)\phantom{{\bf p}_1} =\ \sum_2 C({\bf p}_1,{\bf p}_2,{\bf r},t)
\phantom{{\bf p}_1} =\ 0\ , \\
\sum_1 C({\bf p}_1,{\bf p}_2,{\bf r},t){\bf p}_1 =\ \sum_2 C({\bf p}_1,{\bf p}_2,{\bf r},t){\bf p}_2
=\ \bold{0}\ ,\\
\sum_1 C({\bf p}_1,{\bf p}_2,{\bf r},t)\epsilon _{1}\, =\ \sum_2 C({\bf p}_1,{\bf p}_2,{\bf r},t)
\epsilon_2\,\, =\ 0\ ,
\end{eqnarray}
where $\epsilon_i=p_i^2/2m$ is the kinetic energy for a specified momentum.
These sum rules express the fact that each of the elementary binary
collisions conserves particle number, momentum, and energy, respectively.

\section{Methods to solve the BL equation}

\subsection{Lattice calculations}
The numerical implementation of the Boltzmann-Langevin equation
 is accomplished by
correctly simulating the basic stochastic process, i.e. the
stochastic transition
rate among phase-space cells.
One can realize that this point
is very delicate in finite systems, when we work with a relatively small
number of nucleons, but still need to build a smooth distribution function.
There are different ways to overcome this difficulty.

As mentioned in the Introduction, one possibility,
proposed in Ref.\cite{fluct_Chomaz}, is to solve the BL equation on a lattice;
phase space is therefore divided into a number of cells, each one having
volume $\Delta s = \Delta \mathbf{r} \Delta \mathbf{p}$
($s$ denotes a point in phase space: $s \equiv ({\bf r},{\bf p})$).

Each transition
involves four locations, and the collision integral arises from the sum
of the transitions evaluated for all possible combinations of the cells:
from this consideration one recognizes that the practicality of the method is
limited by the huge computational effort required.
In fact, only two-dimensional implementations exist \cite{fluct_Chomaz,fluct_burgio}.
Each transition
represents a basic stochastic process, and, following the BL theory,
the actual number of such transitions in a time step $\Delta t$
is dispersed around the mean value
according to a Poisson distribution, so its variance amounts to:
\begin{equation}\label{sigma_ni}
\sigma^2_{\nu}=\bar{\nu},
\end{equation}
where the average is given by:
$$
\bar{\nu}(12;34)=\frac{\Delta s_1}{h^D}\frac{\Delta s_2}{h^D}
\frac{\Delta s_3}{h^D}\frac{\Delta s_4}{h^D}f_1 f_2\bar{f}_{3}\bar{f}_{4}
W(12;34)
$$
\begin{equation}
\cdot \delta(\mathbf{r}_{1}-\mathbf{r}_{2})
\delta(\mathbf{r}_{1}-\mathbf{r}_{3)}\delta(\mathbf{r}_{2}-\mathbf{r}_{4)}\Delta t.
\end{equation}\label{ni_med}
\\
This is a compact form of the fluctuation-dissipation theorem, since
it predicts that the fluctuations are simply related to the mean number
of dissipative processes.
Then, a noise term $\delta\nu(12;34)$ is added to the mean
number $\bar{\nu}(12;34)$, and thus the actual number of transitions
is a random number picked from a normal distribution, the center
of which is the mean value, and where the width is given by Eq.(\ref{sigma_ni}).
In refs. \cite{fluct_burgio,chapelle} this procedure is proven to yield
correct results for a 2D system of fermions interacting through hard two-body scattering.
Starting from a non-equilibrium situation (two touching Fermi spheres)
the method
is successful in reproducing the expected
fluctuations, preserving also the average trajectory of the system.
In particular,
 the fluctuations introduced build the expected
statistical value for
the equilibrium one-body density variance:
\begin{equation}
\sigma^2_f(s) = \prec \delta f(s)\delta f(s) \succ = f_{eq}(s)
\bar{f}_{eq}(s) h^D/(g\Delta s),
\label{variance}
\end{equation}
evaluated considering phase-space cells of volume $\Delta s = h^D$, for which
the fluctuating transition rate is implemented.
Here $f_{eq}(s)$ denotes the equilibrium value of the one-body distribution
function and $\delta f(s) = f(s) - \prec f(s) \succ $.
Of course, Eq.(\ref{variance}) holds also for volumes larger than $h^D$
and fluctuations are scaled accordingly.
Also the co-variance, i.e. the correlation between density
fluctuations in two different phase-space points, $s$ and $s'$,
is well reproduced \cite{fluct_burgio}.

The same authors of Ref.\cite{fluct_burgio}
also verified that, for the same idealized system,
but prepared at low spatial density,
early fluctuations developed in momentum space are subsequently
transmitted into density fluctuations and amplified by the nuclear mean field
\cite{fluct_burgio1},
leading to large instabilities and a statistical population of fragments.

\subsection{The pseudo-particle correlation method}

The method of Ref.\cite{fluct_Chomaz} was originally developed as a way to overcome
the problems arising from the solution of transport equations with
the test particle method, which we now discuss briefly.
Within the test particle method, each nucleon is represented by a
collection of $N_{test}$ test particles, that are propagated according
to the mean-field interaction and random two-body scattering.
The use of a large number of test particles allows
to have a smooth distribution function and a good
coverage of the phase space for the Pauli blocking.
On the other side, since collisions are treated stochastically for the
single test particles, the dispersion around the average number of nucleon
collisions (Eq.(\ref{sigma_ni})) is automatically divided by $N_{test}$.
Hence fluctuations
which are introduced by the test particle algorithm can also be seen as the
correct ones expected from the BL approach, divided by the factor $N_{test}$.

Various paths have been followed to overcome this problem. For instance,
one can choose $N_{test}=1$, i.e. work with whole nucleons
\cite{md}, but in this case one has to solve the problem of the
induced numerical errors on the smooth path of the dynamics,
i.e. the mean-field and the Pauli-blocking factors. Alternatively, Bauer et al.
proposed a method to introduce a correlation between close particles in
phase-space \cite{fluct_bauer}.

The method follows the idea, applied in extended Time-Dependent-Hartree-Fock
(TDHF) calculations,
 of evolving the single-particle density
including a statistical mixing of Slater determinants \cite{tdhf}.
The jump from
one configuration to the other is possible under the important assumption
that any coherence between determinants can be ignored
(decorrelation approximation).  The choice of the single-particle basis
for the determinant is however somewhat arbitrary, apart from the requirement
of momentum-energy conservation. These ideas can be
translated into the semi-classical approximation using the test particle
method. The mixing between Slater determinants is realized in the collision
integral by means of the following procedure.
\begin{enumerate}
\item First of all, the nucleon-nucleon (NN) cross
section is reduced by a factor $N_{test}$.

\item Then two test particles $i_1$ and $i_2$ are chosen as
colliding partners, and will be moved from  their positions ${\bf p}_1$
and ${\bf p}_2$, to new positions ${\bf p}_3$
and ${\bf p}_4$,
according to the corresponding transition probability
including the Pauli blocking of the final states
\cite{Aldo_rep}.

\item If the two particles can collide, the
scattering actually involves two ``clouds'' of neighbouring test particles,
corresponding to two entire nucleons ($2\times N_{test}$ particles).
The contiguity criterion is based on the following definition of
phase-space distance:
\begin{equation}\label{distance1}
d_{ij}^2=(\mathbf{p}_i-\mathbf{p}_j)^2/p_F^2 +
(\mathbf{r}_i-\mathbf{r}_j)^2/R^2
\end{equation}
where $p_F$ is the Fermi momentum and $R$ is the radius of the
considered system.
(It should be noticed that  the choice of this phase-space metric
is rather arbitrary).

\item In order to ensure energy and momentum conservation, the final states
are adjusted
using the average momenta of the two clouds as initial states.
\end{enumerate}

Two important remarks have to be done:
\begin{itemize}
\item the Pauli blocking is checked only for the collision of the test
particles $i_1$ and $i_2$ and not for each particle of the cloud.
We will see in the following that this
choice induces a strong violation of the Pauli principle;
\item the effect of the collision is a mere translation (in momentum space)
of the two clouds to final positions, as sketched in Fig. \ref{trasla}.
\end{itemize}
\begin{figure}[htb]
\centering
\includegraphics[scale=0.5]{fig1_flu.eps}
\caption{Illustrative picture of nucleon-nucleon scattering in the
pseudo-particle correlation method.}
\label{trasla}
\end{figure}

This method is conceptually simple, and moreover the implementation into
existing BUU or BNV codes is really straightforward. Indeed it
was immediately used to test the effects of the fluctuations on fragment
mass spectra in a $^{20}Ne+ \,^{20}Ne$ reaction
at $100\,AMeV$ \cite{fluct_bauer}.

Its validity has been checked later by Chapelle et al. \cite{chapelle},
for the same idealized 2D system for which
the BL lattice simulation method was implemented.
Accordingly, a different metric,
involving only distances in momentum space,
was  used.
As shown in \cite{chapelle},
the final equilibrium value obtained for the fluctuation
variance does not reproduce the expected one.
The variance profile, as a function of the energy,  appears to be proportional to $f_{eq}$ rather than
to $f_{eq}\bar{f}_{eq}$.
Actually, this would be the result for a classical
system, obeying to the Boltzmann statistics.
This can be considered as a consequence of the fact that Pauli
blocking is checked only for the first two colliding partners, and the
system progressively evolves towards an equilibrium configuration which
is more consistent with the one of a classical system.

From this discussion it is evident that a proper treatment of the Pauli
blocking is a fundamental requirement for any numerical implementation
 of the  Boltzmann-Langevin theory for fermionic systems.
However, apart from these problems, the pseudo-particle correlation method
represents a simple and practical way to implement fluctuations into
transport codes.
Therefore, it is worthwhile to think about some improvements to
make it more accurate.
This is the aim of the present work, as we will discuss in the
next Section.

\section{The improved pseudo-particle correlation method}\label{improved}

The method devised by Bauer et al. is able to agitate the phase space
function and
build fluctuations, although their strength does not reproduce the
expected value for fermionic systems.
However, a careful modification of the original
procedure can considerably improve the results.

The new procedure can be summarized in the following steps:

\begin{enumerate}

\item The choice of the two colliding partners closely traces the
standard recipe. If the two test particles are allowed to collide,
two clouds of $N_{test}$ particles will be moved,
with conditions specified below.

\item
 Only particles within a sphere, in coordinate
space,
around the
center of mass of the two partners $i_1$ and $i_2$ can belong to the
clouds. The distance criterion is
$|\mathbf{r}_i-\mathbf{r}_{CM}(i_1,i_2)|<d_r$, where $d_r$ is a free
parameter (see later discussion).
Then, for
this considered space sphere,
 a grid is introduced in momentum space,
around $i_1$ and $i_2$,
the size of each cell being $V_{cell}$.

\item Given the momentum space cells $I$ and $J$ containing the partners $i_1$ and $i_2$,
we consider the cells $I'$ and $J'$, corresponding to the final
positions,  in a rotated frame,
as indicated in Fig. \ref{ruota}. For a given set of initial and final cells,
the number of test particles that will be actually moved to final states is
equal to the minimum between the occupation of the initial cells, $n$, and
the availability of the final ones ${\bar n}' = (N_{max}-n')$:
$$
n_t(I,I';J,J')= \min(n_I,n_J,\bar{n}_{I'},\bar{n}_{J'}),
$$
where $N_{max} = g V_r V_{cell}/h^D $ is the maximum number of test particles that can stay
inside a cell.
It should be noticed that this choice corresponds to the maximum of possible moves.

\item
Surrounding momentum space cells
are searched with the same prescriptions, until two entire nucleons
are found. The search procedure is symmetric with respect to the
center of momentum of the two partners
and random in the vicinity of the original cells (I,J).

\item Finally, the two clouds are moved to the new states of the rotated frame
(see Fig.2),  so that
the conservation of energy and momentum is automatically obtained.
However, a further check is performed.
In fact, due to the finite number of test particles, the one-body density is
not perfectly homogeneous inside each cell.
This causes slight violations
of the conservations laws.
Hence the origin of the momentum space grid is eventually slightly displaced in order to
have a perfect energy and momentum conservation.

\end{enumerate}
This method involves two parameters, namely the radius of the sphere in
$\mathbf{r}$ space, $d_r$, and the size of the momentum cells, $V_{cell}$.
\begin{figure}[htb]
\centering
\includegraphics[scale=0.35]{fig2_flu.eps}
\caption{Schematic picture of the collisional procedure in the improved pseudo-particle
correlation method.}
\label{ruota}
\end{figure}
The radius $d_r$ fixes the spatial extension of the nucleon, that in turn
influences its spreading in momentum space.
Hence this parameter fixes somehow the extension of the nucleon wave packet
in phase-space and could affect the transport
dynamics.
It can be constrained by physical
arguments, and in general it should depend on the physical properties of the
system, such as its dilution.

\subsection{General comments}\label{comments}

The method proposed here
has several advantages with respect to the original one of Ref.
\cite{fluct_bauer}. It is
still very flexible, and can be easily implemented into existing
transport codes;
moreover, it is quite fast, and does not require huge amounts
of resources.
The Pauli blocking
is locally checked in each collision, so in principle this procedure
satisfies the basic requirements to preserve the average evolution and
to yield the expected fluctuations. Besides, differently from the original
method, where the nucleon is essentially spherical in phase space,
here it is allowed
to assume any shape in
momentum space, since the only restriction is given by the cell
size $V_{cell}$. In
other words, it is possible to take into account momentum
deformations of the nucleon
wave packet,
due to the momentum distribution of the system in particular conditions.
For instance, when equilibrium
is approached,
the Pauli principle allows collisions essentially for nucleons
lying on the surface of the Fermi sphere;
therefore
such nucleons should not have a spherical shape in momentum space.

It is also interesting
to make a connection between the proposed method and the
lattice simulation method to build fluctuations.

The numerical procedure described in Ref. \cite{fluct_Chomaz},
where the phase space is discretized by a fixed grid of
cells,
allows to construct the correct value of the fluctuation amplitude
in each
cell of the grid, generally taken of volume $h^D$
(but in principle it
can be even smaller), on the basis of the fluctuating gain and
loss terms.

On the other hand, our procedure is based on a different
philosophy: Instead of constructing fluctuations in fixed cells,
when a collision occurs
we move all together the phase space volume containing
one nucleon.  The shape of this
volume is not fixed a priori and depends on the location of the four
points (${\mathbf p}_1,{\mathbf p}_2,{\mathbf p}_3,{\mathbf p}_4$)
that have been chosen as initial and final centroids of the considered
two-body collision. This means that the ``correlation volume'',
i.e. the volume of the sphere that envelops all test particles that move all
together in each collision, is not fixed, but has to be considered as a dynamical
variable that generally exceeds $h^D$ and
could even reach rather large values.
In other words,
compared to what is done in Ref. \cite{fluct_Chomaz}, we build
fluctuations in volumes containing one nucleon, but these volumes
do not correspond to single cells of a fixed grid: instead, their contour changes
from one collision to another, according to the prescription adopted to
construct the clouds around the two first colliding test particles
$i_1$ and $i_2$.

But how can the two proposed methods be connected to the original BL equation ?
It should be noticed that the Boltzmann-Langevin equation (\ref{BLE}) is a semi-classical equation,
derived in the ``continuum limit'' ($h \rightarrow 0$).
This implies that the evolution
of the system under consideration is completely determined by a smooth distribution function,
i.e. it behaves as a fluid in phase space.
In other words, the volume $\Delta s$ of the phase space cells for which
the fluctuating transition rate is evaluated, is supposed to be much larger
than $h^D$.
The Boltzmann-Langevin equation is
essentially based on the idea of a fluctuating collision rate among
phase space cells
containing many nucleons. The formalism does
not provide any additional information about the ``structure'', i.e.
the shape of the single nucleon wave packet.

In both numerical implementations described above,
two-body collisions among nucleons are treated as a stochastic process.
Only the definition of nucleon wave packet and/or of the elementary cells where
fluctuations are built change from one method to the other.
As discussed above, this
ingredient is not contained in the BL equation and such information
is beyond its derivation.
Hence, both procedures can be considered as correct
implementations of the Boltzmann-Langevin theory in the ``continuum limit'',
the definition of ``nucleon wave packet'' being rather arbitrary.

However,
it is also very interesting to investigate the differences between
the
various possible methods used to build the nucleon clouds,
that in turn influence the correlation volume.
In particular, one can also
try to force the
procedure in order to reduce the wave packet smearing and to correlate
particles inside smaller volumes (close to $h^D$)

\subsection{Details of the model}
In the present first implementation of the method,
similarly to what was done for the lattice simulations of Ref.\cite{fluct_burgio},
we neglect
the evolution in coordinate space,
extending $d_r$ to the whole coordinate space.
The cell size $V_{cell}$ is then constrained by the following arguments: it should be
small enough to allow an accurate check of the Pauli blocking, but large
enough to contain a sufficient number of test particles to reduce numerical
uncertainties.
In any case, $V_{cell}$  has not to exceed the volume where, at
most, one nucleon can be accommodated:
$V_p =
h^D/(g V_r)$,
being $V_r$
the volume
in the coordinate space ($V_r =  \frac{4}{3} \pi d_r^3$).
In the following we take $g = 1$, unless expressly indicated.
Since we neglect the coordinate space,
we will refer only to volumes $V$ in momentum space.
The expected equilibrium
value of the variance $\sigma^2_f$ in such volumes is:
\begin{equation}
\label{scaling}
\sigma^2_f= f_{eq}\bar{f}_{eq}\frac{V_p}{V}=\frac{f_{eq}\bar{f_{eq}}}{N_{V}}
\end{equation}
where $N_{V}$ is the number of nucleons that can be accommodated, at most,
in the considered volume $V$.
It should be noticed that,
according to our procedure to build the nucleon wave packet,
the expected value of fluctuations, (Eq.\ref{scaling}),
can be reproduced only in volumes $V$ larger than
the average
correlation volume,
while it is underestimated
in smaller cells,
since fluctuations are built on a larger scale.

\subsection{Illustrative results}
In order to illustrate how our  procedure works,
we will first consider the very simple case
of particles moving only along one direction.
We take a large system, containing 1000 nucleons with average occupancy
$\prec f \succ = 0.5$.  So the total extension in momentum space
of the system is $2000\times V_{p}$.
Then the system is divided into $N = 4000$ smaller cells that can contain at most
half nucleon, whose occupancy can be either $0$ or $1$. Hence $V_{cell} = V_p/2$.
The occupancy is randomly chosen at the beginning and collisions are performed
until equilibrium has been reached.
The initial ($p_1$,$p_2$) and final  ($p_3$,$p_4$)
states of a collision
are chosen randomly among the $N$ cells. We neglect energy and momentum conservation in this
simple example.
 Moreover, for simplicity, we take ${p}_2 = -{p}_1 $ and  ${p}_4 = -{p}_3$.
Hence only half of the considered space is independent and, in this situation,
the probability for a collision to happen is proportional to the product $f({p}_1)(1-f({p}_3))$.
According to Eq. (\ref{scaling}), the expected fluctuation value, in a given
momentum space volume $V$, will be
$\sigma^2_0 = 0.25~V_{p}/V$.
The general prescription (see Sect. \ref{improved}) states that when a
collision happens, a full nucleon has to collide.
Therefore, we move together the cell $i$ and the closest occupied
cell that finds an empty cell in the final position.
This example can be seen as a simple implementation of the general
procedure when considering nucleons represented by only two test
particles.
One can easily realize that the  ``correlation radius'',
i.e. the distance between the two cells that move together,
changes from one
event (collision) to another: sometimes the two co-moving cells are
neighbours, sometimes they can be rather far.
Of course,
for practical reasons, we have to restrict our search within a given radius.
For instance, considering 10 cells on each side (with
respect to $i$), we are able to reconstruct the nucleon and to
perform the collision in more than 98$\%$ of the cases.

\begin{figure}[htb]
\vspace{0.5cm}
\centering
\includegraphics[scale=0.25]{fig3_flu.eps}
\caption{Variance of the distribution function
as a function of the number $N_V$
(see text),
in the idealized case of particles moving only along one direction.}
\label{flu_joseph}
\end{figure}
In Fig. \ref{flu_joseph} we present the results for the
fluctuation variance
as a function of the number $N_V$ of ``nucleon'' cells
(of volume $V_{p}$) that are contained inside the volume $V$ where the function $f$ is evaluated.
One thousand events have been considered.
The fluctuation value presented in the Figure is multiplied by $N_V$.

As one can see in the Figure,
the correct value of the fluctuations, i.e. $N_V\sigma^2_0 = 0.25$, is
approximately recovered (within $\sim5\%$) for volumes
$V \approx 20~V_{p}$.
The slight deviation is due to the upper bound
considered for the  ``correlation radius'' and to the constraints
 imposed by the total mass conservation, that are also responsible for the
decrease of fluctuations observed
for large $N_V$ values.

From this simple example, it already appears that
the exact fluctuation value $(N_V\sigma^2_0 = 0.25)$
is reached only asymptotically.
This is due to the arbitrariness of the nucleon shape and of the
centroid position.
Hence, when testing a given procedure to solve the stochastic
collision integral, the resulting
fluctuations have to be
compared to the expected values
in large volumes, in which all different possibilities
(nucleon configurations) may be accommodated.

These considerations would apply also to the original procedure
proposed by Bauer et al.  Thus, in Ref.\cite{chapelle}, the comparison
between the obtained fluctuation variances and the expected ones
should have been done considering also volumes containing many nucleons,
and not only cells of volume $V_p$
where, as seen in Fig.\ref{flu_joseph}, fluctuations are
naturally suppressed by the fact that the nucleon wave packet may
take any shape.
However, this does not affect the conclusions drawn in \cite{chapelle},
mainly based on the wrong profile of the density variance
(and not on its absolute value).

In the next Section we will move to describe more complex geometries
and different procedures to build the nucleon wave packets.

\section{Fluctuations of a Fermi gas: Analysis on the Fermi surface}

Here we will discuss our procedure to build fluctuations in the context of
a system of fermions having
a given density and temperature, i.e. initialized according to a Fermi-Dirac
distribution. This is a situation that is easily reached in the course
of a dissipative nuclear reaction, after the initial collisional
shock \cite{rep}.
Since collisions happen mostly on the Fermi surface, let us first
consider only particles within a thin stripe $dp$ around the
Fermi momentum, where
$\prec f \succ = 0.5$.

For those colliding particles
the relative velocity $v_{12}$ is constant,
since only particles with opposite
momenta can collide, in order to stay on the Fermi surface. So
the transition rate is essentially governed by
the product
$f_1f_2\bar{f}_3\bar{f}_4$.

This simplified system will allow us to
study the sensitivity of the results to the
ingredients of the procedure used to build, step by step,
the nucleon clouds.
As also discussed before,
we will show that
the correct
fluctuation amplitude is recovered in the ``continuum limit'', i.e. in large
volumes.
However, in the following
we will mostly concentrate
on the fluctuation amplitude in volumes of
size $V_{p}$, in order to probe
the extension
of the nucleon wave packet.
The latter could influence significantly
the evolution of
actual nuclear collisions, where it could be important to work with
more compact nucleon configurations.

Using this simple model we will find that building the maximum fluctuation
value (e.g. $\sigma^2_0=0.25$) in such a volume, i.e. correlating
particles inside $V_p$,
is a quite difficult task and additional efforts are required in order to
reduce  the nucleon cloud smearing.

\subsection{Details of the simplified model}
\label{details}
In order to keep the particles on the Fermi surface,
the original procedure described in Sect. \ref{improved} was slightly
modified so that the modulus of the momentum of the particles does not
change: only particles belonging to cells with opposite momenta can collide.
Then, for symmetry reasons it will result
$f(\mathbf{p}_1)=f(\mathbf{p}_2 = \mathbf{-p}_1)$.
The
transition probability finally depends only on two values of the distribution
function, $f(\mathbf{p}_1)$ and  $f(\mathbf{p}_3)$.
We use
spherical coordinates $\big(\cos (\theta),\phi \big)$
as independent variables. The collision
partners of particles in the cell $\big( \cos (\theta),\phi\big)$
will be the ones in the cell $\big( \cos (\pi-\theta),\pi+\phi\big)$.
The choice of these coordinates implies the grids used in the procedure to be
fixed along the $\cos(\theta)$ axis. They can slide along the $\phi$ axis,
where periodic
conditions
are taken into account.

The procedure is
as follows: a test particle $i_1$ is chosen at random
(together with the corresponding partner $i_2$) and a grid is constructed,
centered around $i_1$ and $i_2$ (see Fig.2).
The rotation angles, also taken from
a flat distribution, fix the final states and the new frame
corresponding to these final positions.

We consider a shell with $40\times40$ cells, whose size
is chosen to accommodate one nucleon at most,
so that
$V_{cell}
=V_p$.
We have checked that the results do not depend on the number of cells,
i.e. the number of nucleons considered,
by performing the same calculations on a $64\times64$ system.
We use $500$ test particles per nucleon.

The growth of fluctuations is described by the variance $\sigma^2_f$
of the distribution function, as discussed before.
Usually it is evaluated performing an ensemble average over a large
number of events having the same initialization. However, in this scheme
the mean value $\prec f\succ$ does not depend on the position ${\bf p}$ and
is not changing in time.
Considering the large number of cells employed, we can also
follow the evolution of the system by calculating in a single event
the following quantity:
\begin{equation}\label{sigmaf}
\sigma^2_f=\langle(f_i-\langle f\rangle)^2\rangle _{cells},
\end{equation}
where $\langle\cdot\rangle $ denotes the average with respect to
the cells and $\langle f\rangle = 0.5$.

\subsection{``Fixed grid'' calculations}
\label{fixed}

In order to understand how the recognition of the nucleon wave packet can be influenced by
the details of our procedure,
we first further simplify the problem
by using a single
fixed grid to make the collisions, i.e. test particles are placed at (and can only move to) the center
of the cells of the grid.
In this case we have a single reference
frame, so the effect of a collision (in the space
$\big(\cos(\theta),\phi\big)$ ) is a mere translation to other states of the same grid.
Apart from the use of test particles, and related differences in the method
followed to solve the collision integral,
in this way the procedure becomes
similar to the one used in Ref. \cite{fluct_Chomaz}.

We will take,
as initial conditions, $f$ strictly equal to $0.5$
in each cell of
the fixed grid,
corresponding to $800$ nucleons.
By construction,
once collisions are performed, the only possible values for $f$ are $0$,
$0.5$ and $1$.

The maximum
fluctuations in $V_p$ ($\sigma^2_0=0.25$) would be reached only if, in each cell, $f$
is either zero or one.
Achieving this result is not  a trivial task:
Cells starting with $f=1$ can be
partially emptied, and cells with $f=0$ partially filled.

The dot-dashed
curve in Fig.\ref{full_esa}, upper panel, shows the variance
of the distribution function
calculated
in the cells of the fixed grid, with increasing number of collisions
(i.e. as a function of time).
It saturates
approximately when almost all nucleons have suffered one collision.
In the equilibrium configuration the number of
cells with the three possible values of $f$ is roughly the same, as shown
 in Fig. \ref{full_esa} (lower panel), so the fluctuation variance is
$\sigma^2_f\approx 0.175$, which is about $2/3$ of the maximum.
\begin{figure}[htb]
\begin{center}
\includegraphics[scale=0.31]{fig4a_flu.eps}
\end{center}
\vspace{0.5cm}
\begin{center}
\includegraphics[scale=0.31]{fig4b_flu.eps}
\end{center}
\caption{Upper panel: variance of the distribution function in $V_p$, as a function of the
number of collisions,  for three
different cases: fixed grid with $f(t=0) = 0.5$ (dot-dashed), fixed grid with randomly
distributed test particles (full), moving grid (dashed). Lower panel:
Distribution of the number of cells as a function of their occupation $f$ at equilibrium,
for the three cases indicated above.}
\label{full_esa}
\end{figure}

The fluctuations reached in $V_p$
also depend on the variance of the initial configuration.
If we initialize the system distributing the test particles
randomly, the occupation in the cells has an average value
$\langle f\rangle=0.5$ and a width given by
$\sigma_{f}(t=0)/ \langle f \rangle=\bar{N}^{-1/2}$,
where $\bar{N}$ is the
average number of test particles per cell (${\bar N} = 250$).
The distribution function can
assume any value between $0$ and $1$,
as shown by the full curve in Fig. \ref{full_esa}, lower panel.
The initial continuous distribution of values of $f$ causes a reduction
of the fluctuations in $V_p$
(full curve in Fig.~ (\ref{full_esa}), upper panel) with respect to the
previous case.
As expected, the fact that $f$ varies continuously between $0$ and $1$
contributes to
enhance the smearing of the nucleon wave packet and hence to reduce
fluctuations in $V_p$.

\subsection{The centroid degree of freedom}
\label{centroid}
Now we turn to a more general situation, relaxing the condition
of the fixed grid.
Hence here we would like to test the full procedure, described in Section 4,  that we will finally
follow in the full 3D case.
The two reference frames used to make
the collision (the one concerning the initial states and the one for the
final states)
are now shifted along the $\phi$ axis, since the angle $\phi$ of the initial
and final states can assume any value between $0$ and $2\pi$,
according to the initial and final positions of the
colliding test particles, $i_1$ and $i_2$.

This carries two
important consequences. First of all, even being the Pauli blocking
checked for each collision in a particular frame,
the prescription $f\leq 1$ can be violated in another frame. However, we
checked that this violation concerns only about $5\%$ of the cells. Second,
the correlation between the cells representing the reconstructed nucleon,
which is maximum in the search frame, will be systematically broken
in another (shifted) reference frame,
increasing the  nucleon delocalization
and thus reducing fluctuations in $V_p$.

The dashed curve in
Fig. \ref{full_esa} shows the fluctuations
obtained in $V_p$
in this case.
Since this feature is due to the arbitrary
position of the centroid of the constructed nucleon with respect to a fixed
reference frame, we will name it \emph{centroid effect}.
We notice
a fluctuation reduction of $\sim 60\%$
with respect to the fixed grid case (dash-dotted line in the same
Figure).

\subsection{Sensitivity to the dimension of the grid cell}

Up to now, the step of the grid has been chosen
so that $V_{cell} = V_p$.
However, in general,
the size of the search cell
 has to be smaller than $V_p$, especially when the
nucleon wave packet is rather delocalized in ${\mathbf p}$ space (being more
compact in coordinate space).
In fact, small search cells are required for the
Pauli blocking to be accurate.
Hence we want
to study the sensitivity of the results to
the size of the search cell.
Let us start with the fixed grid scheme.
In Fig. \ref{smallphi_full} we plot the
fluctuation variance obtained
dividing the step in the $\phi$ direction, $\phi_{step}$,
by two or four;
the volume of the search cell is then
$V_{cell} = V_p/2, V_p/4$, respectively.
This is equivalent to shift, for each collision,
the origin of the grid by $0.5$ or $0.25$ $\phi_{step}$
respectively, and acts also
as a sort of centroid effect, thus reducing
the variance in $V_p$ (see solid and
dashed lines in Fig. \ref{smallphi_full}).
\begin{figure}[htbp]
\centering
\includegraphics[scale=0.3]{fig5_flu.eps}
\caption{Variance of the distribution function in the fixed grid scheme and in the
general case,
for two sizes of the search cell. The lines on the right indicate the values
for $V_{cell}=V_p$, also shown in Fig. \ref{full_esa}.}
\label{smallphi_full}
\end{figure}
The value obtained with $V_{cell} = V_p/2$ is
lower than the corresponding case with  $V_{cell} = V_p$,
indicated as ``fixed grid'' in the Figure.
Moreover,
it further reduces for $V_{cell} = V_p/4$.
Nevertheless, this smearing effect is partly compensated by the fact that the
procedure used to select the cells belonging to the nucleon clouds
favours compact, symmetric configurations.

It is interesting to evaluate the effect
of the use of a small search cell on the fluctuations
also in the general case, i.e. when the fixed grid
condition is relaxed.
Intuitively, we expect
a small effect on the variance, as long as $V_ {cell}$ is not too small:
the main reduction of fluctuations comes
already from the {\it centroid effect} (see previous Subsection),
which is equivalent to
an infinitely small step.
In fact,
the differences between calculations with $V_{cell} = V_p/2$
and $V_{cell} = V_p/4$  (Fig. \ref{smallphi_full}, dotted and
dash-dotted curves respectively) are smaller than in the
previous calculation with fixed grid.
The reduction with respect to the case with $V_{cell} = V_p$ (denoted as ``centroid effect'' in the Figure)
is also smaller.

\subsection{Conclusions on the search procedure}
In conclusion, it appears from our calculations that
the smearing of the nucleon wave packet
is rather dependent on the procedure adopted to build nucleons and is
particularly affected by the so-called {\it centroid effect}, related to
the use of randomly distributed test particles for the phase-space mapping.
 As a consequence, fluctuations
calculated in volumes $V_p$
appear reduced.
As shown above,
the use of a fixed grid to locate the initial and final cells
of the colliding test particles
would help
to increase the value of the variance in $V_p$.
However,
energy and momentum conservation would not be exact, due to the
phase-space discretization.
On the other hand, the standard test particle method allows a good mapping of the phase space and a
good resolution of the average collision integral, with perfect energy and
momentum conservation. However, when constructing fluctuations,
the {\it centroid effect}
may destroy the correlations already built and
increase the nucleon cloud smearing.
However, the expected statistical value
of the fluctuation variance
is recovered
in large volumes.

This last point has been checked
for all the simplified situations considered here (Fermi surface).
For instance, this analysis is shown in Fig.\ref{fluct_new} in the case
of the fixed grid with $V_{cell} = V_p$ and test particles randomly distributed
inside the system (last case of Section 5.2).
\begin{figure}[htbp]
\vspace{0.8cm}
\begin{center}
\includegraphics[scale=0.3]{fig6_flu.eps}
\end{center}
\caption{Variance of the distribution function (rescaled by $N_V$) as a
function of the number $N_V$ of nucleon cells $V_p$, obtained
in the fixed grid case with $V_{cell}=V_p$ and test particles randomly distributed.
\label{fluct_new}
\end{figure}
Here fluctuations are evaluated in volumes $V$, by considering
an ensemble of 800 events.
As expected, the fluctuation variance, rescaled by $N_V$,  almost approaches
the expected value as the volume $V$ increases.
As explained above for the idealized case presented in Fig.3,
the maximum value, $N_V\sigma^2_0 = 0.25$, is not exactly reached
because of the upper bound imposed to the correlation
radius in the nucleon reconstruction procedure and of mass conservation constraints.

\section{Non-equilibrium situations}\label{init_distr}

In the Fermi energy domain, two colliding
nuclei may be schematically represented
by touching Fermi spheres in $\mathbf{p}$ space, so
it is interesting to study the
development of fluctuations starting from
this initial distribution.

In our simplified 2D scheme, we can mimic the two spheres using  a ``chess-board'' initialization, that is
subdividing our momentum space in four regions where $f$ is alternatively zero or one.
Although the system is in a non-equilibrium situation, the variance,
defined as in Eq.(\ref{sigmaf}), starts
from its maximum value
(since $\langle f \rangle =0.5$).
It can be shown that in the fixed grid case
the equilibrium configuration consists of a random
distribution of full and empty cells, and the variance keeps constant in time.
However, a very different result comes out either changing
the search step or removing the fixed grid constraint. In both simulations
the variance rapidly decreases because the initial correlation is destroyed
by the {\it centroid effect}.

As an illustrative example, in Fig. \ref{f1_full}
we plot $\sigma^2_f$ for the choice $V_{cell}=V_p/2$, fixed grid.
Comparison with the corresponding curve
for the configuration $f(\mathbf{p})=0.5 \pm\sigma_{f}(t=0)$
(dashed, see Section 5.4)
shows that the two calculations converge to the same value.
\begin{figure}[htbp]
\vspace{0.8cm}
\begin{center}
\includegraphics[scale=0.3]{fig7_flu.eps}
\end{center}
\caption{Comparison between variances obtained with the ``chess-board''
initialization (solid line) and with the case
$f = 0.5 \pm\sigma_{f}(t=0)$
(dashed curve) in the fixed grid case for $V_{cell}=V_p/2$.
\label{f1_full}
\end{figure}
We get a confirmation of these findings
by removing the fixed grid constraint.
Also in this case, at the end we find no difference with  the case described
before (initialization with randomly distributed test particles, Section 5.3), see Fig. \ref{f1_rand}.
Moreover, in this case the convergence between the two calculations is rather quick.
Then we can conclude that
the final outcome does not depend on the initial
conditions of the system, but only on the final equilibrium situation, as it should be.

\begin{figure}[htbp]
\vspace{0.8cm}
\begin{center}
\includegraphics[scale=0.3]{fig8_flu.eps}
\end{center}
\caption{Comparison between variances obtained with the same initializations
of Fig. \ref{f1_full} but for $V_{cell}=V_p$ and
without the fixed grid constraint.
\label{f1_rand}
\end{figure}

\section{Results for full 3D simulations}

Guided by the simplified schematic cases illustrated above,
let us now discuss the full case of a nuclear Fermi gas of particles at a given
temperature and density interacting through hard two-body scattering.

As in the previous examples,
since we focus on fluctuations in momentum space, only one large cubic cell is
present in coordinate space, and all particles can be chosen to collide
(no restrictions in $\mathbf{r}$ space).
The size of the box is $L=26\,fm$, and we consider $2820$ nucleons, so that
the density has the saturation value $\varrho_0=0.16\,fm^{-3}$; each nucleon
is represented by a collection of $500$ test particles.
 Besides, we do not consider
any distinction between neutrons and protons, so that one nucleon
occupies at least a phase-space volume $h^3/4$ ($g = 4$).
We initialize the momenta so as to
reproduce a Fermi-Dirac profile corresponding to a temperature of $5\,MeV$.
Finally, we consider a constant cross section $\sigma=160\,mb$.
In these calculations the volume $V_{cell}$ corresponds to a cube of side
$l_s=30\,MeV/c$ (and coincides with $V_p = h^3/(4L^3)$).
We notice that, for the full Fermi gas case, a grid
in Cartesian coordinates is easier to use, with respect to
spherical coordinates, in the nucleon construction procedure.
However, spherical coordinates will be adopted to analyze the resulting
fluctuations.
Calculations are stopped when the fluctuation variance saturates.

First of all, we checked
the effect of our fluctuating collision integral on
the average evolution of the system.
In Fig. \ref{effe} we plot the energy
profile of the distribution function at the initial time (solid) and
at the time when the fluctuation variance saturates and
we stop the calculation (dashed line). Slight changes are
observed, probably due to
the finite extension of the nucleon packet,
that induces a discretization of the phase space.

\begin{figure}[htb]
\vspace{0.4cm}
\centering
\includegraphics[scale=0.25]{fig9_flu.eps}
\caption{Distribution function profile at the initial time $t=0$ (solid)
and at the time $t=100\,fm/c$ (dashed line).}
\label{effe}
\end{figure}
We also compared the average collision rate with analytical expectations:
we find that the number of collisions is slightly larger than the
analytical one, but differences are within $10\%$; on the other
hand, similar deviations are also found in usual BUU implementations,
due to the approximate mapping of phase space
induced by the finite number of test particles.

In the following  we analyze in more detail the capability of the procedure to build the
fluctuations in the considered 3D scenario.
We investigated the region in ${\bf p}$ space
where collisions
are most probable and the extension of the nucleon cloud.
From the left panel of Fig. \ref{polar} it is evident
that most collisions involve nucleons lying on the surface of the Fermi
sphere, as expected. Besides, the right panel illustrates the extension of
the clouds in the radial direction, calculated according to the definition:
\begin{equation}\label{delta_p}
\Delta p =\sqrt{\langle(p_i-\langle p\rangle)^2\rangle_{cloud}}, \quad
\end{equation}
where $\langle p\rangle$ denotes the center of momentum of the cloud and
$p_i$ is the modulus of the momentum of the test particles contained in the cloud.

From this distribution one can deduce that the clouds
correspond to a relatively narrow stripe, $2\Delta p \approx 35 MeV/c$,
in the region where $f$ is different
from either zero or one.
\begin{figure}[htb]
\centering
\includegraphics[scale=0.34]{fig10a_flu.eps}
\includegraphics[scale=0.34]{fig10b_flu.eps}
\caption{Left panel: distribution of the position of the centroid of the
nucleon ``clouds'' in the $|p|$ direction. Right panel: distribution of
the size of the nucleons on the same axis.}
\label{polar}
\end{figure}
Moreover, we have checked that the extension of the constructed wave packet
is similar in the three directions $p_x$,$p_y$,$p_z$,
as expected.
We notice also that the tail  of the distribution (Fig.10, right)
extends up to large $\Delta p$ values, $\approx 50 MeV/c$.
Although the nucleon clouds extends, on average,
over a volume $(2\Delta p)^3 \sim~(35\,MeV/c)^3$,
that is not much larger than $V_p$, the {\it centroid effect}, i.e. the fact that
the nucleon centroid may be located anywhere,
causes a significant smearing of fluctuations, as already discussed in
Section 5.
In fact, we find that, for our Fermi gas calculations, fluctuations
in $V_p$ are considerably underestimated, as shown in Fig.\ref{flu_vp},
where the variance $\sigma^2_f$ is reported as a function of the energy
$E  = p^2 / (2m)$.
\begin{figure}[htb]
\vspace{1.2cm}
\centering
\includegraphics[scale=0.3]{fig11_flu.eps
\caption{Profile of the fluctuation variance, as a function of the
energy $E = p^2/(2m)$, in cubic cells
of size $V_p$, for a Fermi gas at equilibrium.}
\label{flu_vp}
\end{figure}
In the following, to test our procedure,
we will evaluate fluctuations considering big volumes,
that contain many nucleons.
However, it should be noticed that, in the geometry considered,
it is not trivial to predict which is the size and shape of
the most appropriate volume to recover the expected value of fluctuations.
While, from one side, it is desirable to consider volumes containing
a large number of nucleons, from the other side one should keep in mind that
spurious fluctuation reduction, due to the
conservation of the total number of particles,
may be a problem.
Hence we have performed a systematic analysis of fluctuations
as a function
of the coordinates
$(p,\theta)$, integrating over the angle $\phi$. The
corresponding volume employed for the calculation of the variance is
given by:
$$
V=2\pi \frac{\Delta p^3}{3}\Delta \cos\theta
$$
We consider volumes corresponding to fixed angular spreads
$\theta$
and/or to a given step $\Delta p^3$.
In Fig. \ref{spicchi_p} we fix the angular spread, $\theta=30^{\circ}$,
and we consider four
different steps, $\Delta p^3$,
corresponding to different volumes, i.e. to a
different number of nucleons $N_{V}$ that can be contained inside.
\begin{figure}[htb]
\vspace{0.4cm}
\centering
\includegraphics[scale=0.31]{fig12_flu.eps}
\caption{Fluctuation variance as a function of the energy, for a fixed angular spread
$\theta=30^{\circ}$ and different steps $\Delta p_{step} = (\Delta p^3)^{1/3}$.}
\label{spicchi_p}
\end{figure}
The considered volumes can
contain up to a rather large number of nucleons (about 90).
The variance, multiplied by $N_V$,  is plotted
as a function of the energy $E$.
We find
that the peak of the variance is around the Fermi energy, as expected.
The rescaled variance gets larger when we increase the step $\Delta p^3$.
However, even in the best situation, the expected equilibrium
value at the peak ($N_V\sigma_0^2 = 0.25$) is underestimated by a factor
$\sim 2.3$.
Therefore, for a fixed value of $\Delta p_{step} = (\Delta p^3)^{1/3}
=190\,MeV/c$, we tested different
angular spreads.
\begin{figure}[htb]

\centering
\includegraphics[scale=0.3]{fig13_flu.eps
\caption{Fluctuation variance as a function of the  energy, for a fixed $\Delta p_{step}=190\,MeV/c$
and different angular spreads (curves with circles). The thick full and dashed
lines are explained in the text.}
\label{spicchi_teta}
\end{figure}
In Fig.\ref{spicchi_teta} we present the obtained results, rescaled by
$N_{V}$.
Now we can observe that, increasing the angular spread, the
variance first increases, then decreases again, due to the mass
conservation constraints.
The variance gets larger
when considering the angular spread $\theta=20 ^{\circ}$, but the expected value is still
underestimated by about a factor $2$, at the peak.
This result can be understood on the basis of
the results obtained in $V_p$ and displayed
in Fig.\ref{flu_vp}.
In fact,
in order to recover  the expected fluctuation value,
the volume employed for the calculation of the variance must be much
larger than $V_p^{1/3}$
in all the directions $(p,\cos (\theta),\phi)$.
Due to the considered geometry,
while this
can be achieved in the angular coordinates, this cannot be easily obtained
along the $p$ direction without losing locality for the value
of $f$ and related fluctuations.
Actually,  in the case of $\Delta p_{step} = 190\,MeV/c$,
the corresponding extension along the $p$ direction,
$\Delta p_{mod}$,
of the volume considered to evaluate the function
$f$ is approximately equal to $30 MeV/c = V_p^{1/3}$ at $p=p_F$.
We thus
expect fluctuations to be reduced, with respect to $f_{eq}{\bar f}_{eq}$,
by a factor
$\alpha^{1/3}\approx 0.48$, being $\alpha\approx 0.11$ the
reduction factor in $V_p$, as evaluated from Fig.\ref{flu_vp}.
The expected fluctuation reduction can be extracted also from
the schematic calculations described in Sect.5.
In fact,
from the inspection of Fig.\ref{full_esa}, it appears that, in the
fixed grid scheme, fluctuations in $V_p$ are reduced by a factor $0.56$.
Relaxing the fixed grid constraint in one
dimension, fluctuations are further reduced
by a factor 0.64 with respect to the fixed grid case.
We thus
expect a fluctuation reduction
by a factor $0.56^{1/3}\times 0.64 \approx 0.5$ (at the Fermi surface).
This is close to our result at $E=E_F = p_F^2 / (2m)$
(see the dashed line with circles in Fig.13, for instance).

One may try to introduce this correction in the results of Fig. 13.
In the considered geometry,
the reduction factor $\alpha$ of the fluctuations in $V_p$ depends
on the energy $E$, i.e. one can write, for the variance in $V_p$:
\begin{equation}
\label{sig1}
\sigma^2_f(E,N_V=1) = {F}(E)~\alpha(E),
\end{equation}
where  ${F}(E)$ indicates the expected fluctuation value, that
in our case should coincide with $f_{eq}{\bar f}_{eq}$.

By approximating the volume dependence of the rescaled variance
by a linear behaviour for small volumes (see Fig.\ref{flu_joseph} for instance),
the fluctuation variance presented in Figs. 12-13 can be rewritten as:
\begin{equation}
\label{sig2}
N_V \sigma^2_f(E,N_V) \approx {F}(E)~ \alpha^{1/3}(E)
\frac{\Delta p_{mod}}{V_p^{1/3}}.
\end{equation}
We notice that $\Delta p_{mod}$ is also depending on $E$.

Eqs.(\ref{sig1})(\ref{sig2}) allow to extract the suppression
factor $\alpha(E)$, as well as the function ${F}(E)$.
The latter is displayed in Fig.13
(thick dashed line).  One can see that
the expected fluctuation value, $f_{eq}{\bar f}_{eq}$ (thick solid line),
is well reproduced at the Fermi surface, while it is underestimated, within $30\%$,
expecially for outer regions.
This underestimation indicates that
in our procedure,
due to the finite extension of the nucleon cloud, the regions far
from the Fermi surface
are less involved in the building of fluctuations than they should be
and can be considered as an intrinsic limitation of the method.

However,
the overall shape of the fluctuation variance is reasonably reproduced, relative to
$f{\bar f}$:
This is a nice indication that the Pauli blocking is
effective, and the equilibrium distribution is almost unchanged by the
procedure. Therefore, it really represents a remarkable improvement with
respect to the original method by Bauer et al. (see Ref.\cite{chapelle}),
which makes us confident
about applications to nuclear reactions.

We can directly visualize the
effect of the fluctuating term by selecting a thin region in momentum space
around the Fermi momentum.
For this purpose, we adopt
a new set of coordinates $(dp,\, pd\theta,\, p \sin(\theta)d\phi )$;
for fixed $p$ and $dp$, the distribution function depends only on
two coordinates, namely:
$$
f(p,\theta,\phi)\rightarrow f(\theta,\eta)
$$
with $\eta=\sin(\theta)\phi$.
This set of coordinates correspond to displacements of equal length
along the three directions: $p,\theta,\phi$.
We choose $p=260\,MeV/c$ (approximately equal to the Fermi momentum) and
$dp=10\,MeV/c$, and in Fig. \ref{tetaeta} we plot the distribution
function $f(\theta,\eta)$ at two different times (initial and final times).
We recognize the sinusoidal
profile on the $\eta$ axis due to the modulation given by the $\sin(\theta)$
factor. At the initial time, $f$ is nearly uniform, and its fluctuations are
simply due to the numerical noise induced by the finite number of test
particles. At later times, we observe
the growth of fluctuations, evidenced by the typical structure with ``peaks
and holes''.
\begin{figure}[htb]
\centering
\includegraphics[scale=0.36]{fig14a_flu.eps}
\includegraphics[scale=0.36]{fig14b_flu.eps}
\caption{Distribution function $f(\theta,\eta)$ at time $t=0$ (left)
and $t=100\,fm/c$ (right).}
\label{tetaeta}
\end{figure}

\section{Methods to increase fluctuations in $V_p$}

\label{optimize}

In actual nuclear collisions the extension of the nucleon wave packet
in phase-space may affect the transport dynamics and it may be important
to correlate the test particles belonging to the nucleon clouds
inside smaller volumes.
In the context of the simplified model (Fermi surface) discussed in Sect.5
we attempted to find a method to
reduce the nucleon cloud smearing and
enhance the fluctuations in $V_p$ without losing
the consistency of the original procedure.

First of all, it should be remembered that the collision probability
is given by $f(\mathbf{p}_1) \bar{f}(\mathbf{p}_3)$
(see Sect. \ref{details}): therefore, we need to consider only
two locations in
$\mathbf{p}$ space.
The main idea is to force the system to fully empty the cells associated with the initial
state or fully fill the final state when a collision occurs.
This condition is not automatically achieved for all cells contributing to
a given collision.
In fact,
the number of particles, $n_t(i)$,
taken from a cell $i$ and employed to build the
nucleon obey to the constraint:
$$
\sum_i n_t(i)=N_{test}
$$
being $N_{test}$ the number of test particles per nucleon.
Hence the number of test particles that are really taken from a cell
is not always
equal to $n_t=min(n_1,\bar{n}_3)$, but instead to $n'_t = min(n_t,n_{res})$,
being
\begin{equation}
n_{res}(i)=N_{test}-\sum_{i'=1}^{i-1}n_t(i')
\label{prob2}
\end{equation}
The idea is to try  to
avoid the situation when the number of particles moved
is less than $n_t$.
In this way
fluctuations get to
a larger value, since
we favour the transitions where cells are completely
filled or emptied.
In order to implement this idea, when constructing the nucleon clouds
 we choose, among the possible ones, the cells for which the
probability $p_t = \frac{n'_t}{n_t}$ is large.

\begin{figure}[htb]
\begin{center}
\includegraphics[scale=0.25]{fig15_flu.eps}
\end{center}
\caption{Variance of the distribution function in $V_p$,
as a function of the number of collisions,
in the fixed grid scheme.
Solid line: no optimization;
dashed line: optimized procedure.
}
\label{opt_full_esa}
\end{figure}
The outcome of this method is plotted
in Fig. \ref{opt_full_esa} for the fixed grid case, where the result
of Fig. \ref{full_esa}
is also shown for comparison (solid line).
As the dashed curve in Fig. \ref{opt_full_esa} clearly shows,
this optimization procedure is successful in
enhancing the fluctuations, which eventually
reach almost their upper limit $\sigma^2_0=0.25$. Actually,
this maximum value is not obtained
due to the fact that  only
adjacent cells are considered in the nucleon construction. In this way, some
half-filled cells remain ``isolated'' and
are not involved in a
collision anymore. Indeed, at the time
when we stop the calculation, isolated cells
amount to about $5\%$ of the total.

The method works also in more general situations, for instance
in the case corresponding to fluctuating initial conditions (see Sect.5.2).
\begin{figure}[thbp]
\begin{center}
\includegraphics[scale=0.25]{fig16_flu.eps}
\end{center}
\caption{Variance of the distribution function in $V_p$, as a function of the
number of collisions, in the fixed grid scheme
with fluctuating initial conditions.
Solid line: no optimization; dashed line: optimized procedure.
}
\label{opt_full_rand}
\end{figure}
The comparison between variances with and without
the optimization procedure is shown in Fig. \ref{opt_full_rand}. As in the
previous case, the optimization favours a faster growth of the fluctuations.
Their saturation value for $V=V_{p}$ is
now enhanced, by $\sim 20\%$ with respect to the standard procedure.
This is an indication that the procedure is able to
increase the number of empty and full cells. This feature is also evident
from the distribution of values of $f$ plotted,  for the same
two simulations of Fig.16, in Fig. \ref{opt_dndf_full}

In the end, this procedure helps the
nucleon packets to keep a similar, more compact shape in $\mathbf{p}$ space during
the time evolution of the system, from one collisional event to another.
As a consequence,
nucleons are more localized and
fluctuations in $V_p$ appear enhanced.

\begin{figure}[tbp]
\vspace{0.3cm}
\centering
\includegraphics[scale=0.25]{fig17_flu.eps}
\caption{Distribution of the number of cells as a function of their occupation $f$ for the same
two calculations of Fig. \ref{opt_full_rand}.
\label{opt_dndf_full}
\end{figure}

\section{Conclusions}
We have proposed a new numerical implementation of the full Boltzmann-Langevin
equation in 3D. The stochastic character of the two-body collision integral,
that in standard transport codes applies only to the single test
particles, leading to a strong reduction of fluctuations,
is recovered, at  the nucleon level, by a careful reconstruction
of the nucleon wave packet once collisions occur.
This is achieved by moving, in each collisional event, two entire test particle clouds,
corresponding to one nucleon each, as proposed in Ref.\cite{fluct_bauer}.
However, while in the original procedure the Pauli blocking was checked
only for the centroids of the nucleon clouds, here we propose a new method that
allows to construct nucleon
wave packets by carefully checking
the Pauli blocking for the entire swarm of test particles.
The last point is essential in the construction of the correct fluctuation
value in fermionic systems.
The nucleon wave packet may take, in principle, any shape. Thus the correlation volume, i.e. the volume of the
sphere that contains all test particles that move together
is generally larger than $h^3$. Moreover, the nucleon centroid is chosen
randomly among the test particle distribution, leading to a partial
overlap of nucleon configurations from one collision to another.
All these effects lead to a reduction of fluctuations in the cells
of volume $h^3$.
However, the correct fluctuation value is recovered in larger volumes,
where all possible nucleon configurations may be accommodated.

The procedure, and its numerical ingredients, has been tested carefully
in several idealized situations, studying in particular the relation
between the results obtained for the fluctuations in cells of volume $h^3$
and the ingredients
of the model.
Considering a fermionic system at equilibrium, at a given density and
temperature, we find that
our procedure builds, within a good approximation,  the expected profile of the fluctuation
variance,  $f_{eq}{\bar f_{eq}}$,  as
a function of the nucleon energy.
This can be considered as a significant improvement with respect
to the original procedure \cite{fluct_bauer}, that opens interesting possibilities
of applications to nuclear collisions.
In fact,
the method proposed can be easily implemented into existing transport
codes, and could allow to treat, in a more complete way,
phenomena where large fluctuations or bifurcations occur and
fluctuations are important.

Finally we stress that the results presented here  are  relevant not only for nuclear
fragmentation studies, but in general for the dynamical description of
quantum many-body systems.

~\\[2ex]

\noindent
{\bf Acknowledgments}
\noindent
We warmly thank H.H.Wolter for the reading of the manuscript and
stimulating discussions.