Cluster formation in nuclear reactions from mean-field inhomogeneities

Perturbing fluids of neutrons and protons (nuclear matter) may lead, as the most catastrophic effect, to the rearrangement of the fluid into clusters of nucleons. A similar process may occur in a single atomic nucleus undergoing a violent perturbation, like in heavy-ion collisions tracked in particle accelerators at around 30 to 50 MeV per nucleon: in this conditions, after the initial collision shock, the nucleus expands and then clusterises into several smaller nuclear fragments. Microscopically, when violent perturbation are applied to nuclear matter, a process of clusterisation arises from the combination of several fluctuation modes of large-amplitude where neutrons and protons may oscillate in phase or out of phase. The imposed perturbation leads to conditions of instability, the wavelengths which are the most amplified have sizes comparable to small atomic nuclei. We found that these conditions, explored in heavy-ion collisions, correspond to the splitting of a nucleus into fragments ranging from Oxygen to Neon in a time interval shorter than one zeptosecond (10$^{-21}$s). From the out-of-phase oscillations of neutrons and protons another property arises, the smaller fragments belonging to a more volatile phase get more neutron enriched: in the heavy-ion collision case this process, called distillation, reflects in the isotopic distributions of the fragments. The resulting dynamical description of heavy-ion collisions is an improvement with respect to more usual statistical approaches, based on the equilibrium assumption. It allows in fact to characterise also the very fast early stages of the collision process which are out of equilibrium. Such dynamical description is the core of the Boltzmann-Langevin One Body (BLOB) model, which unifies in a common description fluctuations in nuclear matter and the disintegration of nuclei into nuclear fragments.


Introduction
The most catastrophic process which can occur in a nuclear complex is its splitting into clusters and fragments when undergoing a violent external action. We want to address this process, which can be probed in a dissipative heavy-ion collision, from the point of view of dynamics, moving from nuclear matter to nuclei, which are finite open self-bound systems. Nuclear clusterisation may appear in various forms. Exotic topologies in nuclear astrophysics (stellar matter), nuclear states explored in low-energy reactions from the recombination and vibrations of existing cluster structures, low-density regions of the equation of state landscape, where clusterisation in intermediate-mass fragments and clusters may arise and can be explained in violent nuclear reactions. Finally, clusterisation characterises in general Fermi liquids as arising from ripples produced by phase-space fluctuations [1].
We constructed a microscopic dynamical framework from applying the theory of Fermi liquids to clusterisation in nuclei (see ref. [2] for a more extended and detailed discussion); within this framework, we explore how the clusterisation progresses from zero-sound propagation. Two avenues are mostly followed to describe nuclear processes, molecular dynamics and mean-field approaches [3]. In mean-field approaches particles evolve independently in their own selfconsistent mean-field. Such treatment, while well suited to describe the collective behaviour, can not handle clusterisation in its pure mean-field form. On the other hand, molecular-dynamics approaches rely on a description in terms of product wave functions. The independent-particle scheme (i.e. mean field) is applied to nucleonic single-particle wavefunctions ψ i , so that manybody correlations in both mean field and scattering can be achieved from the localisation of ψ i (a coherent-state subspace is used in FMD [6] and an even stronger localisation is imposed in AMD [4,5] by also fixing the widths of ψ i ). Such treatment is successful for final-state correlations, but it may approximate collective behaviour and 0-sound propagation.
To profit from a well suited description of 0-sound propagation, we follow thereafter a meanfield approach where, in order to describe clusterisation, extensions to handle large-amplitude dynamics should be explicitly introduced. In particular, correlations beyond the level of kinetic equations are needed or, in terms of BBGKY hierarchy [7], upper orders beyond two-body correlations should be recovered in order to access highly non-linear regimes. A technique to obtain such result in an approximated form consists in applying two treatments in parallel: first, introducing nucleon-nucleon correlations continuously, second, handling a stochastic ensemble of several mean-field trajectories to exploit the introduced correlations. In this case, a mean-field trajectory ρ 1 taken at a given time, would split at a successive time into subensembles ρ (n) 1 : In practice, at the level of kinetic equations, not only collisional correlations should be introduced, but also a term of vanishing mean which injects fluctuations around the collision integral intermittently in time. Such term can be exploited as a stochastic source to revive fluctuations all along a dissipative process. Adding the collisional correlations introduced above and fluctuations to the mean-field produces a scheme which resembles stochastic TDHF [8,9]: whereĪ (n) coll and δI (n) coll are the average collision contribution and a continuous source of fluctuation seeds, respectively. The corresponding Wigner transform yields the Boltzmann-Langevin (BL) equation [10], in terms of an ensemble of distribution function f n , where the ensemble f (n) replaces the ensemble of Slater determinants in Eq. (2) and corresponds to a Fermi statistics at equilibrium, h (n) is the effective Hamiltonian acting on f (n) , and the residual contributions of Eq. (2) are replaced by corresponding Uehling-Uhlenbeck (UU) terms. We focus on the BL equation and we solve it in full phase space through the BLOB approach [11,12,2], where nucleon-nucleon correlations are introduced by constructing a fluctuating collision term which acts on extended equal-isospin phase-space portions for each single in-medium collision: where g is the degeneracy factor, W is the transition rate in terms of relative velocity between the two colliding phase-space portions, and F handles the Pauli blocking of initial and final states over their full phase-space extensions. The extended size of phase-space portions involved in the scattering in Eq. (4) should be large enough so that the occupancy variance in h 3 cells is equal to f (1 − f ), which corresponds to the scattering of two nucleons [13]. For comparison, we consider a second simplified approach to solve the BL equation, Eq. (3), based on the SMF [14] treatment, where fluctuations are injected from an external stochastic contribution U ext and projected on spacial density. In both cases, in all following calculations, a simplified SKM* effective interaction, with momentum dependence omitted, is used [15,16]. A soft isoscalar equation of state with a compressibility modulus k = 200MeV is used. For test, both asy-stiff (linear) and asy-soft (quadratic) parametrisations are used for the symmetry energy. All nuclear-matter calculations employ a periodic box of edge L = 39fm, and the system is initialised with a Fermi-Dirac distribution at a temperature T = 3MeV. In the collision term, either a free or a constant nucleon-nucleon cross section σ NN is used.
Before carrying on a study on heavy-ion collisions, we require that fluctuation amplitudes are consistent with analytic expectations from Fermi liquids: for this purpose, we study nuclear matter in initially homogeneous conditions.

Fluctuations in two-component nuclear matter
While unperturbed, the dynamics of a periodic portion of uniform nuclear matter at low temperature (so that the collision term I UU in the analytic description can be neglected) would evolve along a mean trajectory f 0 . Let us introduce perturbations δf q f 0 in two forms: either neutrons and protons move in phase (isoscalar perturbation, indexed with q = s), so that In a periodic box, this action introduces fluctuations ρ q k associated with plane waves of wave number k, with a corresponding equilibrium variance (σ q k ) 2 (or intensity of response), and with an equilibrium variance of spacial density correlations (σ ρ q ) 2 (this latter variance is obtained from the former through an inverse Fourier transform).
If the Boltzmann-Langevin equation is applied to the disturbance δf q , the following scheme can be established.
• First, if fluctuations are of stable nature, the two equilibrium variances (σ q k ) 2 and (σ ρ q ) 2 can be related to the free-energy density curvature F q (k) through the fluctuation-dissipation theorem: where T is the temperature and ∆V a volume cell in configuration space where (σ ρ q ) 2 is calculated.
An application of such scheme to nuclear matter is well suited to investigate the equilibrium variance of spacial isovector density correlations (σ ρ v ) 2 as a function of the symmetry energy. From an analogous application to open system, the isotopic distributions of clusters emerging from density ripples can be analysed. • Second, in unstable conditions, both initial fluctuation seeds and intermittent fluctuation seeds injected at later times can yield an exponential growth of the intensity of response (σ q k ) 2 over a growth time τ k [17,18]. An application of such scheme to nuclear matter is well suited to sample zero-sound propagation and investigate the instability growth rates Γ k = 1/τ k . Analogous conditions in open system are suited to investigate the arising of clusterisation.
In the following, these steps are attended to.

Relating isovector fluctuations to the symmetry energy
To undertake a simulation on the propagation of isovector fluctuations, it is convenient to isolate them from the overwhelming effect of isoscalar fluctuations; those latter would in fact completely dominate the dynamics with their larger amplitude. Scalar terms are therefore suppressed in the potential so that stable conditions are imposed at any density ρ 0 .
In this case, the fluctuation-dissipation relation Eq. (5) takes a form where the isovectordensity variance is related to an effective symmetry energy, proportional to the symmetry energy E sym at zero temperature [19].
This relation can be solved numerically by calculating the isovector variance (σ ρ v ) 2 in cells of edge size l, where l should be chosen large enough to minimise surface effects and better sample the volume symmetry energy (this value is found around l = 2fm).
asy-soft asy-stiff Figure 1. SMF calculation of the effective symmetry energy, scaled by N test , compared to the analytic expression of the symmetry energy for different system densities. The collision term is either activated (and a free σ NN is used) or suppressed (σ NN = 0). Two parametrisations of the symmetry energy, either asy-stiff (left panel), or asy-soft (right panel) are used.
In a first calculation, shown in Fig. 1, where the simplified SMF approach to solve the BL equation is used, even though the evolution of the effective symmetry energy F v eff as a function of density reproduces correctly the shape of the analytic expression for the symmetry energy E sym , the former and the latter differ of a factor equal to the number of test particles per nucleon N test employed in the numerical sampling of the mean field: The same result is obtained with or without the contribution of the residual term I UU , indicating that explicit isovector terms are missing in the isovector channel in building fluctuations. This is not surprising because the residual term is not defined to build nucleon-nucleon correlation and the corresponding fluctuations. Those latter, even though reflecting the potential employed in the calculation, develop from initial fluctuation seeds which are related to the numerical noise in the sampling of the mean-field, which, on its turn, is related to N test . A second calculation, shown in Fig. 2, employs the BLOB approach, where the residual term is designed to explicitly introduce phase-space fluctuations from nucleon-nucleon correlations. In this case, the isovector variance is enhanced, but equilibrated nuclear matter is not a favourable condition to completely recover the large isovector variance of the analytic expectation. In particular, at small density ρ 0 collisional correlations are ineffective, at larger density collisions become rare and can hardly revive fluctuations and, in addition, the same noise issue already 0 0.04 0.08 0.12 0 Stable symmetric matter, asy-stiff, Figure 2. Isovector variance at equilibrium as a function of the system density for an asy-stiff parametrisation of E sym . Calculated with BLOB, with different values of σ NN , and compared with SMF, with the collision term either activated (with a free σ NN ) or suppressed. For the system density 0.08fm −3 , the inset shows the effect of varying N test (see text).
observed in the SMF calculation is acting as a smearing contribution [17] with a dependence on N test .
For test purposes, it is instructive to appreciate the effect of increasing the nucleon-nucleon cross section and, as a consequence, intensifying the nucleon-nucleon collision rate. This action enhances the effect of the BLOB fluctuating term in reviving nucleon-nucleon collisions, so as to prevail over the smearing effect of the mean-field noise.
A more consistent (but numerically costly) approach is reducing the numerical noise by refining the mean-field paving: this action, consisting in simply increasing N test , also improves the results with respect to the SMF approach.
From these results and trends, we can expect that the limits imposed by the conditions of equilibrated nuclear matter could be actually overcame in out-of-equilibrium conditions, where nucleon-nucleon collision rates are definitely higher: this situation corresponds to the early stages of heavy-ion collisions, which are investigated in Sec. 3.1.

Relating zero-sound propagation to instability growth
To undertake a simulation of zero-sound propagation of collective modes in nuclear matter, the system should be prepared as initially homogeneous and at low temperature. If the temperature increases, a transition from zero to first sound may in fact occur [20,21]. If unstable conditions are investigated, an amplification of the fluctuation amplitude is expected to trigger a catastrophic process. Early times are therefore better suited to compare to analytic expectations. If the zero-sound propagation is correctly sampled, isoscalar fluctuations should develop spontaneously with the correct amplitude, inducing the arising of a mottling pattern.
The analytic expectation can be defined in a linear-response approximation, assuming small deviations from a mean trajectory f 0 ; in this case, from the linearised Vlasov equation (where residual terms are suppressed) expressed as a function of the disturbance wave number k and frequencies ω k , the corresponding dispersion relation [22,23] can be extracted by applying selfconsistency, as follows: of the arrow) is restricted to the Fermi surface so that, introducing the Landau parameter F 0 (k) = (3/2)(ρ 0 / F )∂ ρ U k , the dispersion relation can be written in terms of sound velocity s = ω k /(kv F ) [18,24] 1 As shown in Fig. 3, instabilities in zero-sound conditions can be tracked by placing the system in a location of the equation-of-state landscape where the incompressibility χ −1 ≡ ρ ∂P ∂ρ is negative. In these conditions, corresponding to the spinodal instability and to F 0 (k = 0) < −1, the dispersion relation yields imaginary roots [25] which, by replacing s → iγ, can be put in the Eq. (10) shows that disturbances of wave number k get amplified with a growth time τ k and a corresponding growth rate Γ k = 1/τ k . As shown in the analytic calculations of Fig. 4 (left panel) for zero temperature, the response intensity at zero sound should present the following evolution of the growth rate Γ k as a function of the k number, or with the corresponding wavelength λ.
• For small k: Γ k tends to increase (decrease) linearly with k (λ) because the more matter has to be relocated, the longer time it takes. • For large k: small wavelengths λ are excluded as a function of the interaction range. This ultraviolet cutoff can be described as a Gaussian smearing [27,26] σ of the mean-field potential in configuration space U ⊗ g(k) with g(k) = e − 1 2 (kσ) 2 which reduces Γ k with k. • The combination of these opposite behaviours produces a maximum which corresponds to the fastest growing disturbance, i.e. the leading k mode and leading wavelength.
To explore the equation of state (especially in applications to heavy-ion collisions), it is convenient to have a finite value of T . Still, as mentioned, T should be low in order to avoid transitions to first-sound. It is possible to introduce a finite temperature T rather than working at T = 0 through a low-temperature expansion of the chemical potential ν. A reduction of Γ k with T is the result of this action.
These modifications to take into account the ultraviolet cutoff and to introduce a finite temperature impose to replace F 0 (k) by an effective Landau parameter [27] F 0 (k, T ) = (µ(T )/ F )F 0 g(k) in eq. (10) in order to obtain the analytic response in Fig. 4 for a finite temperature.
The corresponding numerical calculations employ the BLOB approach to extract the response intensity σ 2 k (t), i.e. the amplitude of the isoscalar fluctuation of a mode k from the Fourier transform of the space density (see ref. [2] for details). Fig. 5 studies the evolution of the ratiõ σ 2 k (t) = σ 2 k (t)/σ 2 k (t = 0). From the analysis of the response intensity we infer some general results which can be transposed from nuclear matter to heavy-ion collisions. Average response intensity as a function of k (left column) and time t (right column), calculated with BLOB in conditions of spinodal-instability (top row) and Landau damping (bottom row). The different curves in the left column correspond to different time instants, while the different curves in the right column correspond to different k intervals (their colours corresponds to the coloured segments in the left panel).
• in unstable (spinodal) conditions, calculated at ρ 0 = 0.053 fm − 3 fluctuations get rapidly amplified as expected. • Even outside of the spinodal instability, in calculations at ρ 0 = 0.14 fm − 3, which correspond to Landau-damping conditions, the response intensity reaches significant amplitudes. • The response intensity saturates rather early due to the combination of large k into small k. In heavy-ion collisions, a similar effect manifests in the recombination of small emerging clusters into larger fragments and in an overall mean-field resilience effect (investigated in detail in ref. [28]).
To the analytic expectation for the growth-rate as a function of k, Fig. 4 (right panel) overlaps the corresponding numerical calculation obtained with BLOB for the same conditions and interaction properties, when taking the average over several stochastic dynamical paths, The achieving of a close correspondence between the BLOB calculation and the analytic expectation, as shown in Fig. 4 ensures that the model succeeds in describing the fluctuation phenomenology consistently and, in particular, fluctuations develop with the correct amplitude.
After showing that this correspondence is achieved, some quantitative results can be extracted from the calculation. The leading wavelength ranges between 8 to 9 fm. By analogy to nuclear matter, we can expect for heavy-ion collisions that fragments and clusters should arise in the region of Neon. The corresponding separation time in heavy-ion collisions should not only take into account the growth time of the leading k modes, calculated in Fig. 4, but also the time needed in the collision to generate a low-density phase where instabilities may develop, and the overall effect of the kinematics in an open system.

Dynamics of clusterisation in open systems, applications to heavy-ion collisions
The composition of the two schemes presented at the beginning of Sec. 2 and applied to open systems provides an analysis of the isovector and isoscalar properties of clustering in stable and unstable condition, in connection with analytic expectations from Fermi fluids.

Spinodal clusters from density ripples at Fermi energies
An analysis which directly corresponds to the above calculations in nuclear matter is dedicated to the properties of emerging clusters in an open system at low density, produced in a heavy-ion collision. Even before that clusters separate into well identified blobs of matter in configuration space, we can track the evolution of potential ripples developing in the bulk since early times. In analogy to independent fragments, we can attribute to those inhomogeneities a mean radius and a density averaged over the corresponding potential well in order to obtain mass and element numbers A , Z . An indicative expectation for isotopic yield distributions may be obtained by analogy to eq. (5) as These yields of forming fragments reflect the isovector fluctuations. As already argued in Sec. 2.1, while in equilibrated nuclear matter fluctuations are hard to revive and entertain due to low collision rates, in open systems fluctuations are initially built out of equilibrium, so that collision rates are larger. As a consequence, the isovector fluctuations variance results more consistent with analytic expectations, as illustrated in Fig. 6 (in the frame) for the system 136 Xe+ 124 Sn at 32 AMeV (central collisions), where the isotopic distribution of forming fragments is calculated with BLOB and compared to the expectation of eq. (12).  Figure 7. The isotopic properties of fragments produced in central collisions in 136 Xe+ 124 Sn at 32 AMeV is shown to evolve according to a process of isospin distillation. The average isospin content in potential ripples as a function of ρ well is calculated with BLOB before the onset of the fragmentation process at t = 50fm/c, or when the process is almost achieved at t = 130fm/c. Fig. 7 extends the study of the isospin content of forming clusters and fragments to the whole range of sizes produced in the process, and illustrates the corresponding average isospin content at a very early time (t = 50fm/c), when the system is still rather homogeneous, and at a later time (t = 130fm/c), when inhomogeneities have been built, as a function of the local density ρ well associated to the density ripple. The larger neutron enrichment of the more volatile phase signs the onset of a distillation process [24].
The right panel of Fig. 6 extends the survey of fragment production to a large incident-energy range, examining the multiplicity of fragments with Z > 4. In this case, fragments are studied at t = 300fm/c, when they are completely separated but still highly excited, and at the end of the full decay sequence (where the Simon transition-state model [29] was employed to complete the decay with sequential evaporation). The BLOB+Simon calculation yields mean values and widths which are in agreement with experimental data from INDRA [31,30]. The transition from fusion to multifragmentation occurs below 30AMeV and the dominant spinodal mechanism is gradually replaced by vaporisation beyond 45AMeV [11]. On the other hand, the simplified BL approach to introduce fluctuations used in SMF results less efficient in building fluctuations, and does not succeed in describing the transition from fusion to multifragmentation at the correct incident energy.

Light systems
heavy-ion collisions and nuclear matter involve anditional clustering processes, where alpha and light charge particles are formed. Such mechanism, different from dynamical instabilities, would require the explicit inclusion of further correlations in the present formalism. Light charged particles related to nuclear clustering have in fact too small size, exceeding the ultraviolet cutoff of the dispersion relation, so that they can not belong to the unstable multipole modes which characterise spinodal fragmentation. Solutions for an explicit treatment of cluster formation are proposed in refs. [32,5].
The ability of the stochastic approach of BLOB in handling fragment formation and introducing correlations bejond the level of kinetic equations can be considered a first approximation to light clusters, which are essentially emerging from dynamical fluctuations. Blob already takes charge of a fraction of the light cluster production and a sequential-decay afterburner (Simon [29] was used) adds the missing fraction from the decay of heavier elements. A comparison with data from ref. [33] is presented in Fig. 8 for the system 12 C+ 12 C at 62 AMeV.
At this level of comparison, the production rates of the light masses in the isobaric distribution seems consistently described. Such production is fed by the most central impact parameters and fades with periferal configurations, which are related to the heavier side of the mass distribution.   Figure 8. Clusters and fragment yields in 12 C+ 12 C at 62 AMeV, divided by the geometric cross section, are shown as calculated with BLOB till the last registered split in a time window of 80 to 140 fm/c (BLOB, scaled by 1000) and after the full decay sequence, for which the model Simon [29] was employed (BLOB+Simon). A free σ NN and an asystiff form for the symmetry energy are used. Experimental yields from ref. [33] for some isotopes named on the plot, are indicated with dots; data are normalised so that the sum of d and t yields equals the corresponding calculated quantity. Squares are calculations corresponding to the isotopes that have been measured.

Fluctuations and transparency at intermediate energy
The connection between fragment production and hydrodynamic properties like the characterisation of the flow [34] at intermediate energies is strictly related to the treatment of fluctuations. At intermediate energy the inclusion of fluctuations has two antagonist effects: on the one hand, it enhances the fragmentation of the system, on the other hand it reduces the directed flow. This effect can be studied in the comparisons of Fig. 9, for the collision 197 Au+ 197 Au at 100 AMeV for an impact parameter of 7fm. The simulation is performed with three approaches, SMF without and with a collision term (constant σ NN = 40mb) and BLOB (also withσ NN = 40mb) , using identical parameters for the mean field as defined in ref. [35]. The SMF approach describes the outward deflection of the trajectory imparted by the directed flow, which is absent in the SMF description without collision term. The BLOB approach exhibits a reduced directed flow with respect to SMF, because it competes with the production of fragments and clusters. This latter, due to the Langevin fluctuations, results in a large variety of very different fragment configurations; two of those are shown, one where the fragmentation of the quasi-target and the quasi-projectile is observed (bottom row in the density map, left), the other where the emitting source is situated at midrapidity (bottom row in the density map, right).
A quantitative study of the flow is illustrated in the right panel of Fig.9. For a simulation where the collision rate is identical (due to using the same constant σ NN ), the larger fragmentation rate is reflected in a smaller slope for the directed flow as a function of reduced rapidity. This is a rather general example to spots the main difference between a simplified description of fluctuations and a BL approach solved in full phase space, when applied to intermediate energies.

Conclusions
The process leading a fermionic system, like nuclear matter or violent nuclear reactions, to separate into clusters and fragments is reviewed in some significant steps Among more possible strategies to address the problem (see introduction), we have chosen to extend a mean-field description to include fluctuations in a Boltzmann-Langevin framework [2] and suggested some scheme for testing the transport approach on analytic expectations.
The approach is particularly adapted to describe mean-field inhomogeneities, like the spinodal decomposition, but also conditions like Landau damping and other collective behaviours. It reduces to an approximation when the formation of the light clusters should be described because no explicit cluster contributions are so far treated. For larger clusters and intermediate-mass fragments in general, the approach yields quantitative predictions for the production of fragments in a broad range of situations and mechanisms encountered at Fermi energies [11,36,37]. To complete the survey, a further application to relativistic spallation has been undertaken elsewhere [28].