M Cojonna

Ph Ciiomaz

M Colonna

Ph Chomaz

Spinodal decomposition in nuclear molecular dynamics *

A

We-etud-y Unstable zero-sound waves and cluster formation' in the framework of molecular dynamics approaches based on gaussian single particle wave packets.

We observMhat zero-sound properties are significantly affected when the gaussian width takes large values. In approaches considering the widths as dynamical variables, cluster formation by spinodal decomposition is inhibited, due to the spreading of the nucleon wave packet.

During the past decade numerous theoretical efforts have been devoted to build new approaches able to deal with the (nuclear) many-body problem. Indeed, the advent of more exhaustive experimental data has put in evidence the inadequacy of average descriptions based on mean-field approaches, such as the Time Dependent Hartree-Fock (TDHF) theory [1,2], to explain the large variety of final configurations observed in the data. In the 80's numerical simulations based on molecular dynamics have been introduced [3,4,5,6,7]. The basic idea of these approaches is that the diversity of final states originates from small differences in the initial conditions which are magnified by the chaotic nature of the molecular propagation. These molecular dynamics models are usually of classical nature [7]. However new approaches incorporating some quantum aspects have been recently developed [3,4,5,6]. This is in particular the case of the so-called Antisymmetrized Molecular Dynamics (AMD) [3,6] and Fermionic Molecular Dynamics (FMD) [5].

From a formal point of view, these approaches are, first of all, based on the parameterization of the many body wave function in terms of a Slater determinant, i.e. as an independent particle wave function. In this sense these models are deeply related to TDHF. However, while in TDHF the single particle wave functions are not constrained, in AMD or FMD they are approximated by gaussian wave packets and the parameters of the gaussians are supposed to evolve smoothly with time. In AMD, only the gaussian centroids are considered as dynamical variables, while in FMD also the gaussian widths can evolve with time.

It should be noticed that this additional gaussian ansatz induces a localization of the particles and violates the indistinguishability principle. In fact the indistinguishability of N independent identical particles is mathematically enforced by a U(N) gauge invariance among the N occupied single particle states. This local symmetry in time is explicitly broken by the Gaussian assumption since a 1 unitary mixing of gaussian functions is not gaussian. This is the reason why in these simulations one can always trace back the origin of any wave packet. For example, during a collision, it is possible to recognize if an evaporated nucleon comes from the target or from the projectile. Conversely, in a TDHF description several pieces of different single particle wave functions can always co-operate to build a new orbital.

From this point of view, nuclear molecular dynamics approaches can be seen as an approximation to TDHF. However it should be noticed that breaking symmetries is a general way to introduce correlations, so that the molecular dynamics can be able to afford more physics than the original TDHF. For example, when the width is kept constant one might recover the exact classical many body dynamics that goes much beyond the mean-field approach. This is the reason why several attempts to study the multifragmentation problem within these approaches have appeared in the literature [4,5,6,7,8].

In this letter we will focus on dynamical properties and fragment formation in dilute systems. Mean-field dynamical simulations are known to exhibit zerosound waves that, in the case of systems prepared inside the spinodal region of the nuclear matter phase diagram, become unstable and may lead to fragmentation [9,10,11,12]. The existence of unstable zero-sound is deeply related to the possibility to observe first order liquid-gas phase transitions. We will study sound properties and investigate realistic scenarios of spinodal decomposition in the framework of molecular dynamics approaches. Then the zero-sound properties can be compared to the response of TDHF calculations to local density perturbations.

We consider a quantal variational derivation of the equations of motion for a many-body system. Let us assume that the many-body state is parameterized by a set of time-dependent complex parameters {Z;}: %p(t) = V>(Z(t)). The equations that determine the time evolution of the system are obtained by extremizing the action /, defined as:

/= /dt < il>\ifij -H\il> >= Idt{C 0 -H), (1)
where £ 0 =< V'N^IV' >• We denote by C = Co -% the "Lagrange" function of the system. Ti. is the expectation value of the total Hamiltonian H, 7i =< tp\H\ip >. The resulting equations of motion for the parameters \Z{\ are the Euler-Lagrange equations: ddCodC dtdZi that, for the most general choice of the complex parameters {Zi}, lead to the Schroedinger equation. We notice that the function £ 0 depends linearly on Z/, so it is possible to expand £ 0 as £ 0 = Yli ff 1 ^') an ^ Eq.( 2) can be rewritten in the form:

T,Au>Zi> = j£, (3)
where the matrix Ai t i< is defined as:

d 2 c 0 d 2 c 0 '•'' dZdzp dZdzp
For a system of identical fermions, the most simple N-body state \ij) > is a Slater determinant built from N different single-particle wave packets. This approximation leads automatically to the concept of mean field and to the use of effective two-body interactions, like in TDHF [1,2].

Introducing the single-particle density matrix />, it is easy to demonstrate that, using a Slater determinant for the state \xj) >, eq.( 3) becomes:

ift^-[*,,] = <), ( 5
)
where h is the mean-field Hamiltonian. This approximation has been particularly successful to study nuclear systems, where, due to the delocalization of the nu-cleon wave packets, the use of mean-field concepts seem to be quite appropriate [1,2].

In FMD or AMD the single-particle wave functions \q n > are chosen to be gaussians: 9n(r, s n , r n ) = e-^-x ^^a "\ s (s n ) X r(r n ), (6) with n = 1,2, ...,./V. s n and r n denote the spin and the isospin respectively.

This implies that, as parameters Z/, we shall consider the time evolution of four complex parameters (the three components of the centroid x n and the width a n )

for each gaussian function \q n >. So we will use the notation Z\ = z n ,i, where i = 1, .,4 indicates one of the complex parameters associated with a gaussian n.

As for the effective interaction, in order to simplify the treatment of the exchange term, we use a local interaction of the Skyrme type:

t>i,2 = {t o + -n + c An)6(ri -r 2 ) o
where n(r) is the diagonal part of the one-body density in the coordinate representation. For spin-isospin saturated nuclear matter this leads to the following expression for the expectation value V of the potential part of the Hamiltonian :

(7) r4i o » + <3nc), o ID n
In numerical simulations the effect of the surface term c An is simulated through the folding of the mean-field potential with a function g, that is taken to be gaussian [2].

As stressed before, the development of instabilities in dilute nuclear systems can be considered as one of the possible mechanisms responsible for fragment formation. It should be noticed that important information about the early dynamical evolution of the system considered can be obtained already in the framework of the linear response approach. Complete calculations, not based on the linear response analysis, will be also discussed below.

In order to establish a connection with the unstable zero-sound properties observed in mean-field calculations of the TDHF type, let us first briefly recall the derivation of zero-sound frequencies in a quantal RPA. According to this theory, the response of the system to the small perturbations 8p of the single particle density matrix around a stationary solution characterized by p 0 , are determined by the linearized equation,

ih^-[ho,6p]-[6U, Po ] = 0 (8)
where h 0 is the mean-field Hamiltonian associated with p 0 , and 8U denotes the fluctuating part of /i[ll]. For infinite nuclear matter, ho is uniform. Hence, the equilibrium single particle density matrix is diagonal in the plane wave representation k. Within this representation, carrying out a Fourier transform with respect to time, we obtain, from eq.( 8), the following quantal dispersion relation for the frequency of the collective mode corresponding to the wave number k [11],

i _ dU(k) I dk' p o (k' -k/2) -p o (k' + k/2) dp J (2TT) 3 u k -e k ' +k/2 + £k'-k/2
where e k = k 2 /(2m) and U(k) is the Fourier transform of the mean-field potential.

The characteristic zero-sound frequencies occur in pairs ±u>;t, and they are real for the stable modes and imaginary for the unstable ones, u>t = i/r^ where ±rî s the characteristic growth or decay time of the mode. Unstable solutions are found in the case of dilute systems, initialized inside the spinodal region of the nuclear matter phase diagram.

Let us perform the same study in the framework of FMD. By linearizing eq.( 3) around a stationary solution (z n< i -z Oni + <5z ni ,) one obtains:

d 2 H d 2/ H n.i.n'.i'lo £<<,,•' = Y^ *" a, lo Sz n',i' + R ~ » |o <$<<,,< (10)
n',t' uz n,i ( J z n',i' uz n,i {JZ n',\' functions are placed on a lattice of a given mesh d, so that the sum over n' in eq.( 10) can be considered to run on the lattice sites. Eq.(lO) does not depend on n anymore, because of the spatial symmetry. The imaginary part of the centroid position and of the width is equal to zero in the stationary state because of the time-reversal invariance. If we agitate the complex centroid position only along one direction, the x direction for instance, indicating 8z n j by the vector (8x n ,6a n ), eq.( 10) leads to the following two coupled equations :

f*<> + K?*<-= £ B^fSxn. + B*; a 8a n , + C;?6xl. + C x n ; a 8a' n , n' n' A a n f6x* n , + A a n f6a' n , = J2 B a n f8x n , + B a n f8a n , + C a n ?8x* n , + C n a ; a 6<,, (11)
where A^ = A n ,, B .,|o, W = ^g^lo and C^" = a^fc

Taking the Fourier transform with respect to time and considering the plane wave representation: 8z n < y i = J2k^zk,ie lk ^d Xn '\ we finally obtain:

i(A x k ' x 8x*_ k + A x k ' a 8a*_ k )u k = B x k x 8x k + B x k ' a 8a k + C x k x 8x*_ k + C x k ' a 8a*_ k i(A a k x 8x*_ k + A a k ' a 8a*_ k )u k = B a k ' x 8x k + B a k a 8a k + C k ' x 8x*_ k + C k ' a 8a'_ k , ( 12
)
where, for instance, A x ' x is the Fourier transform of A^f.

To get information about the evolution of the system, we have to look for eigenvalues and eigenvectors of this system of two complex equations. The four eigenvalues come by pairs of solutions with opposite sign, so that there are two independent solutions for o>|.

In order to provide some quantitative application of the previous analytical derivation let us first consider, for the sake of simplicity, two-dimensional infinite nuclear matter at zero temperature. The parameters of the potential are: to -243.5

MeV fm 2 , t 3 = 846 MeV fm 4 [13].

Since the effective interaction considered does not depend on spin or isospin, it is possible, at zero temperature, to group the single particle wave functions four by four in the initial conditions. Then, this symmetry is preserved along the entire time evolution. In the unperturbed state, for spatial symmetry reasons, all gaussian functions have the same width. In Fig. l we represent the behaviour of the energy per nucleon as a function of the width parameter in the case of a 2D square lattice of mesh d = 5 fm (full line). For a large value of the width, the system really resembles uniform nuclear matter (with homogenous density)

and indeed the energy value approaches the value expected for 2D nuclear matter at the density n = A/d 2 = 0.16 frn~2 (the saturation density for the considered force in 2D is n 0 = 0.55 fm~2). However, the minimum value of the energy is obtained for a finite value of the width: a = 2 fm 2 . This means that the system gains energy by grouping the nucleons into a-like structures. For this value of the width the density exhibits bumps centered on the lattice nodes. The presence of this minimum is observed for all dilute systems.

Let us first consider the case where the width of the gaussian function is kept fixed, like in AMD calculations. In this case the width value is arbitrarilly fixed once for ever. The only time-dependent parameter is then the gaussian centroid.

The electric isoscalar modes are then defined by the unique complex equation for the displacement of the centroid positions of the wave functions:

iA x k x 6x'_ k u k = B x k ' x 8x k + C k ' x Sx*_ k ( 13
)
Then we have only one solution u>l for each wave number k.

Let us study a situation where quantal RPA calculations predict instabilities.

For example, we can consider the mesh d = 3.8 fm, that corresponds to an average density, n = n o /2. The eigenvalues a;*'obtained in this case are negative in the specific range of k compatible with the lattice mesh. Therefore, the system is unstable, in agreement with the RPA predictions. However, the instability time Tjt = i/u k , that is displayed in Fig. 2 as a function of k, is affected very much by the value considered for the initial width. For values of a around 2 fm 2 the instability times are close to the ones predicted by the RPA. However, in such a case, due to the small width value, the observed instabilities can be understood as spinodal instabilities of a gas of a-particles. Indeed, the a-a effective interaction, which can be built by calculating the energy per particle as a function of the lattice mesh, resembles to a Van der Waals force.

Increasing the width of the gaussian functions we observe an important quenching of instabilities. In particular, in the case of homogeneous initial density (large values of the width), instabilities are almost suppressed, in contradiction with the quantal RPA calculations (see Fig. 2). This originates from the fact that in AMD the wave-functions are localized so that the zero-sound can only propagate by an explicit motion of the centroids. However, because of the smoothing of the interaction due to the large width of the wave packets, very weak forces act on the centroids. Then, the motion of the centroids is inhibited. It should be noticed also that waves with wave-lengths shorter than twice the lattice mesh cannot propagate. Conversely, in a full TDHF calculation sound waves can be propagated by the co-operative contribution of several single-particle wave functions.

Moreover, the RPA dispersion relation presents a cut at larger wave numbers k > 0.2 /m" 1 , due to quantal effects.

When the width is kept fixed, like in AMD, it is well-known that particles can experience a spurious dynamics, due to conflicts between the fixed width of the wave packets and the need to have anti-symmetric Slater determinants. In particular, particles of the same type exhibit a spurious scattering. Recently, it has been shown that this may lead to chaotic behaviours also in absence of interaction [START_REF] Drozdz | Chaos induced by Pauli blocking in molecular dynamics approach[END_REF]. This problem, that could affect the results presented above, is removed in calculations of the FMD type, by letting also the widths of the wave packets evolve. In this case, in order to have a stationary state as initial condition for the typical mesh size of d = 3.8/ra, we shall consider the value a -2 fm 2 since it corresponds to the minimum of the energy. However, we will also study the case a = 1 fm 2 that corresponds to a quasi-stationary state since the derivative of the energy is small for large values of the width (see Fig. l).

Large widths mimic also stationary states for more disordered systems (where the spin/isospin symmetry is explicitly broken), as we will see later.

With the width as a time-dependent parameter, we get two solutions for u>l associated with different motions in which the centroids and the widths are coupled. In Fig. 2 we display the solution associated with the eigenfunction that corresponds mostly to a displacement of the centroid positions (thick lines). We observe that the inclusion of the width as a dynamical parameter slightly modifies the AMD dispersion relation. We conclude that, also in the framework of FMD calculations, instabilities are suppressed when the stationary width is large, in disagreement with the predictions of quantal RPA calculations.

As illustrated in Fig. s 1-2, when a spin/isospin symmetry is considered, the system condensates into a gas of a particles, the nucleon wave packets taking a rather small width, compatible with the a ground state. In such a case the system does present zero-sound waves. However, this corresponds to a rather peculiar situation. From a physical point of view one would expect that in a nuclear collision, even if in the ground state the partners present an a-cluster structure, this should be partly destroyed during the evolution by both the residual interaction (such as the spin-orbit and the Coulomb interaction) and the strong agitation expected during the dense phase of the collision.

In the case of a nucleon gas, the minimum energy is reached asymptotically for large values of the width a (see the dashed line in Fig. l). Therefore, one can conclude that in an excited situation the natural tendency of the nucleons will be to present large width. This is in fact the well-known effect related to the the wave packet spreading. Then, according to our discussion about dispersion relations, one should expect an important quenching of the spinodal instabilities when the width is considered as a dynamical variable (FMD).

These conclusions are confirmed by direct 3D simulations. We have performed calculations for dilute systems at different excitation energies. We present here the results of these simulations for a spherical source of 160 nucleons. The parameters of the force are t 0 = -1033. MeV fm 3 , t 3 = 14687.5 MeV fm 6 . The potential has been folded with a gaussian function of width g = 0.4 fm. The positions of the various nucleon centroids are randomly distributed inside a sphere of radius R -8 fm so that the average density is about n o /2 (no = 0.16 /m~3). In the results presented here we have given to the system a small self-similar expansion velocity. We consider, as initial width of the gaussian functions, two cases: a = 2 fm 2 , that corresponds to an energy per nucleon of 2 MeV and a = 4 fm 2 , that corresponds to E/A = -2 MeV. In Fig. 3 we show the time evolution of the sources considered in these two examples.

In both cases we observe that the early dynamics is dominated by an increase of the width. This is the behaviour expected for a disordered system since the minimum energy is obtained for large widths (see Fig. l). The more excited system undergoes a vaporization into individual nucleons, while the other system collapses into a single source, evaporating few nucleons from its surface . Almost no clusterization is observed. This absence of spinodal instabilities in the simulations confirms the conclusion drawn from the linear response treatment. Because of the appearance of large widths the residual interaction between centroids is smoothed out so that spinodal instabilities are quenched.

We have also performed 3D calculations considering the a clustering in the initial conditions. In this case the dynamical evolution is dominated by the oscillation of the gaussian width around the equilibrium value (2 /m 2 ) for a particles. When, depending on the initial conditions, these oscillations are not too large, so that the widths remain small during the entire evolution, the formation of "fragments" (that are nothing but clusters of a particles) is observed (see Fig. 4).

As discussed before, this mechanism resembles the spinodal decomposition of a gas of a particles. Finally it should be noticed that, if the widths are kept fixed (as in AMD) at a rather small value (like a = 2 fm 2 ) we also observe the occurrence of instabilities and cluster formation.

In conclusion, we have investigated the dynamical properties of widely used quantum molecular approaches such as AMD and FMD. In particular, we have carefully studied zero-sound properties and compared to the well-known features of mean-field approaches. This is of particular interest since unstable zero-sound waves may play a role in multifragmentation scenarios (like in spinodal decomposition). Performing a linear response treatment, as well as complete 3D simulations, we show that zero-sound properties and fragment formation can be significantly affected by the gaussian ansatz, especially when the widths of the wave packets are large. Indeed, in this case the single particle wave functions might be so broad that all the dynamical features expected at the TDHF level are washed-out.

In the case of a a-particle gas we observe the occurrence of spinodal decomposition essentially because the widths of the nucleon wave packet remain small, as nucleons are bound in the a's. Conversely, in the case of a disordered nucleon gas, if the width is considered as a dynamical parameter, we observe that large widths are generated, during the evolution, in the bulk of nuclear matter, leading to the suppression of unstable zero-sound waves and spinodal decomposition.

Fig.l:

The energy per nucleon as a function of the gaussian width a calculated for a 2D lattice of a-particles (full line) and individual nucleons (dashed line) with an average density n = 0.16 fm~2. The energy associated with a homogeneus medium is equal to -6 MeV per nucleon (dotted line).

Fig. 2 :

 2 Fig.2:The dispersion relation obtained in AMD calculations for a 2D lattice of a-particles. The mesh size of the lattice is equal to 3.8 fm. Results associated with the width value a -2 fm 2 (dot-dashed line), 4 fm 2 (dashed line), 7 fm 2 (dotted line) are presented. The thick lines correspond to the results obtained in FMD. The full line represents the results of a quantal RPA.

Fig. 3 :

 3 Fig.3: Time evolution of a dilute spherical source of 160 nucleons, as obtained in 3D FMD calculations, having an initial width value of 2 fm 2 (top figure) and 4 fm 2 (bottom part). The dots correspond to the centroid positions of the individual nucleons, projected on a plane, and the circles display a typical gaussian width.

Fig. 4 :

 4 Fig.4: Time evolution of a dilute source of 160 nucleons initialized considering a-clusters. The initial value of the widths is taken equal to 4 fm 2 .

Acknowledgements: We warmly thank H.Feldmeier, A.Ohnishi, A.Ono, J.Randrup and J.Schnack for helpful discussions. This work was supported in part by the Commission of the European Community, under Contract No. ERBFMB I-CT-960654.

""This work was supported by the Commission of the European Community, under Contract No. ERBFMB I-CT-960654.