# A soluble statistical model for nuclear fragmentation

## I INTRODUCTION

Nuclear fragmentation resulting from heavy ion collsions is a complex phenomenon. The role of equilibration and dynamics has not yet been determined as a plethora of approaches have been investigated. Examples of approaches are evaporative pictures, percolation models, lattice gas models, and dynamical models based on Boltzmann simulations. In this paper we consider the statistical approach where one considers sampling all configurations of non-interacting clusters. Recently, Chase and Mekjian derived relations which allow the exact calculation of the canonical partition function for such a system. By eliminating the need for computationally intensive Monte Carlo procedures and associated approximations, this technique allows a deeper insight into the thermodynamic principles which drive the statistics of fragmentation.

In the next section we present the recursive technique of Chase and Mekjian and review the thermodynamic properties, some of which have already been presented in the literature. We emphasize that the surface energy is the most important parameter in determining the fragmentation and phase transition properties of the model. In the three subsequent sections, we present extensions of the model which are necessary for serious modeling of nuclear systems: excluded volume, Coulomb effects, and isospin degrees of freedom. In section VI we show how a microcanonical distribution may be generated from the canonical distribution.

## II The Model

For completeness, we present an outline of the model, which is based on the work of Chase and Mekjian. The expressions used here are based on a picture of non-interacting liquid drops. Mekjian and Lee had also applied similar recursion relations to a more algebraically motivated fragmentation model that was not based on a liquid-drop picture.

We consider that there are $`A`$ nucleons which thermalize in a volume $`V`$ much larger than $`V_0`$ where $`V_0=A/\rho _0`$ is the ground state volume of a nucleus of $`A`$ nucleons. These nucleons can appear as monomers but also as composites of $`a`$ nucleons. The canonical partition function of this system can be written as

$`\mathrm{\Omega }_A={\displaystyle \underset{\mathrm{\Sigma }n_ka_k=A}{}}\mathrm{\Pi }_k{\displaystyle \frac{\omega _k^{n_k}}{n_k!}}`$ (1)   
where $`\omega _k`$ is the partition function of a single composite of size $`a_k`$, $`n_k`$ is the number of such composites and the sum goes over all the partitions which satisfy $`n_ka_k=A`$. A priori this appears to be a horrendously complicated problem but $`\mathrm{\Omega }_A`$ can be computed recursively via the formula,

$$\mathrm{\Omega }_A=\frac{1}{A}\underset{k}{}\omega _k\mathrm{\Omega }_{Aa_k}$$
(2)   
Here $`\mathrm{\Omega }_0`$ is 1. It is this formula and the generalisation of this to more realistic case (see later) that makes this model so readily soluble.

All properties of the system are determined by the partition functions of indepedent particles. The recursive formula above allows a great deal of freedom in the choice of partition functions for individual fragments, $`\omega _k`$. Any function of temperature, density and $`A`$ is allowed. However, explicit dependence on the configuration of the remainder of the system is outside the scope of this treatment.

For the illustrative purposes of this section, we assume the form,

$$\omega _k=\frac{V}{\mathrm{}^3}\left(\frac{a_kmT}{2\pi }\right)^{3/2}\times e^{F_{k,\mathrm{int}}/T}$$
(3)   
The first part is due to the kinetic motion of the center of mass of the composite in the volume $`V`$ and the second part is due to the internal structure. Following the choice of reference we assume the form

$$F_{k,\mathrm{int}}=W_0a_kS(T)a_k^{2/3}T^2a_k/ϵ_0$$
(4)   
Here $`W_0`$ is the volume energy per nucleon(=16 MeV), $`S(T)`$ is the surface tension which is a function of the temperature $`T`$. The origin of the different terms in Eq. (4) is the following: $`W_0k+Sk^{2/3}`$ is the ground state energy of the composite of $`k`$ nucleons, and the last term in the exponential arises because the composite can be not only in the ground state but also in excited states which are included here in the Fermi-gas approximation. Following reference  the value of $`ϵ_0`$ is taken to be 16 MeV. Lastly the temperature dependence of $`S(T)`$ in ref is $`S(T)=S(0)[(T_c^2T^2)/(T_C^2+T^2)]^{5/4}`$ with $`S(0)=18`$ MeV and $`T_c=18`$ MeV. Any other dependence could be used including a dependence on the average density.

Upon calculation, the model described above reveals a first order phase transition. In Figure 1 the specific heat at constant volume, $`C_V=(1/A)dE/dT`$, is displayed as a function of temperature for systems of size, $`A=700`$, $`A=1400`$ and $`A=2800`$. The sharp peak represents a discontinuity in the energy density, which sharpens for increasingly large systems. The usual picture of a liquid-gas phase transition gives a discontinuity in the energy density when pressure is kept constant rather than when the volume is kept constant. To understand this result we consider a system divided into one large cluster and many small clusters. The pressure and free energy may then be approximated as

$`E/A`$ $``$ $`ϵ_{\mathrm{bulk}}+{\displaystyle \frac{3}{2}}{\displaystyle \frac{N_{\mathrm{cl}.}}{A}}T,`$ (5)   
$`P`$ $`=`$ $`{\displaystyle \frac{N_{\mathrm{cl}.}}{V}}T,`$ (6)   
where $`N_{cl}`$ is the number of clusters. The bulk term depends only on the temperature and not on the way in which the nucleons are partioned into fragments. We have neglected the surface energy term which is proportional to $`A^{1/3}`$. In this limit, $`C_v`$ and $`C_p`$ become

$`C_V`$ $`=`$ $`{\displaystyle \frac{ϵ_{\mathrm{bulk}}}{T}}+{\displaystyle \frac{3}{2}}{\displaystyle \frac{N_{\mathrm{cl}.}}{A}}`$ (7)   
$`C_p`$ $`=`$ $`C_V+{\displaystyle \frac{N_{\mathrm{cl}.}}{A}}.`$ (8)   
The bulk term depends only on the temperature and is therefore continuous across the phase transition. Thus, a spike in $`C_p`$ is equivalent to a spike in $`C_V`$ since both are proportional to $`N_{\mathrm{cl}.}`$. It is difficult to make a connection between this approach and the standard Maxwell construction, since here interactions between particles enter only through the surface term.

Intrinsic thermodynamic quantities may be calculated in a straightforward manner. For instance the pressure and chemical potentials may be calculated through the relations,

$`\mu `$ $`=`$ $`T\left(\mathrm{ln}\mathrm{\Omega }_A\mathrm{ln}\mathrm{\Omega }_{A1}\right)`$ (9)   
$`P`$ $`=`$ $`T{\displaystyle \frac{\mathrm{ln}(\mathrm{\Omega }_A)}{V}}`$ (10)   
Calculations of $`\mu `$ and $`P`$ are displayed in Figure 2 as a function of density for a system of size $`A=200`$. Both the pressure and chemical potential remain roughly constant throughout the region of phase coexistence. Of particular note is that the pressure actually falls in the coexistence region due to finite size effects.

We now make some comments about influences of various factors in Eq. (4). The bulk terms, $`W_0+T^2/ϵ_0`$, are not affected by the free energy, thus they may be ignored when calculating fragmentation observables. Their influence with respect to intrinsic thermodynamic quantities is of a trivial character. The surface term $`S(T)`$ is completely responsible for determining all observables related to fragmentation and therefore all aspects of the phase transition. Aside from the system size $`A`$, fragmentation is determined by two dimensionless parameters. The first is the specific entropy, $`(V/A)(mT/(2\pi \mathrm{}^2))^{3/2}`$ and the second is the surface term $`S(T)/T`$.

At a given temperature the free energy $`F=ETS`$ of $`A`$ nucleons should be minimized. With the surface tension term, $`E`$ is minmised if the whole system appears as one composite of $`A`$ nucleons but the entropy term encourages break up into clusters. At low temperatures the surface term dominates while at high temperatures entropy prevails and the system breaks into small clusters. The mass distribution may be calculated given the partition function.

$$n_k=\frac{\omega _k\mathrm{\Omega }_{Aa_k}}{\mathrm{\Omega }_A}$$
(11)   
The mass distribution is displayed in Figure 3 for three temperatures, 6.0, 6.25 and 6.5 MeV which are centered about the transition temperature of 6.25 MeV. The distributions have been multiplied by $`a_k`$ to emphasize the decomposition of the system. The mass distribution changes dramatically in this small temperature range. The behavior is reminiscent of that seen in the percolation or lattice gas models.

## III Excluded Volume

The volume used to define to the partition functions of individual fragments, $`\omega _k`$ given in Eq. (4), should reflect only that volume in which the fragments are free to move. Hahn and Stöcker suggested using $`VVA/\rho _0`$ to incorporate the volume taken up by the nuclei. By inspecting Eq. (4) on can see that this affects the partion function by simply mapping the density or volume used to plot observables. More realistically, the excluded volume could depend upon the multiplicity. Nonetheless, in rather complicated calculations not reported here, it was found that for the purpose of obtaining $`pV`$ diagrams in the domain of interest in this paper, it is an acceptable approximation to ignore the multiplicity dependence of the excluded volume.

Incorporating a multiplicity dependence would be outside the scope of the present model, as it would represent an explicit interaction between fragments. However, one could add an $`a`$-dependence to the volume term to account for the difficulty of fitting fragments of various sizes into a tight volume. This might affect the model in a non-trivial fashion.

We like to remind the reader that the parameter $`b`$ in the Van der Waals EOS: $`(p+a/V^2)(Vb)=RT`$ also has its roots in the excluded volume. But there $`b`$ plays a crucial role. We could not for example set $`b`$=0 without creating an instability at high density. Furthermore, the phase transition disappears when $`a`$ is set to zero.

## IV Coulomb effects

It has been understood that the Coulomb effects alter the phase structure of nuclear matter. Although explicit Coulomb interactions are outside the scope of this treatment, they may be approximated by considering a screened liquid drop formula for the Coulomb energy as has been used by Bondorf and Donangelo. The addition to the internal free energy given in Eq. (4) is

$$F_{\mathrm{coul}}=0.70\left(1\left(\frac{\rho }{\rho _0}\right)^{1/3}\right)\frac{a_k^{5/3}}{4}\mathrm{MeV}.$$
(12)   
This form implies a jellium of uniform density that cancels the nucleons positive charge when averaged over a large volume. This may be more physically motivated for the modeling of stellar interiors where the electrons play the role of the jellium.

We display $`C_v`$, both with and without Coulomb terms for an $`A=100`$ system in Figure 4. Coulomb forces clearly reduce the temperature at which the transition occurs. For sufficiently large systems, Coulomb destroys the transition as large drops become unstable to the Coulomb force.

## V Conservation of Isospin

The recursive approach employed here is easily generalized to incorporate multiple species of particles. If there exist a variety of particles with conserved charges $`Q_1`$, $`Q_2\mathrm{}`$, one can write a recursion relation for each charge.

$$\mathrm{\Omega }_{Q_1,Q_2\mathrm{}}=\underset{k}{}\frac{q_{i,k}}{Q_i}\omega (k)\mathrm{\Omega }_{Q_1q_{1,k},\mathrm{}Q_iq_{i,k}\mathrm{}},$$
(13)   
where $`Q_i`$ is the net conserved charge of type $`i`$ and $`q_{i,k}`$ is the charge of type $`i`$ carried by the fragment noted by $`k`$.

For the nuclear physics example, one would wish to calculate $`\mathrm{\Omega }_{N,Z}`$ where $`N`$ and $`Z`$ were the conserved neutron and proton numbers. To find $`\mathrm{\Omega }_{Z,N}`$ one must know $`\mathrm{\Omega }_{N^{},Z^{}}`$ for all $`N^{}<N`$ or $`Z^{}<Z`$. To accomplish this one must use both recursion relations.

## VI Obtaining the microcanonical distribution

In nuclear collisions, one does not have access to a heat bath, but one can vary the excitation energy. A microcanonical treatment is therefore more relevant for practical calculations, particularly given the existence of a first order phase transition which occupies an infinitesimal (in the limit of large $`A`$) range of temperatures in a canonical calculation, but a finite range of energies in a microcanonical ensemble.

The relevant partition function for a microcanonical ensemble is the density of states,

$`\rho (E)`$ $`=`$ $`{\displaystyle \underset{i}{}}\delta (E_iE)`$ (14)   
$`=`$ $`{\displaystyle \frac{1}{2\pi }}{\displaystyle 𝑑\beta \underset{i}{}e^{i\beta (EE_i)}}`$ (15)   
$`=`$ $`{\displaystyle \frac{1}{2\pi }}{\displaystyle 𝑑\beta \mathrm{\Omega }(i\beta )e^{i\beta E}},`$ (16)   
where the sum over $`i`$ represents the sum over all many-body states. Although $`\mathrm{\Omega }(i\beta )`$ is easily calculable given the recursion relations discussed in the previous sections, one must perform the integral over $`\beta `$ numerically.

The true solution for the density of states would be ill-defined given the discreet nature of quantum spectra which can not be combined with a delta function. However, if one defines the density of states in a finite region of size $`\eta `$, the density of states becomes well-behaved even for discreet spectra. For that reason we more pragmatically define the density of states as

$`\rho _\eta (E)`$ $``$ $`{\displaystyle \underset{i}{}}{\displaystyle \frac{1}{\sqrt{2\pi }\eta }}\mathrm{exp}{\displaystyle \frac{(EE_i)^2}{2\eta ^2}}`$ (17)   
$`=`$ $`{\displaystyle \frac{1}{2\pi }}{\displaystyle 𝑑\beta \mathrm{\Omega }(i\beta )e^{i\beta E\eta ^2\beta ^2/2}}`$ (18)   
One might have considered replacing the delta function by a Lorentzian rather than by a Gaussian, but this would be dangerous given that the density of states usually rises exponentially for a many-body system. The finite range $`\eta `$ used to sample the density of states might correspond to the range of excitation energies sampled in an experimental binning. In the limit $`\eta 0`$, $`\rho _\eta `$ approaches the density of states.

As an example of a quantity one may wish to calculate with a microcanonical approach, we consider the average multiplicity of a fragment of type $`k`$ in a system whose total energy is within $`\eta `$ of $`E`$.

$`n_k_\eta (E)`$ $`=`$ $`{\displaystyle \frac{_in_{i,k}\frac{1}{\sqrt{2\pi }\eta }\mathrm{exp}\frac{(EE_i)^2}{2\eta ^2}}{\rho _\eta (E)}}`$ (19)   
$`=`$ $`{\displaystyle \frac{\frac{1}{2\pi }𝑑\beta \omega _k(i\beta )\mathrm{\Omega }_{Aa_k}(i\beta )e^{i\beta E\eta ^2\beta ^2/2}}{\rho _\eta (E)}},`$ (20)   
where $`n_{i,k}`$ is the number of particles of species $`k`$ within the fragment $`i`$.

The integration over $`\beta `$ clearly provides an added numerical challenge that increases for small $`\eta `$. For the purposes of generating a mass distribution, one must perform this integration for every species. It might be worthwile to consider estimating the integrals over $`\beta `$ with the saddle point method, although one should be wary of taking derivatives of $`\mathrm{\Omega }`$ with respect to $`\beta `$ in the phase transition region.

Microcanonical quantities might also be calculated in a completely different manner by discreetizing the energy. For instance one might measure energies in units of 0.1 MeV. One might then treat energy on the same footing as any other conserved charge. One may then write recursion relations for $`N_{A,E}`$, the number of ways to arrange $`A`$ nucleons with net energy $`E`$, where $`E`$ is an integer.

$`N_{A,E}`$ $`=`$ $`{\displaystyle \underset{k,e_k}{}}{\displaystyle \frac{a_k}{A}}\omega _{k,e_k}N_{Aa_k,Ee_k}`$ (21)   
$`=`$ $`{\displaystyle \underset{k,e_k}{}}{\displaystyle \frac{e_k}{E}}\omega _{k,e_k}N_{Aa_k,Ee_k}`$ (22)   
Here, $`\omega _{k,e_k}`$ is the number of ways of arranging a fragment of type $`k`$ with net energy $`e_k`$. All other relavant microcanonical quantities may be calculated in a similar manner.

Since one needs to calculate $`N`$ at all energies $`E^{}`$ less than the targeted energy $`E`$, and must sum over all energies less than $`E^{}`$ to obtain $`N_{A^{},E^{}}`$, the length of the calculation is proportional to $`E^2`$. Typically, nuclear decays occur with on the order of a GeV of energy deposited in a nucleus. Therefore, these calculations may become numerically cumbersome unless the energy is discreetized rather coarsely.

## VII summary

The recursive techniques discussed here have several attractive features. They are easy to work with, incorporate characteristics of nuclear composites and appear to have standard features of liquid-gas phase transitions. In the present forms these models are resricted to low densities. For modeling nuclear disintegration this is not a serious problem, although for completeness it would be nice to be able to modify the model so that it can be extended to higher density.

In this paper we have studied thermal properties of the model, and we emphasize the importance of the surface term in determining these properties. We can associate the discontinuity in the energy density with temperature to the discontinuity in the number of clusters. In addition, we have seen that including Coulomb effects lowers the temperature at which the fragmentation transition occurs and reduces the sharpness of the phase transition. We have also presented an extension of the formalism for the inclusion of isospin degrees of freedom.

For comparing to nuclear physics experiments, development of the microcanonical approaches presented here is of greatest importance. It remains to be seen whether the microcanonical formalisms are tenable, as they have yet to be implemented.

###### Acknowledgements.

This work is supported in part by the Natural Sciences and Engineering Research Council of Canada and by le fonds pour la Formation de Chercheurs et l’aide á la Recherche du Québec, by the US Department of Energy, Grant No. DE FG02-96ER 40987, and by the US National Science Foundation, grant 96-05207.