Production of light nuclei at colliders -- coalescence vs. thermal model

The production of light nuclei in relativistic heavy-ion collisions is well described by both the thermal model, where light nuclei are in equilibrium with hadrons of all species present in a fireball, and by the coalescence model, where light nuclei are formed due to final-state interactions after the fireball decays. We present and critically discuss the two models and further on we consider two proposals to falsify one of the models. The first proposal is to measure a yield of exotic nuclide $^4{\rm Li}$ and compare it to that of $^4{\rm He}$. The ratio of yields of the nuclides is quite different in the thermal and coalescence models. The second proposal is to measure a hadron-deuteron correlation function which carries information whether a deuteron is emitted from a fireball together with all other hadrons, as assumed in the thermal model, or a deuteron is formed only after nucleons are emitted, as in the coalescence model. The $p\! -\!^3{\rm He}$ correlation function is of interest in context of both proposals: it is needed to obtain the yield of $^4{\rm Li}$ which decays into $p$ and $^3{\rm He}$, but the correlation function can also tell us about an origin of $^3{\rm He}$.


Introduction
Production of light nuclei in nucleus-nucleus collisions has been studied for decades but plethora of experimental results from Relativistic Heavy Ion Collider (RHIC) [1,2] and Large Hadron Collider (LHC) [3,4,5] have revived an interest in the problem and attracted a lot of attention. At high-energy collisions light nuclei occur as fragments of incoming nuclei, as at low-energy collisions, but we also deal with a genuine production process -the energy released in a collision is converted into masses of baryons and antibaryons which form nuclei and antinuclei.
The remnants of initial nuclei occur at rapidities of projectile and target nuclei while the genuinely produced nuclides populate a midrapidity domain. Therefore, products of the two mechanisms are kinematically well separated. Further on we are interested only in the midrapidity domain where numbers of the nuclei and antinuclei are approximately equal to each other at RHIC and are exactly equal at LHC. The baryon-antibaryon symmetry clearly shows that the matter created in the collisions a e-mail: stanislaw.mrowczynski@ncbj.gov.pl arXiv:2004.07029v1 [nucl-th] 15 Apr 2020 is (almost) baryonless -there is no baryon charge. Together with light nuclei and antinuclei up 4 He and 4 He there are also produced hypertritons and antihypertritons both at RHIC and LHC [6,7].
According to the coalescence model [8,9] production of light nuclei is a two-step process: production of nucleons and formation of nuclei due to final-state interaction among nucleons which are close phase-space neighbors. The energy scale of the first step, which is a double nucleon mass, is much bigger than that of the second one which is a nuclear binding energy. Consequently, a probability to produce a given nuclide factorizes into the probabilities to produce the nucleons and to form the nuclide. The latter probability takes into account an internal structure of a produced nucleus. Although the coalescence model is known to work well in a broad range of collision energies, the model is well justified when nucleons are truly produced because of the energy scale separation. So, it is not surprising that the model properly describes production of light nuclei and antinuclei at RHIC and LHC [10,11,12,13,14].
The thermodynamical model, see the review [15], is also more reliable and simpler at the highest available collision energies than at lower ones. Since thousands of hadrons are produced at colliders, it is easier to justify the statistical assumption of equipartition of energy. The model is also simpler because the matter produced at midrapidity is baryonless and consequently the baryon chemical potential vanishes. Therefore, particle's yields are determined solely by two thermodynamic parameters: the temperature and system's volume. Nevertheless the model predicts very well not only the yields of all hadron species measured at RHIC and LHC but also of light nuclei and hypernuclei [16,17,18]. The predictions depend on masses and numbers of internal degrees of freedom of the light nuclei but are independent of their internal structures.
The evident success of the thermal model, which has attracted a lot of interest [14,19,20,23,24,25,26,27,28,30], is very puzzling, as it is hard to expect that nuclei exist in the hot and dense fireball environment. The temperature is much bigger than a binding energy and the system is so dense that the inter-particle spacing is smaller than the typical inter-nucleon distance in a nucleus. Therefore, proponents of the thermal model speculate [18] that nuclei are produced as colorless droplets of quarks and gluons with quantum numbers that match those of the final-state nuclei.
The thermal and coalescence models are physically different but its was observed long ago that the models give rather similar yields of light nuclei [31]. The observation has been recently confirmed [12,19], using a more refined version of the coalescence model [32,33,34,35,36], which properly treats a quantum-mechanical character of the formation process of light nuclei.
The aim of this short but rather pedagogical review is to critically discuss the thermal and coalescence models and to analyze why the models predict similar yields of light nuclei. Subsequently we consider a possibility to falsify one of the models. We present two proposals. The first one is to compare the yield of 4 He to that of exotic nuclide 4 Li [19,20,22]. Since the masses of both nuclei are close to each other, the yield of 4 Li is according to the thermal model about five times bigger than that of 4 He due to five spin states of 4 Li and only one of 4 He. The coalescence model predicts instead a significantly bigger yield of 4 He than that of 4 Li because the former nuclide is well bound and compact while the latter one is weakly bound and loose. The model also predicts that the ratio of yields strongly depends, in contrast to the thermal model, on the collision centrality [20,22].
The second proposal relies on the observation [21] that a hadron-deuteron correlation function can tell us whether deuterons are directly emitted from a fireball or they are formed later on due to final-state interactions. The radii of the source of deuterons differ from each other by the factor 4/3 in the two cases. Therefore, knowing the radius of a nucleon source from the proton-proton correlation function we can quantitatively distinguish the emission of a deuteron from the fireball, as in the thermal approach, from the formation of a deuteron afterwards, as in the coalescence model.
We also discus the p− 3 He correlation function which is important in the context of both proposals. The correlation function needs to be measured to obtain the yield of 4 Li which is unstable and decays into the p − 3 He pair. The p − 3 He correlation function also carries information analogous to that of the hadron-deuteron correlation function. If one assumes that 3 He is emitted directly from the fireball the source radius inferred from the p− 3 He correlation function is smaller by the factor 3/2 than that corresponding to the scenario where nucleons emitted from the fireball form the nuclide 3 He due to final-state interactions.
Throughout the article we use the natural units with c = = k B = 1.

Coalescence and thermal models
Let us introduce the coalescence and thermal models. We stress again that we do not consider light nuclei which are fragments of colliding nuclei but only those genuinely produced at midrapidity in collider experiments at RHIC or LHC.

Coalescence model
As already mentioned in the introduction, production of light nuclei is a two-step process in the coalescence model [8,9]. At first nucleons are produced and later on a formation of nuclei proceeds due to final-state interactions of nucleons which are close to each other in momentum and coordinate space. The fact that the energy scale of the first step, which is the double nucleon mass, is much higher than that of the second one, which is a typical nuclear binding energy, is important for two reasons. First of all, the two steps of the process can be distinguished because the temporal scales are roughly inverse of the energy scales. Consequently, the production of nucleons, which occurs at first, is a fast process while the formation of nuclei is a slow process which occurs subsequently. Secondly, a probability to produce a nucleus of A nucleons can be factorized into the probability to produce (independently) A nucleons and the probability that nucleons fuse into the nucleus. Therefore, the number of nuclei with where p A = Ap and p is assumed to be much bigger than the characteristic internal momentum of a nucleon in the nucleus of interest, A A is the formation rate of a nucleus under consideration. One often assumes, as suggested long ago in [9], that nucleons form a nucleus if they occur in a momentum sphere of a radius p 0 . Then, where g S and g I are spin and isospin factors which take care of probability that quantum numbers of A nucleons match to those of the nuclide of interest. The nucleons are assumed to be unpolarized. The momentum distributions of protons and of neutrons are assumed to be of the same shape but the numbers of protons and of neutrons can differ. The parameter p 0 , which is roughly a momentum of a nucleon in a nucleus, is a free parameter of the model to be inferred from experimental data.
It is also often required that nucleons, which fuse into a nucleus, must be close to each other not only in the momentum space but in the coordinate space as well, see e.g. [31]. The formula (2) is then modified.
To formulate a relativistically covariant coalescence model one uses, see e.g. [32], the Lorentz invariant momentum distributions and writes down the relation analogous to (1) as where E A and E N are energies of the nucleus and nucleons under consideration. Demanding that the relations (1) and (3) are identical in the center-of-mass frame of the nucleus, which is formed, the parameter B A is found to be where m is the nucleon mass. The form of the relation (3) implies that the parameter B A is a Lorentz scalar. However, one should realize that the formula (4) does not have a solid foundation. We return to this point after a quantum-mechanical approach to the formation rate is introduced.
The phenomenological approaches to production of light nuclei, which are based of the formula (2) or their variations, do not take into account a quantum-mechanical character of the process of a bound state formation. However, it was discovered by Sato and Yazaki [32] and discussed later on by several authors, see e.g. [33,34,35], that the formation of a nucleus driven by final-state interactions is fully analogous to the process responsible for short range correlations observed among final-state hadrons with small relative velocities. Therefore, the quantum-mechanical formula which gives the deuteron formation rate is almost identical to that of neutron-proton correlation function [36]. The two quantities are actually related to each other due to completeness of quantum states [37,38,39].
The formation rate of a nucleus of A nucleons A A is given as where g S and g I are, as previously, the spin and isospin factors; the multiplier (2π) 3(A−1) results from our choice of natural units where = 1; V is the normalization volume which disappears from the final formula; the source function D(r) is the normalized to unity position distribution of a single nucleon at the kinetic freeze-out and Ψ (r 1 , r 2 , . . . r A ) is the wave function of the nucleus of interest. The formula (1) does not assume, as one might think, that the nucleons are emitted simultaneously. The vectors r i with i = 1, 2, . . . A denote the nucleon positions at the moment when the last nucleon is emitted from the fireball. For this reason, the function D(r i ) actually gives the space-time distribution. It is often chosen in the isotropic Gaussian form where √ 3R s is the root-mean-square (RMS) radius of the nucleon source. The Gaussian parametrization of the source function (6) is not only convenient for analytical calculations but there is an empirical argument in favor of this choice. The imaging technique [40] allows one to infer the source function from a two-particle correlation function provided the inter-particle interaction is known. The technique applied to experimental data from relativistic heavy-ion collisions showed that non-Gaussian contributions to the source functions are rather small and do not much influence the correlation functions [41].
As already mentioned, Lorentz invariant momentum distributions are used in the relativistically covariant coalescence model and the coalescence rate formula (5) is modified, see e.g. [32,34]. However, the modifications are actually heuristic as a theory of strongly interacting bound states faces serious difficulties. In particular, there is no factorization of a center-of-mass and relative motion. To avoid complications one considers the formation process in the center-of-mass frame of the nucleus to be formed where the process can be treated nonrelativistically even so momenta of nucleons are relativistic in both the rest frame of the source and in the laboratory frame. The point is that the formation rate is non-negligible only for small relative momenta of the nucleons. Therefore, the relative motion can be treated as nonrelativistic and the corresponding wave function is a solution of the Schrödinger equation. The source function, which is usually defined in the source rest frame, needs to be transformed to the center-of-mass frame of the pair as discussed in great detail in [43].
The practical calculations of the formation rate A A for A = 2, 3 and 4, which require a separation of a center-of-mass and relative motion, are presented in Secs. 3.2, 4.3 and 4.1, respectively.
It was repeatedly stated in the literature -starting from the very first paper on the coalescence model [8] -that a neutron and proton must interact with a third body to form a deuteron as otherwise the energy and momentum cannot be conserved simultaneously. The statement, which was also extended to nuclei heavier than a deuteron, is indeed correct if the neutron and proton are both on mass shell. However, as observed long ago [34], nucleons, which are emitted from a fireball, are not on the mass shell due to the finite space-time size of a fireball. The space-time localization of a nucleon within the fireball washes out its four-momentum due to the uncertainty principle. Using a more formal language of scattering theory, the nucleons are not in an asymptotic state in the remote past or remote future which indeed requires the mass-shell condition. Instead the nucleons are in an intermediate scattering state. Therefore, there is no reason to impose the mass-shell constraint. Because the spacetime size of the fireball is of the same order as that of the nucleus, which is formed, the mismatch of the energy-momentum is washed out by the uncertainty of energy and momentum of the nucleons.
We close the presentation of the coalescence model by saying that whenever we refer to the model we keep in mind the expression (1) with the formation rate given by Eq. (5).

Thermal model
The fundamental postulate of the thermal model is the equipartition of fireball's energy among all degrees of freedom active in the system. Light nuclei are assumed to be populated as all other hadrons and when the fireball decays the nuclei show up in a collision final state. Their yield reflects a thermodynamic state of the fireball at the moment of chemical freeze-out when inelastic collisions of fireball's constituents become no longer operative.
In the fireball rest frame a momentum distribution of hadrons h at the moment of chemical freeze-out is, see e.g. [15], where g h is the number of internal degrees of freedom of the hadron h, m h is its mass and E p ≡ m 2 h + p 2 is the energy, V chem is the system's volume at the chemical freeze-out, β chem ≡ 1/T chem is the inverse temperature and µ is the chemical potential related to conserved charges carried by the hadron. The baryon chemical potential usually plays an important role but, as already mentioned, the matter created in midrapidity at collider energies is baryonless and consequently the chemical potential vanishes. We note that the formula (7) is classical that is it neglects effects of quantum statistics. The effects are usually minor because of many hadron species which populate many quantum states. The formula (7) also neglects effects of inter-hadron interactions.
The yield of hadrons h is where K 2 (x) is the so-called Macdonald function which for x 1 can be approximated as A microscopic mechanism responsible for production of light nuclei in the fireball is unspecified and may be even unknown. As the formula (8) shows, the yields of hadrons are determined by their masses and internal degrees of freedom, and by two thermodynamical parameters: V chem and T chem . The temperature, however, is much bigger than a typical nuclear binding energy and the inter-particle spacing is smaller than a characteristic size of light nuclei. So, the nuclei cannot exist in the fireballthey are as 'snowflakes in hell' [42]. Proponents of the thermal model argue [18] that there are colorless droplets of quarks and gluons present in the fireball and those with appropriate quantum numbers are converted later on into light nuclei.

Do the models differ?
The coalescence is a microscopic picture while the direct thermal production is a macroscopic description. So, one wonders whether the production mechanisms of light nuclei behind the coalescence and thermal models are physically different from each other.
One can argue that instead of the two models we should rather consider, as in the study [24], hadron-hadron and hadron-deuteron interactions which are responsible for a deuteron production and disintegration in a fireball before its decay. Such an approach is physically sound if a particle source and an average inter-hadron spacing in the source are both much bigger than a deuteron size. Additionally the lifetime of the source should be much longer than the characteristic time of deuteron formation.
However, the assumptions are rather far from reality of relativistic heavy-ion collisions. The deuteron radius is about 2 fm and the time of deuteron formation, which is of the order of the inverse binding energy, is roughly 100 fm/c. Consequently, the size of the particle source is of the same order as the deuteron radius, the inter-hadron spacing in the source is smaller than a deuteron, and the lifetime of the source is significantly shorter than the deuteron formation time.
The coalescence mechanisms of deuteron formation and direct thermal production are physically different in relativistic heavy-ion collisions because the particle source is small and dense when compared to a deuteron and the source lifetime is shorter than the deuteron formation time. According to the coalescence model, light nuclei are formed long after nucleons are emitted from the source. The thermal model assumes that light nuclei are emitted directly from the source.

Yield of deuterons
As already mentioned in the introduction, its was observed long ago that the coalescence and thermal models give rather similar yields of light nuclei [31]. The results of the thermal model can be easily reproduced by means of simple formulas but those of the coalescence model are usually obtained using Monte Carlo generators, see e.g. [12]. Therefore, it is hard to see how it happens and why its happens that the models predict similar yields of light nuclei. For this reason we derived [19] simple analytical formulas which give the ratio of yields of deuterons -the simplest nuclei -in the two models. The model parameters were inferred from experimental data. In this section we recapitulate and improve the analysis presented in [19].

D/p in thermal model
The yield of protons is given by the formula (8). Since β chem m, where m is the proton mass, equals 6 for T chem = 156 MeV, we use the approximation (9) and write the proton yield as where except the spin degeneracy factor 2 we have included the parameter λ p which takes into account a sizable contribution of protons coming form decays of baryon resonances [15]. The parameter will be estimated later on. Since the number of deuterons is given by the formula analogous to (10), the ratio of the deuteron to proton yield is where the spin degeneracy factor of a deuteron is 3 and its mass is approximated as 2m. We note that the parameter λ D analogous to λ p is not introduced because the contribution of deuterons, which originate from decays of excited light fragments, is negligible [29].

D/p in coalescence model
The momentum distribution of the final-state deuterons is given by the formulas (1) and (5) both with A = 2. Introducing the center-of-mass variables and writing down the deuteron wave function as ψ(r 1 , r 2 ) = e iPR φ D (r), the deuteron formation rate equals where the normalized to unity 'relative' source function is The latter equality holds for the Gaussian single-particle source function (6). In this section the source radius carries the index 'kin' not 's' to stress that we deal with the kinetic freeze-out. The factor 3 4 reflects the fact the deuterons come from the neutron-proton pairs in three spin triplet states out of four possible spin states of unpolarized two nucleons.
To compute the deuteron yield, the nucleon momentum distribution needs to be specified. We write down the proton distribution in terms of the transverse momentum (p T ), transverse mass m T ≡ m 2 + p 2 T , and rapidity (y) as and we choose the distribution at midrapidity which is flat in rapidity and azimuthal angle and it exponentially decays with the transverse mass that is where ∆y is a small rapidity interval centered at y = 0 and T kin is the effective temperature at the kinetic freeze-out which takes into account a sizable radial expansion of the fireball. One checks that the distribution (16) is normalized to N p . The number of deuterons is found as To obtain the formation rate A 2 in a simple analytic form, we do not use the Hulthén wave function of a deuteron, as we did in [36], but we choose not only the Gaussian parameterizations of the source function but also of |φ D (r)| 2 . Thus, we get where R D is the root-mean-square radius of a deuteron. In our original paper [19] the factor 2/3 in front of R 2 d in Eq. (17) was missing which influenced though insignificantly our numerical results.
Using the formula (17), the ratio of the deuteron to proton yields equals where the number of protons N p is assumed to be the same as in the thermal model and it is given by Eq. (10). The ratio of the ratios (18) and (11), which is denoted as Q, equals the ratio of deuteron yields in the coalescence and thermal models because the proton yield is assumed to be the same in both approaches. The ratio Q equals In the next section, after estimating the parameters which enter Eq. (19), a magnitude of the ratio Q is computed.

Discussion of D/p ratio
The D/p ratio found within the thermal model (11) is determined by the proton mass m, the temperature of chemical freeze-out T chem and the parameter λ p which we choose in such a way that at T chem = 156 MeV the ratio (11) reproduces the experimental value 3.6 × 10 −3 measured in Pb-Pb collisions at √ s NN = 2.76 TeV [3].
Thus, one finds λ p = 2.51 which is used further on.
To obtain the D/p ratio within the coalescence model (18), one needs, except m, T chem and λ p , the values of ∆y, R d , V chem , R kin and T kin . The measurement [3] was performed in the rapidity window ∆y = 1. The root-mean-square radius of the deuteron is R d = 2 fm [44]. V chem can be found from the deuteron analog of the formula (10), using the measured number of deuterons at different collision centralities given in [3].
In our original paper [19] we used the femtoscopic pion data [45] to get a value of R kin . However, the experimental analysis [46] shows that the radii of pion sources are significantly bigger than those of proton sources. Therefore, we use here the radii of (anti-)proton sources inferred from proton-proton and antiproton-antiproton correlation functions at the smallest transverse momentum [46]. Since the data presented in [46] are for different centrality bins than those in [3] we have performed a linear interpolation or extrapolation to get values of R kin for centrality bins given in the first column of the Table 1.
The parameter T kin from the formula (19) is the effective temperature at kinetic freeze-out which takes into account a sizable radial expansion of the fireball. To determine T kin we express it through the mean transverse momentum of deuterons p T given in [3] using the formula Since the effective kinetic temperature is comparable to the nucleon mass, the approximation (9) cannot be applied to the Macdonald function in Eq. (20).
In Table 1 we list the values of the ratio Q defined by Eq. (19) for the four collision centralities together with the parameters of Pb-Pb collisions at √ s NN = 2.76 TeV.
The predictions of the models are seen to differ by the factor smaller than 2 for all Table 1. The ratio Q and the centrality dependent parameters of Pb-Pb collisions at √ sNN = 2.76 TeV. The numbers in the first three columns are taken from the experimental study [3]. The parameters V chem , T kin and R kin are estimated as explained in the text. The ratio Q is given by Eq. (19). four centralities where number of deuterons grows five times from the peripheral to central collisions. We conclude that the two models indeed predict very similar yields of deuterons. We do not see any deeper reason for the similarity but it is also no accidental, this is a game of numbers -the parameters, which characterize the produced matter at the chemical and kinetic freeze-out, are correlated with each other in a specific way.
4 How to falsify one of the models?
As we discussed in the previous section, the thermal and coalescence model predict similar yields of light nuclei. One asks whether one of the models can be falsified. For this purpose we need a situation that the models give quantitatively different predictions. Below we discuss two such situations and two proposals to distinguish the models. The first one is to measure the yield of exotic nuclide 4 Li and compare it to that of 4 He. The ratio of yields of 4 Li to 4 He is different in the thermal and coalescence models [19,20,22]. The second proposal is to measure a hadron-deuteron correlation function which appears to carry information [21] whether a deuteron is emitted from a fireball together with all other hadrons, as assumed in the thermal model, or a deuteron is formed only after nucleons are emitted, as in the coalescence model. Another version of the second proposal is to measure a hadron− 3 He correlation function which can tell us about an origin of 3 He [22].

4 Li vs. 4 He
The mass of exotic nuclide 4 Li is close to the mass of 4 He. However, there are five spin states of 4 Li, which has spin 2, and only one spin state of 4 He, which has zero spin. Consequently, the thermal model predicts five times bigger yield of 4 Li than that of 4 He. The nuclide 4 Li is weakly bound and loose while 4 He is well bound and compact. The coalescence model is thus expected to predict a significantly smaller yield of 4 Li than that of 4 He. So, let us derive the ratio of yields of 4 Li to 4 He, which in the thermal model simply equals 5, and in the coalescence model it is given by the ratio of formation rates to be computed according to the formula (5) with A = 4.
The modulus squared of the wave function of 4 He is chosen as |Ψ He (r 1 , r 2 , r 3 , r 4 )| 2 = C α e −α(r 2 12 +r 2 13 +r 2 14 +r 2 23 +r 2 24 +r 2 34 ) , where C α is the normalization constant, r ij ≡ r i − r j and α is the parameter related to the RMS radius R α of 4 He. Further calculations are performed using the Jacobi variables defined as which have the nice property that the sum of squares of particles' positions and the sum of squares of differences of the positions are expressed with no mixed terms of the Jacobi variables that is With the help of relations (23) and (24), one easily finds where V is the normalization volume of the plane wave describing a free motion of the center of mass. Substituting the formulas (6) and (21) with the parameters (25) into Eq. (5), one finds the coalescence rate of 4 He as where the spin and isospin factors have been included. Since 4 He is the state of zero spin and zero isospin, the factors are g S = g I = 1/2 3 because there are 2 4 spin and 2 4 isospin states of four nucleons and there are two zero spin and two zero isospin states. The stable isotope 6 Li is a mixture of two cluster configurations 4 He − 2 H and 3 He− 3 H [47]. Since 4 Li decays into p+ 3 He, we assume that it has the cluster structure p− 3 He . Following [47] we parametrize the modulus squared of the wave function of 4 Li, which is treated here as a stable nucleus, as where the nucleons number 1, 2 and 3 form the 3 He cluster while the nucleon number 4 is the proton; z is the Jacobi variable (22); Y lm (Ω z ) is the spherical harmonics related to the rotation of the vector z with quantum numbers l, m. The summation over m is included in the spin factor g S . Using the Jacobi variables, one computes the constant C Li together with the parameters β and γ as where R c and R Li are the root-mean-square radii of the cluster 3 He and of the nuclide 4 Li, respectively. Substituting the formulas (6) and (27) where the isospin and spin factors are computed as follows. The nuclide has the isospin I = 1, I z = 1 and thus the isospin factor is g I = 3/2 4 because there are three isospin states I = 1, I z = 1 of four nucleons. The spin of the ground state of 4 Li is 2 which can be arranged with the orbital angular momentum l = 1 and l = 2. We assume here that the cluster 3 He has spin 1/2 as the free nuclide 3 He. When the spins of 3 He and p are parallel or antiparallel, the orbital number is l = 1 or l = 2, respectively. However, the ground state of 4 Li is of negative parity which suggests that l = 1. We note that the parity of a two-particle system is P = η 1 η 2 (−1) l and internal parities η 1 , η 2 of 1 H and 3 He are both positive. When l = 1, the total spin of 3 He and p has to be one and there are 3 2 such spin states of four nucleons. Consequently, there are 3 2 angular momentum states with 5 states corresponding to spin 2 of 4 Li and thus g S = 5/2 4 . The ratio of the formation rates A Li 4 and A He 4 depends on four parameters: R s , R α , R Li and R c . The fireball radius at the kinetic freeze-out R s can be inferred from the proton-proton correlation functions which have been precisely measured at LHC [46,48]. The RMS radius of 4 He is R α = 1.68 fm [44] and the RMS radius of the cluster 3 He, which is identified with the radius of a free nucleus 3 He, is R c = 1.97 fm [44]. The radius R Li is unknown but it is expected to be 2.5-3.5 fm. The ratio of A Li 4 to A He 4 is shown in Fig. 1 as a function of R s for four values of R Li = 2.0, 2.5, 3.0 and 3.5 fm.
As already mentioned, the ratio of yields of 4 Li to 4 He in the thermal model equals 5 and it is independent of the size of particle's source. Fig. 1 shows that the ratio is significantly smaller in the coalescence model and it significantly depends on R s that is it depends on collision centrality. Therefore, performing measurements at several centralities it should be possible to quantitatively distinguish the coalescence mechanism of light nuclei production from the creation in a fireball.
The nuclide 4 Li is unstable and it decays into p + 3 He with the width of 6 MeV [49], see also [50]. Therefore, the yield of 4 Li can be experimentally obtained through a measurement of the p− 3 He correlation function which is discussed in Sec. 4.4.

Hadron-deuteron correlations
In this section we show that a hadron-deuteron correlation function carries information about the source of deuterons and allows one to determine whether a deuteron is directly emitted from the fireball or if it is formed afterwards. At first we derive the hadron-deuteron correlation function treating a deuteron as in the thermal model, that is as an elementary particle emitted from a source together with all other hadrons. Further on a deuteron is treated as a neutron-proton bound state formed at the same time when the hadron-deuteron correlation is generated.

Deuteron as an elementary particle
The h−D correlation function R is defined as where dN h d 3 p h , dN D d 3 p D and dN hD d 3 p h d 3 p D are number densities of h, D and h − D pairs with momenta p h , p D and (p h , p D ); q is the relative momentum of h and D in their center-of-mass frame. If the correlation results from quantum statistics and/or finalstate interactions, the correlation function is known to be [51,52] where the source function D(r) is, as previously, the probability distribution of emission points and ψ q (r h , r D ) is the wave function of the hadron and deuteron in a scattering state. Let us eliminate the center-of-mass motion of the h−D pair in a non-relativistic manner. We introduce the center-of-mass variables where M ≡ m h +m D , and we write down the wave function as ψ q (r h , r D ) = e iRP φ q (r) with P being the momentum of the center of mass. The correlation function (31) is then found to be where the 'relative' source is With the Gaussian single-particle source function (6), the relative source function equals which is independent of particle masses even so the variable R given by Eq. (32) depends on m h and m D . The single particle source function (6) is assumed to be independent of particle's momentum and particle's mass. This is not quite right as, in general, a source radius depends on both particle's mass m and momentum. More precisely, it scales with the particle's transverse mass m ⊥ ≡ m 2 + p 2 ⊥ . For the case of one-dimensional analysis relevant for our study, the effect is well seen in Fig. 8. of Ref. [46] where experimental data on Pb-Pb collisions at LHC, which are of particular interest for us, are presented. The dependence of the source radius on m ⊥ is evident when we deal with pions and m ⊥ 0.9 GeV. However, the dependence becomes much weaker for protons when m ⊥ 1.0 GeV. The figure shows that the radius of proton source tends to decrease in central Pb-Pb collisions when m ⊥ grows from 1.1 GeV to 1.7 GeV but the decrease is not seen for the collision centrality 10 − 30% nor 30 − 50%. The behavior is well understood as the decrease of the source radius with growing m ⊥ is caused by the collective radial flow which is stronger in central than in peripheral collisions.
In case of proton-deuteron correlations, which are of our particular interest, the interval of m ⊥ from 1 to 2 GeV is of crucial importance. The experimental data from non-central collisions, which are presented in Fig. 8 of Ref. [46], show no dependence of the source radius on m ⊥ in the interval. Since we are interested in rather peripheral collisions, where the source radii are sufficiently small and the effect we suggest to measure is significant, it is legitimate to assume that the source radius is independent of particle's transverse mass.
If the Coulomb interaction is absent but there is a short-range strong interaction, the wave function can be chosen, as proposed in [52], in the asymptotic scattering form where q ≡ |q| and f (q) is the s−wave (isotropic) scattering amplitude.
With the source function (35) and the wave function (36), the correlation function (33) equals The remaining integral needs to be taken numerically. The formula (37) has been repeatedly used to compute correlation functions of various two-particle systems. When one deals with charged particles, the formula (36) needs to be modified because the long-range electrostatic interaction influences both the incoming and outgoing waves. However, the Coulomb effect can be approximately taken into account [53] by multiplying the correlation function by the Gamow factor that equals where the sing + (−) is for the repelling (attracting) particles and a B is the Bohr radius of the pair.

Deuteron as a bound state
Let us now derive the h−D correlation function treating the deuteron as a neutronproton bound state created due to final-state interactions similarly to the h−D correlation. Then, the correlation function is defined as where p n = p p = p D /2 and q is the relative momentum of h and D. The deuteron formation rate A 2 is given by the formula (5) with A = 2. The correlation function multiplied by the deuteron formation rate equals where ψ hnp (r h , r n , r p ) is the wave function of a h−D system. The spin factor 3/4 has the same origin as that in Eq. (13).
To compute the integral in Eq. (40), we introduce the Jacobi variables of a three- with M ≡ m n + m p + m h , m D ≡ m n + m p and we write down the wave function as ψ hnp (r h , r n , r p ) = e iPR ψ q hD (r hD ) ϕ D (r np ), (42) where ψ q hD (r hD ) and ϕ D (r np ) are the wave functions of the relative motion of p and D and of internal motion of D, respectively.
Using the Gaussian source (6), the integral over the center-of-mass position R in Eq. (40) gives where D r (r) is again given by Eq. (35) and the normalized function D 3r (r) equals As a result of the integration over R in the right-hand-side of Eq. (40), the formation rate, which is given by Eq. (13) factors out. Consequently, the rate, which is also present in the left-hand-side of Eq. (40), drops out and the correlation function equals The formula (45) has the same form as (33) but the source function differs. When deuterons are directly emitted from the fireball as 'elementary' particles the radius of deuteron source is the same as the radius of proton source. When deuterons are formed only after emission of nucleons from the fireball, the source becomes bigger because the deuteron formation is a process of spatial extent. More quantitatively, the source radius of deuterons treated as bound states is bigger by the factor 4/3 ≈ 1.15 than that of 'elementary' deuterons.

p−D correlation function
We have first considered the K − −D correlation function which has appeared not very sensitive to a source radius because the strong interaction of K − and D is rather weak that is the scattering length is short. The p−D system is better suited for our purpose.
Since the p−D pair can have spin 1/2 or 3/2 there are two interaction channels. The s−wave scattering lengths of p−D scattering in the 1/2 and 3/2 channels are, respectively, 4.0 fm and 11.0 fm [54]. As nucleons are assumed unpolarized, the p−D correlation function is computed as the average where the weights factors 1/3 and 2/3 reflect the numbers of spin states in the two channels. The average p−D correlation function is shown in Fig. 2 for three values of the source radius which are chosen in such a way that R s = 2.00 fm = 4/3 · 1.73 fm = The dependence of the p − D correlation function on R s becomes weaker as R s grows. Consequently, the analysis of higher p T particles from non-central events, when the sources are relatively small, is preferred.
Our proposal to distinguish the scenario of deuterons directly emitted from the fireball from that of deuterons formed due to final-state interactions does not relay on an absolute value of the source size inferred from the p−D correlation function but on a comparison of source size parameters inferred from the p − D and p − p correlation functions. Therefore, systematic uncertainties of the femtoscopic method, both experimental and theoretical, are not of crucial importance here, as they are expected to influence in a similar way the source parameters inferred from the p−D and p−p correlation functions.
We note that the size of the proton source in pp collisions at LHC was measured with an experimental accuracy of 7% where the statistical error is only 2% [48]. Our proposal requires an accuracy better than 15% which, however, does not include systematic experimental and theoretical uncertainties. Therefore, the required accuracy of the measurement seems achievable.

p− 3 He correlations
The p− 3 He correlation function is interesting for two different reasons. It is needed to obtain the yield of 4 Li but at the same time the correlation function carries information about a source of 3 He. It allows one to determine whether 3 He is directly emitted from the fireball, as in the thermal model, or it is formed afterwards, as in the coalescence approach.

General formula of p− 3 He correlation function
If the nucleus 3 He is treated as an elementary particle emitted from a source, the p− 3 He correlation function is defined as and it is given by Eq. (33) where the source function (35) enters. The wave function φ q (r), however, describes not the relative motion of h and D but of p and 3 He. Taking into account that the nucleus 3 He is a bound state of (p, p, n) formed due to final-state interactions at the same time when the correlation among 3 He and p is generated, the correlation function is defined as where A 3 is the formation rate of a nucleus 3 He given by the formula (5) with A = 3. The product of the formation rate and correlation function is where D(r i ) with i = p, 1, 2, 3 is again the source function while ψ p 3 He (r p , r 1 , r 2 , r 3 ) is the wave function of p and 3 He. Let us compute A 3 . Using the Jacobi variables for a system of three particles with equal masses and writing down the wave function of 3 He as with φ3 He (x, y) being the wave function of relative motion, the formation rate equals where the normalized two-particle relative source function D r (x, y) is The latter equality holds for the Gaussian parametrization (6).
To derive the p− 3 He correlation function we use the Jacobi variables for a system of four particles defined by Eqs. (22) and we write down the wave function as where ϕ q (z) is the wave function of relative motion in the center of mass of p− 3 He system. Computing the integral (49), one finds that A 3 factors out and the correlation function equals where E 0 and Γ are the resonance energy and its width and the parameter λ R , which is assumed to be real, controls a strength of the resonance. In case of the 4 Li resonance, the energy difference, which enters the amplitude, is and the momentum q 0 , which corresponds to the resonance peak, is q 0 = 75.7 MeV. When compared to the original formula (134.12) from the textbook [56], we have introduced the parameter λ R and have replaced the factor (2l + 1)/q by (2l + 1)/q 0 , where q 0 corresponds to the energy E 0 , to avoid the divergence of the amplitude at q = 0. The modification is legitimate as, strictly speaking, the amplitude is valid only in the vicinity of the resonance. In our numerical calculations we assume that λ R = 1 but, as we discuss further on, the parameter λ R can be and should be inferred from experimental data.
The correlation function, which is computed with the source function (57) and takes into account the resonance interaction, is where R 0 (q) is the correlation function due to s−wave scattering only. For l = 1 the coefficient J l and the function K l (q) are In Fig. 3 we show the p− 3 He correlation function which takes into account only the resonance interaction. The peak at q = q 0 = 75.7 MeV is well seen. Fig. 4 shows the spin-average correlation function which takes into account the s−wave scattering, the resonance interaction and the Coulomb repulsion included by means of the Gamow factor (38). The resonance contributes only to the triplet correlation function because, as explained in Sec. 4.1, the spin of p and 3 He, which are constituents of 4 Li, equals unity. One observes that the s−wave scattering and Coulomb repulsion strongly deform the resonance peak.
A measurement of the p − 3 He correlation function in relativistic heavy-ion collisions at LHC is difficult but possible. The correlation function was measured in 40 Ar−induced reactions on 197 Au at the collision energy per nucleon of 60 MeV [57]. It was also measured in relativistic Au−Pt collisions at AGS [58].
One asks whether the p− 3 He correlation function is sensitive enough to the change of the source radius from R s to 3/2 R s that allows one to judge about the origin of 3 He: whether the nuclides are directly emitted from the source or they are formed due to final-state interactions. The answer obviously depends on accuracy of the p− 3 He correlation function to be experimentally measured but Fig. 4 suggests that it will be a difficult task. The problem is that the source radius R s must be inferred from experimental data together with the parameter λ R which controls the resonance strength.

Yield of 4 Li
As discussed in the previous section, the resonance peak of the p − 3 He correlation function is distorted by the Coulomb repulsion and s−wave scattering. So, it is not obvious how to infer the resonance yield from the distribution of the p− 3 He pairs.
To derive the yield of 4 Li, we start with the formula (47), which is written as where P ≡ p p + p3 He . The correlation function strongly depends on q but the dependence of the product dNp d 3 pp dN3 He d 3 p3 He on q is rather weak in the momentum domain of the resonance. So, it is taken at q = 0. To get the yield of 4 Li, one should sum up the number of correlated p− 3 He pairs within the resonance peak. However, the peak is deformed by the Coulomb repulsion and s−wave scattering. So, we suggest to fit an experimentally obtained correlation function with the theoretical formula (63) where λ R , which enters the amplitude (61) to control a strength of the resonance, is treated as a free parameter. Then, the contribution from the resonance can be disentangled.
Denoting the correlation function shown in Fig. 3, which differs from unity solely due to the resonance interaction, as R R (q), the yield of 4 Li of the momentum P equals dN4 Li where the factor 3/4 takes into account that 4 Li is produced only in the triplet channel and Since the source function is assumed to be isotropic and the correlation function depends on q only through q, the trivial angular integration has been performed in Eq. (69). The upper limit of the integral (69) should be chosen in a such way that the integral covers the resonance peak centered at q 0 ≈ 76 MeV/c. To get the ratio of the yields of 4 Li to 4 He one has to express the yields of 3 He and of protons through the yields of nucleons. Keeping in mind the coalescence formula (1), Eq. (68) is written as where the additional factor 1/2 takes into account that half of nucleons are protons. Consequently, the ratio of 4 Li to 4 He yields equals In Fig. 5 we show the ratio as a function of q max for λ R = 1 and four values of R s . There are indicted the values of relative momenta of 3 He and p (in the center-ofmass frame) which correspond to the energy of the resonance peak E 0 , to E 0 + Γ , to E 0 + 3Γ , etc. The integral (69) is seen to change rather slowly for q max bigger than, say, 150 MeV but it is not clear whether the integral saturates when q max → ∞. As observed in Ref. [38] and further studied in [39], the analogous integrals of correlation functions usually diverge as q max → ∞ because the correlation functions tend to unity as q −3 or slower. However, it is not physically reasonable to extend the integral (69) to a value of q max higher than, say, q max = 177 MeV which corresponds to E 0 + 3Γ . The value of S R does not change very much when q max is increased from 177 MeV to 286 MeV with the latter value corresponding to E 0 + 9Γ . Finally we note that for λ R = 1 and q max ≈ 150 MeV the ratio varies between 2 and 3.

Closing remarks
It is truly surprising result that production of light nuclei at the highest accessible energies of heavy-ion collisions is equally well described by two completely different models. Our objective was to broadly present the problem and to consider proposals to falsify one of the models. The measurements which are proposed are challenging but possible. Hopefully, we will learn weather the ideas discussed here are actually useful not in a remote future.