# Virial corrections to simulations of heavy ion reactions

\[

## Abstract

Within QMD simulations we demonstrate the effect of virial corrections on heavy ion reactions. Unlike in standard codes, the binary collisions are treated as non-local so that the contribution of the collision flux to the reaction dynamics is covered. A comparison with standard QMD simulations shows that the virial corrections lead to a broader proton distribution bringing theoretical spectra closer towards experimental values. Complementary BUU simulations reveal that the non-locality enhances the collision rate in the early stage of the reaction. It suggests that the broader distribution appears due to an enhanced pre-equilibrium emission of particles.

\]

The Boltzmann equation including the Pauli blocking (the BUU equation ) and the closely related method of quantum molecular dynamics  (QMD) are extensively used to interpret experimental data from heavy ion reactions. Due to their quasi-classical character, they offer a transparent picture of the internal dynamics of reactions and allow one to link observed particle spectra with individual stages of reactions.

An ambition to cover the heavy ion reactions within experimental errors has been recently cooled down by a failure of BUU simulations to describe the energy and angular distribution of neutrons and protons in low and mid energy domain . Indeed, the Boltzmann equation is not the full story. As noticed in numerical studies of hard sphere cascade by Halbert  and more general by Malfliet , it is disturbing that all dynamical models rely more or less on the use of the space- and time-local approximation of binary collisions inherited from the Boltzmann equation. This approximation neglects a contribution of the collision flux to the compressibility and the share viscosity which control the hydrodynamic motion during the reaction. To include the collision flux and other virial corrections, the non-local character of binary collisions has to be accounted for. Malfliet also demonstrated that non-local collisions can be easily incorporated into simulation BUU codes.

In absence of a first principle theory, Malfliet in his pioneering study, and more recently Kortemeyer, Daffin and Bauer , had to use classical hard-sphere-like non-local collisions which naturally do not result in a full quantitative agreement with experimental data. This ad hoc approximation reflects a gap in former quasi-classical theories of quantum transport: authors either cared about non-local collisions leaving aside quasiparticle features or vice versa. Moreover, quantum theories of binary collisions treated non-local collisions via gradient contributions to the scattering integral  which are numerically inconvenient and thus have never been employed in demanding simulations of heavy ion reactions. Recent theoretical studies have filled this gap. Danielewicz and Pratt  pointed out that the collision delay can be used as a convenient tool to describe the virial corrections to the equation of state for the gas of quasiparticles. Although their discussion is limited to the equilibrium, it marks a way how to introduce virial corrections also to dynamical processes. This task was approached in , where a quasi-classical kinetic equation was derived by a systematic quasi-classical approximation of non-equilibrium Green’s function in the Galitskii–Feynman approximation. The derived kinetic equation is suitable for simulations having the quasiparticle renormalizations in the standard form of Landau’s theory and non-local collisions reminding classical hard spheres.

In this letter we follow the line initiated by Malfliet, however, with more advanced theoretical and numerical tools. We start from the kinetic equation from  and evaluate non-local collisions from the Paris potential solving the two-particle T-matrix. These non-local corrections are then incorporated into the QMD simulation. Results show that the non-local corrections bring the energy distribution of protons closer to experimental values. On a complementary BUU simulation we enlighten the microscopic mechanism leading to this improvement.

The scattering integral of the non-local kinetic equation derived in  corresponds to a following picture of a collision. Assume that two particles, $`a`$ and $`b`$, of initial momenta $`k`$ and $`p`$ start to collide at time instant $`t`$ being at coordinates $`r_a`$ and $`r_b`$. Due to a finite range of the interaction, at the beginning of collision particles are displaced by $`r_br_a=\mathrm{\Delta }^{\mathrm{be}}`$. The collision has a finite duration $`\mathrm{\Delta }_t`$, i.e., it ends at $`t+\mathrm{\Delta }_t`$. During the collision, both particles move so that their end coordinates differ from those at the beginning, $`r_a^{}r_a=\mathrm{\Delta }_a`$ and $`r_b^{}r_b=\mathrm{\Delta }_b`$. The particle $`a`$ transfers a momentum $`q`$ to the particle $`b`$, therefore their relative momentum changes from $`\kappa =\frac{1}{2}(kp)`$ to $`\kappa ^{}=\frac{1}{2}(kp)q`$. Their sum momentum is modified by an external field acting on the colliding particles during the collision going from $`K=k+p`$ to $`K^{}=k+p+\mathrm{\Delta }_K`$. The same field changes the sum energy of colliding particles from $`E=ϵ_a+ϵ_b`$ to $`E^{}=ϵ_a^{}+ϵ_b^{}=ϵ_a+ϵ_b+\mathrm{\Delta }_E`$.

The values of $`\mathrm{\Delta }`$’s are given by derivatives of the total scattering phase shift $`\varphi =\mathrm{Im}\mathrm{ln}T_R(\mathrm{\Omega },k,p,q,t,r)`$ of the two-particle T-matrix $`T_R`$,

$$\begin{array}{cccccc}\mathrm{\Delta }_t\hfill & =& \frac{\varphi }{\mathrm{\Omega }}|_E^{}\hfill & \hfill \mathrm{\Delta }^{\mathrm{be}}& =& \left(\frac{\varphi }{p}\frac{\varphi }{q}\frac{\varphi }{k}\right)_E^{}\hfill \\ & & & & & \\ \mathrm{\Delta }_E\hfill & =& \frac{\varphi }{t}|_E^{}\hfill & \hfill \mathrm{\Delta }_a& =& \frac{\varphi }{k}|_E^{}\hfill \\ & & & & & \\ \mathrm{\Delta }_K\hfill & =& \frac{\varphi }{r}|_E^{}\hfill & \hfill \mathrm{\Delta }_b& =& \frac{\varphi }{p}|_E^{}\hfill \end{array}.$$
(1)   
Note that energy $`\mathrm{\Omega }`$ enters as an independent quantity so that one needs to know the scattering phase shift out of the energy shell. The on-shell energy, $`\mathrm{\Omega }=E^{}`$, is substituted after derivatives are taken. Figure 1 illustrates the nonlocal concept derived in . We like to point out that this concept leads to a continuous trajectory in the kinetic picture replacing real potential scattering. Consequently the energy, momentum, density and angular momentum conservation are conserved including second order quantum virial corrections .

It is our intention to incorporate these features of collisions into the QMD and BUU simulation codes. The selfconsistent evaluation of all $`\mathrm{\Delta }`$’s for all collisions would be too demanding. We employ two kinds of additional approximations. First, following approximations used within the BUU equation, we neglect the medium effect on binary collision, i.e., use the well known free-space T-matrix. Second, we rearrange the scattering integral into an instant but non-local form. This instant form parallels hard-sphere-like collisions what allow us to employ computational methods developed within the theory of gases  similarly as it has been done in .

In the instant approximation we let particles to make a sudden jump at time $`t`$ from $`r_a`$ and $`r_b`$ to effective final coordinates $`\stackrel{~}{r}_a`$ and $`\stackrel{~}{r}_b`$. These effective coordinates and momenta $`\stackrel{~}{\kappa }`$ and $`\stackrel{~}{K}`$ are selected so that at time $`t+\mathrm{\Delta }_t`$ particles arrive at the correct coordinates, $`r_a^{}`$ and $`r_b^{}`$, with the correct momenta, $`\kappa ^{}`$ and $`K^{}`$. Accordingly, in the asymptotic region, after $`t+\mathrm{\Delta }_t`$, there is no distinction between the non-instant and instant pictures, which is shown as solid line in figure 2. This asymptotic condition is naturally met if one extrapolates the out-going trajectories from known coordinates and momenta at $`t+\mathrm{\Delta }_t`$ back to the time $`t`$. Doing so one finds that the effective coordinates read

$`\stackrel{~}{r}_a`$ $`=`$ $`r_a^{}{\displaystyle \frac{kq}{m}}\mathrm{\Delta }_t=r_a+\mathrm{\Delta }_a{\displaystyle \frac{kq}{m}}\mathrm{\Delta }_t,`$ (2)   
$`\stackrel{~}{r}_b`$ $`=`$ $`r_b^{}{\displaystyle \frac{p+q}{m}}\mathrm{\Delta }_t=r_b+\mathrm{\Delta }_b{\displaystyle \frac{p+q}{m}}\mathrm{\Delta }_t.`$ (3)   
The change of the relative momentum is insensitive to the instant approximation, $`\stackrel{~}{\kappa }=\kappa ^{}=\kappa q`$. The sum momentum and energy, however, get simplified because during the instant process the mean-field has no time to pass any momentum and energy to the colliding pair. Accordingly, $`\stackrel{~}{K}=K`$ and $`\stackrel{~}{E}=E`$. In other words, the momentum and energy gains are naturally covered by the effect of the mean field on particles during the time interval $`(t,t+\mathrm{\Delta }_t)`$ which, in the instant picture, is already covered by a free motion. Similarly, in agreement with the continuity of the center of mass motion, one finds that $`\stackrel{~}{r}_a+\stackrel{~}{r}_b=r_a+r_b`$.

When incorporating the displacements into the QMD simulation code, we have to face the fact that two particles are selected for a collision if they meet at the point of closest approach. This distance is different from the distance $`\mathrm{\Delta }^{\mathrm{be}}`$ required from the equivalent scattering scenario presented in figure 2 as solid line.We consider now the time required to travel from $`\mathrm{\Delta }_{\mathrm{be}}`$ to the distance of closest approach $`\stackrel{~}{\mathrm{\Delta }}_t=\frac{m}{2\kappa ^2}\kappa \mathrm{\Delta }_{\mathrm{be}}`$ in analogy to . Within this scenario we are allowed to jump at the point of closest approach to the final asymptotics (2) and (3)whith the additional distance the particle travel during $`\stackrel{~}{\mathrm{\Delta }}_t`$. The effective final coordinates thus have to be evaluated as

$$\stackrel{~}{r}_{a,b}=\frac{R_a+R_b}{2}\mathrm{\Delta }^\mathrm{f},$$
(4)   
with the effective displacement

$$\mathrm{\Delta }^\mathrm{f}=\frac{1}{2}\mathrm{\Delta }^{\mathrm{be}}\mathrm{\Delta }_a+\frac{kq}{m}(\mathrm{\Delta }_t\stackrel{~}{\mathrm{\Delta }}_t).$$
(5)   
Since the center of mass does not jump in the collision, the final displacement can be also written in an alternative way, $`\mathrm{\Delta }^\mathrm{f}=\frac{1}{2}\mathrm{\Delta }^{\mathrm{be}}+\mathrm{\Delta }_b\frac{p+q}{m}(\mathrm{\Delta }_t\stackrel{~}{\mathrm{\Delta }}_t)`$. The non-local corrections are thus performed as follows. When the collision is selected, we evaluate $`\mathrm{\Delta }^\mathrm{f}`$ from (5) and (1), redisplay particles into $`\stackrel{~}{r}_a`$ and $`\stackrel{~}{r}_b`$ and continue with the simulation.

At this point it is possible to establish a connection of the present theory to the hard-sphere-like corrections used by Malfliet  and Kortemeyer, Daffin and Bauer . For hard spheres of the diameter $`d`$, the phase shift has a classical limit $`\varphi =\pi |q|d`$ which gives $`\mathrm{\Delta }_{a,b}=0`$ and $`\mathrm{\Delta }^\mathrm{f}=\mathrm{\Delta }^{\mathrm{be}}=\frac{q}{|q|}d`$. The displacement thus has the same amplitude $`d`$ for all binary collisions and points in the direction of the transferred momentum, as it is known from the Enskog equation . In the present theory, an amplitude of displacement (5) depends on the selected channel, $`\mathrm{\Delta }^\mathrm{f}(k,p,q)`$, and the direction does not coincide with the transferred momentum,

$$\mathrm{\Delta }^\mathrm{f}=\frac{q}{|\kappa |}H_2+\frac{kpq}{|\kappa |}H_1.$$
(6)   
The second term is perpendicular to the transferred momentum and stays in the collision plane. The notation of coefficients $`H_{1,2}`$ is identical to the corresponding interface subroutines which can be obtained from authors.

In order to investigate the effect of non-local shifts on realistic simulations of a heavy ion reaction, we have evaluated $`\mathrm{\Delta }^\mathrm{f}`$ from the two-particle scattering T-matrix $`T^R`$ in the Bethe-Goldstone approximation  using the separable Paris potential . The comparison of the shifts calculated for different potentials concerning partial wave coupling up to D-waves can be found in . We have incorporated these shifts into a QMD code for the central collision of <sup>129</sup>Xe$``$<sup>119</sup>Sn at $`50`$ MeV/A. Figure 3a shows the exclusive proton spectra subtracting the protons bound in clusters. This procedure is performed within a spanning tree model which is known to describe a production of light charged cluster in a reasonable agreement with the experimental data, Figs. 3b-f. Within the local approximation, however, the remaining distribution of high-energy protons is too low to meet the experimental values. As one can see, the inclusion of non-local collisions corrects this shortage of the QMD simulation. As demonstrated in Fig. 3, productions of light clusters are rather insensitive to the non-local corrections. This also shows that the improvement of the proton production is not on cost of worse results in other spectra.

A microscopic mechanism leading to the increase in the high-energy part of the particle spectrum can be traced down to an enhancement of the number of collisions at the pre-equilibrium stage of the heavy ion reaction demonstrated in Fig. 4 for the BUU simulation of the same reaction. This enhancement gives rise to an immediate proton production which itself translates into a high energetic spectra. In other words, the strong production of the high-energy protons follows from the pre-equilibrium emission of particles. The BUU simulation also shows that non-local corrections are important namely in the early stage of reaction well before most of light clusters are formed. It explains why the production of protons is affected while the formation of light clusters is nearly untouched by the non-local corrections.

In summary, as documented by the improvement of the high-energy proton production, the non-local treatment of the binary collisions brings a desirable contribution to the dynamics of heavy ion reactions. According to an experience from the theory of gases, one can also expect a vital role of non-localities in the search for the equation of state of the nuclear matter. It is encouraging that the non-local corrections are easily incorporated into the BUU and QMD simulation codes and do not increase computational time.

We thank the INDRA collaboration for the use of data prior to publication. J.Aichelin and W. Bauer are thanked for the well documented QMD and BUU codes. This work was supported from Czech Republic, the GACR Nos. 202960098 and 202960021, and GAAS Nr. A1010806, and Germany, the BMBF, Nr. 06R0884, and the Max-Planck-Society.