Vol.3, No.7, 594-599 (2011) Natural Science
http://dx.doi.org/10.4236/ns.2011.37082
Quantum molecular dynamics simulation for
multifragmentation resulting from an expanding
nuclear matter
Atef Abdel-Hafiez1*, Mohamed. Abdo Khalifa2, Albiomy Abd El-Daiem3
1Department of Experimental Nuclear Physics, Nuclear Research Center, Atomic Energy Authority, Cairo, Egypt;
*Corresponding Author:
2Mathematics Department, Faculty of Science, Tanta University, Egypt.
3Physics Department, Faculty of Science, Sohag University, Sohag, Egypt.
Received 20 December, 2010; revised 20 April, 2011; accepted 3 May, 2011.
ABSTRACT
Quantum molecular dynamics (QMD) is used to
investigate multifragmentation resulting from an
expanding nuclear matter. Equation of state, the
structure of nuclear matter and symmetric nu-
clear matter is discussed. Also, the dependence
of the fragment mass distribution on the initial
temperature (Tinit) and the radial flow velocity (h)
is studied. When h is large, the distribution
shows exponential shape, whereas for small h,
it obeys the exponentially falling distribution
with mass number. The cluster formation in an
expanding system is found to be different from
the one in a thermally equilibrated system. The
used Hamiltonian has a classical kinetic energy
term and an effective potential term composed
of four parts.
Keywords: Quantum Molecular Dynamics;
Multifragmentation
1. INTRODUCTION
Multifragmentation has been a long-standing topic in
nuclear physics. In particular, a fragment mass distribu-
tion has been extensively investigated both theoretically
and experimentally as a phenomenon which might have
some connection with nuclear phase transition [1-6].
Among the theoretical approaches, the statistical model
[7] is a familiar one which predicts the fragment mass
distribution.
According to this model, the exponentially falling dis-
tribution with mass number is derived under the condi-
tion that each fragment can be regarded as a thermally
equilibrated droplet. This exponentially falling distribu-
tion with mass number behavior has attracted many au-
thors’ interest because the exponentially falling distribu-
tion with mass number is a signature of a second order
phase transition [8-13].
In order to investigate the relation between the cluster
formation and the fragment mass distribution, many
heavy ion collision experiments have been conducted
[3-6]. In these experiments, two types of fragment mass
distribution have been observed. The exponentially fal-
ling distribution with mass number is one of them and is
originating from a spectator region where thermal equi-
librium is often assumed [3,4]. On the other hand, in a
participant region, the source of fragments undergoes
strong compression and each nucleon has large collec-
tive momentum at its decompression stage. Thus, we
cannot regard this source of fragments as thermal equili-
brated object that the statistical model assumes. In this
region, the mass distribution does not obey the power
law but obeys an exponential distribution and the shape
of which strongly depends on the collision energy [5,6].
There are several factors to prevent the direct under-
standing of cluster formation in heavy ion collision ex-
periments. The finiteness of heavy ion collision is a
source of complication. For instance, the surface effect
during fragment formation might affect the distribution
significantly. Moreover, it is doubtful whether thermal
equilibrium is achieved even in a spectator region. To
check the basic mechanism of fragmentation, therefore,
it is desirable to adopt a model which does not assume a
thermal equilibrium and is free from the complication
arising from the finiteness of the system.
The analysis of nuclear matter is useful to avoid the
complication coming from finiteness. Instability of nu-
clear matter against multifragmentation has been studied
based on the Boltzmann-Langevin approach [14], clas-
sical molecular dynamics [11], kinetic equation [1], and
linear response theory [15]. These are simulations of an
infinite system studying under what condition the insta-

A. Abdel-Hafiez et al. / Natural Science 3 (2011) 594-599 595
bility against the multifragmentation occurs at high ex-
citation. These studies, however, have not taken into
account the collective expansion of the system explicitly.
In [16], the effect of collective expansion has been
treated in the framework of Boltzmann-Langevin ap-
proach and the fragmentation is discussed for a finite
system. The instability to multifragmentation is also
treated in the quantal RPA approach [17].
The aim of this work is to investigate multifragmenta-
tion by using an expanding nuclear matter. When the
matter expands slowly, it simulates a spectator region
where each nucleon moves slowly. On the other hand,
the rapid expansion simulates a participant region where
the system has a rapid collective motion. The rate of
expansion is controlled by only one parameter (h). By
changing h we can investigate both spectator and par-
ticipant regions. In the case of h = 0, the system is just in
a static thermal equilibrium state with temperature (T).
Such a system is similar to what is assumed in the statis-
tical model. Our main purpose is to know how the frag-
ment mass distribution depends on the expansion veloc-
ity and the temperature.
The used quantum molecular dynamics (QMD) and
the nuclear matter at the saturation density

0 and a
temperature Tinit.
On the basis of QMD, a series of simulations is car-
ried out. There are four steps in the calculation.
1) The nuclear matter at the saturation density

0 and
a temperature Tinit is prepared by Metropolis sampling
method.
2) To this thermal matter with Tinit , a radial flow ve-
locity (h) is imposed.
3) The time evolution of the system is calculated from

0 to 0.001

0 by QMD.
4) A fragment mass distribution is calculated at
0.001

0.
The Hamiltonian has a classical kinetic energy term
and an effective potential term composed of four parts:
2
12
A
i
nuclsurface Pauli coulmb
i
p
HVVVV
M


 (1)



11
00
1
0
2
,,
00
1
2
2
,
0
2
2(1 )
1
12
22
1
1
2
1
AA
nucl ii
ii
AA
ex
s
ij ijij
ij iij iij
A
ex
ij
ij iij
V
C
Ccc
C

















 














pp
pp
(2)
 
3
ij
5/3
,
0
d
2
A
SF
surface
ij i
V
Vrrr




 (3)

3
00
22
,i,j
22
.00
2
exp
22
p
pauli
Aij ij
ij
ij i
C
V
qp c
qp








 




pp pp(4)
 
2
,
33
e11
222
1
dd
coulombi j
ij i
ij
V 


 

 
 
 


 rrr r
rr
(5)
where

i is the nucleon spin and

i is the isospin [

i =1/2
for protons and –1/2 for neutrons], 

i is an overlap of
one body density

i for ith nucleon with other nucleons
defined as
 
3
ij
d
AA
iij
ji ji
r  



rr (6)
and the explicit form of 
i is given by
 

2
32
11
exp 2
2π
ii
L
L





rrR
(7)
where L is the squared width of the wave packet, L =
1.95 fm2. This wave packet with fixed width is charac-
teristic of the QMD model. It has the advantage of deal-
ing with the many body correlation directly but imposes
considerable restrictions on the instability of the nuclear
matter.
We approximate the fermionic nature of the system by
introducing the Pauli potential between two identical
nucleons [18,19], instead of antisymmetrization. The
nuclear potential Vnucl has several components: The
first two terms are Skyrme type nuclear interactions. The
term with Cs0 is a symmetry potential where Ci is 1 for
protons and 0 for neutrons. The last two terms with Cex(1)
and Cex(2) are momentum dependent potentials originat-
ing from the exchange term of the Yukawa interaction.
The term Vsurface stands for a surface interaction to de-
scribe the fragment property. The values of the parame-
ters in the interaction are fitted to give good agreements
with the nuclear matter saturation properties, the real
part of optical potential, the rms radii of nuclei, and their
binding energies.
The incompressibility constant K at the saturation
point is simulated to the following parabolic form

2
0
2
0
16
18
K
EA 



Equation of state and the structure of nuclear matter.

A. Abdel-Hafiez et al. / Natural Science 3 (2011) 594-599 596
In this section we study properties of nuclear matter at
several conditions. It is desirable to use a cell large
enough to include several periods of structure and to
avoid the spurious effects of the boundary condition on
the structure of matter. Though our calculation with typi-
cally 1024 particles in a cell is not fully satisfactory in
this respect, we consider it is enough for semi qualitative
discussions at the beginning of this study. Actually, the
global quantities, e.g., the ground state energy of the
system, are well saturated at this number of particles in a
cell.
In Figure 1, we show the energy per nucleon of the
infinite system as a function of the particle number (NA)
in a cell for four average densities from  = 0.4 
0 to 2.2

0. For all densities, the energy of the system has al-
ready approached the asymptotic value above 256 parti-
cles within 100 KeV/ nucleon. In this study we simulate
an infinite system by the periodic boundary condition
mainly with 1024 or 2048 particles in a cell, and inves-
tigate the ground state properties of nuclear matter.
Figure 2 shows the energy per nucleon as a function
of the average density. The solid squares indicate the
energy of ‘‘uniform’’ nuclear matter while open squares
are the results of energy-minimum configurations.
The ‘‘uniform’’ matter energy is calculated as follows:
First we distribute nucleons randomly and cool the sys-
tem only with the Pauli potential. The Pauli potential is
repulsive and does not spoil the uniformity of the system.
Then we impose the other effective interactions and cool
only in the momentum space, fixing the positions of
particles. The system turns out to be approximately uni-
form with this procedure. Note that simulated ‘‘uniform’’
matter is not exactly the same as ideal nuclear matter
since the latter is continuous and completely uniform.
Both cases of uniform and energy-minimum configura-
tions have almost the same energy per nucleon for the
higher densities as is seen in this figure. Below the satu-
Figure 1. Particle number (NA) dependence of the en-
ergy per nucleon (E/A) of the infinite system for four
average densities from  = 0.4 
0 to 2.2

0.
Figure 2. The energy per nucleon of symmetric nuclear
matter [E/A = 0.5] at zero temperature as a function of the
average densities.
ration density

0, the energy per nucleon of the energy-
minimum configuration is lower than the uniform case.
The deviation amounts to about 5 MeV. As we see in the
following, this is due to the structure change of matter
from uniform to nonuniform such as a clusterized one.
From the left, the open squares are the results with
soft [K = 210 MeV] and hard EOS [K = 380 MeV] ob-
tained by the damping equation of motion searching the
energy-minimum configuration in the full phase space.
The solid squares indicate results obtained from the spa-
tially uniform distribution. The kinetic energy of the
electron is not included. We use 1024 particles in a cell
for all cases.
If a system is finite and has high temperature and/or
pressure, it expands spontaneously. However, an infinite
matter cannot expand without an initial collective mo-
tion even if it has high temperature and pressure. There-
fore, we give a collective momentum to each constituent
nucleon so that the matter could expand homogeneously.
This collective momentum is expressed in terms of one
parameter h as
13
0
coll i
iF
ph p


R (8)
where i
R is the position of i th particle [the center of
wave packet] measured from the center of the primitive
cell and pF is the Fermi momentum at

0. The distance is
normalized by using 
0
–1/3. The nearest neighbor pair of
nucleons, on average, have a relative momentum hpF
under this condition. Time evolution from 
0 to 0.001

0
by QMD.
The time evolution of the nuclear matter with a given
temperature Tinit and radial flow velocity with [h,Tinit]
by QMD has been followed by Shinpei Chikazumi, et al.
[20]. Figure 3(a) shows a snapshot of an initial state at

= 
0. Figure 3(b) shows a part of the whole system at

= 0.05 
0 during the expansion. The expansion is stopped
when the average density reaches  = 0.001 
0 shown in
Figure 3(c). In this procedure, the normal periodic

A. Abdel-Hafiez et al. / Natural Science 3 (2011) 594-599 597
Table 1. Effective interaction parameter set.
K = 210 MeV[soft] K = 380 MeV[hard]
 [MeV] –121.9 –21.21
 [MeV] 197.3 97.93
 1.33333 1.66667
Cs0 [MeV] 25.0 25.0
Cex(1) [MeV] –258.5 –258.54
Cex(2) [MeV] 375.6 375.6

1 [MeV] 2.35 2.35

2 [MeV] 0.4 0.4
VSF [MeV] 20.68 20.68
CP [MeV] 115.0 115.0
p0 [MeV] 120.0 120.0
q0 [fm] 2.5 2.5
L [fm2] 1.95 2.05
(a)  = 1.000 
0 (b)
 = 0.05 
0 (c)
 = 0.001 
0
Figure 3. Snapshots of expanding nuclear matter at different
densities. (a) An initial state with
 = 1.0 
0, (b) An intermedi-
ate state at  = 0.05 
0 during the expansion. (c).
boundary condition is not applicable to the expanding
matter. Therefore, we introduce an extended periodic
boundary condition.
2. FRAGMENT MASS DISTRIBUTION
The matter leave to expand until the average density
becomes small enough so that each fragment can be
identified. All the fragments are isolated when the aver-
age distance between nearest neighbor pair reaches
about 5 fm at 0.05

0. However, we noticed that the in-
trinsic expansion of fragments is not yet stopped at
 =
0.05

0 especially when h is large. Therefore, we set the
final density equal to 0.001

0. At this density, all the
fragments are stabilized so that we can identify them
only by their positions.
We show our main result of the fragment mass dis-
tribution. The distribution strongly depends on the radial
flow velocity h whereas the initial temperature has little
effect on the distribution. First, we concentrate on the
distribution for Tinit = 30 MeV and investigate how the
radial flow affects the distribution.
Figure 4 shows the distributions calculated for h = 0.1,
0.5, 1.0, 2.0 with the initial temperature Tinit = 30 MeV.
In this figure, all the distributions follow straight lines
except for h = 0.1. These straight lines in the semi loga-
rithmic scale mean that the number of large fragments
decreases exponentially as fragment size increases. This
exponential shape can be understood as the manifesta-
tion of random distribution[21]. When h is large enough
to suppress the interaction, the final distribution is sim-
ply an enlarged copy of the initial configuration. Under
this limit, the particles are randomly distributed irrespec-
tive of the interaction because the initial configuration is
prepared at saturation density and the particles are ran-
domly distributed under any temperature. The value of h
directly corresponds to the radial flow velocity in the
participant region. Actually, this kind of heavy ion colli-
sion is investigated by a simulation of a finite system
with a radial flow velocity [22]. Figure 5 shows the
same quantity as Figure 4 in double-logarithmic scale.
In our method, we can discuss most clearly the effect of
radial flow because our system is an infinite expanding
nuclear matter which is free from any finite size effect of
the parent nucleus.
Figure 4. The fragment mass distribution obtained in our
QMD simulation in case of semi logarithmic scale.
Figure 5. The fragment mass distribution obtained in our
QMD simulation in case of double logarithmic scale.

A. Abdel-Hafiez et al. / Natural Science 3 (2011) 594-599 598
When h is small, the interaction must play an impor-
tant role in fragmentation. In this case, h = 0.1 follows a
straight line whereas the others do not. The straight line
in double-logarithmic scale refers to the so-called expo-
nentially falling distribution with mass number. The ex-
ponentially falling distribution with mass number has
been discussed in many works in connection with critical
phenomena. Usually, the exponentially falling distribu-
tion with mass number is thought to be a manifestation
that the system undergoes a second order phase transi-
tion. In the model assuming thermal equilibrium, the
power law appears only when the second order phase
transition happens. However, we must remember that the
expanding system is completely in nonequilibrium.
Even when h = 0.1, the system does not have enough
time to reach thermal equilibrium where the temperature
can be determined. Therefore this exponentially falling
distribution with mass number we obtained must be ex-
plained in a different way.
3. EXPANDING FRAGMENTATION VS.
ISOTHERMAL FRAGMENTATION
It is significant to compare two types of the power law,
i.e., dynamical and isothermal ones. The isothermal frag-
mentation can be investigated by Metropolis sampling
method. Like creating initial configuration, we prepare
1000 samples with a fixed temperature T at
 = 0.05 
0
instead of  = 
0. From investigation of isothermal pres-
sure, we have already known that the critical tempera-
ture is around 8 MeV. Figure 6 shows how fragment
mass distribution depends on the temperature T = 5, 8,
18 MeV at  = 0.05 
0. For T = 5 MeV, the distribution
shows the U-shape where large fragments and small
fragments exist but there is no middle size fragments. T
= 8 MeV is just a critical temperature in which the
power law appears. The exponent
 , i.e., the slope of
distribution, is around 2.5 which is consistent with the
Fisher’s droplet model. As the temperature increases
Figure 6. The comparison between isothermal fragmentation and expanding fragmentation.
The three figures [[i], [ii], and [iii]] are isothermal distributions [t = 5, 8 and 18 MeV] at
 =
0.05 
0. The two figures [[iv] and [v]] are the distributions obtained by the expansion [h =
0.10].

A. Abdel-Hafiez et al. / Natural Science 3 (2011) 594-599 599
beyond the critical temperature Tc = 8 MeV, the distribu-
tion deviates from the exponentially falling distribution
with mass number and becomes an exponential shape. In
the bottom figure [T = 18 MeV], thermal motion com-
pletely overcomes attractive interaction. Therefore, a
random distribution appears like the rapid expansion
case.
4. CONCLUTIONS
In this paper we have presented the QMD approach
for studying multifragmentation resulting from an ex-
panding nuclear matter. We can reproduce well the finite
nuclear properties for various mass ranges by inclusion
of the Pauli and surface potential. We have investigated
the EOS of nuclear matter by simulating an infinite sys-
tem with the used QMD. The fragmentation during ex-
pansion can be classified according to the speed of ex-
pansion h. Also, the symmetric nuclear matter is dis-
cussed.
Our QMD model contains a further possibility for the
simulation of the dynamical evolution of infinite nuclear
matter such as supernova explosions, the glitch of the
neutron star, and the initial stage of the universe. An
intensive and systematic study of nuclear matter with the
present model will be important since it contains fewer
assumptions than the foregoing models as to the struc-
ture of matter.
REFERENCES
[1] Schmelzer, J., Ro¨pke, G. and Ludwig, F.-P. (1997) Nu-
clear multifragmentation processes and nucleation theory.
Physical Review C, 55, 1917-1927.
doi:10.1103/PhysRevC.55.1917
[2] Strachan, A. and Dorso, C.O. (1999) Temperature and
energy partition in fragmentation. Physical Review C, 59,
285-294. doi:10.1103/PhysRevC.59.285
[3] Gilkes, M.L. et al., (1994) Determination of critical
exponents from the multifragmentation of gold nuclei.
Physical Review Letters, 73, 1590-1593.
doi:10.1103/PhysRevLett.73.1590
[4] Mastinu, P.F. et al., (1996) Circumstantial evidence for
critical behavior in peripheral Au+Au collisions at 35
MeV/nucleon. Physical Review Letters, 76, 2646-2649.
doi:/10.1103/PhysRevLett.76.2646
[5] Petrovici, M. et al., (1995) Cluster formation during
expansion of hot and compressed nuclear matter pro-
duced in central collisions of Au on Au at 250 A MeV.
Physical Review Letter ,74, 5001-5004.
http://dx.doi.org/10.1103/PhysRevLett.76.2646
[6] Reisdorf, W. et al., (1997) Central collisions of Au on Au
at 150-A/Mev, 250-A/Mev and 400-A/Mev, Nuclear
Physics A, 612, 493-556.
doi:10.1016/S0375-9474(96)00388-0
[7] Das Gupta, S., Mekjian, A.Z. and Tsang, B. (2001)
Symmetry and surface symmetry energies in finite nuclei.
Advances in Nuclear Physics, 26, 89.
[8] Mishustin, I.N. (1998) Fluctuations near the deconfine-
ment phase transition boundary. Nuclear Physics A, 630,
111-125. doi:10.1016/S0375-9474(97)00748-3
[9] Belkacem, M., Latora, V. and Bonasera, A. (1995) Criti-
cal evolution of a finite system. Physical Review C, 52,
271-285. doi:10.1103/PhysRevC.52.271
[10] Colonna, M. and Chomaz, Ph. (1998) Quantal effects on
growth of instabilities in nuclear matter. Physical Review
B, 436, 1-9. doi:10.1016/S0370-2693(98)00882-X
[11] Finocchiaro, P., Belkacem, M., Kubo, T., Latora, V. and
Bonasera, A. (1996) The nuclear liquid-gas phase transi-
tion within fermionic molecular dynamics. Nuclear
Physics A, 600, 236-250.
doi:10.1016/0375-9474(96)00040-1
[12] Ohnishi, A. and Randrup, J. (1997) Quantum fluctuation
effects on nuclear fragment formation. Physical Review B,
394, 260-268. doi:10.1016/S0370-2693(97)00032-4
[13] Sugawa, Y. and Horiuchi, H. (1999) Nucleon-induced
fragment formation with antisymmetrized molecular dy-
namics. Physical Review C, 60, 064613-054621.
doi:10.1103/PhysRevC.60.064613
[14] Colonna, M. andChomaz, Ph. (1994) Unstable infinite
nuclear matter in stochastic mean field approach. Physi-
cal Review C, 49, 1908.
[15] Baran, V., Colonna, M., Di Toro, M. and Larionov, A.B.
(1998) Fragment isotope distributions and the isospin
dependent equation. Nuclear Physics A, 632, 287-303.
doi:10.1016/S0375-9474(98)00006-2
[16] Guarnera, A., Colonna, M. and Chomaz, Ph. (1996) Iso-
spin effects in nuclear fragmentation. Physics Letters B,
373, 267-274. doi:10.1016/0370-2693(96)00152-9
[17] Ayik, S. and Colonna, M. (1995) Quantal effects on
growth of instabilities in nuclear matter. Physical Review
B, 353, 417-421. doi:10.1016/0370-2693(95)00596-D
[18] Dorso, C.O., Duarte, S. and Randrup, J. (1987) Classical
simulation of the fermi gas. Physics Letters B, 188, 287-294.
doi:10.1016/0370-2693(87)91382-7
[19] Maruyama, T. Niita, K., Oyamatsu, K. Maruyama, T.,
Chiba, S. and Iwamoto, A. (1998) Quantum molecular
dynamics approach to the nuclear matter below the satu-
ration density. Physical Review C, 57, 655-665.
doi:10.1103/PhysRevC.57.655
[20] Chikazumi, S., Maruyama, T., Chiba, S., Niita, K. and
Iwamoto, A. (2001) Quantum molecular dynamics simu-
lation of expanding nuclear matter and nuclear multi-
fragmentation. Physical Review C, 63, 024602.
doi:10.1103/PhysRevC.63.024602
[21] Holian B.L. and Grady, D.E. (1988) Fragmentation by
molecular dynamics: The microscopic ‘‘big bang”. Phy-
sical Review Letters, 60, 1355-1358.
doi:10.1103/PhysRevLett.60.1355
[22] Reisdorf, W.et al., (1997) Central collisions of Au on Au
at 150-A/Mev, 250-A/Mev and 400-A/Mev. Nuclear
Physics A, 612, 493-556.
doi:10.1016/S0375-9474(96)00388-0