A new approach based on relativistic kinetic equations is proposed to solve the long-standing puzzle of light cluster formation, also called nucleosynthesis, in high-energy heavy-ion collisions. This method addresses the tension between STAR data and previous studies relying on either statistical equilibrium or coalescence approaches.

The creation of the first light nuclei approximately 20 minutes after the Big Bang is commonly known as nucleosynthesis [1], which formed deuterium, tritium (short-lived), helium, and lithium. The remnants of this process can still be observed today in the yields of deuterium and helium and minuscule traces of lithium and beryllium in the solar system. In line with the small binding energy of deuterons (the nucleus of deuterium), which is of the order of a few MeV, Big Bang nucleosynthesis proceeds at a temperature below 1 MeV. Heavier elements are then created via the fusion processes in stars, supernovae, and neutron star mergers [2, 3].

Today’s heavy-ion accelerators recreate conditions similar to those shortly after the Big Bang, enabling detailed studies on the formation of light elements. Over the past decade, two seemingly contradicting ideas have emerged to understand the production of light clusters during the collision of heavy ions. The first one is direct thermal statistical production and emission of clusters from a hot fireball at a temperature of approximately 155 MeV [4, 5]. The second one is the formation of clusters via coalescence [6-8] during the expansion stage by interactions of relatively cold protons and neutrons.

A recent study by Sun, Wang, Ko, Ma, and Shen [9] has significantly advanced our understanding of light cluster formation in relativistic heavy-ion collisions. The study focused on the production of deuterons (²H), tritons (³H), and helium (³He) in collisions of gold nuclei over an energy range of $\sqrt{s_{\text{NN}}}=7.7-200$ GeV (as explored in the STAR experiment at the Brookhaven National Lab’s RHIC collider [10, 11]) and in the ALICE experiment at the LHC with collisions at $\sqrt{s_{\text{NN}}}=2.76$ and 5.02 TeV [12, 13].

To solve the puzzle of light cluster formation, Sun et al. modeled the most violent stages of the reaction using a state-of-the-art (3+1)-dimensional relativistic viscous fluid dynamics model (MUSIC) [14], coupled with a network of relativistic kinetic equations [15, 16] to examine the formation of light clusters while maintaining partial chemical equilibrium [17, 18, 5]. The starting time of the kinetic approach is ascertained by a local and time-dependent particlization criterion [19], chosen here as the fixed energy density of hadronization, which provides the initial values based on local Cooper–Frye sampling from the above-mentioned constant energy density hypersurface [20].

After converting the viscous fluid to hadrons, Sun et al. solved the kinetic equations for cluster formation and dissociation in the expanding hadronic matter, which in the case of the deuteron read $\frac{\partial f_{\text{d}}}{\partial t}+\frac{\text{P}}{E_{\text{d}}}\frac{\partial f_{\text{d}}}{\partial \text{R}}=-L{f}_{\text{d}}+G\left(1+f_{\text{d}}\right)\mathrm{,}$ (1) where f_d denotes the deuteron phase-space distribution, P the deuteron momentum, and R the deuteron coordinate. The functions L and G on the left-hand sides of Eq. (1) are, respectively, the loss and gain rates given by the relativistic scattering kernel. Because the loss and gain rates may generally involve more than two particles, the geometric collision criterion applied in most dynamical approaches cannot be used. To overcome this difficulty, a stochastic scattering term needs to be implemented. For example, the gain rate ${\pi }^{+}+n+p\to {\pi }^{+}+d$ can be expressed as a scattering probability in the time interval Δt per volume ΔV, $P_{3\to 2}\propto v_{{\pi }^{+}p}^{\text{rel}}{\sigma }_{{\pi }^{+}p\to {\pi }^{+}p}\frac{W_{\text{d}}}{N_{\text{test}}^2}\frac{\Delta t}{\Delta V}+\left(p\leftrightarrow n\right)\mathrm{,}$ (2) where $v_{{\pi }^{+}p}^{\text{rel}}$ represents the relative velocity of the pion and proton, and the deuteron Wigner function is denoted as W_d. (The deuteron Wigner function is a phase-space representation of the deuteron wave function [21].) Similar kinetic equations allow to calculate the time-dependent gain and loss rates for all considered clusters without assuming equilibrium in the system.

The main idea is illustrated in Fig. 1, which depicts the time evolution of a relativistic heavy-ion collision. The evolution proceeds from left to right, starting with the initial colliding nuclei experiencing a pre-equilibrium phase, followed by the quark–gluon plasma phase (modeled by relativistic hydrodynamics) until hadronic freeze-out with cluster formation (modeled by relativistic kinetic equations).

Fig. 1Full-screen View

(Color online) Time evolution of a relativistic heavy-ion collision. The evolution proceeds from left to right, starting with the initial colliding nuclei, via a pre-equilibrium phase, followed by the quark–gluon plasma phase (modeled by relativistic hydrodynamics) until hadronic freeze-out with cluster formation (modeled by the relativistic kinetic equations). The arrows denote the momenta of the particles, and the colors encode the particle type. The figure is taken from Ref. [9]

After careful testings of their novel approach in static box calculations to validate the numerical implementation and demonstrate that the kinetic equations yield the correct thermodynamic asymptotic state, the recent STAR data on cluster production can be analyzed.

The strong effect of the hadronic stage, as modeled by the kinetic equations, is illustrated in Fig. 2, which shows the time evolutions of the deuteron (d) and triton (³H) numbers for gold-gold (Au+Au) collisions at $\sqrt{s_{\text{NN}}}=200$ GeV for calculations with (yellow band) and without (blue band) hadronic rescattering after thermal cluster production at the end of the hydrodynamic phase.

Fig. 2Full-screen View

(Color online) Calculated time evolution of the deuteron (d) and triton (³H) numbers for calculations with (yellow band) and without (blue band) hadronic rescattering for gold-gold (Au+Au) collisions at $\sqrt{s_{\text{NN}}}=200$ GeV. The red symbol indicates the experimental data from the STAR collaboration. The figure is taken from Ref. [9]

One clearly observes that especially for the triton (³H) a substantial depletion emerges due to the hadronic dynamics which is absent in the case of direct emission from the hydrodynamic stage (blue band). The main reason for this depletion is that the kinetic equation allows for the evolution of the triton yield through the hadronic stage either in partial equilibrium or even out of equilibrium.

Physically, the suppression of the triton yield is in line with other indications of a long-lasting hadronic rescattering stage after the quark-gluon-plasma phase transition. For example, the yields of hadronic resonances (states such as K^*, ρ, Λ^*, ∑^*) are suppressed relative to their ground states. This is usually interpreted as the absorption of their decay daughter particles which are rescattered in the hadronic medium and destroy the resonance signal in an invariant mass analysis [22]. In addition, these analyses reveal a hadronic lifetime of the order of 10 fm/c [23].

The energy dependence of cluster production is shown in Fig. 3 (left) for the triton-to-proton ratio in central Au+Au collisions over an energy range of $\sqrt{s_{\text{NN}}}=$7.7–200 GeV. The red symbols represent STAR data [11], whereas the kinetic equation results are depicted as yellow bands. The blue line indicates the statistical hadronization model (SHM) (taken from [11]). One clearly observes that the inclusion of further hadronic interactions after the hadronization from the hydrodynamic stage leads to a reduction of the triton production by a factor of 1.8 compared to the statistical model. As shown in Fig. 3 (right), this effect is even more pronounced by studying the ratio of the product of the number of tritons and protons to the square of the number of deuterons. The inclusion of hadronic re-scatterings led again to a substantial reduction and brings the calculations more closely with experimental data. One should note that there have been intense speculations of a potential local maximum of this ratio, which might indicate the existence of enhanced neutron number fluctuations, potentially related to a critical point [24].

Fig. 3Full-screen View

(Color online) Left: Energy dependence of the triton (³H)-to-proton ratio from the kinetic equation approach (yellow band) compared with the statistical hadronization model, SHM (blue line and dashed blue line, SHM divided by 1.8) and the data from the STAR experiment (symbols). Right: The ratio of the product of the number of tritons and protons to the square of the number of deuterons as a function of energy (yellow band) compared with the statistical hadronization model, SHM (blue line), and the STAR data (symbols) for gold-gold (Au+Au) collisions from $\sqrt{s_{\text{NN}}}=7.7-200$ GeV. The figure is taken from Ref. [9]

In summary, the production of light clusters in relativistic heavy-ion collisions resembles the process of nucleosynthesis as it was proceeding around 20 minutes after the Big Bang. During the Big Bang nucleosynthesis, the formation of deuteron, triton and helium proceeded only after the temperature dropped below a few MeV. In heavy-ion collisions, two conflicting interpretations were put forward: statistical production directly from the phase boundary and late stage coalescence in kinetic models. The current work by Sun et al. has made important progress to clarify this question and to lead to a unified picture of light cluster production. By solving a network of relativistic kinetic equations, following the hadronization they quantified the deviations from the statistical hadronization model with high precision, while still benchmarking their rate equations to yield the statistical result for asymptotic times. Additionally, the study provides the first quantitative description of the full RHIC data on the ratio $N_{{}^3\text{H}}\times N_{\text{p}}/N_{\text{d}}^2$ and indicates that the previous speculations about enhanced neutron number fluctuations are not well justified and require further investigations.

Finally, one should not underestimate the potential of the developed method for the further exploration of more exotic cluster state, e.g. multi-strange hypernuclei, which will be in the focus of new experiments at FAIR [25, 26] and HIAF [27].