Revisiting angular momentum conservation in transport simulations of intermediate-energy heavy-ion collisions

We use simplified but reasonable setups in the IBUU transport model to illustrate the effect of the angular momentum conservation on the dynamics of intermediate-energy heavy-ion collisions. The simulations are mainly focused on non-central Au+Au collisions at different beam energies. The initial coordinates of neutrons and protons are sampled according to their density distributions in Au nucleus obtained based on the Skyrme-Hartree-Fock model with the MSL0 force Chen et al. (2010), and their initial momenta are sampled isotropically within the local isospin-dependent Fermi sphere. The nucleon momenta are then boosted according to the beam energy, and the simulation is performed in the C.M. frame of Au+Au collisions. We use the Skyrme-like momentum-independent mean-field potential including a symmetry potential linear to the density, which reproduces the empirical nuclear matter properties, and their detailed forms are the same as those used in Ref. Xu et al. (2016). To implement the mean-field potential, we employ the lattice Hamiltonian method Lenk and Pandharipande (1989). We use point particles in the implementation of the Coulomb force, for which the cut-off distance is set to be 1 fm to avoid divergence.

We use the modified Bertsch’s prescription Bertsch and Das Gupta (1988); Zhang et al. (2018) for nucleon-nucleon collisions. The minimum distance of two colliding particles in their C.M. frame perpendicular to their relative velocity is

where $\vec{r}_{1}^{\star}$ and $\vec{r}_{2}^{\star}$ are positions of the two particles, and $\vec{v}_{12}^{\star}=\vec{v}_{1}^{\star}-\vec{v}_{2}^{\star}$ is their relative velocity, with the asterisk representing the quantity in the C.M. frame of the colliding particles. The collision can happen if the condition

is satisfied, and we use a constant and isotropic baryon-baryon elastic cross section $\sigma=40$ mb in the present study. Whether the collision happens in this time step is determined by the condition of the closest approach, i.e.,

We set $\delta t=\Delta t/\gamma$, where $\gamma=1/\sqrt{1-\vec{\beta}^{2}}$ is the Lorentz factor with $\vec{\beta}$ being the average velocity of the colliding pair in the computational frame. We have also removed the spurious collisions by setting that the two particles, which have collided once, can not collide again unless one of them has collided with a third particle. We use point nucleons in the implementation of the Pauli blocking, where the cubic local phase-space occupation probability with the dimension $\Delta r=2$ fm and $\Delta p=0.1$ GeV/c is calculated to evaluate the Pauli blocking probability together with the interpretation method. If the collision is Pauli blocked, both the momenta and coordinates of the colliding particles are retained. Besides elastic $N+N\rightarrow N+N$ collisions, we have also incorporated the channels for the pion production, including $N+N\leftrightarrow N+\Delta$ and $\Delta\leftrightarrow N+\pi$ channels for different isospin states of nucleons, $\Delta$ resonances, and pions. The detailed inelastic cross sections, decay width, and $\Delta$ and pion masses are exactly the same as those in the homework setups of Ref. Ono et al. (2019).

For the ease of discussion on the collision prescription for particle 1 and particle 2, we first give the basic quantities and their relations in their C.M. (collision) frame and computational (heavy-ion simulation) frame. The momentum $\vec{p}_{1}^{\star}$, the coordinate $\vec{r}_{1}^{\star}$, the time $t_{1}^{\star}$, and the energy $E_{1}^{\star}$ of particle $1$ before collision in the C.M. frame can be expressed in terms of those in the computational frame as follows

with $a_{1}=\gamma(\frac{\gamma}{\gamma+1}\vec{p}_{1}\cdot\vec{\beta}-E_{1})$ and $b_{1}=\gamma(\frac{\gamma}{\gamma+1}\vec{r}_{1}\cdot\vec{\beta}-t_{1})$. For particle 2, the same relations hold with the subscript $1\rightarrow 2$. In the C.M. frame, the momenta satisfy $\vec{p}_{1}^{\star}=-\vec{p}_{2}^{\star}\equiv\vec{p}^{\star}$. We assume that the collision occurs at the same time ($t_{1}=t_{2}=0$) in the computational frame, but generally at different times in the C.M. frame ($t_{1}^{\star}\neq t_{2}^{\star}$). Similarly, the corresponding quantities in the computational frame after collision can be expressed in terms of those in the C.M. frame as

with ${a^{\prime}}_{1}^{\star}=\gamma(\frac{\gamma}{\gamma+1}\vec{p^{\prime}}_{1}^{\star}\cdot\vec{\beta}+{E^{\prime}}_{1}^{\star})$ and ${b^{\prime}}_{1}^{\star}=\gamma(\frac{\gamma}{\gamma+1}\vec{r^{\prime}}_{1}^{\star}\cdot\vec{\beta}+{t^{\prime}}_{1}^{\star})$. As is seen, quantities with a prime represent those after collisions. Similarly, $\vec{p^{\prime}}_{1}^{\star}=-\vec{p^{\prime}}_{2}^{\star}\equiv\vec{p^{\prime}}^{\star}$ is satisfied. In order to set the same time $t^{\prime}_{1}=t^{\prime}_{2}=0$ for particles after collision in the computational frame, we set ${t^{\prime}}_{1}^{\star}=-\vec{r^{\prime}}^{\star}_{1}\cdot\vec{\beta}$ and ${t^{\prime}}_{2}^{\star}=-\vec{r^{\prime}}_{2}^{\star}\cdot\vec{\beta}$.

Defining the average and relative coordinates before collision in the C.M. frame as

the angular momentum before collision in the computational frame can be expressed as

and the angular momentum after collision in the computational frame can be expressed similarly as

In the above, $E^{\star}={E^{\prime}}^{\star}=E_{1}^{\star}+E_{2}^{\star}={E^{\prime}}_{1}^{\star}+{E^{\prime}}_{2}^{\star}$ is satisfied from the energy conservation condition. $\vec{J}^{\star}=2\vec{r}^{\star}\times\vec{p}^{\star}$ and $\vec{J^{\prime}}^{\star}=2\vec{r^{\prime}}^{\star}\times\vec{p^{\prime}}^{\star}$ are the angular momentum in the C.M. frame before and after collision, respectively. To assure the conservation of the angular momentum in the computational frame ($\vec{J}=\vec{J^{\prime}}$), we require the angular momentum conservation in the C.M. frame ($\vec{J}^{\star}=\vec{J^{\prime}}^{\star}$) as well as that from the C.M. motion, i.e., the rest parts in Eqs. (4) and (5) should be equal. The latter requires a shift of the C.M. coordinate

Therefore, the coordinates of particle 1 and particle 2 after collision in their C.M. frame can be respectively written as

representing the so-called “side jump”.

We illustrate with Fig. 1 how to conserve the angular momentum ($\vec{J}^{\star}=\vec{J^{\prime}}^{\star}$) in the C.M. frame of an elastic collision. In the original collision prescription, as shown in Fig. 1 (a), the coordinates of colliding particles are unchanged after collision, and their momenta are sampled in the $4\pi$ solid angle according to the differential cross section. In this way, both the direction and the magnitude of the angular momentum in the C.M. frame of the collision are generally changed. To conserve the angular momentum, the coordinates of the particles after collision should generally be changed, once the directions of their momenta after collision are determined. Given $|\vec{p}_{1}^{\star}|=|\vec{p^{\prime}}_{1}^{\star}|=|\vec{p}_{2}^{\star}|=|\vec{p^{\prime}}_{2}^{\star}|$, the constraint of the angular momentum conservation requires an in-plane collision, i.e., without changing the azimuthal angle, and the distance between the horizontal blue dashed lines, representing the direction of the initial momenta, should be the same as that between the inclined blue dashed lines, representing the directions of the final momenta. Then, there are infinite options to achieve that, and in the present study we consider two typical options. Figure 1 (b) shows a similar collision prescription as in Ref. Liu et al. (2020), where both the coordinates and the momenta of the colliding particles rotate around their C.M. with a certain angle, and we dub this prescription as “rotation” in the following discussion. Figure 1 (c) shows a collision prescription that the coordinates of particles after collision move vertically to the blue dashed line, and we dub it as “vertical”, which is identical to the prescription in Ref. Gale and Das Gupta (1990). In the “vertical” prescription, we have the least change of the coordinates after collision in the C.M. frame, so in principle this prescription should be most similar to the original one, except that an in-plane collision is required.

The above prescription for elastic collisions can be easily generalized to inelastic collisions. For the $N+N\leftrightarrow N+\Delta$ channel in the present study, we have $|\vec{p}^{\star}|\neq|\vec{p^{\prime}}^{\star}|$, and the distance between the inclined blue dashed lines should be modified to satisfy $\vec{r}^{\star}\times\vec{p}^{\star}=\vec{r^{\prime}}^{\star}\times\vec{p^{\prime}}^{\star}$. In the $\Delta\rightarrow N+\pi$ process, $N$ and $\pi$ are produced at the same position as $\Delta$, and the angular momentum is always conserved. We are unable to conserve the angular momentum in the $N+\pi\rightarrow\Delta$ process without further incorporating spin degree of freedom, and leave this channel as it is in the original IBUU model.

Before simulating heavy-ion collisions, we first check with the elastic nucleon-nucleon collision rate in a box with the periodic boundary condition. The box has a size of $20\times 20\times 20$ fm³. Initial coordinates of 640 neutrons and 640 protons are sampled uniformly within the box, and their initial momenta are sampled within the Fermi sphere. Thus, the system is initialized at zero temperature and at the saturation density. Without incorporating the mean-field potential, Pauli blocking, or Coulomb potential, the system evolves with only elastic nucleon-nucleon collisions, and the momentum distribution gradually changes to a Boltzmann distribution at the temperature of about $T=14.28$ MeV as a result of energy conservation, as in Ref. Zhang et al. (2018). Figures 2 (a) compares the collision rates from the three collision prescriptions as in Fig. 1 with the theoretical limit $dN_{coll}/dt=111.4$ c/fm. The increase of the collision rate at early stage represents the slightly higher collision rate from the Boltzmann distribution than the initial Fermi-Dirac distribution. The “original” collision prescription overestimates the theoretical limit, as shown in Ref. Zhang et al. (2018), as a result of higher-order correlations. The “rotation” prescription leads to a collision rate consistent with the theoretical limit, while the “vertical” prescription underestimates the theoretical limit. On the other hand, it is remarkable to see that the new prescriptions reproduce the theoretical limit within $5\%$. We have checked that the three prescriptions lead to exactly the same momentum distribution after reaching thermal equilibrium. To understand the different collision rate, we have further checked the distribution of nucleon coordinates. While the coordinates are generally uniformly distributed within the box, there are always local density fluctuations, and different collision prescriptions may lead to different local density fluctuations and thus different collision rates. We use the event average value of the neighboring distance for each nucleon $\Delta r_{min}=\sum_{i=1}^{N}\min\{|\vec{r_{i}}-\vec{r_{j}}|\}/N$, where particle $j$ is in the same event of particle $i$, and $N=1280$ is the total nucleon number in each event, to qualify the local density fluctuation. As shown in Fig. 2 (b), the initial value of $\Delta r_{min}$ is about 1.02 fm, while it soon becomes stable at different values for different collision prescriptions. Prescriptions of in-plane collisions, especially the “vertical” prescription which has the minimum change of the coordinate after collision, lead to larger $\Delta r_{min}$, thus weaker local density fluctuations and lower collision rates. We note that for a simple and face-centered cubic lattice, $\Delta r_{min}$ can be as large as 1.842 and 2.067 fm at the saturation density in a box system, respectively.

We then move to simulations of non-central Au+Au collisions, and compare the total angular momentum from three collision prescriptions at different collision energies in Fig. 3. In the Cascade mode with only elastic and inelastic collisions as well as Pauli blockings, the angular momentum $L_{y}$ perpendicular to the reaction plane is exactly conserved for the “rotation” and “vertical” prescriptions at 100 AMeV, and slightly violated at higher collision energies. The latter is due to the difficulty of conserving the angular momentum in the $N+\pi\rightarrow\Delta$ channel without incorporating spin degree of freedom, as mentioned before. The violation becomes even larger after incorporating the mean-field potential in the Full-mode calculation, but after all the violations are at the $0.1\%$ level. Incorporating the Coulomb potential between point charged particles in the present study always conserves the angular momentum. For the “original” collision prescription, although the angular momentum violation in one collision may cancel with that in another collision, the violation of the total angular momentum is much larger ($3\%$ at higher collision energies) compared to the “rotation” and “vertical” prescriptions. The violation of the total angular momentum in the “original” collision prescription is mostly due to nucleon-nucleon collisions at 400 and 1000 AMeV, but is seen to be much enhanced by the mean-field interaction at 100 AMeV.

The constraint of the angular momentum conservation generally requires the change of the coordinates of colliding particles after collision. Figure 4 compares the distribution of the side-jump distance $\Delta r$ in the computational frame from the “rotation” and “vertical” prescriptions in non-central Au+Au collisions at $E_{lab}=100$ AMeV, and the results are found to be insensitive to the collision energy. One sees that $\Delta r$ is smaller than 1.2 fm in all collisions, and is smaller in the “vertical” prescription than in the “rotation” prescription, even if the distance is now calculated in the computational frame.

We now compare the nucleon dynamics and nucleon observables in detail from three collision prescriptions in non-central Au+Au collisions. The comparisons are only in the Full-mode calculation with elastic and inelastic collisions as well as Pauli blockings, mean-field potential, and Coulomb potential.

Figure 5 compares the central density evolution from three collision prescriptions in non-central Au+Au collisions at different collision energies. As shown in Fig. 2, in-plane collisions, i.e., without changing the azimuthal angle after collisions, generally lead to a larger average distance between nucleons, thus a faster expansion of the system and a lower density. For the “rotation” prescription, the coordinates of the colliding particles are effectively pulled back a little after collisions, so this prescription leads to a slower expansion and a higher density compared to the “vertical” prescription, which happens to result in a central density evolution closer to that from the “original” prescription.

Figure 6 compares the C.M. energy dependence of the successful elastic nucleon-nucleon collision rate from three collision prescriptions. With the same Pauli blocking implementation, the successful collision rate is dominated by the density evolution as shown in Fig. 5. One sees that the “vertical” prescription leads to fewest collisions compared with other two prescriptions. The relative successful collision rates from the three prescriptions are qualitatively consistent with those in Fig. 2 (a) where Pauli blockings are not implemented. We note that the correlation between the successful collision number and the dynamics becomes different in different collision prescriptions.

Figures 7 and 8 compares the rapidity dependence of the directed flow $v_{1}=\langle p_{x}/p_{T}\rangle$ and the transverse momentum dependence of the elliptic flow $v_{2}=\langle(p_{x}^{2}-p_{y}^{2})/p_{T}^{2}\rangle$ for free nucleons, which are determined by a cut-off density $\rho<0.15\rho_{0}$ at their freeze-out. $v_{1}$ from the three prescriptions are very similar, consistent with the smaller effect on the transverse flow found in Ref. Gale and Das Gupta (1990), and the slightly stronger $v_{1}$ from the “vertical” prescription is likely due to the fast expansion of the system as mentioned above. The negative $v_{2}$ at large transverse momenta manifests the expansion of the participant matter blocked by the spectator matter and the squeeze-out of energetic nucleons perpendicular to the reaction plane. It seems that such squeeze-out effect is strongest from isotropic nucleon-nucleon collisions in the “original” prescription, since in non-central collisions nucleons are more likely to stay in the reaction plane with in-plane collisions. The “rotation” prescription has an effectively pulling-back movement, so nucleons are more likely to stay in the participant region. This further weakens the squeeze-out effect and makes the $v_{2}$ less negative at high transverse momenta, especially at higher collision energy.

For all results from Fig. 5 to Fig. 8, the difference among results from different collision prescriptions increases with increasing collision energy. This is expected since the Pauli blocking becomes less important and there are more successful nucleon-nucleon collisions at higher collision energies. This also shows the importance of incorporating properly the angular momentum conservation in the simulation of heavy-ion collisions at higher energies.

We now compare with the pion productions from three collision prescriptions in non-central Au+Au collisions, by turning on the $N+N\leftrightarrow N+\Delta$ and $\Delta\leftrightarrow N+\pi$ channels for different isospin states. We neglect potentials for $\Delta$s and pions as well as the possible threshold effect in the present calculation, and only illustrate the effect of the angular momentum conservation in the collision prescription on the pion multiplicity and $\pi^{-}/\pi^{+}$ yield ratio.

Figure 9 compares the time evolutions of the $\Delta$ and pion multiplicities from three collision prescriptions at different collision energies. The peak values of the $\Delta$ multiplicities in the three collision prescriptions are consistent with their maximum densities as shown in Fig. 5. Compared with the “original” prescription, both channels of $N+N\rightarrow N+\Delta$ and $N+\Delta\rightarrow N+N$ are suppressed in the “vertical” prescription due to the larger $\Delta r_{min}$ as shown in Fig. 2 (b). Since the $\Delta$ in the $N+\Delta\rightarrow N+N$ channel is farther from the nucleon in the “vertical” prescription compared to the “original” prescription, the $\Delta$ has a larger chance to decay. Therefore, the “vertical” prescription leads to a pion multiplicity that is about $20\%$ larger than that from the “original” prescription. Such effect is stronger at 400 AMeV but weaker at 1000 AMeV, due to the higher density and thus the weaker $\Delta r_{min}$ effect at higher collision energies. For the “rotation” prescription, there are more $\Delta$s produced compared with the “vertical” prescription, and the larger pion multiplicity than the “original” prescription can also be explained by the above $\Delta r_{min}$ effect.

Figures 10 compares the time evolutions of the ratio $(\pi^{-}/\pi^{+})_{like}=(\pi^{-}+\Delta^{-}+\frac{1}{3}\Delta^{0})/(\pi^{+}+\Delta^{++}+\frac{1}{3}\Delta^{+})$ from three collision prescriptions at different collision energies. While the $(\pi^{-}/\pi^{+})_{like}$ ratios are similar at $E_{lab}=400$ AMeV from different prescriptions within statistical error, the “vertical” prescription leads to a ratio which is about $1.5\%$ larger than those from the other two prescriptions at $E_{lab}=1000$ AMeV. The effect is smaller than that from the symmetry energy Li (2002a, b) but not negligible. The larger $(\pi^{-}/\pi^{+})_{like}$ ratio from the “vertical” prescription is likely due to the stronger suppression on the $p+\Delta^{-}\rightarrow n+n$ channel as a result of less protons in the high-density region, compared with other prescriptions with smaller $\Delta r_{min}$ as shown in Fig. 2 (b).