³H and ³He nuclei production in a combined thermal and coalescence framework for heavy-ion collisions in the few-GeV energy regime

In this study, we use a statistical hadronization model to describe hadron production in heavy-ion collisions in the few-GeV energy regime [1−3]. This approach effectively describes the transverse momentum and rapidity spectra of protons and pions. Furthermore, it accurately reproduces the experimentally measured deuteron yield [4], provided that, among the two options, the freeze-out corresponding to the higher temperature is selected. We consider a version of the model that incorporates spheroidal hydrodynamic expansion [2], which is more realistic than the one that assumes spherical symmetry. The latter was used in the original blast-wave model by Siemens and Rasmussen [5]; see also Ref. [6]. It is important to note that our assumptions differ from the commonly used boost-invariant blast-wave models applied at ultrarelativistic energies [7], which are not suitable at lower beam energies.

Our analysis is based on data collected by the HADES Collaboration for Au-Au collisions at a beam energy $ \sqrt{s_{NN}}= $ 2.4 GeV and a centrality class of 10% [8−10]. Fits to the particle abundances performed by us and other groups [11] suggest two distinct sets of possible freeze-out thermodynamic parameters. In particular, they correspond to two different freeze-out temperatures: T = 49.6 MeV and T = 70.3 MeV. Since the analysis of deuteron production favors the higher temperature scenario [4], we adopt T = 70.3 MeV throughout this study. In this case, the system size R is relatively small, approximately 6 fm. As a result, the probability of nucleus formation via nucleon triplet is relatively large (compared to the low-temperature case where $ R \approx 16 $ fm). The thermodynamic parameters used in our calculations, i.e., temperature, baryon chemical potential, and isospin chemical potential, are listed in Table 1. The thermodynamic parameters related to strangeness production are not shown, as they are irrelevant to the present analysis.

T/MeV	$ \mu_{{\rm{B}}} $/MeV	$ \mu_{{\rm{I}}_3} $/MeV	R/fm	H/MeV	δ	$ v_R $	$ \gamma_R $
70.3	876	−21.5	6.06	22.5	0.4	0.60	1.25

Table 1.  Thermodynamic parameters (temperature T, baryon chemical potential $ \mu_{{\rm{B}}} $, isospin chemical potential $ \mu_{{\rm{I}}_3} $), system radius R, Hubble expansion parameter H, longitudinal eccentricity δ, radial flow $ v_R $, and Lorentz gamma factor at the system's boundary ($ r=R $) $ \gamma_R $, extracted from experimental data.

In the collisions of gold nuclei considered, one detects on average per event about 28.7 nuclei of ²H, 8.7 nuclei of ³H, and 4.6 nuclei of ³He [8, 9]. ¹ Hence, about 46.5 protons are found in the bound states. The remaining average number of directly produced protons is 77.6. In the framework defined in Refs. [1, 2] we assume that all protons eventually detected in bound states and those measured as free particles originally form a thermal system. Thus, while comparing the predictions of our model with the data, the final results for the proton spectra are rescaled by the factor 77.6/(77.6+46.5). This implicitly assumes the existence of a certain physical mechanism that combines the remaining protons with neutrons into the light nuclei mentioned above. In the present work, we assume that this process can be interpreted as coalescence and use the standard expressions of the coalescence framework (see, for example, Refs. [12, 13]) to predict the light nuclei spectra. Continuing the approach from Ref. [4], we do not use Gaussian ad hoc parametrizations of the thermal source but rather consider initial particle production as independent thermal production within a sphere of the radius R (which is consistent with thermal model assumptions adopted in Refs. [1, 2]). Our present work does not include contributions from decaying resonances. We have verified using THERMINATOR [14, 15] that for the low-temperature case they are negligible, while for the high-temperature case only the Delta resonance plays a noticeable, though very small, role.

The topic of producing light nuclei has recently attracted considerable interest (see, for example, Refs. [16−28]). In relativistic heavy-ion collisions, the production of such systems is well reproduced by both the thermal approach and the coalescence model. These two approaches are usually treated as exclusive alternatives (for a short and critical review of this issue, see Ref. [29]). In this study, we present a scenario for lower energies where the initial thermal production of particles is combined with a subsequent coalescence mechanism in a consistent way, showing that the two pictures may coexist and supplement each other.

The remainder of this paper has the following structure. In Sec. II, we introduce the coalescence model and discuss the nuclei formation rate. In Sec. III, the freeze-out model is defined. The spheroidal form of expansion is discussed in Sec. IV, where we introduce a compact description of the distribution functions with given symmetries, which becomes very useful for implementing the coalescence model. Our main numerical results are presented and discussed in Sec. V. We conclude the paper in Sec. VI. Throughout the paper, we use natural units: $ c=\hbar=k_B=1 $. The metric tensor has the signature $ (+1,-1,-1,-1) $.

In the coalescence model for light nuclei production, the particle spectrum is obtained as a product of proton and neutron spectra taken at a fraction of the nucleus momentum [30, 31]. Specifically, with the proton and neutron three-momentum distributions defined by

the three-momentum distribution of a nucleus X with a mass number A and an atomic number Z is the product

where $ A_{{\rm{FR}}}^{{\rm{X}}} $ is the particle formation rate coefficient discussed in more detail below, and the subscripts p and n refer to protons and neutrons, respectively.

Equation (2) assumes the additivity of the three-momenta of the nuclei that form the considered nucleus, $ {\boldsymbol{p}}_{{\rm{X}}} = A {\boldsymbol{p}} $. Classically, this leads to a small violation of energy conservation: for example, owing to the finite nucleus binding energy, $ m_{\rm \,^3 He} = 2809 $ MeV, while $ m_{{\rm{p}}}+2m_{{\rm{n}}} = (938+2\times940) $ MeV = 2818 MeV. This and other related problems connected with the conservation laws in the coalescence model are usually circumvented with reference to the quantum nature of the real coalescence process, which introduces the natural uncertainty of the energies and momenta of interacting particles [12]. In any case, small differences due to the binding energy and the inequality of $ m_{{\rm{p}}} $ and $ m_{{\rm{n}}} $ do not affect our present analysis. Therefore, in the following, we use the approximation $ m_{{\rm{p}}} \approx m_{{\rm{n}}} \approx m $, where m is the mean nucleon mass, and $ m_{{\rm{X}}} \approx Am $. Hence, while switching from the variable $ {\boldsymbol{p}}_{{\rm{X}}} $ to the particle transverse momentum $ p_{\perp{\rm{X}}} $ and rapidity $ y_{{\rm{X}}} $ in Eq. (2), we use the simple rules

As in both theory and experiment, one usually deals with invariant momentum distributions, $ E \, {\rm{d}} N/( {\rm{d}}^3 p) $, rather than distributions of forms presented in Eq. (1). It is convenient to recall that, for systems with cylindrical symmetry (with respect to the beam axis z) considered in this study, we have

where E is the on-mass-shell energy of a particle, $ E=\sqrt{m^2+{\boldsymbol{p}}^2} $, and $ m_\perp $ is its transverse mass, $ m_\perp=\sqrt{m^2+p_\perp^2} $. Therefore, from Eq. (2), we obtain

Note that, since the momenta and energies of protons and neutrons are equal, the proton and neutron distributions differ only on account of the different thermodynamic parameters used in equilibrium distribution functions. Note also that, within our approximations $ y_{{\rm{X}}}=y_{{\rm{p}}}=y_{{\rm{n}}}=y $ and at zero rapidity, Eq. (5) reduces to the formula

Previous works introduce two models that reproduce well the spectra of protons and neutrons [1, 2]. Here, we use one of these models to predict, based on Eq. (6), the spectra of ³H and ³He nuclei. For finite values of rapidity, Eq. (5) takes the form

By integrating this equation over the transverse mass of the particle, we obtain the particle rapidity distribution.

A popular form of the coefficient $ A_{{\rm{FR}}}^{{\rm{X}}} $ used in literature is [12]

where the function $ D({\boldsymbol{r}}_1,\dots,{\boldsymbol{r}}_A) $ represents the distribution of the space positions of neutrons and protons at freeze-out, normalized to unity; $ \phi_X({\boldsymbol{r}}_1,\dots,{\boldsymbol{r}}_A) $ is the light nucleus wave function; and s and $ s_{{\rm{X}}} $ stand for the nucleon spin and spin of X, respectively. Many studies on the coalescence model assume the Gaussian profile [32, 33]

where $ R_{{\rm{kin}}} $ is the radius of the system at freeze-out. However, this formula may be considered inconsistent with the physical assumptions used in our model, where original particles are produced independently inside a sphere of radius R at a fixed laboratory time t [1, 2]. The expression (9) gives the root-mean-squared value $ r_{\rm rms} = \sqrt{6} R \approx 2.45 R $, which implies hadron production far away from the original thermal system and its long formation time ². Thus, as an alternative to the Gaussian distribution (9), we use the position distribution of particles that are produced independently within a sphere of radius R

with the normalization condition

In this study, since we are interested in the production of ³H and ³He nuclei, we consider the case $ A = 3 $ and use the form

Here, $ \theta_{{\rm{H}}}(x) $ is the Heaviside step function. In this case, it is convenient to introduce a set of normalized Jacobi coordinates

The coordinates $ {\boldsymbol{r}}_{{\rm{cm}}} $, $ {\boldsymbol{r}}_{12} $, and $ {\boldsymbol{r}}_{312} $ can be interpreted as the center of mass of the X system, the relative distance between two particles, and the distance between this pair's center of mass and a third particle, respectively. Changing to these variables, we write

For the wave functions of ³H and ³He, we use the Gaussian approximation in the same coordinates [19]

where $ R_{\rm \,^3 H} = 1.76\; $fm and $ R_{\rm \,^3 He} = 1.96\; $fm are the radii of ³H and ³He nuclei, respectively. The function $ \left| \phi_{{\rm{X}}}(r_{12},r_{312})\right|^2 $ is normalized to unity. In Ref. [4], it was argued that relativistic corrections can be neglected, and we continue this approach in this study. We note that the wave function's dependence on its arguments can be factorized, allowing us to write

The inner integral is invariant under the rotation of $ {\boldsymbol{r}}_{{\rm{cm}}} $; therefore, we can set

and change the variables to

which gives

The inner integral is equal to

with

In the outer integral, we change the variables again to eliminate the step functions

where $ B(1) $ denotes a ball of radius 1 centered at the origin. The values of the formation rate coefficient for different freeze-out scenarios are listed in Table 2. The wave functions of ³H and ³He are shown in Figs. 1 and 2.

Parameter	Value/MeV6
$ A_{{\rm{FR}}}^{^{3}{\rm{H}}} $	$ 5.51\times 10^{11} $
$ A_{{\rm{FR}}}^{^{3}{\rm{He}}} $	$ 5.03\times 10^{11} $

Table 2.  Values of the formation rate parameter $ A_{{\rm{FR}}} $ for ³H and ³He nuclei.

Figure 1.  (color online) Nucleon space distribution integrated over the center of mass coordinate in terms of the moduli of the $ {\boldsymbol{r}}_{12} $ and $ {\boldsymbol{r}}_{312} $ coordinates, $ D(r_{12},r_{312}) = \left(4\pi\right)^2 r_{12}^2 r_{312}^2 \int {\rm{d}} r^3_{{\rm{cm}}} $$ D(r_{{\rm{cm}}},r_{12},r_{312}) $.

Figure 2.  (color online) Wave function space distribution in terms of the moduli of the $ {\boldsymbol{r}}_{12} $ and $ {\boldsymbol{r}}_{312} $ coordinates $ \left(4\pi\right)^2 r_{12}^2 r_{312}^2 \left| \phi_X(r_{12},r_{312})\right|^2 $ in the Jacobi coordinates of ³H (left) and ³He (right).

The freeze-out models specify the hydrodynamic conditions for particle production at the latest stages of the system's spacetime evolution. In the single-freeze-out scenario [34, 35] adopted here, the freeze-out stage is defined by a set of thermodynamic variables such as temperature T, baryon chemical potential $ \mu_{{\rm{B}}} $, and shape of the freeze-out hypersurface Σ. In addition, one defines the form of the hydrodynamic flow $ u^\mu(x) $ on Σ.

The standard starting point for quantitative calculations is the Cooper-Frye formula [36], which describes the invariant momentum spectrum of particles

Here, $ f(x,p) $ is the phase-space distribution function of particles, and $ p^\mu = \left( E, {\boldsymbol{p}} \right) $ is their four-momentum with the mass-shell energy $ E = \sqrt{m^2 + {\boldsymbol{p}}^2} $.

The infinitesimal element of a three-dimensional freeze-out hypersurface from which particles are emitted $ {\rm{d}}^3\Sigma_\mu(x) $ may be obtained from the formula (see, for example, Ref. [37])

where $ \epsilon_{\mu \alpha \beta \gamma} $ is the Levi-Civita tensor with the convention $ \epsilon_{0123}=-1 $, and $ a, b, c $ are the three independent coordinates introduced to parametrize the hypersurface. This allows us to construct a six-dimensional, Lorentz-invariant density of the produced particles

The independent variables in such a general parametrization are the three components of three-momentum and the variables a, b, and c.

In local equilibrium, the distribution function $ f(x,p) $ has the following general form

where T is the freeze-out temperature, $ \chi=-1 $ ($ \chi=+1 $) for Fermi-Dirac (Bose-Einstein) statistics, and $ g_s = 2s +1 $ is the degeneracy factor connected with spin. In this study, we use such distributions of protons and neutrons only (as input for the light nuclei coalescence formula); hence, we take $ \chi=-1 $ and $ s=1/2 $. Since we use the same mass for protons and neutrons, their equilibrium distributions differ only by the values of thermodynamic potentials that determine the fugacity factor

Here, $ \mu_{{\rm{B}}} $ and $ \mu_{{\rm{I}}_3} $ are the baryon and isospin chemical potentials, respectively. The values of T, $ \mu_{{\rm{B}}} $, and $ \mu_{{\rm{I}}_3} $ for the considered expansions scenario are listed in Table 1. In the following, the distribution in Eq. (27) is used in the coalescence formula (Eq. (2)).

In the case of freeze-outs that are spheroidally symmetric with respect to the beam axis, it is expedient to use the following parametrization of the spacetime points on the freeze-out hypersurface:

Here, $ {{\boldsymbol{e}}}_{r\perp} = \left( \cos\phi\sin\theta, \sin\phi \sin\theta\right) $, and the parameter $ \epsilon\; $ controls deformation from a spherical shape. For $ \epsilon\; >\; 0 $, the hypersurface is stretched in the (beam) z-direction. The resulting infinitesimal element of the spheroidally symmetric hypersurface has the form

where the prime denotes the derivatives taken with respect to ζ. We also introduce the spheroidally symmetric hydrodynamic flow

where $ \gamma(\zeta,\theta) $ is the Lorentz factor, which is given by the formula

resulting from the normalization condition $ u\cdot u=1 $. The four-momentum of a hadron is parameterized as

with $ E=m_\perp \cosh(y) $, $ p_\perp = \sqrt{m_\perp^2-m^2} $, and $ p_\parallel=m_\perp \sinh(y) $. Thus, we obtain

and

where $ \kappa (\xi)= \sqrt{1+\xi}\, p_\parallel\cos\theta + \sqrt{1-\xi}\, p_\perp\sin\theta \cos(\phi-\phi_p) $.

Prior analyses of the spectra indicated that a satisfactory description of the data can be obtained by assuming $ t^\prime=0 $ (which allows considering $ \zeta = r $), $ \epsilon=0 $, and $ \delta \neq 0 $. The flow profile is represented by $ v(r)=\tanh(H r) $ [38], where the value of the Hubble constant H is as given in Table 1. In this case,

and the Cooper-Frye formula for fermions takes the form

where

The scalar product $ u\cdot p $ is given by Eq. (34), and spherical symmetry allows us to consider $ \phi_p = 0 $. Changing the variables from p and $ \theta_p $ to rapidity and transverse mass, we write

where

Within our approximations, the angle in Eq. (40) is the same for nucleons and nuclei. In the zero-rapidity case,

The model calculations of the spectra of ³H ($ A=3, Z=1 $) and ³He ($ A=3, Z=2 $) nuclei are based on Eq. (5) with the nucleon distributions defined by Eq. (39), where we use different values of the fugacity factor Υ for protons and neutrons. Our results for the rapidity distributions and transverse-mass spectra are shown in Figs. 3 and 4.

Figure 3.  (color online) Predicted transverse-mass spectra at zero rapidity for ³H (red solid line) and ³He (blue dashed line).

Figure 4.  (color online) Predicted rapidity spectra of ³H (red solid line) and ³He (blue dashed line).

Integrating the spectra, we find the total yields, which we compare with the experimental results in Table 3. We observe that the model calculations underestimate the experimental results by factors of 2.7 and 2.0 (for ³H and ³He, respectively). The main uncertainty of the theoretical results comes from the formation rate factors that are calculated in a relatively simple manner, as described in Sec. II B. More realistic calculations of the formation rate may improve the agreement with the data (here, we consider, for example, more realistic forms of the wave functions). We note that the values of the coefficients $ A_{{\rm{FR}}} $ are very large (see Table 2), which means that our final results are products of very large and small numbers (all in units of MeV). Therefore, it is quite intriguing that we correctly reproduce the order of magnitude of the yields. This is achieved by the assumption that the coalescence takes place just after freeze-out within a sphere of radius R. We note that other theoretical frameworks, for example, IQMD + MST, also underestimate the yields; see Refs. [39, 40].

Nucleus	Quantity	Model	Experiment [8, 9]
3H	N	3.16	$ 8.65 \pm (0.01)_{{\rm{stat}}} \pm (1.05)_{{\rm{sys}}} $
	$ ( {\rm{d}} N/ {\rm{d}} y)_{y=0} $	4.06	–
3He	N	2.26	$ 4.55 \pm (0.01)_{{\rm{stat}}} \pm (0.29)_{{\rm{sys}}} $
	$ ( {\rm{d}} N/ {\rm{d}} y)_{y=0} $	2.89	–

Table 3.  Model predictions for the total yield N and particle yield at zero rapidity $ ( {\rm{d}} N/ {\rm{d}} y)_{y=0} $ compared with the experimental results [8, 9].

Additional verification of the model can be obtained from a comparison of the rapidity and transverse-momentum spectra. In fact, very often only the correct description of the shapes of the spectra of composite particles is regarded as the indication that the coalescence mechanism is at work. To check whether we can reproduce the rapidity profiles of the preliminary HADES data [8, 9], we have varied the values of $ A_{{\rm{FR}}} $ by rescaling them by a factor α (compared to the original values derived above). In the case of the rapidity distribution of ³H, the best description is obtained for $ \alpha = 2.3 $, which is slightly smaller than the factor 2.7 needed to reproduce the experimental yield. For ³He, the best adjustment is obtained for $ \alpha = 1.7 $, compared to the factor 2.0 found above to fit the yield, see Fig. 5. ³ The corresponding total multiplicities are given in Table 4, where we also show the dependence on the integration range in rapidity.

Figure 5.  (color online) Rapidity spectra prediction for ³H (left) and ³He (right) for different scaling factors: exact model value (red solid line, $ \alpha =1.0 $), fit to the experimental yield (blue dashed line, $ \alpha =2.7 $ for ³H and $ \alpha =2.0 $ for ³He), and best fit to data (teal dash-dotted line, $ \alpha =2.4 $ for ³H and $ \alpha =1.9 $ for ³He). The data shown are preliminary HADES results [8, 9].

Nucleus	Scaling factor α	$ \int^{0.65}_{-0.65} $	$ \int^{\infty}_{-\infty} $	Experiment [8, 9]
3H	1.0	3.10	3.16	$ 8.65 \pm (0.01)_{{\rm{stat}}} \pm (1.05)_{{\rm{sys}}} $
	2.7	8.50	8.67	–
	2.4	7.02	7.16	–
3He	1.0	2.21	2.26	$ 4.55 \pm (0.01)_{{\rm{stat}}} \pm (0.29)_{{\rm{sys}}} $
	2.0	4.44	4.55	–
	1.9	3.75	3.85	–

Table 4.  Model predictions for the total yield N calculated for different scaling factors α and integration ranges compared with the experimental results [8, 9].

On account of the uncertainty of the normalization of the experimental transverse-mass spectra, we have refrained from making comparisons between theory and experiment in this case. However, we observe that the slopes of our model for the $ m_{\perp} $-distributions are steeper than those shown as preliminary results.

In this study, we use the previously constructed model of hadron production in heavy-ion collisions in the few-GeV energy regime [1, 2, 4] to determine the yields and spectra of ³H and ³He. Our approach combines the thermal framework with single freeze-out with the concept of coalescence of emitted nucleons. This approach proves to be very successful in determining the yield of ²H [4]. The present study extends the previous framework to the next lightest nuclei. Our predictions for the yields of ³H and ³He are smaller by factors of 2.7 and 2.0 than the experimental values, respectively. Despite this discrepancy, it is noteworthy that our simple model correctly describes the order of magnitude of the multiplicities studied. We also note that the model yield values are controlled by the formation rate factors $ A_{{\rm{FR}}} $, whose calculation can be improved in the future by using more realistic wave functions of ³H and ³He. The comparison of the preliminary HADES results on the spectra shows that we are missing a factor of approximately 2 to reproduce the rapidity spectra.

Since a similar factor is missing to reproduce the yields of both ³H and ³He, our approach effectively reproduces the experimental ratio of these two yields. The importance of the studies of such ratios has been recently stressed in Ref. [41]. Finally, we note that there might be other sources of light nuclei production that enhance their yields, such as the break-up of spectator nuclei [28, 42]. However, such a mechanism is not included in our model, which describes the production of particles in the participant zone only.

We thank Małgorzata Gumberidze, Szymon Harabasz, and Piotr Salabura for useful discussions concerning the experimental data.