A realistic coalescence model for deuteron production

A microscopic understanding of (anti)deuteron production in hadron–hadron collisions is the subject of many experimental and theoretical efforts in nuclear physics. This topic is also very relevant for astrophysics, since the rare production of antinuclei in our Universe could be a doorway to discover new physics. In this work, we describe a new coalescence afterburner for event generators based on the Wigner function formalism and we apply it to the (anti)deuteron case, taking into account a realistic particle emitting source. The model performance is validated using the EPOS and PYTHIA event generators applied to proton–proton collisions at the centre-of-mass energy s=13\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{s}= 13$$\end{document} TeV, triggered for high multiplicity events, and the experimental data measured by ALICE in the same collision system. The model relies on the direct measurement of the particle emitting source carried out by means of nucleon–nucleon femtoscopic correlations in the same collision system and energy. The resulting model is used to predict deuteron differential spectra assuming different deuteron wavefunctions within the Wigner function formalism. The predicted deuteron spectra show a clear sensitivity to the choice of the deuteron wavefunction. The Argonne v18\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{18}$$\end{document} wavefunction provides the best description of the experimental data. This model can now be used to study the production of (anti)deuterons over a wide range of collision energies and be extended to heavier nuclei.


Introduction
The microscopic understanding of the formation of light nuclei in high-energy collisions is a fundamental open problem that is being addressed since several decades both experimentally and theoretically.The main question is how nuclear bound states are formed and their structure emerges from the properties of the strong interaction and the laws of Quantum Chromodynamics.This problem is also relevant in astroparticle physics, where insights on the production mechanism are necessary to interpret future measurements of antinuclei in cosmic rays by spaceborne experiments searching for dark matter [1,2,3,4,5].Indeed, it has been suggested that, among the products of dark matter particle annihilation and decay, light antinuclei could be a promising signal to search for, since the only known background source are cosmic rays interactions with the interstellar medium.As cosmic rays and the interstellar medium mostly consist of hydrogen (90%) and helium (8%), and only in small percentage of heavier nuclei, most of the relevant interactions for the produc-tion of antinuclei in the Galaxy are proton-proton (pp) and proton-helium collisions.Studying these processes in the laboratory provides a unique opportunity to emulate under controlled conditions the production of antinuclei in the Galaxy.
From the theory side, the production of light (anti)nuclei is mainly described by the statistical hadronisation or the coalescence approach.Statistical hadronisation models (SHMs) assume that nuclei are emitted from a source at hadron-chemical and thermal equilibrium and their abundances are fixed at the so-called chemical freeze-out [25,26,27,28,29,30].Yields depend on the hadron mass, on the temperature and the volume of the system at freeze-out, and on the conservation of quantum numbers, which is imposed to hold on average if the volume of the produced system is large (e.g. in Pb-Pb collisions), or exactly if the volume is small (e.g. in pp collisions).The SHM is a macroscopic model: it predicts hadron yields requiring the equilibrium conditions, without any control on the microscopic production mechanism.Also, no detailed information on the structure of bound objects is contained in the model [31] and the predicted yields are integrated over transverse momentum.
On the contrary, coalescence aims at a microscopic description of the formation of (anti)nuclei, being the result of final state interactions between (anti)protons and (anti)neutrons.Coalescence was first introduced in the 1960s to explain the production of deuterons in protonnucleus collisions at the CERN SPS [32,33] and has been employed since then to describe the production of light nuclei in relativistic heavy-ion as well as pp collisions [34,35,36,2,37,38,39].In this approach, nucleons that are close to each other in space and momentum can coalesce to form a nucleus [32,33] when the system produced in high-energy hadronic or nuclear interactions undergoes freeze-out.The early works of Pearson and Butler [32], and Kapusta [33] emphasized the momentum dependence of the coalescence probability, by identifying the coalescence momentum p 0 as the main parameter governing the probability of cluster formation.Sato and Yazaki [34] introduced the density matrix formalism to calculate the nucleus formation probability, given by the integral over coordinate space of the projection of the proton and neutron wavefunctions over the deuteron one.They assumed no correlation between spatial and momentum distributions, nor any dynamical correlation between the proton and the neutron emitted from the particle source created after the hadron-hadron collision.
Following works [35,36] drew attention to the constraint in coordinate space to be considered in heavyion collisions, where the source is a spatially-extended and collectively-expanding system.In their seminal work, Scheibl and Heinz [36] included a detailed modelling of the system and pointed out that the source could be identified with the volume out of which particles are emitted with similar momenta, the size of which can be extracted by measuring two-particle momentum correlations with femtoscopy techniques [40,41,42].In recent years, a similar approach was employed to explain the production of (anti)nuclei in pp and p-Pb collisions [2,37], where the final-state particles are expected to occupy a small volume, of the order of the size of the proton (r p ≈ 1 fm) and even smaller than the typical size of a light nucleus (for deuteron, r d ≈ 2 fm).A key aspect of these works is having identified that the coalescence process depends on the size of the hadron emission region and more specifically, on the size of the nucleus relative to the size of the source.This relation provides the relevant length scale for the process [2,43,39].
Using a Wigner function representation of nucleons and nuclei, Blum et al. [37,39] and Mrowczynski et al. [38] obtained a relation between femtoscopic correlations and coalescence, being the continuum and discrete result of final state interactions among nucleons, respectively.Notably, they employed a Gaussian wavefunction for the deuteron, which allows for a fully analytic calculation of the coalescence parameter.All the mentioned models are based on the Wigner formalism and differ by the details of their implementation, including how the source size is accounted for.However, they have in common that they provide an analytical or numerical solution for the coalescence probability to be directly compared with experimental data.
Motivated by the need to model the production of light antinuclei in cosmic ray interactions for the estimate of their flux near Earth, recent developments of the coalescence model focus on pp collisions, where small particle sources are expected.In particular, the work by Kachelriess et al. [44,45] provides a solution to apply coalescence as an afterburner to particle production.In their approach, Monte Carlo (MC) event generators are used to obtain two-nucleon momentum correlations from the simulated pp collision events.To form a deuteron, the Wigner function-based coalescence is employed on an event-byevent basis in the WiFunC framework [45].The size of the formation region, a free parameter of the model, is taken as collision process-dependent.This parameter is fitted to ALICE data in [44], whereas it is predicted in [45] within WiFunC from two-particle momentum correlations, thus being limited to the accuracy of the description of the twonucleon correlations native to the generators.It is shown that the predicted size parameter is consistent with AL-ICE femtoscopy data for the baryon emitting source [46] and is close to 1 fm, from e + e − to pp collisions.
In our work, inspired by [45], we provide a coalescence afterburner that takes into account realistic particle emission and correlation, to be used to simulate event-by-event deuteron and antideuteron 1 production in high-energy hadronic collisions.We focus on pp collisions due to their relevance for the searches for antinuclei in cosmic rays as potential signatures for dark matter particles.We chose to test and validate our model by simulating deuteron production in high multiplicity pp collisions at √ s = 13 TeV because this corresponds to the only data sample for which simultaneous measurements of deuteron yields and baryon source radius are available [46,24].As it will be evident later, the size of the emitting source is a key ingredient of our implementation, allowing us to provide a realistic model for deuteron production.
First, we generate pp collision events using two distinct MC generators, EPOS 3 [48,49] and PYTHIA 8.3 [50], chosen for their known capability to reproduce most features of the LHC pp data.Then, the momentum distributions of the final state nucleons are tuned to data 1 We assume that the formation mechanism of antinuclei and nuclei in high-energy collisions is the same.This is justified by the present evidence that the force responsible for the nuclear binding, a residual of the strong force among partons, acts in the same way among antinucleons inside antinuclei as it does among nucleons inside matter nuclei [47].In the following, we omit the prefix "anti" for brevity.and are used as inputs for the afterburner.The latter is implemented based on a state-of-the-art Wigner-function coalescence approach, described in Sec. 2. Compared to previous approaches, we improve the generator-borne description of the nucleon-emitting source, as described in Sec. 3, by including a) a realistic account of its size derived from the measurement of two-baryon correlations, and b) a modelling of the cocktail of short-lived hadron resonances, which lead to a delayed nucleon emission.Using this improved framework, we predict the differential deuteron spectra with the Hulthén [41], Argonne v 18 [51] and Chiral Effective Field Theory (N 4 LO) [52] wavefunctions, moving past the traditional Gaussian wavefunction approach, and compare the results to the measured spectra.We discuss our results for the deuteron p T distributions as well as for the coalescence parameter B 2 in comparison to ALICE data in Sec. 4.

Wigner function formalism
As mentioned in Sec. 1, with the use of Wigner functions, it is possible to describe the production of nuclei via coalescence.In the process of deuteron formation via coalescence, the interactions of the nucleons of the pair with the rest of the particles are assumed to be subdominant due to the low particle density [42].Hence, the Lorentz-invariant yield of deuterons with momentum ⃗ P can be written as where Ψ d,P (x 1 , x 2 ) is the bound-state Bethe-Salpeter amplitude describing the deuteron, ρ 1,2 is the density matrix of the two nucleons, and S d = 3/8 is a factor that takes into account spin-isospin statistics.
We assume that the two-nucleon density matrix can be factorised into single-nucleon densities The single nucleon density ρ 1 can be written in terms of the single particle Wigner function f W 1 [37] Moreover, following the procedure described in [39], the deuteron Bethe-Salpeter amplitude can be written factoring out the motion of the deuteron where r d is the space-time position of the deuteron, P its four-momentum, and φ d (r) is the deuteron spatio-temporal wavefunction.
Hence, the deuteron spectrum takes the form where we define the relativistic internal Wigner density as Adapting the Wigner approach used in [45] to a fourdimensional space leads to 2π) 4 D(q, r)× × W np (P/2 + q, P/2 − q, r, r d ) , (7) where we defined Using the smoothness and equal-time approximations, as done in [39], in the pair rest frame (PRF) we obtain P = (M, ⃗ 0), q = (0, ⃗ q) and r = (t * , ⃗ r * ), and the Bethe-Salpeter amplitude becomes independent of time in the non-relativistic limit Ψ (q, r) → Ψ (⃗ q, ⃗ r * ) .
Defining t * 1 and t * 2 the emission time of the two nucleons in the PRF, in Eq. 7 one can write the relative distance of the nucleons as r = (t * 1 − t * 2 , ⃗ r ) and the center-of-mass coordinate of the nucleons as r d = (r 0 d , ⃗ r d ).Therefore, the integral over the four-momentum q becomes an integral over the momentum ⃗ q γ dN d where D(⃗ q, ⃗ r) is the Wigner density in a three-dimensional space.The time of kinetic freeze-out, r 0 d = t f , represents the moment in which the momentum of the final-state particles is fixed.Separating the space-and time-integrals, Eq. 10 becomes Assuming that t f is fixed and is the same for all particles, and considering that the particle yield is fixed at a common time (chemical freeze-out), the integral over t f can be omitted.In addition, due to the time equalisation in the PRF 2πδ(t * − t eq ), choosing arbitrarily t eq = 0 one obtains t * = 0. Hence, the integration over t * in Eq. 11 removes the dependence on t * , giving, as a result, a genuine three-dimensional equation The three-dimensional Wigner function of the deuteron D(⃗ q, ⃗ r) is defined as Notably, the choice of the deuteron wavefunction φ d affects only the term D(⃗ q, ⃗ r) in Eq. 12, while the other terms remain the same.The starting point in this theoretical derivation are the single free nucleon momentum distributions.In principle, there is no overlap between the deuteron Wigner function (D(⃗ q, ⃗ r)) and the free nucleon one (W np ).This problem arises from the conservation of energy and it can be resolved by introducing a third particle (usually a pion), see e.g.Ref. [36].However, as done in previous works [3,45,36], we make a semi-classical approximation where we assume the binding energy to be only a negligible correction, as it is much smaller than the mass scale of the nucleons (2.2 MeV ∼ E B ≪ m p = 0.938 GeV).
In [45] a factorisation of space and momentum dependence of the proton-neutron Wigner function is assumed where G np is the two-particle momentum distribution, taken from MC generators, containing the nucleon singleparticle momentum distributions and their initial-state correlation.Furthermore, for the space term H np , the spatial correlation is neglected, and the spatial single-particle distributions h(⃗ r p,n ) are assumed to be Gaussian, hence Here, ⃗ r d ≡ ⃗ r p + ⃗ r n , r 0 is the size of the two-particle emitting source, and ⃗ r ≡ ⃗ r p − ⃗ r n as before.In our work, possible space-momentum correlations and two-particle correlations in momentum and space coordinates at the hadron production point are considered in the source model discussed in Sec. 3. Finally, the coalescence probability4 P(r 0 , q) can be obtained by folding the spatial distribution of nucleons with the deuteron Wigner function Equations 16 and 17 are based on the assumption of a Gaussian source, for which the mean value (r µ ) of the two-particle distance distribution is related to the source size r 0 through r µ = (4/ √ π) r 0 .However, in particle generators the source has a shape that is generally not Gaussian [46].For this reason, we evaluate r µ from a fit to the distribution and then we obtain the source size r 0 through their relation.After this derivation, the deuteron momentum distribution in Eq. 12 assumes the final form 6  . ( 3 Source To account for realistic particle emission and correlation, in our implementation of the coalescence afterburner we rely on the ALICE measurement of the nucleon-emitting source [46], which has been performed differentially in transverse mass only in pp collisions at √ s = 13 TeV.As mentioned before, a measurement of the deuteron production spectra [24], which we use to validate our model, is available in the same system as discussed in Sec. 4.
To simulate pp events we employ two different MC generators, EPOS 3.117 and PYTHIA 8.3 with the Monash 2013 tune configuration.EPOS 3 [49] is a hybrid Monte-Carlo event generator in which the reaction volume is divided into "core" and "corona", depending on the local density and transverse momentum of the string segments.The core represents a thermalised bulk of matter evolving according to 3+1D viscous hydrodynamics and hadronises according to a Cooper-Frye mechanism [53,54].The particles in the corona originate from string fragmentation.As a function of final-state particle multiplicity, the relative fraction of core and corona evolves dynamically.On the other hand, PYTHIA 8.3 [50] is a parton-based microscopic event generator, where the primary process of a pp collision is represented with hard parton scatterings via 2 → 2 matrix elements defined at leading order.It is complemented by parton showers that include initial-and finalstate radiation via the leading-logarithm approximation.The hadronisation from partons is performed using the Lund string fragmentation model [55].The Monash 2013 tune [56] is chosen because it improves the description of minimum-bias and underlying-event observables in pp collisions at LHC energies and it includes a multi-parton interaction-based color-reconnection scheme.The underlying event in the model consists of particles originating from multi-parton interactions as well as from beam remnants.In the color reconnection picture, the strings between partons can be rearranged in a way that the total string length is reduced.
A proper modelling of the nucleon-emitting source, needs to take into account a) the overall final-state particle multiplicity, b) the possible contribution of feed-down from strongly-decaying resonances, and c) the phase-space correlations among the particles of interest.
Starting from point a), the MC simulations are required to reproduce the average charged-particle multiplicity density at midrapidity as measured.The considered ALICE data were collected with a high-multiplicity trigger, that corresponds to 0-0.01% of the total inelastic pp cross section.Following the nomenclature used in [24], this multiplicity class is referred to as HM I .In these events, the average multiplicity density is ⟨dN ch /dη⟩ |η|<0.5 = 35.8±0.5 [24].Therefore, the simulated events of interest are selected (triggered) based on the number of charged particles produced at backward and forward rapidity, mimicking the selection method of ALICE [24].The trigger acceptance factors are 10 −2 for EPOS 3 and 2 × 10 −5 for PYTHIA 8.3.
Regarding b), as described in detail in [46], short-lived resonances that decay into protons and neutrons, the particles of interest, play a crucial role in the determination of the source size.Indeed, the emitting source can be modeled as the convolution of a Gaussian core, which is the same for all baryons [46], and an exponential tail related to the decay of resonances.This means that a wrong resonance cocktail will influence the overall source size.For this reason, the relative fractions of the resonance cocktail of the MC simulations are tuned using a statistical hadronisation model, namely ThermalFIST [57].As input to the model, the measured temperature of chemical freezeout and the correlation volume [58] corresponding to the HM I class are provided.
The source size r 0 obtained directly from the event generators after the tuning of the resonance cocktail does not match the ALICE measurement [46] and hence it needs to be calibrated event-by-event.We reproduce the measured source size by acting at the particle-propagation level in the event generators, allowing for the conservation of the phase-space correlations provided by the event generator.In this regard, our approach differs from that used in [45], which employs an analytical form for r 0 , fitted to the radii measured by ALICE.
The starting point of our source model are the spacetime coordinates of the nucleons produced by the event generator.Figure 1 illustrates an example of the implementation for a primordial neutron paired with a proton stemming from a ∆ + resonance.To take into account that the particles in this pair are not necessarily created at the same time, the particle created earlier is propagated along its momentum for the time difference between the two.In Fig. 1, the neutron and the ∆ + resonance are depicted at time t 1 at distances r n (t 1 ) and r ∆ (t 1 ) from the production point, respectively, after the aforementioned time equalisation.The resonance decays after a time interval ∆t ∆ , during which the neutron moves.To estimate the distance between the neutron and the final-state proton, the neutron is propagated along its momentum for the time ∆t ∆ .The distance between the proton and the neutron d native pn is evaluated at the time t 1 +∆t ∆ .In the case in which both nucleons come from resonances, the one with the smallest time component is propagated along its momentum until the time equalisation is achieved.Finally, if a resonance decays into the particle of interest in a multi-step process (e.g. as it is the case for ∆(1900) ++ → N * (1440) + → p), ∆t ∆ is defined as the time difference between the last decay and the production of the first resonance.After the resonances decay and the final state particles have equal time (t 1 + ∆t ∆ ), the distance d native pn , in the protonneutron pair rest frame, and the average transverse mass ⟨m T ⟩ of the resulting pair 5   In summary, with this model, we succeed in preserving the space-momentum correlations among the nucleons, which were explicitly broken in the factorisation shown in Eq. 15, by reproducing the measured source size.This is shown in Fig. 2 for both EPOS 3 (red line) and PYTHIA 8.3 (blue line).Finally, the angular correlations of the nucleon pairs in the event generators do not match the ∆ϕ-∆η correlations measured by ALICE [59].However, reweighting them would destroy the nucleon space-momentum correlations provided by the event generators.More details on the angular correlations are given in Appendix B.

Results
The source modelling and the Wigner function formalism described in the previous sections are employed to obtain the deuteron spectra starting from nucleons simulated with EPOS 3 and PYTHIA 8.3.The proton and neutron 6 spectra are reweighted such to match the pro- 6 Neutrons and protons are assumed to have the same production spectra, since both belong to the same isospin doublet.Model / Data Fig. 3. Proton spectra generated by EPOS 3 and PYTHIA 8.3, compared with the ALICE measurement [24].In the bottom panel, the data-to-model ratios are shown.ton measurement by ALICE [24] (see Fig. 3) in order to start from realistic transverse momentum distributions.
Starting from the generated protons and neutrons, the event-by-event coalescence is implemented as a statistical rejection method.For each pair, we apply the source modelling described in Sec. 3 and calculate the coalescence probability for the single proton-neutron pair P √ π 4 d scaled pn , q using Eq. 17.The exact form of P depends on the deuteron wavefunction.For this study, we have considered three different wavefunctions: a Gaussian, Hulthén [41], and Argonne v 18 [51].Whereas the Gaussian is the simplest functional form of the wavefunction, the Hulthén and Argonne v 18 are based on physical properties of the deuteron.While the Argonne v 18 is expected to yield the best results, as it reproduces modern scattering data with a χ 2 ∼ 1, the Hulthén and the Gaussian are merely considered for comparison with previous studies.The shape of the deuteron wavefunction in the asymptotic region (r ≳ 1.5 fm) is determined by the solution of the Schrödinger equation for an interaction potential that reproduces the deuteron binding energy correctly (E B ∼ 2.2 MeV).In this range, Hulthén and Argonne v 18 are very similar, as shown in Fig. 8.In the range r ≲ 1.5 fm, the Hulthén and Argonne v 18 wavefunctions drastically differ, as a consequence of the different p-n potentials.Indeed, the Hulthén potential corresponds to an attractive interaction, while the Argonne v 18 contains a repulsive core in the interaction.A difference in the predicted yields computed with the two wavefunctions suggests that the Model / Data Fig. 4. Deuteron spectra obtained from EPOS 3 and PYTHIA 8.3, applying the coalescence model with different hypotheses for the deuteron wavefunction, compared with the ALICE measurement [24].The width of the bands represents the statistical uncertainty of the models.The systematic uncertainty of the spectra (6%) is not shown.In the bottom panel, the data-tomodel ratios are shown.
nuclear production mechanism is sensitive to the shortrange strong interaction between nucleons.In this work, the Wigner functions of the Hulthén and Argonne v 18 wavefunctions are computed for the first time, with details given in Appendix A. Figure 4 shows the deuteron spectra obtained with the different wavefunctions and different event generators.The results using the Argonne v 18 wavefunction are in excellent agreement with the data measured by ALICE, regardless of the event generator used.This proves that with our model, given the correct source size, nucleon spectra and average charged-particle multiplicity density, and a realistic wavefunction, it is possible to predict deuteron yields accurately.The systematic uncertainties of the model on the final deuteron spectra are estimated to be around 6%, independent of p T .The first source of systematic uncertainties is related to the source size.For this, we varied the source size used in the model by ±7%, based on the uncertainties reported in [46].The resulting systematic uncertainty is obtained from the relative deviation in the final spectra between the default source size and the varied one.The second source of systematic uncertainties is related to the fraction of primordial nucleons.To account for this, the primordial nucleon fraction is varied by ±10% and the relative deviation is the final spectra is considered.The two sources of uncertainties are summed in quadrature.In order to further test the impact of the core part of the strong-interaction potential on the deuteron yield, the deuteron wavefunction obtained from ab initio Chiral Effective Field Theory (χEFT) [52] (N 4 LO) calculations are employed to compute the deuteron yield and the results are compared with those obtained using the Argonne v 18 wavefunction (see Fig. 5).As for the Hulthén and Argonne v 18 wavefunctions, also the χEFT one is computed for the first time in this work, with details reported in Appendix A. On the one hand, the Argonne v 18 potential constructs the core part of the interaction using a combination of central, tensor, and spin-orbit interactions.On the other hand, the χEFT NN potential is derived systematically from the underlying theory of QCD and involves a perturbative expansion in powers of a small parameter related to the pion mass.The two-body interactions in χEFT are calculated using Feynman diagrams that involve pions and nucleons.The leading order (LO) twobody interactions involve only one-pion exchange, while higher-order (NLO, NNLO, etc.) interactions involve multiple pion exchanges and nucleon self-interactions.While Argonne v 18 and χEFT are two different approaches, the deuteron wavefunctions obtained from these approaches do not differ significantly, except for very short ranges.Indeed, the deuteron wavefunction from χEFT shows less repulsion than the one obtained from the Argonne v 18 potential.Nevertheless, the predicted deuteron yields using Argonne v 18 and χEFT are compatible, indicating that the production of the deuteron is not affected by extremely short-range interactions, as shown in Fig. 5.
Figure 6 shows the deuteron spectrum using the Argonne v 18 wavefunction, with and without the re-modelling of the source.Only the modeled source is able to describe the data.The difference between the two predictions is up to a factor of two, depending on p T .
The impact of the correlation between the relative momentum q and the distance between particles r taken from the event generator, shown in Fig. 10, on the final results is evaluated to be around 10%.For this, we compare the deuteron spectrum shown in Fig. 4 with one obtained by randomly sampling distances from the source size measurement shown in Fig. 2. Using the m T -dependent parameterisation of the source size, the final spectra change by around 30% with respect to the ones obtained with the native source provided by the generator, as shown in Fig. 6.Lastly, the difference between the results obtained using a realistic wavefunction and a Gaussian one is around 50%.
Using the spectra of protons and deuterons shown in Figs. 3 and 4, it is possible to compute the coalescence parameter B 2 , as The labels d and p indicate the deuteron and the proton, respectively, and the transverse momentum of protons is half of that of deuterons (p p T = p d T /2).The comparison among the B 2 measured by ALICE [24] and those obtained using EPOS 3 and PYTHIA 8.3, is shown in Fig. 7.A similar comparison in terms of B 2 was done in [24], where the coalescence predictions were obtained using the formalism described in [37].In that case, the Gaussian wavefunction provided the best description of data, while Hulthén overestimated the measurement by a factor of two.In comparison, in this work, the predictions using a Gaussian wavefunction underestimate the coalescence parameter by about 50 to 70%, while the predictions using the Hulthén wavefunction overestimate the measured B 2 by about 20 to 50%.Among different assumptions, the main difference between the coalescence predictions shown in [24] and in the present work is that in the formalism described in [37] the difference in momentum between the nucleons (⃗ q = (⃗ p p −⃗ p n )/2) is neglected in the Wigner function of the p-n state.The authors of [37] state that such approximation, motivated by the ease of calculations, is valid to an accuracy of around 10% in Pb-Pb collisions, while the accuracy in pp collisions is not estimated and could potentially be much larger.

Conclusions
In this paper, we show the implementation of a coalescence afterburner based on a state-of-the-art Wigner function formalism and use it to reproduce the (anti)deuteron spectra measured by ALICE in pp collisions at √ s = 13 TeV, collected with a high-multiplicity trigger.The novelty of our work consists in the preservation of the space-momentum correlation of nucleons, obtained by correcting the m T scaling of the source in the event generators with a parameterisation anchored to experimental measurements.At the same time, the Wigner function formalism with realistic deuteron wave functions is employed.The constraint to the measured source size allows for an accurate prediction of the (anti)deuteron spectra.In this work, three different hypotheses for the internal wavefunction of the deuteron are tested: a simple Gaussian, the Hulthén and the Argonne v 18 wavefunction.The Wigner function formalism is applied for the first time to the Hulthén and Argonne v 18 wavefunctions.The Argonne v 18 wavefunction, which is anchored to a realistic description of the nucleon-nucleon interactions, provides the best agreement with the deuteron spectra measured by AL-ICE.The predictions obtained with the Argonne v 18 wavefunction are compared with those obtained with a χEFT (N 4 LO) one and they are found to be in excellent agreement.This suggests that the production of the deuteron is not affected by extremely short-range interactions, where the two approaches differ.The good agreement between data and predictions shows that this model, is able to predict deuteron spectra if the correct proton spectra and source sizes are provided as input.Our work shows how important it is to measure the proton production spectra and the size of the emitting source simultaneously for different collision energies.This would allow for a prediction of the production yields of (anti)deuterons for different energies and multiplicities, providing a reliable estimation of antideuterons produced in the collisions of cosmic rays with the interstellar medium, which constitutes the background for the search for dark-matter annihilation with antinuclei in the final state.This work has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant Agreement No 950692).This work has been supported by the Deutsche Forschungsgemeinschaft through grant SFB 1258 "Neutrinos and Dark Matter in Astro-and Particle Physics".This research was supported by the Munich Institute for Astro-and Particle Physics (MIAPP) of the DFG cluster of excellence "Origin and Structure of the Universe".

A Deuteron Wavefunctions
In this Appendix, we report the calculation of Wigner densities D(⃗ r, ⃗ q) for different hypotheses of the deuteron wavefunction.Namely, we will consider a simple Gaussian, the Hulthén , the Argonne v 18 wavefunctions, and χEFT as shown in Fig. 8.All wavefunctions are normalised according to d 3 r |ϕ(⃗ r)| 2 = 1, thus the effect on the deuteron yields arises from the different shapes of the wavefunction.

A.1 Gaussian wavefunction
The most simple assumption is a single Gaussian wavefunction where d is a parameter related to the nucleus radius.For this calculation, d = 3.2 fm, as in [60].Using Eq. 13, the corresponding Wigner density is A.

Hulthén wavefunction
The Hulthén wavefunction represents a more realistic hypothesis for the deuteron wavefunction and it is based on the Yukawa theory of interaction.The wavefunction has the form where α = 0.23 fm −1 and β = 1.61 fm −1 are parameters taken from [41].For convenience, the Wigner density is calculated starting from the Fourier transform ψ( ⃗ k) of the wavefunction In Fourier space, Eq. 13 has the form Using the substitutions ⃗ k 2 = 2⃗ q + ⃗ k 1 and ⃗ k 1 = ⃗ k + ⃗ q, and integrating over ⃗ ξ and ⃗ k 2 , one obtains the following expression The integral depends on the angle θ between ⃗ r and ⃗ k.To eliminate this dependence on the angle θ, the angular average over θ is performed using sin(θ).With these simplifications, the Wigner density of the Hulthén wavefunction becomes A.3 Argonne v 18 wavefunction The Argonne v 18 potential is a phenomenological potential constrained to proton-neutron scattering measurements [51].
In such a potential, the deuteron wavefunction has the form where is the spin tensor, χ 1m is a spinor, and u(r) and w(r) are radial S and D wavefunctions, respectively.We define ⃗ r 1 the coordinates of the proton, ⃗ r 2 the coordinates of the neutron, ⃗ r = ⃗ r1−⃗ r2 2 , and ⃗ R = ⃗ r1+⃗ r2 2 .The spin averaged density for the deuteron is and the wavefunction is normalised as (29) In the previous integral, the contribution of the first addend is dominant, as the first part of the integral is equal to 0.9424 [51].Since the Argonne v 18 potential has only a numerical evaluation and no analytical form, an analytical form of its Wigner density is obtained through a fit to the numerical values of u(r)/r and w(r)/r.The fit is performed using the function where N 1 , N 2 , N 3 , a, b, c, and f are fit parameters.F (r) can describe both u(r)/r and w(r)/r individually.Therefore, u(r)/r and w(r)/r are fitted separately and two different sets of fit parameters are obtained for the S and D wave components, respectively.F (r) describes the shape of u(r)/r in the range 0 < r < 15 fm with a χ 2 ndf ∼ 6.83 • 10 −8 and w(r)/r with χ 2 ndf ∼ 1.3 • 10 −10 in the same range.The fit parameters for the S and D waves are reported in Tab. 1.The corresponding Wigner density for F (r) has the form where the terms T 1 , T 2 , T 3 , T 4 , T 5 , and T 6 are defined as In the previous equations, Ei is an exponential integral defined as Ei(x) =

∞
x dt e −t /t.

A.4 Chiral Effective Field Theory wavefunction
The Chiral Effective Field Theory (Chiral EFT or χEFT) is a theoretical framework used to study the low-energy QCD, such as atomic nuclei or hadrons (protons and neutrons).The Chiral EFT approach is systematic in the sense that the various contributions to a particular dynamical process can be arranged as an expansion in terms of a suitable expansion parameter.The expansion parameter is chosen to be the ratio of a typical low momentum (soft scale) to the chiral symmetry breaking scale (Λ χ ∼ 1 GeV, hard scale).The systematic expansion allows for the derivation of low-energy observables with controlled uncertainties.The deuteron wavefunction is obtained using the N N potentials through five orders of Chiral EFT, ranging from leading order (LO) to next-to-next-to-nextto-next-to-leading order (N 4 LO) [52], with the normalisation defined in [61].A cutoff at Λ c = 500 MeV is used.As in the case of Argonne v 18 , the wavefunction is composed by two components u(r) and w(r), which correspond to the radial S wave and to the radial D wave, respectively.The two components are shown separately in Fig. 8. Also in this case, only the numerical values of u(r)/r and w(r)/r are available and an analytic expression of the wavefunction is obtained with a fit, using the function where N 0 , N i , α i , β i , a, b, c, and r 0 are fit parameters.For w(r)/r, only the second term of F(r) is used since it is sufficient to describe the numerical values properly.u(r)/r and w(r)/r are fitted separately and the two sets of fit parameters are reported in Tab. 2. The extracted values of χ 2 ndf of the fits in the range 0 < r < 15 fm are ∼ 3.31 • 10 −8 and ∼ 3.23 • 10 −3 for u(r)/r and w(r)/r, respectively.The Wigner density of F (r) has the form dζ dγ sin(qζ) ζ× × (κ 0 + κ 1 + κ 2 + κ 3 ) , (39) where the terms κ 0 , κ 1 , κ 2 , and κ 3 are defined as

B Angular correlations
In this Appendix, we discuss the effect of angular correlations, namely ∆φ∆η, on the deuteron predictions.In Fig. 9 the ∆η-integrated ∆φ correlation function C(∆φ) of sign-like (anti)proton pairs measured by ALICE is shown and compared to predictions by EPOS 3 and PYTHIA 8.3 with the Monash 2013 tune [59].Note that, while the EPOS 3 prediction was obtained from pp collisions at √ s = 13 TeV, the ALICE measurement and PYTHIA prediction were obtained from pp collisions at √ s = 7 TeV.However, no qualitative difference between the predictions at these different energies is expected.It is noteworthy that at ∆φ close to zero ALICE measures a depletion in the correlation function, while both models predict a peak.On the contrary, the away-side peak around ∆φ = π is much more pronounced in the ALICE measurement compared to Monte Carlo studies.These discrepancies show a fundamental flaw in the production mechanism of baryons in these models.This means that it is impossible to properly correct these effects a posteriori.Instead, a rework of the hadron production mechanism inside the event generators would be required.A worst-case estimation of the effect on deuteron spectra can be performed by reweighting the p-n pairs according to the ratio of measured and predicted C(∆φ).Such a conservative estimate would lower the overall deuteron yield by no more than 20%.This effect is neglected in this study.C correlation of q and r Fig. 10 shows the distribution of relative momentum q and distance r for proton-neutron pairs, evaluated in the pair rest frame.Clearly visible is a positive correlations between q and r, where small relative momenta are preferred for pairs with small distances.This effect enhances the deuteron yield by roughly 10% since the phase-space region interesting for coalescence (q ≲ 0.5 GeV/c, r ≲ 2 fm) is more populated than a sample with no correlation.

Fig. 2 .
Fig.2.Comparison between the source size r0 measured by ALICE[46], the native ones for EPOS 3 (orange) and PYTHIA 8.3 (green) and those obtained after the source modelling (in red and blue for EPOS 3 and PYTHIA 8.3, respectively), as a function of the average transverse mass ⟨mT⟩ of the protonproton (antiproton-antiproton) pair.For the ALICE data, statistical and systematic uncertainties are summed in quadrature and shown as vertical bars, while for EPOS 3 and PYTHIA 8.3 uncertainties are negligible.

Fig. 6 .
Fig.6.Deuteron spectra obtained with EPOS 3, with the source modelling (model) and without (native), for the same wavefunction hypothesis, i.e.Argonne v18.Predictions are compared with the corresponding ALICE measurement[24].

Fig. 10 .
Fig. 10.Distribution of relative momenta and distances of proton-neutron pairs in EPOS 3.

Table 1 .
Fit parameters for F (r) obtained from the numeric values of u(r)/r (2 nd column ) and w(r)/r (3 rd column) in the range 0 < r < 15 fm.

Table 2 .
Fit parameters for F(r) obtained from the numeric values of u(r)/r (2 nd column ) and w(r)/r (3 rd column) in the range 0 < r < 15 fm.