Spin Quantization in Heavy Ion Collision

We analyzed recent experimental data on the disassembly of $^{28}$Si into 7$\alpha$ in terms of a hybrid $\alpha$-cluster model. We calculated the probability of breaking into several $\alpha$-like fragments for high $l$-spin values for identical and non-identical spin zero nuclei. Resonant energies were found for each $l$-value and compared to the data and other theoretical models. Toroidal-like structures were revealed in coordinate and momentum space when averaging over many events at high $l$. The transition from quantum to classical mechanics is highlighted.


I. INTRODUCTION
Recent experimental data [1] have shown evidence of resonances in the disassembly of the 28 Si nucleus into 7α. The data were obtained from the collision of a 28 Si beam at 35 MeV/A on a 12 C target, the experiment was performed at the Cyclotron institute, Texas A&M university. The authors of [1] tentatively associated these structures to the population of toroidal high-spin isomers as predicted by a number of theoretical models [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20]. In particular the experimental analysis concentrated on the disassembly of the projectile nucleus into α-like particles. The data show to a high degree of confidence some structures at excitation energies 114, 126 and 138 MeV, respectively, close to the predicted toroidal state at 143 MeV [12]. These encouraging results call for more experimental and theoretical efforts to uncover these resonances also for different nuclei, different disassembly routes and as a function of excitation energy. Due to the dynamics involved in the disassembly, microscopic models such as the Anti-symmetrized Molecular Dynamics (AMD) model [21] or Constrained Molecular Dynamics (CoMD) model [22] could be used, but they may become numerically difficult to handle when a large number of events is needed. Furthermore, they may not be able to describe in detail the α-like events as selected from the data [1]. Hybrid models may help in overcoming numerical problems at the expense of some physical insights [23].
It was observed already in the 1930s that α-like nuclei ( 12 C, 16 O,...) [24][25][26][27][28][29] display many properties that can be easily explained by assuming that those nuclei are made of α particles with no internal structure. Inspired by these many works we implemented a dynamical model where α particles interact through suitable chosen two body forces. We enforced two main variations with respect to what can be found in the literature. The first is the α-α interaction. For simplicity we used the phenomenological Bass potential (for A = 4) widely used for low energy heavy ion collisions [30]. This potential is very attractive at short distances thus the particles strongly overlap, overcoming the Coulomb repulsion. As a second ingredient of the model, we treat α particles as Gaussian distributions with widths given by their radius. Overlapping particles experience Pauli blocking because of the internal structure of the αs. Thus, we include a repulsive effect due to the increase of the Fermi energy opportunely adjusted to take into account finite size effects [31]. With such simple assumptions we are able to reproduce the binding energies of even-even N = Z nuclei up to mass 104 with less than 5% discrepancy to the experimental data. We do not want to stress much the properties of the model since our main goal is to simulate the dynamics of the disassembly to compare to data at least qualitatively. Furthermore, we would like to confirm or disprove the existence of exotic unstable shapes using a simple and transparent model and hope to be of guidance for future experiments. We dubbed the model the hybrid α-cluster model (hαc). It is a semi-classical model since it includes Pauli blocking effects. In fact the model ground states display strongly overlapping α-particles and a strong repulsion due to the increase of the Fermi energy. It means that ground states cannot be described by α-point particles and the nucleons degrees of freedom are essential. Systems excited by some external probes expand and the α-degrees of freedom may become dominant. Notice that since the repulsion is due to the Pauli blocking and the Coulomb potential, heavy ion collisions using this model can be simulated at energies above the Coulomb barrier up to maybe 80 MeV/A as we discuss below. At higher energies we may need to introduce a suitable collision term, which is a task to be discussed in the future. We introduced another quantum effect in the initial conditions, i.e., we give to the nucleus at time zero a quantized angular momentum l = l z along the z-axis. We assume that the angular momentum is transferred in the collision of 28 Si and 12 C, or 28 Si and 28 Si. The substantial difference between the two systems is that only even-l values in the entrance channel are allowed in the latter case. Changing the initial angular momentum revealed a wealth of model features ranging from a first order phase transition of dynamical origin to the formation of short living toroids when averaging over events. Due to its simplicity and numerical affordability, we can make prediction to be tested in future experiments.

II. THE HYBRID α-CLUSTER MODEL
In our model α-degrees of freedom are treated explicitly while nucleon (protons and neutrons) degrees of freedom are treated implicitly hence the hαc acronym. The interaction between the α-particles is given by the Coulomb repulsion (in the monopole-monopole approximation for simplicity) and the nuclear attraction. The latter is approximated as V αα = V Bass (A = 4), i.e., the Bass potential for mass A = 4 nuclei [30]. Coulomb repulsion is not sufficient to prevent a strong overlap among α-particles. Overlapping nuclei increase the repulsion due to the combined action of the Pauli principle and Heisenberg uncertainty principle, in particular the Fermi energy (per α-particle) is given by: MeV is the average kinetic energy of infinite nuclear matter, the factor of four takes into account the fact that we are dealing with α and not nucleons. For small nuclei corrections are needed to take into account finite size effects, which reduce the Fermi energy thus the parameter x F . In ref. [31], Equation 5.5, such correction was discussed for medium light nuclei resulting in x F = 0.65 for 8 Be. For overlapping α-clusters only one nucleon in one α particle is identical to another nucleon in the other α. This parameter takes into account the fact that the Heisenberg uncertainty principle is at play as well for non-identical nucleons [32]. The overlap between α-particles can be described as Gaussian distributions with standard deviation proportional to the α-radius r α = r 0 4 1/3 : (2) The parameter β = 1.22 is fitted to reproduce the binding energy of 12 C, it is the only free parameter entering the model if we exclude the radius parameter r 0 . The value of r 0 has some consequences regarding the moment of inertia, which we discuss below. For the purpose of this paper we use r 0 = 1.15 fm, unless otherwise noted, similar to the parameters entering the Bass potential [30]. At maximum overlap ρ = 2 and ε F Nα = 86.7 MeV which is the maximum repulsion in the two body channel to compare to the nuclear attraction V αα (r = 0) = −58 MeV. This implies that colliding nuclei will become transparent at beam energies well above the Coulomb barrier, similar to Time Dependent Hartree-Fock calculations [33]. A suitable collision term may remedy this shortcoming but it is outside the purpose of this work [34]. Equations (1) and (2) give the repulsion between particles and we treat it as a classical two-body force.
The classical Hamilton equations of motion for interacting α-particles are solved numerically using the O(dt 5 ) Runge-Kutta method, dt = 1 fm/c is a typical time step used in the calculations. At the highest excitation energies or angular momenta discussed in this paper, the particle velocities become very large thus we implemented relativistic kinematics. This correction is important but we stress that the description in terms of classical interactions is still valid. To obtain the nuclear ground states and their binding energies, the equations of motion were solved adding a friction force until a minimum and stable configuration is reached. The particles position are saved on a file and used as initial positions in dynamical simulations. To generate events, the initial positions are rotated randomly for each event and/or many different ground states are generated.
In Figure 1, we plot the binding energies of α-cluster nuclei as function of the mass number A. The free parameter β of the model was fixed to reproduce the 12 C binding energy. This leads to an overestimation of the binding energy of 8 Be of 2.6 MeV: 59.1 MeV theory vs. 56.5 MeV experiment [35]. This is an important feature since fixing the free parameter to the binding energy of 8 Be would lead to a general underestimation of all the other nuclei. It implies that the α-particles must be more overlapping for heavier nuclei, thus the increase in Fermi energy. It confirms our discussion above that the correct description of nuclear ground states must be in terms of nucleonic degrees of freedom while α-clusters may dominate at lower densities, i.e., in the expansion stage of the nuclear dynamics. Our hybrid model reproduces the binding energies of nuclei up to mass 104 with an error less than 5%. In Figure 1, we have included for comparison the contribution to the binding due to the α binding energy, full diamonds. We notice that changing the value of r 0 to 1.26 fm produces a similar agreement to the binding energies with β = 1.02. Thus, these data are not able to constrain the parameter values to high degree and we will investigate fusion cross sections of even-even N = Z nuclei for further constraints.
Once the initial conditions are found, we can generate many initial ground states to be used as initial conditions in dynamical calculations. We treat each α-particle as a Gaussian distribution normalized to one of radius (variance) r α . In Figure 2, we plot the density averaged over ensembles at two different times. Naturally the system is stable and the central density is rather reasonable. The displayed system is 28 Si and we are going to concentrate on this nucleus for the remainder of this paper since it was carefully investigated in ref. [1]. The calculated binding energy is 236.9 MeV (236.5 MeV from experiments). An important quantity is the moment of inertia I, which can be obtained This result should not surprise since the initial configurations are obtained by randomly oriented ground state initial condition and this procedure produces spherical shapes on average. Notice that increasing the value of r 0 → 1.26 fm gives Ψ = 0.125 MeV for a sphere, a result that we test briefly below.

III. FUSION CROSS SECTIONS OF IDENTICAL SPIN ZERO NUCLEI
The total fusion cross section of the nuclear reaction is where E cm is the reaction energy in the center of mass frame, µ is the reduced mass of the reaction system and Π l is the fusion probability of the reaction at angular momentum l. We simulate the reactions of identical spin zero nuclei, thus only even l-values are allowed, 28 Si + 28 Si and 12 C + 12 C at different angular momenta l for a given E cm to obtain the fusion probability Π l with hαc model. In Figure 3, we plot the fusion cross section of 28 Si + 28 Si and 12 C + 12 C as function of the reaction energy in the center of mass frame. The fusion cross sections calculated from the neck model are also presented for comparison [36,37]. The hαc model can reproduce the experimental cross section data qualitatively. While for 12 C + 12 C at high E cm where there are no experimental data, the difference between neck model and hαc model is quite large. Naturally, the model interaction can be improved for a better description of the data, but of course fusion below the barrier needs the inclusion of tunneling [37].

IV. ROTATIONS AND DYNAMICAL FIRST ORDER PHASE TRANSITION
In this section, we will explore the dynamical properties of a 28 Si nucleus rotating along the z-axis with initial orbital angular momentum l = l z , in units of . We assume that the orbital angular momenta are transferred through the collision with the 28 Si target nuclei. Due to angular momentum and parity conservation, l must be even in the entrance channel. The amount of angular momentum transferred during the interaction of the two nuclei must be simulated microscopically but this is presently outside the validity of this model. For simplicity and numerical convenience we will restrict our investigation to even l. There are many methods theoretically to give the nucleus an initial angular momentum l, a popular one is the cranking model [1,[3][4][5][6][7][8][9][10][11][12][13][14]. Since we are dealing with individual particles (7α) we will use the following ansatz to give the initial momenta K y , K x ( -units) to particle i: In Equation (5) r 2 xy = i (x(i) 2 + y(i) 2 ) and the sum is extended to all the constituents α-particles of the nucleus. Different events are obtained by different ground states initial positions and we give the initial momenta according to Equation (5). Because of the finite number of particles, this will produce classical fluctuations, while the orbital momenta are quantized. This method is more justified when the excitation energy or angular momenta l is larger, i.e., at higher entropies. Of course one should not be surprised that we are utilizing the concept of entropy since our equations of motion are time reversible invariant. Entropy arises from events mixing and averages over phase space. Notice also that for each event there may be some angular momentum along the x and y directions, see Equation (5), but on average we have l x = l y = 0, this contributes to the initial classical fluctuations. In Figure 4, we plot the excitation energy as function of l(l + 1) recovering the familiar linear behavior [38]. The error bars are given by the variances obtained from event averaging and they are quite small at small l values as expected. The slope of this plot is proportional to the inverse of the moment of inertia and agrees with our estimate in the previous section for a rigid sphere. Changing the values of r 0 and β produces the expected variation, see Figure  4. The moment of inertia is a dynamical quantity since the system expands and breaks into fragments at high E * , thus the obtained value refers to time t = 0 fm/c. The calculations give the most probable energies E * for each value of l, the values obtained are reported in Table I.   An interesting physical quantity is the excitation energy as function of the kinetic energy of the particles. If the system would reach thermal equilibrium, the latter quantity could be related to the temperature. In Figure 5, we plot these quantities and notice the peculiar behavior for E k near 25 MeV. The increase in excitation energy for fixed kinetic energy signals the occurrence of a first order phase transition and signals the opening of new channels such as evaporation, 'fission' and fragmentation, i.e., the sudden increase of the degrees of freedom of the system, from one nucleus to many fragments. We notice that we get a finite probability of breaking into 7α for l = 16, and E * (7α) = 54 MeV, which gives the model lowest excitation energy (in 1000 events) when breaking into 7α. The question remains if the transition is of thermal origin. A signature of thermal equilibrium is to observe the same features for each coordinate. We know that the nucleus is rotating along the z-axis, thus we expect the momenta along the same axis to be small, zero on average, see Equation (5). To reach equilibrium, energy must be transferred from the other directions and this may be impossible for high l-values (and angular momentum conservation). From the large 'error bars' (i.e., variances due to the initial conditions) in Figure 5 we may guess that the system becomes more and more chaotic but that does not prove that thermal equilibrium is reached. A definite answer to this problem may be obtained by repeating the plot as function of the kinetic energy along the z-axis, E kz . If the system reaches thermal equilibrium multiplying E kz by a factor of 3 should reproduce Figure 5. In Figure 6, we plot E * as function of 3 × E kz , compare to Figure 5. At low energy we observe an increase of E * up to about 25 MeV where the phase transition occurs. Higher l-values or E * values do not produce an increase in E kz but just an increase in the variances. It means that even if we increase the excitation energy the system does not have enough time to transfer kinetic energy from the x-y plane to the z-direction, i.e., to reach thermalization. This can be interpreted as an apparent maximum temperature that the system may sustain. Such behavior could be compared to the Lyapunov exponent of an expanding system as discussed in ref. [39], it proves that the phase transition is of dynamical origin. These findings could be experimentally tested [1]. The small density increase at large r is due to α-particle evaporation in some events. For this plot, 50 events were generated.
In Figure 7, we plot the density distribution at two different times t = 150 fm/c and 1500 fm/c at E * = 34 MeV, i.e., in the region of the phase transition. The bump that we observe at later times is due to the escape (evaporation) of one α-particle in some events thus doubling the number of degrees of freedom. There have been suggestions that highly rotating nuclei may display toroidal shapes [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19]r and this was also the focus of the experimental investigation in ref. [1]. Our model allows us to study the shape evolution as function of time for given l or E * . In Figure 8, we plot the coordinates of each α-particle in the x-y plane at two different times, 200 fm/c and 600 fm/c. Notice that only the events breaking into 7α are included in the plot. We use a simple algorithm to recognize the fragments, i.e., we assume that two particles belong to the same fragment if their relative distance is less than 5 fm. This may not be the best approach to recognize fragments at earlier times but it is of little importance since we follow the expansion for very long times. Using this algorithm we can easily estimate the probability of decays into all possible channels allowed by dynamics.
In Figure 8, we see that at 200 fm/c matter is missing at the center during the expansion, left panel. At a later time, more α-particles are recognized and a toroidal shape is observed, notice the change of scales. This is due to the combined effect of the angular momentum and averaging over events. We can repeat the same considerations in momentum space, see Figure 9. At later times the system expands under the influence of Coulomb only and this explains why little expansion is seen in momentum space at the two different times, compare to Figure 8.

V. CROSS SECTION ESTIMATE
It is instructive to give an estimate of the cross-section in order to get some deeper insight into the process and be of guidance and stimulus to more experimental and theoretical investigations [1]. The cross section is given in Equation (4). We consider the reaction 28 Si on 12 C at 35 MeV/A [1], thus E cm = 294 MeV. In order to extend the sum to even l-values only, we assume that the angular momenta are transferred in symmetric Si + Si collisions. T l gives the probability that in the collision a certain angular momentum l and/or excitation energy E * is transferred to the 28 Si, see Table I. This is the part missing in the calculation and we will estimate it from the available phase space in the reaction. We assume that the maximum excitation energy E M that can be transferred to the Si is proportional to its mass, thus: To get this excitation energy in Si + Si we estimate a beam energy E/A = 29.4 MeV and this could be an interesting experiment to confirm our findings and widen the results of ref. [1] and also to investigate the angular momentum transfer to each nucleus during the reaction. Our simple ansatz for the available phase space is: T l → 0 if the excitation energy E * > E M . This crude approximation guarantees that the cross section vanishes at large excitation energies or large l-values. In particular from Table I we   The hαc model provides the probability Π l for the system to break into different open channels for a given l-value using the simple fragment recognizing algorithm discussed above. This quantity is plotted in Figure 10 for the channels as indicated in the inset. We have included in the 7α channel events where one or more 8 Be are produced, this is to simulate the experimental data where those events are implicitly included; however, recall that the model binding energy of 8 Be is overestimated. For reference we have included the experimental points from ref. [1] with huge error bars just to indicate the position of the resonances. We have seen in the Figures 8 and 9 that those values correspond indeed to short living toroidal states. These results can be compared to the toroidal shell model results reported in ref. [1], table II. Notice that the two models differ on the way the orbital angular momenta is given to the system. The 7α channel dominates at high E * values, while lower resonances are dominated by events where at least one large fragment is present. The larger is the heaviest fragment the smaller is the resonant energy. This is qualitatively similar to what was observed in the data [1]. We stress the fact that all the possible α-decay channels can be estimated in the present model. To include channels where nucleons or other fragments are produced we could couple the hαc model to a statistical one [23]. These competing channels will decrease the probability that only α-channels are relevant, thus we expect that our cross sections are overestimated.
The experiment [1] provided the reaction cross section per unit energy and different α-decay channel. The hαc model gives the energy distribution for each channel whose central values E * are reported in Table I, and the variances Σ l are reported in the Figures 3-5. We notice that when selecting particular channels (i.e., 7α), the most probable energies E * and their variances change, see Figure 10, and we use these values in the following calculations. We assume that the energy distribution for each l-value is given by the Gaussian distribution normalized to one, g l (E, E * , Σ l ), with units 1/MeV. Thus, we write the differential cross section as: In Figure 11, we plot the differential cross section for the 7α channel only. Different l-values contributions are included and indicated in the inset. The estimated cross-section is generally above the experimental one given by the open crosses. The integrated cross section is a factor 8 above the data (1.9 mb) partly due to the many different open decay channels not included in the model. Notwithstanding these differences there are some interesting features to notice. The model gives distinguishable bumps at low energies and l = 16-22. The lowest l-values are in a region where experimental values have large error bars. The experimental cross section starts from 62 MeV (55 MeV in the hαc model) [40,41]. Another feature worth noticing is the increase of the widths for increasing l-values, see also Figure  4. When those widths become very large, different l-values distributions overlap and they cannot be distinguished anymore. This signals the crossing into classical dominated dynamics where using the impact parameter becomes a good approximation. Our results suggest that repeating the experiment say at higher beam energies may not reveal more values of E * because they overlap with nearby l-values. Lower beam energies may reveal the wealth of distinguishable E * as in Figure 11. Of course a crucial point would be to have a detector with better granularity and coverage. If not just even l-values are admitted in the entrance channel (i.e., Si + C collisions) may open different scenarios. On the same footing experiments using a 29 Si on an identical target would also be interesting to confirm these findings. The 29 Si beam would have the problem that neutrons are emitted and those are usually difficult to detect in coincidence with 7α.
Interesting consequences of our model can be derived from a simple inspection of Equations (5)-(8) and the oscillations displayed in Figure 11. Indeed similar features have been investigated in the fusion of two heavy ions below and above the Coulomb barrier [37,42,43]. Equation (8) suggests some simple scaling of the cross section by defining the dimensionless quantity: In the equation above, the E cm term is important when comparing the same system at two different beam energies. It does not guarantee overall scaling since at different beam energies, different l-values maybe relevant and we expect variations on the tail of the distributions. Similar to fusion reactions [42], Equation (9) gives the 'energy-weighted excitation function' (EWEF). Notice that in low energy fusion reactions it could be convenient to replace Equation (9) with [42,43]: , where σ R is the reaction cross-section. Since many exit channels are available in fragmentation reactions, it is useful to define the dimensionless excitation energy as: where the Q-value depends on the exit channel, in ref. [1] the decay of 28 Si in 7α was analyzed in detail and Q(7α) = 33.6 MeV. In Figure 12, we compare the experimental dimensionless quantities [1], open crosses, to the hαc model (full squares). The model calculations have been divided by a factor 8 to take into account the difference to the data in the total cross-section [1], compare to Figure 11. Standard errors have been included to the data to indicate the region of low statistics, thus of particular interest is the region 2.4 < E < 5.2. Figure 12 does not add much to Figure 11 but it will become a more interesting observable when data at different beam energies and different combinations of projectile and target will be available. Oscillations in the EWEF can be better displayed by defining its first and second derivative: The first derivative of the EWEF is plotted in Figure 13 as function of the dimensionless excitation energy. Interesting structures in the energies of interest can be noticed. In particular the data show peaks at E = 3.4, 3.6 and 4.1, corresponding to E * = 114, 122 and 138 MeV, respectively, very close to the data analysis of ref. [1]. More peaks may be seen at E = 3.1, 4.5 and 4.7 (E * = 104, 151 and 158 MeV) but in a region where statistics is rather low thus the need for higher statistics experiments. As expected, the model calculations display definite peaks especially at low excitation energies where the data is poor thus impossible to compare, see also Table I. The average data trend is rather well reproduced by the model. In Figure 14, we repeat our analysis for the second derivative of the EWEF. In this plot the difference between model and data is more marked. Peaks in the model calculations occur especially for E <3 while for the data E >3. The peaks are consistent to those obtained from the first derivative of the EWEF.
To complete our comparison to fusion reactions we define the logarithmic derivative as [42,43]: This quantity has the obvious property that all constants entering the EWEF and its derivative cancel out, for instance the factor of 8 used to normalize the model to the data. In Figure 15, we plot L(E) as function of E and remark the substantial differences to the previous plots. Peaks in the data remain only for the values reported in ref. [1]. The model calculations show large fluctuations in the low energy region and in analogy to the fusion hindrance we can interpret these as a signature of the lowest resonances for this particular decay channel. Notice however that the model peaks reported in Figure 11 at E = 55 MeV have very low statistics and are not included in this plot (off scale). If this interpretation is correct we can assume that the L(E) becomes monotonic at large energies, i.e., in the classical limit. Unfortunately not much can be inferred from the low energy data since error bars are large, notice however that in the high-energy region of low statistics, the data show a flat behavior thus suggesting that low energy resonances are reachable with more statistics, better detectors and maybe lower beam energy.

VI. CONCLUSIONS
In this work, we have introduced a semi-classical model for nuclei whose constituents are α particles. In the ground state, the α-particles strongly overlap and this gives rise to repulsion due to the increase of the Fermi energy. This means that a proper description of the nuclear ground state must be done in terms of nucleonic degrees of freedom with the possibility that expanding nuclei coalesce into α-clusters. We have applied our hybrid model to the fragmentation results [1] of 28 Si breaking into even-even N = Z nuclei only. We have found preferential values of the excitation energies for each l-value but also larger and larger variances in the energy distribution due to the fluctuations in the initial conditions, which are classical in origin. For high l-values, the fluctuations become very large and different l-distributions overlap signaling the approach to classical mechanics. We have shown that the spin quantization could be determined experimentally in various different experimental situations by changing the beam energies and the masses of the colliding nuclei including radioactive species. These experiments require very well performing 4π detectors, high statistics and one could take advantage of running the experiments in inverse kinematics. A dynamical phase transition was also revealed together with the limiting temperature that the system could sustain. Above the phase transition, toroidal-like shapes are observed when averaging over many events. These findings may open up a new route of research based on the seminal results of ref. [1] and be linked to the oscillations seen in fusion reactions above and below the Coulomb barrier [42,43]. We have discussed scaled energy weighted excitation functions and its derivatives. We have shown how these quantities consistently display interesting features similar to the barrier fluctuations in fusion reactions. We believe this work may open a new route of investigations to link fusion reaction to deep-inelastic, incomplete fusion and fragmentation. Well performing 4π detectors are needed to improve the energy resolution and granularity and high statistics data are essential. More theoretical work is finally needed to link the proposed analysis to nuclear fundamental properties.