source: trunk/documents/UserDoc/UsersGuides/PhysicsReferenceManual/latex/hadronic/theory_driven/Multifragmentation/MultFragSim.tex @ 1211

\section{Multifragmentation process simulation.}
\hspace{1.0em} The GEANT4 multifragmentation model is capable to predict 
final states as result of an highly excited nucleus statistical break-up.

The initial information for calculation of
multifragmentation stage consists from the atomic mass number $A$,
charge $Z$ of excited nucleus and its excitation
energy $U$.  At high excitation energies $U/A > 3$ \ MeV the
multifragmentation mechanism, when nuclear system can eventually breaks
down into fragments, becomes the dominant. Later on the excited primary
fragments propagate independently in the mutual Coulomb field and
undergo de-excitation.  Detailed description of multifragmentation
mechanism and model can be found in review \cite{BBIMS95}.

\subsection{Multifragmentation probability.}

\hspace{1.0em}The probability of a breakup channel $b$ is given by the
expression (in the so-called microcanonical approach \cite{BBIMS95},
\cite{Botvina87}):
\begin{equation}
\label{MFS1}W_b(U,A,Z)=\frac{1}{\sum_{b}\exp [S_b(U,A,Z)]} \exp [S_b(U,A,Z)], 
\end{equation}
where $S_b(U,A,Z)$ is the entropy of a multifragment state corresponding
to the breakup channel $b$. The channels $\{b\}$ can be parametrized by
set of fragment multiplicities $N_{A_f,Z_f}$ for fragment with atomic
number $A_f$ and charge $Z_f$. All partitions $\{b\}$ should satisfy
constraints on the total mass and charge:
\begin{equation}
\label{MFS2}\sum_{f}N_{A_f,Z_f}A_f = A 
\end{equation}
and
\begin{equation}
\label{MFS3}\sum_{f}N_{A_f,Z_f}Z_f = Z. 
\end{equation}
It is assumed \cite{Botvina87} that thermodynamic equilibrium is
established in every channel, which can be characterized by the channel
temperature $T_b$.

The channel temperature $T_b$ is determined by the equation constraining
the average energy $E_b(T_b, V)$ associated with partition $b$:
\begin{equation}
\label{MFS4}E_b(T_b, V)= U+E_{ground} = U+M(A,Z), 
\end{equation}
where $V$ is the system volume, $E_{ground}$ is the ground state (at
$T_b = 0$) energy of system and $M(A,Z)$ is the mass of nucleus.

According to the conventional thermodynamical formulae the average energy
of a partitition $b$ is expressed through the system free energy $F_b$
as follows
\begin{equation}
\label{MFS5}E_b(T_b, V)= F_b(T_b,V) +T_bS_b(T_b,V).
\end{equation} 
Thus, if free energy $F_b$ of a partition $b$ is known, we can find
the channel temperature $T_b$ from Eqs. ($\ref{MFS4}$) and ($\ref{MFS5}$),
then the entropy $S_b = -dF_b/dT_b$ and hence, decay probability $W_b$
defined by Eq. ($\ref{MFS1}$) can be calculated.

Calculation of the free energy is based on the use of the liquid-drop
description of individual fragments \cite{Botvina87}.  The free energy
of a partition $b$ can be splitted into several terms:
\begin{equation}
\label{MFS6}F_b(T_b,V) =  \sum_{f}F_f(T_b,V) + E_{C}(V), 
\end{equation}
where $F_f(T_b,V)$ is the average energy of an individual fragment
including the volume
\begin{equation} 
\label{MFS7} F^V_f = [-E_0 - T^2_b/\epsilon(A_f)]A_f, 
\end{equation}
surface 
\begin{equation}
\label{MFS8}F^{Sur}_f = \beta_0[(T_c^2 - T^2_b)/(T_c^2 + T^2_b)]^{5/4}A_f^{2/3}
= \beta(T_b)A_f^{2/3},
\end{equation}
symmetry 
\begin{equation}
\label{MFS9}F^{Sim}_f = \gamma(A_f - 2Z_f)^2/A_f, 
\end{equation}
Coulomb 
\begin{equation}
\label{MFS10}F^{C}_f = \frac{3}{5}\frac{Z_f^2e^2}{r_0A_f^{1/3}}
[1 - (1+ \kappa_{C})^{-1/3}]
\end{equation}
and translational 
\begin{equation}
\label{MFS11}F^{t}_f = -T_b\ln{(g_fV_f/\lambda^3_{T_b})} + 
T_b\ln{(N_{A_f,Z_f}!)}/ N_{A_f,Z_f}
\end{equation}
terms  and the last term
\begin{equation}
\label{MFS12} E_{C}(V)=\frac{3}{5}\frac{Z^2e^2}{R} 
\end{equation}  
is the Coulomb energy of the uniformly charged sphere with charge $Ze$
and the radius $R = (3V/4\pi)^{1/3}= r_0A^{1/3}(1 + \kappa_{C})^{1/3}$,
where $\kappa_{C} = 2$ \cite{Botvina87}.

Parameters $E_0 = 16$ \ MeV, $\beta_0 = 18$\ MeV, $\gamma = 25$ \ MeV
are the coefficients of the Bethe-Weizsacker mass formula at $T_b = 0$.
$g_f=(2S_f+1)(2I_f+1)$ is a spin $S_f$ and isospin $I_f$ degeneracy
factor for fragment ( fragments with $A_f > 1$ are treated as the
Boltzmann particles), $\lambda_{T_b} = (2\pi h^2/m_N T_b)^{1/2}$ is the
thermal wavelength, $m_N$ is the nucleon mass, $r_0 = 1.17$ \ fm,
$T_c=18$ \ MeV is the critical temperature, which corresponds to the
liquid-gas phase transition. $\epsilon(A_f) = \epsilon_0[1 + 3/(A_f-1)]$
is the inverse level density of the mass $A_f$ fragment and
$\epsilon_0=16$ \ MeV is considered as a variable model parameter, whose
value depends on the fraction of energy transferred to the internal
degrees of freedom of fragments \cite{Botvina87}. The free volume $V_f
=\kappa V=\kappa\frac{4}{3}\pi r_0^4 A$ available to the translational
motion of fragment, where $\kappa \approx 1$ and its dependence on the
multiplicity of fragments was taken from \cite{Botvina87}:
\begin{equation}
\label{MFS13}\kappa =[1 + \frac{1.44}{r_0A^{1/3}}(M^{1/3} - 1)]^{3} - 1.
\end{equation}
For $M = 1$ $\kappa = 0$.

The light fragments with $A_f < 4$, which have no excited
states, are considered as elementary particles characterized by the
empirical masses $M_f$, radii $R_f$, binding energies $B_f$, spin
degeneracy factors $g_f$ of ground states.  They contribute to the
translation free energy and Coulomb energy.


\subsection{ Direct simulation of the low multiplicity multifragment
disintegration.} 

\hspace{1.0em}At comparatively low excitation energy (temperature)
system will disintegrate into a small number of fragments $M \leq 4$ and
number of channel is not huge. For such situation a direct
(microcanonical) sorting of all decay channels can be performed. Then,
using Eq. ($\ref{MFS1}$), the average multiplicity value $<M>$ can be
found.  To check that we really have the situation with the low
excitation energy, the obtained value of $<M>$ is examined to obey the
inequality $<M> \leq M_0$, where $M_0 = 3.3$ and $M_0 = 2.6$ for $A \sim
100$ and for $A \sim 200$, respectively \cite{Botvina87}.  If the
discussed inequality is fulfilled, then the set of channels under
consideration is belived to be able for a correct description of the
break up. Then using calculated according Eq. ($\ref{MFS1}$)
probabilities we can randomly select a specific channel with given
values of $A_f$ and $Z_f$.
 

\subsection{ Fragment multiplicity distribution.}

\hspace{1.0em}The individual fragment multiplicities $N_{A_f,Z_f}$ in
the so-called macrocanonical ensemble \cite{BBIMS95} are distributed
according to the Poisson distribution:
\begin{equation}
\label{MFS14} P(N_{A_f,Z_f}) = \exp{(-\omega_{A_f,Z_f})}
\frac{\omega_{A_f,Z_f}^{N_{A_f,Z_f}}}{N_{A_f,Z_f}!}
\end{equation}
with mean value $<N_{A_f,Z_f}>=\omega_{A_f,Z_f}$ defined as
\begin{equation}
\label{MFS15} <N_{A_f,Z_f}> =g_fA_f^{3/2}\frac{V_f}{\lambda^3_{T_b}}
\exp{[\frac{1}{T_b}(F_f(T_b,V)-F^{t}_f(T_b,V) - \mu A_f - \nu Z_f)]}, 
\end{equation}
where $\mu$ and $\nu$ are chemical potentials. The chemical potential
are found by substituting Eq. ($\ref{MFS15}$) into the system of
constraints:
\begin{equation}
\label{MFS16}\sum_{f}<N_{A_f,Z_f}>A_f = A 
\end{equation}
and
\begin{equation}
\label{MFS17}\sum_{f}<N_{A_f,Z_f}>Z_f = Z
\end{equation}
and solving it by iteration. 
 
 
\subsection{ Atomic number distribution of fragments.}

\hspace{1.0em}Fragment atomic numbers $A_f > 1$ are also distributed
according to the Poisson distribution \cite{BBIMS95} (see
Eq. ($\ref{MFS14}$)) with mean value $<N_{A_f}>$ defined as
\begin{equation}
\label{MFS18} <N_{A_f}> = A_f^{3/2}\frac{V_f}{\lambda^3_{T_b}}
\exp{[\frac{1}{T_b}(F_f(T_b,V)-F^{t}_f(T_f,V) - \mu A_f - \nu <Z_f>)]}, 
\end{equation}
where calculating  the  internal free energy  
$F_f(T_b,V)-F^{t}_f(T_b,V)$  one has to substitute $Z_f
\rightarrow <Z_f>$. The average  charge $<Z_f>$ for fragment having atomic 
number $A_f$ is given by
\begin{equation}
\label{MFS19}<Z_f(A_f)> = \frac{(4\gamma + \nu)A_f}{8\gamma + 2[1 - (1 +
\kappa)^{-1/3}]A_f^{2/3}}.
\end{equation}

\subsection{ Charge distribution of fragments.}

\hspace{1.0em}At given mass of fragment $A_f > 1$ the charge $Z_f$
distribution of fragments are described by Gaussian
\begin{equation}
\label{MFS20} P(Z_f(A_f))\sim \exp{[-\frac{(Z_f(A_f) - 
<Z_f(A_f)>)^2}{2(\sigma_{Z_f}(A_f))^2}]}
\end{equation} 
with dispertion 
\begin{equation}
\label{MFS21}\sigma_{Z_f(A_f)} = \sqrt{\frac{A_fT_b}
{8 \gamma + 2[1 - (1 +\kappa)^{-1/3}]A_f^{2/3}}} 
\approx \sqrt{\frac{A_fT_b}{8 \gamma}}.
\end{equation} 
and the average charge $<Z_f(A_f)>$ defined by Eq. ($\ref{MFS17}$).

\subsection{ Kinetic energy distribution of fragments.}

\hspace{1.0em}It is assumed \cite{Botvina87} that at the instant of the
nucleus break-up the kinetic energy of the fragment $T^{f}_{kin}$ in the
rest of nucleus obeys the Boltzmann distribution at given temperature
$T_b$:
\begin{equation}
\label{MFS22} \frac{dP(T^{f}_{kin})}{dT^{f}_{kin}}\sim \sqrt{T^{f}_{kin}}
\exp{(-T^{f}_{kin}/T_b)}.
\end{equation}
Under assumption of thermodynamic equilibrium the fragment have
isotropic velocities distribution in the rest frame of nucleus. The
total kinetic energy of fragments should be equal $\frac{3}{2}MT_b$,
where $M$ is fragment multiplicity, and the total fragment momentum
should be equal zero. These conditions are fullfilled by choosing
properly the momenta of two last fragments.

The initial conditions for the divergence of the fragment system are
determined by random selection of fragment coordinates distributed with
equal probabilities over the break-up volume $V_f = \kappa V$. It can be
a sphere or prolongated ellipsoid. Then Newton's equations of motion are
solved for all fragments in the self-consistent time-dependent Coulomb
field \cite{Botvina87}.  Thus the asymptotic energies of fragments
determined as result of this procedure differ from the initial values by
the Coulomb repulsion energy.

\subsection{Calculation of the fragment excitation energies.} 

\hspace{1.0em}The temparature $T_b$ determines the average excitation
energy of each fragment:
\begin{equation}
\label{MFS23} U_{f}(T_b) = E_f(T_b) - E_f(0) =  \frac{T_b^2}{\epsilon_0}A_f +
[\beta(T_b) - T_b \frac{d\beta(T_b)}{dT_b} - \beta_0]A^{2/3}_f,  
\end{equation} 
where $E_f(T_b)$ is the average fragment energy at given temperature
$T_b$ and $\beta(T_b)$ is defined in Eq. ($\ref{MFS8}$).  There is no
excitation for fragment with $A_f < 4$, for $^{4}He$ excitation energy
was taken as $U_{^{4}He} = 4T^2_b/\epsilon_o$.