Alternative coalescence model for deuteron, tritium, helium-3 and their antinuclei

TL;DR

This paper introduces a microscopic coalescence model combining Monte Carlo simulations and Wigner function representations to improve predictions of antinuclei production, relevant for dark matter detection, with a universal parameter fitted to collider data.

Contribution

It presents a new coalescence model integrating microscopic details and momentum correlations, fixed by collider data, enhancing the accuracy of antinuclei production predictions.

Findings

Model accurately describes antinuclei production in collider experiments.

Single universal parameter effectively fits multiple reaction types.

Improved microscopic understanding of antinuclei formation processes.

Abstract

Antideuteron and antihelium nuclei have been proposed as a detection channel for dark matter annihilations and decays in the Milky Way, due to the low astrophysical background expected. To estimate both the signal for various dark matter models and the astrophysical background, one employs usually the coalescence model in a Monte Carlo framework. This allows one to treat the production of antinuclei on an event-by-event basis, taking thereby into account momentum correlations between the antinucleons involved in the process. This approach lacks however an underlying microscopic picture, and the numerical value of the coalescence parameter obtained from fits to different reactions varies considerably. Here we propose instead to combine event-by-event Monte Carlo simulations with a microscopic coalescence picture based on the Wigner function representations of the produced antinuclei states. This approach allows us to include in a semi-classical picture both the size of the formation region, which is process dependent, and momentum correlations. The model contains a single, universal parameter which is fixed by fitting the production spectra of antideuterons in proton-proton interactions, measured at the Large Hadron Collider. Using this value, the model describes well the production of various antinuclei both in electron-positron annihilation and in proton-proton collisions.

Tables

Table 1. Table 1: Fit results for the constant ζ 𝜁 \zeta factor, in comparison to the old model.

Experiment	one-Gaussian		two-Gaussian		Old model
	$\sigma$ [$\frac{1}{\mathrm{GeV}}$]	$\frac{{\chi }^2}{N-1}$	$\sigma$ [$\frac{1}{\mathrm{GeV}}$]	$\frac{{\chi }^2}{N-1}$	$p_0$ [MeV]	$\frac{{\chi }^2}{N-1}$
ALICE 0.9 TeV	$3.5\pm 0.7$	7.5/2	$6.2\pm 0.3$	6.0/2	181	7.3/2
ALICE 2.76 TeV	$4.3\pm 0.3$	44/6	$6.6\pm 0.1$	32/6	174	45.6/6
ALICE 7 TeV	$4.1\pm 0.2$	182/19	$6.6\pm 0.1$	133/19	176	177/19
ALICE combined	$4.1\pm 0.1$	235/29	$6.6\pm 0.1$	172/29	176	229/19
ALICE helium-3	$4.5\pm 0.9$	1.7/2	-	-	179	1.2/2
ALEPH	$0_{-0}^{+2.3}$	-	${5.0}_{-0.6}^{+0.9}$	-	${214}_{-26}^{+21}$	-
ALEPH + OPAL	$0_{-0}^{+4.4}$	3.2/1	${5.5}_{-1.1}^{+1.3}$	3.2/1	201	3.2/1

Table 2. Table 2: Fit results obtained taking into account the modification of σ ⟂ subscript 𝜎 perpendicular-to \sigma_{\perp} by transverse boosts.

Experiment	one-Gaussian		two-Gaussian
	$\sigma$ [$\frac{1}{\mathrm{GeV}}$]	$\frac{{\chi }^2}{N-1}$	$\sigma$ [$\frac{1}{\mathrm{GeV}}$]	$\frac{{\chi }^2}{N-1}$
ALICE 0.9 TeV	$3.9\pm 0.8$	6.7/2	$6.9\pm 0.3$	2.6/2
ALICE 2.76 TeV	$4.9\pm 0.3$	35/6	$7.5\pm 0.1$	8.6/6
ALICE 7 TeV	$4.7\pm 0.2$	143/19	$7.6\pm 0.1$	29/19
ALICE combined	$4.7\pm 0.4$	186/29	$7.6\pm 0.1$	45/29
ALICE helium-3	$5.2\pm 1.0$	1.1/2	-	-
ALEPH	$0_{-0}^{+2.4}$	-	${5.3}_{-0.6}^{+1.0}$	-
ALEPH + OPAL	$0_{-0}^{+4.6}$	3.2/1	${5.8}_{-1.1}^{+1.4}$	3.2/1

Equations

$$E_A\text{derivative}[3]N_A{P}_A=B_A{\left(E_p\text{derivative}[3]N_p{P}_p\right)}^Z{{\left(E_n\text{derivative}[3]N_n{P}_n\right)}^N|}_{P_p=P_n=P_A/A}.$$

$$B_A=A{\left(\frac{4\pi }{3}\frac{p_0^3}{m_N}\right)}^{A-1},$$

$$\text{derivative}[3]N_d{P}_d=\text{tr}{\rho }_d{\rho }_{\mathrm{nucl}}.$$

$$⟨{\psi }_n{\psi }_p||{\psi }_p{\psi }_n⟩=N_p{N}_n.$$

$$\text{derivative}[3]N_d{P}_d=S\int \text{differential}[3]x_1\text{differential}[3]x_2\text{differential}[3]x_1^{\prime }\text{differential}[3]x_2^{\prime }{\varphi }_d^{\ast }({\mathbf{x}}_1,{\mathbf{x}}_2){\varphi }_d({\mathbf{x}}_1^{\prime },{\mathbf{x}}_2^{\prime })⟨{\psi }_n^{†}({\mathbf{x}}_2^{\prime }){\psi }_p^{†}({\mathbf{x}}_1^{\prime }){\psi }_p({\mathbf{x}}_1){\psi }_n({\mathbf{x}}_2)⟩,$$

$${\varphi }_d({\mathbf{x}}_1,{\mathbf{x}}_2)=(2\pi {)}^{-3/2}\exp \mathrm{i}{\mathbf{P}}_d\cdot ({\mathbf{x}}_1+{\mathbf{x}}_2)/2{\phi }_d({\mathbf{x}}_1-{\mathbf{x}}_2).$$

$$⟨{\psi }_n({\mathbf{x}}_2^{\prime }{)}^{†}{\psi }_p({\mathbf{x}}_1^{\prime }{)}^{†}{\psi }_p({\mathbf{x}}_1){\psi }_n({\mathbf{x}}_2)⟩=$$

$$\times \exp [\mathrm{i}{\mathbf{p}}_n\cdot ({\mathbf{x}}_2-{\mathbf{x}}_2^{\prime })]\exp [\mathrm{i}{\mathbf{p}}_p\cdot ({\mathbf{x}}_1-{\mathbf{x}}_1^{\prime })].$$

$$\text{derivative}[3]N_d{P}_d=\frac{S}{(2\pi {)}^6}\int \text{differential}[3]q\int \text{differential}[3]r_p\text{differential}[3]r_n\mathcal{D}(\mathbf{r},\mathbf{q})W_{np}({\mathbf{P}}_d/2+\mathbf{q},{\mathbf{P}}_d/2-\mathbf{q},{\mathbf{r}}_n,{\mathbf{r}}_p),$$

$$\mathcal{D}(\mathbf{r},\mathbf{q})=\int \text{differential}[3]\xi \exp -\mathrm{i}\mathbf{q}\cdot \mathbf{\xi }{\phi }_d(\mathbf{r}+\mathbf{\xi }/2){\phi }_d^{\ast }(\mathbf{r}-\mathbf{\xi }/2)$$

$${\phi }_d(\mathbf{r})={\left(\pi d^2\right)}^{-3/4}\exp -\frac{r^2}{2d^2},$$

$$\mathcal{D}(\mathbf{r},\mathbf{q})=8{\mathrm{e}}^{-r^2/d^2}{\mathrm{e}}^{-q^2{d}^2}.$$

$$\int {\mathrm{d}}^3{r}_p{\mathrm{d}}^3{r}_n{W}_{np}({\mathbf{p}}_n,{\mathbf{p}}_p,{\mathbf{r}}_n,{\mathbf{r}}_p)=N_p{N}_n|{\psi }_{np}({\mathbf{p}}_n,{\mathbf{p}}_p){|}^2\equiv G_{np}({\mathbf{p}}_n,{\mathbf{p}}_p),$$

$$W_{np}({\mathbf{P}}_d/2+\mathbf{q},{\mathbf{P}}_d/2-\mathbf{q},{\mathbf{r}}_n,{\mathbf{r}}_p)=H_{np}({\mathbf{r}}_n,{\mathbf{r}}_p)G_{np}({\mathbf{P}}_d/2+\mathbf{q},{\mathbf{P}}_d/2-\mathbf{q}).$$

$$h(\mathbf{r})={\left(2\pi {\sigma }^2\right)}^{-3/2}\exp -\frac{r^2}{2{\sigma }^2}.$$

$$\text{derivative}[3]N_d{P}_d=\frac{3\zeta }{(2\pi {)}^6}\int \text{differential}[3]q{\mathrm{e}}^{-q^2{d}^2}{G}_{np}({\mathbf{P}}_d/2+\mathbf{q},{\mathbf{P}}_d/2-\mathbf{q}),$$

$$\zeta \equiv {\left(\frac{d^2}{d^2+4{\sigma }^2}\right)}^{3/2}\le 1$$

$$h(\mathbf{r})\propto \exp -\frac{{\mathbf{r}}_{\parallel }^2}{2{\sigma }_{\parallel }^2}-\frac{{\mathbf{r}}_{\perp }^2}{2{\sigma }_{\perp }^2}.$$

$$L_{\mathrm{had}}\sim \gamma L_0,$$

$${\sigma }_{\parallel (pp)}^2={\sigma }_{\parallel (e^{\pm })}^2+{\sigma }_{\parallel (\mathrm{geom})}^2\approx 2{\sigma }_{\parallel (e^{\pm })}^2.$$

$${\sigma }_{(\mathrm{geom})}^2=⟨r_1^2{⟩}_O-⟨{\mathbf{r}}_1{⟩}_O^2.$$

$$⟨A⟩=\frac{\int \text{differential}[2]r_1\text{differential}[2]r_{2}A{\rho }_1(r_1){\rho }_2(r_2)w_{\mathrm{int}}(|\mathbf{b}-{\mathbf{r}}_1+{\mathbf{r}}_2|)}{\int \text{differential}[2]r_1\text{differential}[2]r_2{\rho }_1(r_1){\rho }_2(r_2)w_{\mathrm{int}}(|\mathbf{b}-{\mathbf{r}}_1+{\mathbf{r}}_2|)},$$

$$w_{\mathrm{int}}(|\mathbf{b}-{\mathbf{r}}_1+{\mathbf{r}}_2|)\propto {\delta }^{(2)}(\mathbf{b}-{\mathbf{r}}_1+{\mathbf{r}}_2),$$

$${\rho }_i=1/(\pi R_i^2)\exp -r^2/R_i^2,$$

$${\sigma }_{\perp (\mathrm{geom})}^2=\frac{R_1^2{R}_2^2}{R_1^2+R_2^2}.$$

$${\sigma }_{\perp (pp)}^2={\sigma }_{\perp (e^{\pm })}^2+{\sigma }_{\perp (\mathrm{geom})}^2\approx 2{\sigma }_{\perp (e^{\pm })}^2,$$

$${\sigma }_{\parallel (pp)}={\sigma }_{\perp (pp)}={\sigma }_{(pp)}=\sqrt{2}{\sigma }_{(e^{\pm })}.$$

$${\tilde{\sigma }}_{\perp }=\frac{{\sigma }_{\perp }}{\sqrt{{\cos }^2\theta +{\gamma }^2{\sin }^2\theta }}.$$

$$\zeta =\frac{d^2}{d^2+4{\tilde{\sigma }}_{\perp }^2}\sqrt{\frac{d^2}{d^2+4{\sigma }_{\parallel }^2}}.$$

$$w=3\zeta {\mathrm{e}}^{-q^2{d}^2}$$

$$N_d=$$

Full text

††journal: Eur. Phys. J. A

1 Introduction

Antideuteron and antihelium nuclei have been suggested as promising detection channels for dark matter, because of the low astrophysical background expected for such signatures Donato:1999gy : The dominant background source of antideuterons are cosmic ray protons interacting with the interstellar medium. The high threshold energy for this reaction channel implies that the antideuterons produced by cosmic rays have relatively large kinetic energies. Low-velocity antideuterons are therefore an ideal tool to search for exotic sources of antimatter. In the case of antihelium nuclei, the suppression of astrophysical backgrounds at low velocities is even stronger, but the maximal event rates expected in dark matter models are challenging for square-meter sized detectors. At present, the search for antinuclei is performed by the AMS-02 experiment on board of the International Space Station, while the GAPS balloon experiment is planned to fly in the next Solar minimum period around 2020 or 2021 Battiston:2008zza ; Aramaki:2015laa .

The production of light clusters of antinuclei like antideuteron, antihelium or antitritium111Since our discussion applies equally well to the production of particles and of antiparticles in $pp$ and $e^{+}e^{-}$ collisions, the preposition ‘anti’ is dropped further on. is usually described in the context of coalescence Schwarzschild:1963zz ; butler ; Sun:2015jta ; Zhu:2015voa ; Zhu:2017zlb or of statistical-thermal models Acharya:2017bso ; Andronic:2017pug ; Vovchenko:2018fiy ; Bellini:2018epz ; Chen:2018tnh ; Xu:2018jff ; Oliinychenko:2018ugs . In coalescence models, cluster formation has been traditionally parametrised by an invariant coalescence factor $B_{A}$ which relates the invariant yield $E_{A}\derivative*[3]{N_{A}}{P_{A}}$ of nuclei with mass number $A$ formed out of $Z$ protons and $N$ neutrons to the invariant yields $E_{i}\derivative*[3]{N_{i}}{P_{i}}$ of protons ($i=p$) and neutrons ($i=n$) via

[TABLE]

In $e^{+}e^{-}$ and $pp$ collisions, one imposes typically the coalescence condition in momentum space, requiring that the momenta of merging nucleons in their two-body center-of-mass (CoM) system are smaller than some critical value $p_{0}$. In the limit of isotropic and equal proton and neutron yields, the so-called coalescence momentum $p_{0}$ is related to $B_{A}$ via

[TABLE]

where $m_{N}$ denotes the nucleon mass. This scheme can be improved taking into account the momentum correlations between nucleons, which are provided by Monte Carlo (MC) simulations on an event-by-event basis. Such an approach, which was first suggested in Refs. Kadastik:2009ts ; Dal11 , is commonly used for the prediction of the antideuteron yield both from dark matter annihilations or decays and from cosmic rays interactions Cui:2010ud ; Ibarra:2012cc ; Fornengo:2013osa ; Dal:2014nda ; Grefe:2015jva ; Herms:2016vop ; Korsmeier:2017xzj ; Coogan:2017pwt ; Poulin:2018wzu , for a review see Ref. Aramaki:2015pii . The only free parameter of this model is the coalescence momentum $p_{0}$, which should be independent of the reaction type and the center-of-mass energy $\sqrt{s}$ in order to be predictive. However, the numerical value of the coalescence parameter obtained from fits to different reactions varies considerably Ibarra:2012cc ; Aramaki:2015pii .

An alternative scheme was developed to describe the formation of light nuclear clusters in heavy-ion collisions. There, the coalescence condition was imposed in coordinate space, assuming that the coalescence factor $B_{A}$ of a cluster with mass number $A$ is proportional to $V^{A-1}$, where $V$ denotes the volume of the emission region of hadrons from the expanding cloud of partons Csernai:1986qf ; Nagle:1996vp . There have been considerable efforts to combine these two approaches and to develop coalescence models which are based on a microscopical picture using, e.g., Wigner functions Sato:1981ez or a diagrammatic approach Duperray:2002pj . Many of these attempts impose the coalescence condition in phase space, using either a classical or quantum mechanical description as a starting point. Such models have been mainly applied to heavy ion collisions and are reviewed, e.g., in Refs. Danielewicz:1991dh ; Scheibl:1998tk ; Chen:2018tnh . An interesting application of this approach to the prediction of antideuteron and antihelium production by cosmic rays has been made recently in Ref. Blum:2017qnn .

The coalescence process has also been modelled as a dynamical process where the formation probability of a deuteron is proportional to the scattering cross section of the reaction $\bar{N}_{1}\bar{N}_{2}\to\bar{d}X$. The amplitude for such processes has been derived, e.g., from models for the non-relativistic nucleon-nucleon potential. As an alternative, Ref. Dal:2015sha used experimental data to determine the cross sections $\bar{p}\bar{n}\to\bar{d}X$ for $X=\{\gamma,\pi^{0},\ldots\}$. The coalescence probability was then determined as $\sigma_{\rm tot}(\bar{N}_{1}\bar{N}_{2}\to\bar{d}X)/\sigma_{0}$ with $\sigma_{0}$ as a free parameter. As a result, antideuterons were mainly produced with momenta close to the delta resonance, $\sim 1$ GeV, and the fit to the antideuteron production data in $pp$ collision data from the ALICE experiment improved significantly.

In this work, we develop a coalescence model for the formation of light nuclei in $e^{+}e^{-}$ and $pp$ collisions, which can be applied as well to dark matter annihilations and decays. Such a combined approach which describes successfully the production of light nuclei both in “point-like” reactions ($e^{+}e^{-}$ , dark matter) and in $pp$ collisions is especially needed for indirect dark matter searches, where the consistent prediction of “signal” antinuclei produced by dark matter annihilations and of background antinuclei created in cosmic ray interactions is required. Our model is based on the Wigner function representations of the produced antinuclei states and allows us to include in a semi-classical picture both the size of the formation region, which is process dependent, and momentum correlations. The model contains a single, universal parameter which is fixed by fitting the production spectra of antideuterons in $pp$ collisions, measured by the ALICE experiment at the LHC Acharya:2017fvb . The obtained value, $\sigma_{(e^{+}e^{-})}=\sigma_{(pp)}/\sqrt{2}\simeq 5\,{\rm GeV}^{-1}\simeq 1$ fm, agrees with its physical interpretation as the size of the formation region of the light nuclei. Using this value, the model describes well the data on the production of antihelium in $pp$ interactions and of antideuterons in $e^{+}e^{-}$ annihilations at the $Z$-resonance Schael:2006fd ; Akers:1995az .

2 Wigner function based deuteron formation model

We develop our model first for the case of deuteron production. The generalization to helium-3 and tritium is straightforward and will be performed in the next section. In the following, we use the fact that the binding energy $B$ of these nuclei is small, e.g. $B\simeq 2.2$ MeV for the deuteron. Therefore we can assume that a nucleus $A$ is formed through the process $N_{1}+\ldots+N_{n}\to A^{*}$, and that the excitation energy is later released by the emission of a photon.

2.1 Derivation

The starting point for the derivation of our new coalescence model is inspired by the approach using Wigner functions, presented in Ref. Scheibl:1998tk . We consider a system consisting of a proton and a neutron in a frame where the motion of their CoM is nonrelativistic. The number of deuterons with a given momentum $\boldsymbol{P}_{d}$ can be found by projecting the deuteron density matrix $\rho_{d}$ onto the two-nucleon density matrix $\rho_{\mathrm{nucl}}$,

[TABLE]

The deuteron density matrix describes a pure state, $\rho_{d}=\ket{\phi_{d}}\bra{\phi_{d}}$. The spin and isospin values of the two-nucleon state can be taken care of by introducing a statistical factor $S=3/8$ Mattiello:1996gq , such that the two-nucleon density matrix can be written as $\rho_{\mathrm{nucl}}=\ket{\psi_{p}\psi_{n}}\bra{\psi_{n}\psi_{p}}$ and is normalized as

[TABLE]

Here, $N_{p}$ and $N_{n}$ are the average multiplicities of protons and neutrons per event, respectively.222We neglect for the moment the double counting of nucleons involved in different pairs.

By evaluating the trace in the coordinate representation $\ket{\boldsymbol{x}_{1}\boldsymbol{x}_{2}}$, where the two indices refer to the positions of the two nucleons, one finds

[TABLE]

where $\phi_{d}(\boldsymbol{x}_{1},\boldsymbol{x}_{2})$ and $\psi_{i}(\boldsymbol{x})$ are the wave functions of the deuteron and nucleon $i$, respectively. Next we factorise the deuteron wave function into a plane wave describing the CoM motion with momentum $\boldsymbol{P}_{\mathrm{d}}$ and an internal wave function $\varphi_{d}$,

[TABLE]

Then we replace the two-nucleon density matrix by its two-body Wigner function,

[TABLE]

Further, we introduce as new coordinates the “average” positions of the proton and neutron, $\boldsymbol{r}_{p}=(\boldsymbol{x}_{1}+\boldsymbol{x}_{1}^{\prime})/2$ and $\boldsymbol{r}_{n}=(\boldsymbol{x}_{2}+\boldsymbol{x}_{2}^{\prime})/2$, as well as their separation $\boldsymbol{r}=\boldsymbol{r}_{n}-\boldsymbol{r}_{p}$, $\boldsymbol{\xi}=\boldsymbol{x}_{1}-\boldsymbol{x}_{1}^{\prime}-\boldsymbol{x}_{2}+\boldsymbol{x}_{2}^{\prime}$ and $\boldsymbol{\rho}=(\boldsymbol{x}_{1}-\boldsymbol{x}_{1}^{\prime}+\boldsymbol{x}_{2}-\boldsymbol{x}_{2}^{\prime})/2$. Changing also the momentum integration variables to $\boldsymbol{p}=\boldsymbol{p}_{n}+\boldsymbol{p}_{p}$ and $\boldsymbol{q}=(\boldsymbol{p}_{n}-\boldsymbol{p}_{p})/2$, and performing then the $\boldsymbol{\rho}$ and $\boldsymbol{p}$ integrals, we arrive at

[TABLE]

where

[TABLE]

is the Wigner function333Our conventions for the normalisation of the Wigner function are described in the Appendix A. of the internal deuteron wave function $\varphi_{d}$.

Using a Gaussian as ansatz for the deuteron wave function,

[TABLE]

its Wigner function follows as

[TABLE]

The measured deuteron rms charge radius $r_{\mathrm{rms}}=1.96$ fm Zhaba:2017syr is reproduced choosing444Note that the variable $r$ in the deuteron wave function describes the diameter, such that $r_{\mathrm{rms}}^{2}=\int\differential[3]r(r/2)^{2}|\varphi(r)|^{2}$. $d=3.2\;\mathrm{fm}$. To proceed, we have to choose also an ansatz for the Wigner function of the two-nucleon state: While Ref. Scheibl:1998tk choses a thermal equilibrium state motivated by the picture of a “fireball” formed in heavy-ion collisions, we are in this work interested in the scattering of “small” systems as in $e^{+}e^{-}$, dark matter or $pp$ collisions. Therefore we use that Monte Carlo simulations of strong interactions provide the momentum distribution of the produced nucleons, $G_{np}(\boldsymbol{p}_{n},\boldsymbol{p}_{p})$, which includes also relevant momentum correlations. On the other hand, $G_{np}(\boldsymbol{p}_{n},\boldsymbol{p}_{p})$ can be obtained from the Wigner function as

[TABLE]

where $\psi_{np}(\boldsymbol{p}_{n},\boldsymbol{p}_{p})$ is the normalized two-nucleon wave function in momentum space. We assume therefore a factorization of the momentum and coordinate dependences,

[TABLE]

Note that this assumption implies a transition from a full quantum mechanical treatment to a semi-classical picture. Finally, we neglect spatial correlations between the proton and the neutron, $H_{np}(\boldsymbol{r}_{n},\boldsymbol{r}_{p})=h(\boldsymbol{r}_{n})\,h(\boldsymbol{r}_{p})$ and choose a Gaussian ansatz for $h(\boldsymbol{r})$,

[TABLE]

Eq. (8) then takes the form

[TABLE]

where the factor

[TABLE]

depends on the characteristic spatial spread of the nucleons and on the spatial extension of the deuteron wave function. The coalescence probability is also suppressed for large $q^{2}d^{2}$ as a Gaussian in our model.

2.2 Parameter estimation

In order to estimate the characteristic values for the parameter $\sigma$ in the spatial distribution $h(\boldsymbol{r})$ defined in Eq. (14), one generally has to consider separately the longitudinal and transverse directions,

[TABLE]

Let us discuss first the case of $e^{+}e^{-}$ annihilation into hadrons in the center-of-mass frame of the collision, choosing the $z$-axis along the direction of the outgoing quark and antiquark, as shown in Fig. 1. Before we proceed, it is worth remarking that this reaction involves three different time and distance scales Dokshitzer:1991wu : The annihilation of the electron-positron pair into the quark and antiquark happens during the time $t_{\rm ann}\sim 1/\sqrt{s}$. For $s\gg\Lambda_{\rm QCD}^{2}$, the hard process is thus almost point-like in coordinate space. The perturbative cascading of the produced (anti-) quark proceeds via parton branchings with the characteristic momentum transfer squared $\Lambda_{\rm QCD}^{2}\ll|q^{2}|\ll s$. This implies that the corresponding longitudinal proper distance scales are smaller than $\Lambda_{\rm QCD}^{-1}$. Therefore, the third and last step, the nonperturbative conversion of the final partons into hadrons, corresponds to the longest time and distance scales: The so-called hadronisation time or formation length $L_{\rm had}$ required for a hadron to build up its parton “coat” is

[TABLE]

where $\gamma$ is the gamma factor of the hadron in the considered frame and $L_{0}$ equals approximately the nucleon size, $L_{0}\sim R_{p}\sim 1$ fm. The coalescence process involves nucleons which have (almost) completed their formation and the process proceeds on distance scales which are comparable to $L_{\rm had}$. Boosting to the rest frame of the produced deuteron compensates the gamma factor in Eq. (18), hence we expect $\sigma_{\parallel}\sim L_{0}\sim 1$ fm in that frame.

The characteristic transverse spread of a produced hadron can be estimated using the uncertainty relation: The transverse displacement of the hadron is obtained summing over the perpendicular components of the random walk performed by the previous generations of partons during both the perturbative and nonperturbative parton cascading. The contribution of a single branching is inversely proportional to the transverse momentum of the parton, $\Delta b_{i}\sim 1/p_{\perp,i}$. Also here, nonperturbative physics gives the dominating contribution with555The numerical value of $\Lambda_{\rm QCD}$ depends on the renormalisation scheme used and varies between 0.3 and 0.9 GeV for three flavors Deur:2016tte . $p_{\perp}\sim\Lambda_{\rm QCD}$. Since the bulk of deuterons is produced with relatively small transverse momenta, boosting to the rest frame has a small effect on $\sigma_{\perp}\sim\Lambda_{\rm QCD}^{-1}$. In the simplest option we consider, we neglect therefore this boost. Since $\Lambda_{\rm QCD}^{-1}$ is of the same order of magnitude as $L_{0}$, we set in the following $\sigma_{\parallel}=\sigma_{\perp}=\sigma$, to minimize the number of parameters. An alternative set-up which takes into account the effect of the transverse boost on $\sigma_{\perp}$, will be discussed at the end of this subsection.

Let us now move to proton-proton, proton-nucleus and nucleus-nucleus collisions. Here, the picture is modified by multiple scattering processes involving multiparton interactions: The proton and neutron taking part in the coalescence process can thus originate from different parton-parton interactions. Therefore it is necessary to take into account the longitudinal and transverse spread of the initial parton clouds of the projectile and the target. Starting with the former, it is important to keep in mind that the effect of the Lorentz contraction is different for fast and slow partons. At a given rapidity $y$ in the laboratory system, partons from, say, the target proton cloud are distributed over the longitudinal distance $\sim R_{p}/\gamma=R_{p}\,e^{-y}$. Boosting to the deuteron rest-frame compensates again the gamma factor, such that the resulting “geometrical” contribution to $\sigma_{\parallel}$, due to the longitudinal extension of the parton cloud, equals $\sigma_{\parallel{\rm(geom)}}\sim R_{p}\sim 1$ fm. Summing the two contributions in quadrature, we obtain

[TABLE]

Finally, we have to consider the geometrical contribution to $\sigma_{\perp}$. One may naively expect it to depend on the impact parameter for a proton-proton (proton-nucleus) collision. Let us show that this is not the case in the simple geometrical picture of Fig. 2 and derive the geometrical contribution to $\sigma_{\perp}$. We define $\sigma_{\perp(\mathrm{geom})}$ as the transverse spread of the overlapping region ($O$) of the projectile and target parton clouds,

[TABLE]

The expectation value $\langle A\rangle_{O}$ follows then as

[TABLE]

where $\boldsymbol{b}$ is the impact parameter for the collision, $\rho_{i}(r_{i})$ are the transverse parton densities of the projectile ($i=1$) and the target ($i=2$), and $w_{\mathrm{int}}$ is the probability for a parton-parton interaction. Assuming for simplicity that the latter is point-like,

[TABLE]

and approximating the density distributions by Gaussians,

[TABLE]

with $R_{i}$ being the transverse radii of the projectile and the target, respectively, we obtain

[TABLE]

For the particular case of $pp$ collisions, we have $\sigma_{\perp\mathrm{(geom)}}^{2}=R_{p}^{2}/2$. Since this is of the same order of magnitude as $\sigma_{\perp(e^{\pm})}^{2}$, we set

[TABLE]

such that we can also use for proton-proton collisions one universal parameter,

[TABLE]

It is noteworthy that such an assumption would generally be unjustified in the case of proton-nucleus and nucleus-nucleus collisions since the corresponding geometrical contributions to $\sigma_{\parallel}$ and $\sigma_{\perp}$ may differ significantly. Considering, as an example, proton-lead collisions, we have666Note that the geometrical contributions to $\sigma_{\parallel}$ and $\sigma_{\perp}$ should generally vary from event to event, depending on the corresponding rate of multiple scatterings. For proton-proton collisions, such a variation is expected to be relatively weak and, hence may be neglected in a first approximation. This is, however, different for proton-nucleus and nucleus-nucleus interactions: Our estimations are valid for relatively central collisions involving numerous pair-wise inelastic rescatterings between the projectile and target nucleons, which provide the bulk contribution to the formation of (anti-)nuclei. On the other hand, peripheral interactions at large impact parameters are dominated by a single binary collision between a pair of projectile and target nucleons, which gives rise to $\sigma_{\parallel{\rm(geom)}}\sim R_{p}$, like in the proton-proton case. This leads to a pronounced correlation between the size of the source region and the multiplicity of secondary hadrons produced. $\sigma_{\parallel{\rm(geom)}}\sim R_{\rm Pb}$, while Eq. (24) yields for the transverse spread $\sigma_{\perp(\mathrm{geom})}\simeq R_{p}$.

Boosted $\sigma_{\perp}$

In an alternative set-up, we take into account that $\sigma_{\perp}$ is defined in the collider frame, while in the derivation of Eq. (15), all quantities and wave functions were evaluated in the rest frame of the deuteron. This requires that we Lorentz transform $W(\boldsymbol{r},\boldsymbol{q})\propto h(\boldsymbol{r})$ between the two frames. Such a transformation includes both a longitudinal boost, with a Lorentz factor $\gamma_{\parallel}\simeq\gamma\cos\theta\simeq\gamma$, and a transverse one. Here $\theta\simeq(p_{p\perp}+p_{n\perp})/(p_{p\parallel}+p_{n\parallel})$ is the small angle between the direction of motion of the nucleon pair in the CoM and the $z$-axis in the original frame, before the boost. While the former transformation has been accounted for in our definition for $\sigma_{\parallel}$, the effect of the latter is to replace $\sigma_{\perp}$ defined in the original frame by $\tilde{\sigma}_{\perp}$ in the CoM, with

[TABLE]

Thus, the factor $\zeta$ in Eq. (16) changes to

[TABLE]

2.3 Numerical implementation

As one can see in Eq. (15), a given proton-neutron pair with momentum difference $2\boldsymbol{q}$ in its CoM has the probability

[TABLE]

to form a deuteron. Depending on whether we use the simplified approach or take into account the modification of $\sigma_{\perp}$ by the transverse boost, the factor $\zeta$ is defined by Eq. (16) or (28), respectively.

At this point, we have to take some care of potential double (triple, etc.) counting since a given proton may be paired with different neutrons and vice versa. Let us assume that for a given event, the final state contains $N_{p}$ protons and $N_{n}$ neutrons. Denoting by $w_{ij}$ the coalescence probability, Eq. (29), for a pair formed out of the $i$-th proton and the $j$-th neutron, we have the following expression for the average number of deuterons produced in such an event,

[TABLE]

where in the last line we have taken into account the smallness of the coalescence probabilities and have neglected the contributions of triple and higher contributions.

As one can see from Eq. (30), the contribution of a given proton-neutron pair  $ij$ to the binning of the deuteron spectrum should be taken with the weight,

[TABLE]

Since we bin the deuteron distribution ${\rm d}^{3}N_{d}/{\rm d}P_{d}^{3}$ in the reference frame of the detector, no additional factor accounting for the Lorentz transformation of the yield is necessary.

2.4 Improving the deuteron wave function

In the treatment above, a Gaussian which reproduces the measured $r_{\mathrm{rms}}$ value of the deuteron charge distribution was used as wave function. However, it is known that the deuteron wave function is stronger peaked at $r=0$ than a Gaussian. An alternative is the Hulthen wave function,

[TABLE]

with $a=0.23\,\mathrm{fm}^{-1}$ and $b=1.61\,\mathrm{fm}^{-1}$, which gives a good description of the deuteron ground state Zhaba:2017syr . Using this wave function, an analytical derivation of the weights would however be not possible. To obtain a better description of the deuteron, and at the same time to keep the problem analytical solvable, we use instead the sum of two Gaussians as an ansatz for the deuteron wave function,

[TABLE]

where we include a relative phase $\alpha$ between the two terms. Choosing ${\rm e}^{{\rm i}\alpha}={\rm i}$ leads to some simplifications. First, the probability distribution

[TABLE]

contains with this choice no mixed terms. Moreover, we will see below that this choice leads to the same weight function as in the one-Gaussian case.

Next we fit $|\varphi_{d}(\boldsymbol{r})|^{2}$ to the Hulthen wave function (32) in order to fix $\Delta$, $d_{1}$, and $d_{2}$. Two possible methods are to fit either $|\varphi_{d}(0)|^{2}$, $\langle r\rangle$, and $\langle r^{2}\rangle$, or $\langle r\rangle$, $\langle r^{2}\rangle$, and $\langle r^{3}\rangle$. The first method will be called $\varphi_{0}$-fit and the second $r^{3}$-fit. The $\varphi_{0}$-fit yields $\Delta=0.581$, $d_{1}=3.979$ fm, and $d_{2}=0.890$ fm, while the $r^{3}$-fit yields $\Delta=0.247$, $d_{1}=5.343$ fm, and $d_{2}=1.810$ fm. The resulting probability distributions are plotted in Fig. 3 together with the one for the one-Gaussian (Eq. (10)) and the Hulthen (Eq. (32)) wave functions. One can see in the figure that the two-Gaussian ansatz resembles the Hulthen probability distribution more closely than the Gaussian wave function does, in particular, regarding the peak around $r=0$. The $\varphi_{0}$-fit reproduces visually the behaviour around $r=0$ best and will therefore be used when comparing with experimental data later on.

The deuteron Wigner function follows then as

[TABLE]

where the function $A$ is odd in $\boldsymbol{r}$, $A(-\boldsymbol{r}\cdot\boldsymbol{q})=-A(\boldsymbol{r}\cdot\boldsymbol{q})$. Therefore the new term $A$ drops out performing the spatial integrals over $\boldsymbol{r}_{p}$ and $\boldsymbol{r}_{n}$ in Eq. (8) and, thus, does not contribute to the weights in the binning procedure. The weights for the two-Gaussian case are thus

[TABLE]

where the $\zeta_{i}$ are given by Eq. (16).

The weighted $q^{2}$ distributions with the one-Gaussian weight and two-Gaussian weights for $pp$ collisions at $\sqrt{s}=7$ TeV are shown in Fig. 4 for the ALICE setup discussed in appendix B.1, using $\sigma=7$ GeV*-1* and a constant $\zeta$. For the old model, $p_{0}=0.2$ GeV is used and the resulting distribution is rescaled by a factor 0.3 to make the figure clearer. Double counting has in all cases only a minor effect on the resulting distributions. In the two-Gaussian case, the better description of the peak at $r=0$ in the probability distribution significantly enhances the contribution from proton-neutron pairs with a relatively large momentum difference.

3 Formation of helium-3 and tritium

The cases of helium-3 and tritium nuclei are similar to the deuteron case, but the derivation of the weight is more cumbersome. We account for the Coulomb interaction between the two protons in the helium nucleus only insofar as we allow for a different rms radius of the two nuclei in the fitting procedure. In this approach, our model applies in the same way for helium-3 and tritium. Moreover, this assumption is supported by the data of the ALICE experiment which found a comparable yield of helium-3 and tritium nuclei. The binding energies of both nuclei are still low ($\simeq 8$ MeV), and the same approximations as in the deuteron case thus still apply.

The number of helium nuclei with momentum $\boldsymbol{P}_{\mathrm{He}}$ is found by projecting the helium density matrix onto the three-nucleon one, cf. with Eq. (3). As in the deuteron case, the nucleus wave function is factorised into a plane wave describing the CoM motion with momentum $\boldsymbol{P}_{\mathrm{He}}$ and an internal wave function which depends on the relative coordinates,

[TABLE]

Here, the Jacobi coordinates $\boldsymbol{R}$, $\boldsymbol{\rho}$, and $\boldsymbol{\lambda}$ are expressed via $\boldsymbol{x}_{1}$, $\boldsymbol{x}_{2}$ and $\boldsymbol{x}_{3}$ as

[TABLE]

with $\boldsymbol{x}_{1}^{2}+\boldsymbol{x}_{2}^{2}+\boldsymbol{x}_{3}^{2}=3\boldsymbol{R}^{2}+\boldsymbol{\rho}^{2}+\boldsymbol{\lambda}^{2}$, $\boldsymbol{\rho}^{2}+\boldsymbol{\lambda}^{2}=(\boldsymbol{x}_{1}-\boldsymbol{x}_{2})^{2}+(\boldsymbol{x}_{1}-\boldsymbol{x}_{3})^{2}+(\boldsymbol{x}_{3}-\boldsymbol{x}_{2})^{2}$, $\differential[3]{r_{1}}\differential[3]{r_{2}}\differential[3]{r_{3}}=3^{3/2}\differential[3]{R}\differential[3]{\rho}\differential[3]{\lambda}$. The internal wave function is again approximated by a Gaussian in the relative coordinates $\boldsymbol{\rho}$ and $\boldsymbol{\lambda}$,

[TABLE]

with $b$ being the rms radius of the nucleus,

[TABLE]

The ${}^{3}\mathrm{He}$ and ${}^{3}\mathrm{H}$ nuclei have rms radii equal 1.96 fm and 1.76 fm, respectively karshenboim_precision_2008 .

Performing the same steps as in Eqs. (5) and (7) in the deuteron case, we obtain for the momentum spectrum of the produced nuclei

[TABLE]

where $S=1/12$ is the statistical factor accounting for the different isospin and spin states and $W_{N_{1}N_{2}N_{3}}$ is the Wigner function for the three-nucleon state. We approximate again $W_{N_{1}N_{2}N_{3}}$ by a product of momentum and coordinate distributions of the nucleons, neglecting spatial correlations between the latter,

[TABLE]

where $h(\boldsymbol{r})$ is given by Eq. (14).

Expressing further the product $\varphi_{\mathrm{He}}(\boldsymbol{\rho},\boldsymbol{\lambda})^{*}\,\varphi_{\mathrm{He}}(\boldsymbol{\rho}^{\prime},\boldsymbol{\lambda}^{\prime})$ via the Wigner function of the helium nucleus and doing the spatial integrals, changing to the coordinates (38), we finally obtain

[TABLE]

where

[TABLE]

accounts for the overlap of the wave functions and

[TABLE]

is a measure of the relative momentum difference between the nucleons. The procedure for finding the correct Lorentz transformation is similar to the deuteron case. The result is

[TABLE]

where $\tilde{\sigma}_{\perp}$ is again given by (28).

The numerical procedure for treating the formation of tritium and helium-3 nuclei is identical to the one described in section 2.3, apart from the different weight factor; $P^{2}$ is now defined in the CoM frame of the three-particle state. One may argue that it is sufficient to calculate the momentum differences between nucleons, entering Eq. (45), in the rest frames for the corresponding nucleon pairs since, because of the exponential factor in Eq. (43), those practically coincide with the ones defined in the rest frame of the nucleus. This approach was used throughout this work.

4 Comparison with experimental data

The predicted yield of antinuclei depends on the hadronisation scheme Dal:2012my , and a comparison to different experimental data sets should be therefore made using a single MC simulation. In this work, we choose to perform all our simulations of $pp$ and $e^{+}e^{-}$ collisions with Pythia 8.230 Sjostrand:2006za ; Sjostrand:2014zea which describes the antiproton spectra at LHC energies within 20–30% Adam:2015qaa . We set $\Delta\tau=0$, switching thereby decays of long-lived particles in Pythia off777This choice is similar to the approach of Ref. Ibarra:2012cc , as Pythia 8 only stores non-zero lifetimes in the event table that are relevant for displaced vertices in collider experiments.. Thereby we exclude nucleons which are produced mostly outside the source region. In each run, we take into consideration all produced nucleon pairs with $q<0.25$ GeV for the one-Gaussian case and with $q<0.5$ GeV for the two-Gaussian ansatz for the deuteron wave function. There are two obvious methods to generalise the standard per-event coalescence model to helium-3 and tritium: One can require that each of the relative momenta lie within a sphere with radius $p_{0}$ in momentum space, or that the absolute momentum difference for each pair of particles is smaller than $p_{0}$ Carlson:2014ssa ; Cirelli:2014qia . The latter approach was used here.

Let us now test our model with available data on antideuteron production in $e^{+}e^{-}$ annihilation at the $Z$ resonance energy, from the ALEPH and OPAL experiments, and with ALICE data on antideuteron and antihelium production in $pp$ collisions. Details of the experimental setups are described in Appendix B. Both for the parameter $\sigma$ of the new model and for $p_{0}$ of the standard coalescence model, we perform $\chi^{2}$ fits to these data sets. The best-fit values, their 1 $\sigma$ errors and the reduced $\chi^{2}$ of the various fits are given in Table 1 for the case of constant $\sigma_{\perp}$. In turn, the fit results reported in Table 2 take into account the modification of $\sigma_{\perp}$ by transverse boosts, i.e., they have been obtained using $\zeta$ defined in Eq. (28), with $\tilde{\sigma}_{\perp}$ from Eq. (27). We first note that the two-Gaussian cases lead to significantly reduced $\chi^{2}$ values, compared to the one-Gaussian ansatz or to the standard coalescence model. At the same time, they favor larger values for the parameter $\sigma$, which is related to an enhanced contribution from nucleon pairs with relatively large $q^{2}$. Taking into account the modification of $\sigma_{\perp}$ by transverse boosts improves significantly the quality of the fits, as one can see in Table 2. Moreover, in that case, the best-fit values of the parameter $\sigma$ agree well with our estimates in Section 2.2: The ratio of the values of $\sigma$, determined from fits to the ALICE data for $pp$ collisions and to the ALEPH data for $e^{+}e^{-}$ annihilation, agrees well with the expected one, equal $\sqrt{2}$ (c.f. Eq. (26)), and the absolute value, $\sigma_{(e^{+}e^{-})}\simeq 5\,{\rm GeV}^{-1}\simeq 1$ fm, is consistent with its interpretation as the characteristic hadronic length scale ($\sim R_{p}$).

In Fig. 5, we compare the ALICE antideuteron data with the best-fit spectra obtained for the various cases considered. It becomes apparent that the slope of the $p_{\perp}$ distributions is best described using the two-Gaussian ansatz for the deuteron wave function, when the modification of $\sigma_{\perp}$ by transverse boosts is taken into account.

In Fig. 6, we compare the predictions of our model to the ALICE data for light antinuclei, antihelium and a single data point for tritium. Using the best-fit values for $\sigma$ obtained from fitting the antideuteron data in the two-Gaussian case, the experimental data are satisfactorily reproduced. Additionally, we show the best-fit spectra obtained by fitting the antihelium data. Since the difference between the two fits compared to the errorbars is small, we conclude that the two data sets can be consistently described using the same value of $\sigma$. The goodness-of-fit parameter is $\chi^{2}/(N-1)=6.0/2$, when the $\sigma$ determined in the combined antideuteron fit is used for the helium-3 and tritium data. Thus there is a slight tension, and it will be interesting to check it by future antihelium data with reduced errors.

5 Conclusion

We have developed a new coalescence model for the formation of antinuclei, which combines an event-by-event Monte Carlo description of the collision process with a microscopic coalescence treatment based on the Wigner function representations of the produced antinuclei states. This approach has allowed us to include in a semi-classical picture both the size $\sigma$ of the formation region and the momentum correlations of the nucleons forming the nuclei. Since the size $\sigma$ is process dependent, the difference in the observed antideuterons yields in $e^{+}e^{-}$ and $pp$ collisions can be naturally accounted for. Fitting the single, universal parameter contained in our model to experimental data, we obtained a best-fit value, $\sigma\simeq 1$ fm, which corresponds well to its physical interpretation as the size of the formation region of the light nuclei. If in the future, antideuteron and antihelium data sets with reduced errors and for a larger $p_{\perp}$ range will be available, an independent fit of the two parameters $\sigma_{\perp}$ and $\sigma_{\|}$ might be warranted.

We have examined different approximations for the deuteron wave functions as well as two different implementations of the transverse size $\sigma_{\perp}$ of the formation region. The fits to the ALICE antideuteron data prefer the two-Gaussian wave function and the approach where the effect on $\sigma_{\perp}$, due to a Lorentz boost to the deuteron frame, is taken into account. Both correspond to the physically expected choices: The two-Gaussian wave function takes into account that the deuteron wave function is rather peaked at $r=0$, while $\sigma_{\perp}\simeq R_{p}$ is expected to hold in the CoM frame of the collider. Using the best-fit values for antideuteron, we could describe well the data for the production of antihelium in $pp$ interactions and of antideuterons in $e^{+}e^{-}$ annihilation at the $Z$-resonance energy.

Our model is therefore well suited to investigate in detail the antideuteron and antihelium fluxes predicted in models for dark matter annihilations and from cosmic ray interactions. In particular, it will be interesting to see whether and how the tentative antihelium events announced by the AMS-02 collaboration ting16 can be explained within our model.

Acknowledgements

M.K. and J.T. acknowledge partial support from the Research Council of Norway (NFR). S.O. acknowledges support from project OS 481/2-1 of the Deutsche Forschungsgemeinschaft.

Appendix A Wigner function

Our definition (9) of the one-particle Wigner function implies in $d=1$ as normalisation (with $\hbar=1=h/(2\pi)$)

[TABLE]

The corresponding probability distributions for the space and momentum variables are obtained from

[TABLE]

For our ansatz $W(x,p)=h(x)g(p)$, it follows that $h(x)$ describes the probability distribution of the nucleon in coordinate space, while the probability distribution of the nucleon momenta $g(p)$ is normalised as

[TABLE]

Appendix B Experiments

B.1 ALICE

The ALICE Collaboration measured the invariant differential yields of deuterons and antideuterons,

[TABLE]

in inelastic proton-proton collisions at centre of mass energies $\sqrt{s}=0.9,2.76$ and 7 TeV, in the $p_{T}$ range $0.8<p_{T}<3$ GeV and for rapidity888Note that while an additional pseudo-rapidity cut $|\eta|<0.8$ was required in the data selection, the measurements have been corrected to correspond to the $|y|<0.5$ selection, including also the contribution of $|\eta|>0.8$, using a MC simulation Acharya:2017fvb . $|y|<0.5$ Acharya:2017fvb . Here $E$ and $\boldsymbol{p}$ are the deuteron energy and momentum, $N_{\mathrm{inel}}$ is the number of inelastic events, $N$ is the total number of detected deuterons, and $n\equiv N/N_{\mathrm{inel}}$. The experiment included a trigger (V0) consisting of two hodoscopes of 32 scintillators that covered the pseudo-rapidity ranges $2.8<\eta<5.1$ and $-3.7<\eta<-1.7$, used to select Non-Diffractive (ND) inelastic events. An event was triggered by requiring a hit (charged particle) on either side (positive or negative $\eta$) of the V0 triggering setup.

Pythia 8 generates general inelastic collisions, including single-diffractive (SD), double-diffractive (DD) and ND events. The minimum bias events selected by the V0 trigger generally include those that Pythia treats as SD and DD events. While we used Pythia 8 to generate general minimum bias $pp$ collisions, only events satisfying the V0 trigger have been included in our analysis.

B.2 ALEPH and OPAL

The ALEPH collaboration at LEP studied the deuteron and antideuteron production in $e^{+}e^{-}$ collisions at the $Z$ resonance energy. The measured production rate of antideuterons was $(5.9\pm 1.8\pm 0.5)\times 10^{-6}$ per hadronic $Z$ decay, for the antideuteron momentum range from 0.62 to 1.03 GeV and for the production angle $\theta$ satisfying $|\cos\theta|<0.95$ Schael:2006fd .

In a similar experiment performed by the OPAL collaboration Akers:1995az , no antideuteron events were detected. Reference Dal:2015sha noted that the resulting upper limit on the antideuteron yield has previously been neglected, but should also be taken into account. The measurements were performed in the antideuteron momentum range $0.35<p<1.1$ GeV, with an estimated detection efficiency $\epsilon=0.234$, which includes the angular acceptance. The expected total number of antideuterons was

[TABLE]

where $N_{\mathrm{ev}}=1.64\times 10^{6}$ is the number of events in the OPAL analysis and $n_{\bar{d},\mathrm{MC}}$ is the MC prediction for the number of antideuterons per event. We follow Ref. Dal:2015sha and assume a Poissonian uncertainty $\sigma_{\bar{d}}=\sqrt{N_{\bar{d}}}$ for the expected number of antideuterons. The $\chi^{2}$ is in this case given by

[TABLE]