On nuclear coalescence in small interacting systems

The formation of light nuclei can be described as the coalescence of clusters of nucleons into nuclei. In the case of small interacting systems, such as dark matter and e+e-\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$e^+e^-$$\end{document} annihilations or pp collisions, the coalescence condition is often imposed only in momentum space and hence the size of the interaction region is neglected. On the other hand, in most coalescence models used for heavy ion collisions, the coalescence probability is controlled mainly by the size of the interaction region, while two-nucleon momentum correlations are either neglected or treated as collective flow. Recent experimental data from pp collisions at LHC have been interpreted as evidence for such collective behaviour, even in small interacting systems. We argue that these data are naturally explained in the framework of conventional QCD inspired event generators when both two-nucleon momentum correlations and the size of the hadronic emission volume are taken into account. To include both effects, we employ a per-event coalescence model based on the Wigner function representation of the produced nuclei states. This model reproduces well the source size for baryon emission and the coalescence factor B2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B_2$$\end{document} measured recently by the ALICE collaboration in pp collisions.


Introduction
The production mechanism for light clusters of nucleons, such as deuteron, helium-3, tritium and their antiparticles, in particle interactions has recently attracted increased attention from both the astroparticle and heavy ion communities. In heavy ion collisions, their small binding energies make these particles sensitive probes for two-nucleon correlations and density fluctuations, which may shed light on the QCD phase diagram [1]. These particles are also of particular interest in cosmic ray studies, because the expected low astrophysical backgrounds makes them ideal probes for exotic a e-mail: Michael.Kachelriess@ntnu.no (corresponding author) physics [2]. Furthermore, the sensitivities of the AMS-02 and GAPS experiments [3,4] are close to the expected fluxes of antideuterons from secondary production and, for optimistic parameters, from dark matter annihilations [5]. In order to correctly interpret the results of these experiments, a precise description of the production mechanism of light nuclei 1 is important.
In small interacting systems, such as dark matter and e + e − annihilations or pp collisions, the production of light nuclei is usually described by the coalescence model in momentum space [6][7][8], where nucleons originating from a particle collision merge to form a nucleus if their invariant momentum difference is smaller than the coalescence momentum p 0 . Traditionally, the yield of a nucleus with mass number A = Z + N and charge Z has been linked to the yields of protons p and neutrons n via the coalescence factor B A as Here, P A /A = p n = p p is the momentum of the nucleus and nucleons, respectively. In the limit of isotropic nucleon yields, the relation between B A and p 0 is This picture can be improved by taking into account twoparticle correlations provided by Monte Carlo event generators for strong interactions, as proposed in Refs. [9,10]. Such two-particle correlations are especially important in small interacting systems, since there the nucleon yields deviate strongly from isotropy. This approach is commonly used to predict the antinucleus yield in cosmic ray interactions, as well as from dark matter decays or annihilations [11][12][13][14][15][16][17][18][19][20][21][22], for a recent review see Ref. [5]. In order to be predictive, B A and p 0 must be independent of the centre-of-mass (c.m.) energy and the interaction process. However, the latter is not the case if the coalescence condition is only imposed in momentum space, since then the process dependent size of the formation region is neglected. An alternative scheme where the coalescence condition is imposed in position space is often employed for heavy ion collisions [23,24]. Here, the coalescence factor scales with the volume of the emission region of hadrons as B A ∝ V A−1 . Much efforts have been spent on unifying these pictures using, e.g., Wigner functions [25] and imposing the coalescence condition in phase space, see Ref. [26] for a review of early works. Such models differ mainly in the way the Wigner function of the nucleons is determined: The phase-space distributions of nucleons used in the coalescence models may be obtained, e.g., from transport models like the AMPT scheme [27] or hybrid schemes combining a hydrodynamical with a microscopic hadron cascade model [28]. Alternatively, analytical coalescence formula like the COAL-SH scheme [29] or statistical models which relate the phase-space volume at kinetic freeze-out to the entropy per nucleon have been proposed [30]. Finally, Refs. [31,32] have studied the influence of preclustering of baryons due to nucleon interactions on the coalescence process.
A key observation in all approaches relying on the phase space picture is that the coalescence probability depends on the size of the hadronic emission region, which can be measured in femtoscopy (often also called Hanburry-Twiss-Brown) correlation experiments [25]. This connection has recently been applied to pp collisions, both in cosmic ray [33] and LHC studies [34,35]. In particular, it was argued in Ref. [35] that the success of the femtoscopy analysis is strong evidence that coalescence is the major production mechanism of light nuclei. Moreover, these authors suggested that the use of experimental data from femtoscopy correlation experiments allows one to reliably predict the yield of light antinuclei in cosmic ray interactions, thereby avoiding the need of additional theoretical inputs.
The approaches discussed above are all based on the coalescence picture, but differ on how the coalescence condition is implemented and how the two-nucleon states are determined. In a competing approach one employs statistical thermal models [36][37][38][39][40][41][42]. Here one assumes that both hadronisation and the formation of light nuclei occurs as a chemical freeze-out process in a radially expanding "fireball" of a Quark-Gluon Plasma (QGP). These models are motivated by the observation that the spectra of light nuclei are consistent with a thermal distribution, with the same freeze-out temperature as for mesons and nucleons [37]. Intriguingly, experimental data from collisions of small systems, such as pp and pPb, show features characteristic for collective flows, or even for the formation of a QGP, see Ref. [43] for a review. It has therefore been suggested that the thermal production of light nuclei can be applied even to small interacting systems [5,39,44,45]. However, it is difficult to reconcile how the nuclei with their small binding energies survive the chemical freeze-out. Even more, the energy spectrum of the nucleons is in the coalescence picture inherited by the nuclei (up to a quantum mechanical correction factor [23]), and the apparent quasi-thermal spectra of light nuclei can therefore be explained by coalescence as well.
In Refs. [46,47], we developed a coalescence model based on the Wigner function representation of the produced nuclei states, which includes two-nucleon momentum correlations obtained from QCD inspired event generators (we will use the abbreviation WiFunC, i.e. Wigner Functions with Correlations, for this model). In this work, we argue that neither two-particle correlations nor the source size can be neglected when describing the cluster formation in small interacting systems 2 . Furthermore, we will use this model to describe the production of hadrons and nuclei in high energy pp collisions and compare it to recent experimental data by the ALICE collaboration on the size of the baryon emitting source [48] and on the multiplicity and transverse momentum dependence of the coalescence factor B 2 [45,49,50]. Both data sets have been interpreted as evidence of collective flows, but we will show that the same characteristics are described using QCD inspired event generators, like QGSJET II [51,52] and Pythia 8.2 [53,54]. Finally, we comment on the suggestion that femtoscopy data alone are sufficient to predict the yield of light antinuclei for astrophysical applications.
This paper treats several different topics related to the formation of nuclei by the coalescence mechanism in small interacting systems, with a focus on recent experimental data, and is structured as follows. We review the WiFunC model in Sect. 2 and its relation to the femtoscopy framework in Sect. 3. In Sect. 4, we compare our predictions for the size of the baryon emitting source to recent measurements of the ALICE collaboration in a femtoscopy experiment. In Sect. 5, the multiplicity and transverse momentum dependencies of the coalescence factor B 2 in pp collisions at 13 TeV, measured by the ALICE collaboration, are compared to the WiFunC model. In Sect. 6 we make comments on the use of isotropic models in astrophysical applications.

The quantum mechanics of coalescence and the WiFunC model
The WiFunC model is based on the quantum mechanical description of the coalescence process reviewed in, e.g., Refs. [25,39]. Here we will highlight only the main steps.
In this approach, the final state produced in a particle collision is described by a density matrix. Thus, one can find the deuteron spectrum in the sudden approximation by projecting the deuteron density matrix, ρ d = |φ d φ d |, onto the reduced density matrix ρ nucl = |ψ p ψ n ψ p ψ n | describing the coalecsing nucleons, By factoring out the c.m. motion of the deuteron, where the statistical factor 3/8 arises from averaging over spin and isospin and r ≡ r n − r p . Here, is the deuteron Wigner function, W np is the Wigner function of the two-nucleon state, and ϕ d is the internal deuteron wave function. If one approximates the deuteron wave function as a Gaussian, then D(r, However, apart from analytical estimates a more accurate wave function should be used, such as a two-Gaussian fit to the Hulthen wave function, chosen in Ref. [46]. To proceed, one has to specify the Wigner function W np in Eq. (4). One possibility is to use simulations in order to determine the phase-space distribution of nucleons. Both the perturbative and non-perturbative evolution in Monte Carlo generators of strong interactions are based on momentum eigenstates and, hence, they provide only information on momentum correlations of nucleons. The addition of spatial information requires thus the transition to a semi-classical picture. Alternatively, one can neglect two-nucleon correlations and assume an isotropic source, as it is often done when describing heavy ion collisions. Finally, one can derive two-particle correlations from experimental data. This is the approach used in the femtoscopy framework that will be discussed in the next section.
The first case is used in the WiFunC model [46] which combines two-nucleon momentum correlations obtained from QCD inspired event generators, with a simple analytical model for the spatial distribution of nucleons. Assuming a factorisation of the momentum and position dependence in the Wigner function, as well as neglecting spatial correlations, H np (r n , r p ) = h(r n )h(r p ), and choosing a Gaussian ansatz for the spatial distribution, Eq. (4) becomes The function ζ reflects the spatial distribution of the nucleons, and is thus clearly process dependent. It is in general given by Here we distinguished between the longitudinal and transverse spreads σ ,⊥ of the emission volume. The transverse spread is modified when boosting from the c.m. frame of the original particle collision to the deuteron frame. Thus γ is the Lorentz factor of the produced deuteron in the collider frame, while θ is the angle between the deuteron momentum and the beam axis. Note that, in contrast to our earlier treatement in Refs. [46,47], we have included the Lorentz boost in only one of the two transverse components: If the x y cordinates are rotated such that P d is contained in, e.g., the xz plane, then the σ y component will not be affected by the Lorentz boost.
Nucleon momentum correlations are provided by the event generator, while the process dependence is incorporated in the spread σ . The spread will in general have a geometrical contribution due to a finite spatial extension of the colliding particles, and a contribution related to the perturbative cascade and hadronisation, The geometrical contributions can be approximated as where R 1 and R 2 are the radii of the colliding particles, while the point-like contributions are given by σ (e ± ) R p 1 fm and σ ⊥(e ± ) Λ −1 QCD 1 fm. This simple picture is expected to give accurate results for pp interactions, while in the case of p A and A A collisions the geometrical contribution varies from event to event: While peripheral interactions which are dominated by binary collisions between a pair of projectile and target nucleons are characterised by σ (geom) R p , the size may increase to σ (geom) R A for the most central collisions. Consequently, the multiplicity of secondaries and the size of the source region are strongly correlated. Neglecting for the moment this correlation, and approximating the radius of a nucleus by with a 0 1.1 fm, allows us to use only one free parameter, whose physical interpretation is the size of the emission region of nucleons. Ideally, also the position integral in Eq. (4) should be evaluated event-by-event. It is therefore worth pointing out that some event generators like Pythia 8.2 have implemented semi-classical trajectories of the produced hadrons [55]. Thus, using Pythia one can instead directly evaluate relying on the semi-classical description of the spatial correlations provided by the simulation. A simple model applying a hard cut-off in both momentum and position space has been considered using UrQMD in Ref. [56]. The approach of the WiFunC model could be carried over in straight-forward way to these models, replacing the hard cutoffs with Eq. (14).
Because of the generality of Eq. (8), the WiFunC model can in principle be used to describe the production of other nucleus-like systems with small binding energies if the approximate wave function of the produced system is known. One additional application of the WiFunC model could therefore be the production of exotic bound states such as the X (3872) or the Z cs (3985), if they are deuteron-like bound states [57][58][59][60][61][62].

Relation to the femtoscopy framework
The emission volume probed in femtoscopy correlation experiments is directly linked to the distribution of nucleons, and can thus be used to check the validity of the WiFunC model. In a similar fashion, the emission volume can be related to the coalescence factor B A , as was done in Refs. [25,34,35]. However, in order to derive their analytic relationship, the so-called smoothness approximation [63] was applied on top of the sudden approximation used in the previous section. In this approximation, the q dependence in the nucleon Wigner function is assumed to be negligible so that the q integral in Eq. (4) can be evaluated. As remarked in Ref. [35], this may be justified for heavy ion collisions where the size of the produced nuclear clusters can be neglected compared to the size of the emitting source. However, a more careful treatment is warranted for small interacting systems. To see this, we note that applying the sudden approximation to Eq. (8) implies that two-nucleon correlations are neglected, but these correlations should be kept for small interacting systems. The WiFunC model evades these problems because it evaluates the momentum integral using the momentum distributions supplied by an event generator.
Within the smoothness approximation, the deuteron spectrum (4) can be written as while the nucleon spectra are given by 3 Following the authors of Refs. [34,35], we assume for simplicity E d /(E p E n ) = 2/m N in the deutron rest frame. Then the coalescence factor (1) becomes with the source function defined as Measured particles will always be affected by final state interactions. This significantly affects two-particle correlation experiments: Even from initially uncorrelated particles one will measure a correlation where S(r) is the emission source function and the final state interactions are encoded in the wave function Ψ [63]. This is very similar to Eqs. (8) and (17): Coalescence is effectively a final state interaction that affects the two-nucleon correlations.
The authors of Refs. [34,35] used Eq. (17) to derive numerical estimates of the B 2 factor as a function of the source radius r measured in femtoscopy experiments. This approach looks very promising, since it allows one to express the coalescence factor only in terms of measurable quantities. Unfortunately, any numerical evaluation is additionally based on three assumptions on the two-nucleon wave function: i) the spatial distribution has to be prescribed, ii) its characteristic size is assumed to be much larger than the one of the produced antinucleus states, such that the smoothness approximation can be used, iii) the two-nucleon momentum correlations are negligible. Yet, all these assumptions are generally not valid for collisions of small systems, as correctly noted already in Ref. [34]. Furthermore, the correlation function has to be inferred from experimental data, and is thus only available for the central rapidity region. The approximations required in the approach of Refs. [34,35] are avoided in the WiFunC model, since the used Monte Carlo generators provide two-nucleon momentum correlations which in turn leads to a non-trivial source function.

Size of baryon-emitting source
The source radius of the baryon emission in pp collisions at 13 TeV was recently measured by the ALICE collaboration, assuming a Gaussian source profile, using the femtoscopy framework, cf. with Eq. (19) of Ref. [48]. Here, the distance r = r p − r n between the two nucleons is defined in their pair rest frame. This study indicates that protons, antiprotons, Λ andΛ originate from the same source volume. Furthermore, a decrease in the source size with increasing transverse mass was observed. This decrease is often attributed to a collective flow, but is, as we will see next, also naturally described in the WiFunC model. Inserting the Gaussian ansatz for the spatial distribution of nucleons (7) into the expression (18) for the source leads to where we have taken into account that the Wigner functions and their spread, cf. with Eq. (9), are defined in the collider frame. Moreover, we chose the coordinate system such thatẑ is directed along the initial beam direction andŷ is perpendicular to bothẑ and P d . Furthermore, we used the identity m 2 T /m 2 = γ 2 sin 2 θ + cos 2 θ , m T being the transverse mass. Using the polar coordinates r x /r = sin ϕ sin ϑ and r y /r = cos ϕ sin ϑ, we find with and Hence the WiFunC model predicts a non-trivial source function described by a Gaussian source modified by the function I(r, m T , σ , σ ⊥ ).
In order to compare the predicted source function to the measurement by the ALICE collaboration, Eq. (22) must be compared to the Gaussian source profile (20) to fix r 0 (m T ). In order to determine r 0 (m T ), we perform a least-squares fit using as uncertainty μ ∝ 1/ √ S 2 (r ) as the expected Gaussian error. Additionally, we consider also a simple analytical approximation: By comparing the Taylor expansion of Eqs. (20) and (22), one finds In the analysis of the data on the source function in pp collisions at 13 TeV by ALICE only high multiplicity events (0-0.17% INEL > 0) were included [48]. However, the WiFunC model says that there is no (or only a weak) multiplicity dependence of the emission volume in pp collisions. In Fig. 1, we compare the source size r 0 estimated for protonproton pairs 4 , using both the exact source function (22) (blue solid line) and the approximation (25) (orange dashed line). Additionally, we show the source size obtained in the limit σ σ ⊥ (green dashed-dotted line), which corresponds to the steepest slope r 0 (m T ) possible in our model. It is worth noticing that the data tend to give better fits for σ > σ ⊥ , as expected from their physical interpretations. Even so, we find not yet any need to fit them separately due to the relatively large experimental uncertainties.
From Fig. 1 one can infer σ = (0.95 ± 0.1) fm. Intriguingly, the WiFunC model thus describes the data well with values of σ similar to those obtained in Refs. [46,47] by a fit to antideuteron measurements. More importantly, we have shown that the decrease of the source size with increasing transverse momentum, which is often attributed to collective flows, is correctly described by the WiFunC model using QCD inspired MC generators.

Multiplicity dependence of coalescence in small interacting systems
In the previous section, we focused on how the emission region of nucleons is related to the source size measured in femtoscopy experiments. Now we consider the effect of twoparticle correlations on the deuteron yield. To this end, we investigate how the coalescence factor B 2 of antideuterons measured at mid-rapidity (|y| < 0.5) in pp collisions at 13 TeV depends on multiplicity and transverse momentum 5 . The experimental results are reported for a specific event class (INEL > 0) and are divided into different multiplicity classes in terms of the percentage of the inclusive cross section, see Ref. [45] and references therein for their definition. We aim to reproduce the data, generating inelastic pp collisions at 13 TeV with QGSJET II and Pythia 8.2, while describing the coalescence by the WiFunC model with σ = 0.9 fm, using the two-Gaussian wave function for the deutron. We check the trigger condition and classify the multiplicity class on an event-by-event basis. For comparison, we consider the standard per-event coalescence model with a hard cutoff p 0 ∼ 0.2 GeV. This serves as a benchmark on what effects are caused by particle correlations, and what by the source size in the WiFunC model.
The results are compared to the experimental data in Fig. 2. Both QGSJET II and Pythia 8 reproduce well the overall yield in the various multiplicity classes. Furthermore, the qualitative behaviour of an increasing transverse momentum 5 We constrain this discussion to the data obtained at 13 TeV because of their small experimental uncertainties, but the same qualitative features are seen also at 7 TeV [50].

Fig. 2
The coalescence factor B 2 for different multiplicity classes, measured by the ALICE collaboration is compared to the predictions by QGSJET II (above) and Pythia 8.2 (below), using the WiFunC model (solid lines). The results for the standard coalescence model (dashed lines) are shown for comparison. Class I corresponds to largest multiplicities, while the multiplicity decreases with increasing class p T slope of B 2 with increasing multiplicity is also reproduced. This increase is often attributed to a collective flow, but our results indicate that it is also well described by the WiFunC model combined with QCD inspired event generators. While the overall behaviour and trends of the experimental data are reasonably well reproduced, deviations are expected as the event generators are not tuned to two-particle correlations. Comparing the results from the WiFunC model, shown as solid lines, to those of the standard coalescence model (dashed lines), one can notice that the multiplicity dependence of the slope of B 2 is stronger in the WiFunC model. Even so, there is also an increase in the slope of B 2 in the standard coalescence model, which is stronger in the case of Pythia. This indicates that two-particle correlations, although not the only effect responsible for the growing slope of B 2 , are not negligible for pp collisions in the kinematical range considered. In the WiFunC model, the multiplicity dependence emerges due to two-nucleon momentum correlations and the dependence of the emission region of nucleons on the event kinematics. In combination, these effects lead to the non-trivial multiplicity dependence visible in Fig. 2: For increasing multiplicity, the momentum phase space available for single nucleons will on average decrease, which implies an increased coalescence probability according to Eq. (8). The main multiplicity dependence of the emission region in pp collisions comes from the modification of the transverse spread by the Lorentz boost, as it can be seen from Eq. (9). In order to get a sense of this dependence, we plot in Fig. 3 the multiplicity dependence of the transverse spread using Pythia and QGSJET at 13 TeV. In both cases, σ ⊥ = 1 fm is used. Both event generators lead qualitatively to the same multiplicity dependence: The average transverse momentum increases with increasing number of produced particles, leading to a decrease in the transverse spread. Such an increase of the average p T with multiplicity has been observed by all experiments at LHC, being reasonably reproduced by Pythia (see, e.g., Ref. [64]) and leading to a gradual decrease ofσ ⊥ up to the rather high values of dN ch /dη. On the other hand, this effect is not properly described by QGSJET-II, in which case the decrease ofσ ⊥ is saturated already for relatively small values of dN ch /dη.

Astrophysical applications
Thus far, we have considered only particles at central (pseudo-) rapidity, which are accessible experimentally. The bulk of produced particles will, however, in general have large longitudinal momenta. In high energy collisions at the LHC, the use of a constant B 2 as function of p z is a good approximation. Therefore, one may naively expect this assumption to be a good approximation for astrophysical processes as well. This is, as we will discuss in this section, however not the case. Even so, an isotropic model with constant B A ( p z ) is still regularly applied in the literature to antinuclei production in proton-proton collisions [33,65,66].
Cosmic ray antideuterons are expected to originate in secondary production, i.e. in collisions between primary cosmic rays and the interstellar matter. The main contribution comes from protons with energies E prim ∼ 20-100 GeV colliding with protons in the interstellar matter, while the bulk of the produced antideuterons has kinetic energies per nucleon in the range T ∼ 2-20 GeV/n [47]. In order to check the validity of a constant B 2 ( p z ) for astrophysical applications, we therefore plot the coalescence parameter B 2 ( p z ) obtained using QGSJET II for primary energies 6 E prim = 50 and Fig. 3 Spread ofσ ⊥ as a function of the number of charged particles in the central pseudo-rapidity region, for QGSJET II (above) and Pythia 8.2 (below). The mean value ofσ ⊥ at each N ch and its standard deviation are shown in solid and dashed lines, respectively; the colour code shows the probability density of events with a givenσ ⊥ and dN ch /dη (η < 0.5) 100 GeV in Fig. 4 as function of the momentum p z in the lab frame. The range of B 2 determined using the femtoscopy framework in Ref. [33] is shown as a violet band. For comparison, we also show the coalescence factors B 2 obtained for √ s = 50 GeV and 13 TeV. In the case of collider energies, the values obtained agree well with the value inferred by femtoscopy experiments in Ref. [33]. At energies most relevant for astrophysical processes, however, the femtoscopy data at the LHC overestimate the coalescence parameter. More importantly, the coalescence parameter depends strongly on the longitudinal momentum at these energies 7 . In order to obtain the correct energy spectra of the produced antinuclei in astrophysical processes, a careful treatment taking into account two-particle correlations is therefore required. Fig. 4 The coalescence factor ford production, as a function of longitudinal momentum in the lab frame in pp collisions for various energies relevant for astrophysical processes and collider energies

Conclusions
The WiFunC model is a per-event coalescence model based on the Wigner function representation of the produced nuclei states, which allows one to account for both two-nucleon momentum correlations and the size of the hadronic emission volume. We have shown that this model reproduces well the source size for baryon emission and the coalescence factor B 2 measured recently by the ALICE collaboration in pp collisions. While these measurements have characteristics that are often attributed to the collective flow of the Quark-Gluon Plasma, our results show that the same properties are well reproduced describing the underlying physical processes by conventional QCD inspired event generators as QGSJET or Pythia. Finally, we have demonstrated that the coalescence parameter depends strongly on the longitudinal momentum for the energy range most relevant for astrophysical processes. Therefore, the use of a constant B A value in astrophysical applications should be abandoned.