Antideuteron production in $\Upsilon(nS)$ decays and in $e^+e^- \to q\overline{q}$ at $\sqrt{s} \approx 10.58 \mathrm{\,Ge\kern -0.1em V}$

We present measurements of the inclusive production of antideuterons in $e^+e^-$ annihilation into hadrons at $\approx 10.58 \mathrm{\,Ge\kern -0.1em V}$ center-of-mass energy and in $\Upsilon(1S,2S,3S)$ decays. The results are obtained using data collected by the BABAR detector at the PEP-II electron-positron collider. Assuming a fireball spectral shape for the emitted antideuteron momentum, we find $\mathcal{B}(\Upsilon(1S) \to \bar{d}X) = (2.81 \pm 0.49 \mathrm{(stat)} {}^{+0.20}_{-0.24} \mathrm{(syst)})/! \times /! 10^{-5}$, $\mathcal{B}(\Upsilon(2S) \to \bar{d}X) = (2.64 \pm 0.11 \mathrm{(stat)} {}^{+0.26}_{-0.21} \mathrm{(syst)})/! \times /! 10^{-5}$, $\mathcal{B}(\Upsilon(3S) \to \bar{d}X) = (2.33 \pm 0.15 \mathrm{(stat)} {}^{+0.31}_{-0.28} \mathrm{(syst)})/! \times /! 10^{-5}$, and $\sigma (e^+e^- \to \bar{d}X) = (9.63 \pm 0.41 \mathrm{(stat)} {}^{+1.17}_{-1.01} \mathrm{(syst)}) \mbox{\,fb}$.

We present measurements of the inclusive production of antideuterons in e + e − annihilation into hadrons at ≈ 10.58 GeV center-of-mass energy and in Υ (1S, 2S, 3S) decays. The results are obtained using data collected by the BABAR detector at the PEP-II electron-positron collider. Assuming a fireball spectral shape for the emitted antideuteron momentum, we find B(Υ (1S) →dX) = (2.81 ± 0.49(stat) +0. 20 −0.24 (syst))×10 The production of nuclei and anti-nuclei in hadronic collisions and in hadronization processes has recently attracted considerable theoretical and experimental interest [1-3], since cosmic anti-nuclei may provide a sensitive probe of dark matter annihilation. Dark matter particles might annihilate into two colored partons -quarks (q) and gluons (g) -which could hadronize into mesons and baryons, potentially forming bound states such as light (anti)nuclei. The latter process, requiring at least six q or q in close proximity, is poorly understood both theoretically and experimentally, and precise measurements of both total rates and momentum spectra are needed. With no initial-state hadrons, e + e − annihilations provide a clean probe of this process not only for q and q but also for g via decays of Υ and other vector resonances.
Experimental measurements focus on antideuteron (d) production as such studies are not limited, as in the deuteron (d) production case, by the high nuclei background due to interactions with the detector materials. The ARGUS [4] and CLEO [5] experiments observed d production at the level of 3 × 10 −5 per Υ (1S) and Υ (2S) decays, and set limits on production in Υ (4S) decays and e + e − → qq at 10.6 GeV. The ALEPH [6] experiment observed a 3σ evidence for d production in e + e − → qq at 91.2 GeV. In these measurements, the accessible kinematic range, 0.4 -1.7 GeV/c, was representing less than 20% of the phase space where d and d could be identified.
In this Letter, we present studies of d production in e + e − annihilation data taken on and just below the Υ (2S), Υ (3S) and Υ (4S) resonances. We also study the d production from Υ (1S) using the Υ (2S) → Υ (1S)π + π − decay chain. The boost of the CM allows a wide momentum to be accessed, 0.3 -3 GeV/c, which corresponds to 0.5 -1.5 GeV/c in the laboratory frame. We confirm the Υ (1S) rates, improve the coverage and precision for Υ (2S) decays, measure d production in Υ (3S) decays, and also, for the first time, in e + e − → qq near 10.6 GeV.
The results presented here are obtained from the complete BABAR Υ (2S, 3S, 4S) (Υ (nS)) datasets (Onpeak), including data collected at a CM energy 40 MeV below the peak of each resonance (Offpeak). The luminosity [7] collected for each dataset and the corresponding number of Υ (nS) decays are reported in Table I. We also use Monte Carlo (MC) simulated data samples generated using JetSet [8] for e + e − → qq (q = u, d, s) events and EvtGen [9] for Υ (nS) decays. The interaction of simulated particles with the BABAR detector is modeled using Geant4 [10]. Neither d nor d production is implemented in JetSet, and d cannot be simulated in the version of Geant4 that is used. Therefore we use the EvtGen phase-space generator and Geant4 to simulate Υ (2S, 3S) → dNN (5h) decays (where h indicates a K ± , π ± , or π 0 , and N, N = p, n) and e + e − → qq → dNN (5h) events for studying reconstruction efficiencies. The additional five hadrons are a representative average of additional particles in the decay and restrict the phasespace, so we have d's with CM momentum lower than 3 GeV/c. These samples will be referred to as "signal MC" throughout this Letter. To account for differences between data and MC samples, corrections are applied and systematic uncertainties are assigned, as discussed below.
The BABAR detector, trigger, and the coordinate system used throughout, are described in detail in Refs. [11,12]. The most relevant part of the detector for this analysis is the tracking system, composed of a 5-layer inner silicon strip tracker, the Silicon Vertex Tracker (SVT), and the 40-layer small-cell Drift Chamber (DCH) inside a 1.5 Tesla axial magnetic field. The SVT provides information on track parameters near the interaction point (IP), while the DCH has a 98% efficiency for detecting charged particles with p T > 500 MeV/c. The p T resolution is σ p T /p T = (0.13 (GeV) −1 · p T + 0.45)%. The ionisation energy loss (dE/dx) is measured by the two systems, with a resolution of approximately 14% and 7%, for the SVT and DCH, respectively. Additional particle identification information is provided by a Detector of Internally-Reflected Cherenkov light (DIRC), which, as described later, is employed in this analysis to provide a veto.
Hadronic events are selected by a filter which requires greater than two reconstructed tracks and a ratio of second to zeroth Fox-Wolfram moments [13] less than 0.98. The reconstructed momentum of the candidate tracks is corrected by 0.019( GeV/c) 3 /p 2 to account for the underestimation of energy loss due to the pion mass assumption employed in the track fit. Candidates are retained only if they are within the full polar angle acceptance of the DCH (−0.80 ≤ cos θ LAB ≤ 0.92) and within 0.5 ≤ p LAB ≤ 1.5 GeV/c, where the most probable dE/dx of d and d is well separated from that of other particle species. Here and throughout this Letter "LAB" denotes observables in the laboratory frame. To reject candidates with poorly-measured dE/dx, we require that the ionization along the track trajectory be sampled at least 24 times by the DCH. A relevant background contribution to the observed d signal comes from "secondary" d's produced in nuclear interactions with the detector material that travel inward toward the IP and are wrongly reconstructed as outward-travelling d's. To suppress this contribution we require that the transverse distance of closest approach (DOCA) of the reconstructed trajectory to the beamspot be less than 400 µm. The effect of underestimated energy loss on the measured DOCA of tracks is found to be well reproduced in the simulation, and we do not apply any correction to this quantity. Finally, d's and d's in the considered momentum range are below the threshold for radiating Cherenkov light in the DIRC quartz-glass bars, so we reject all the candidates with more than 10 associated Cherenkov photons for the the best-fit DIRC mass hypothesis, π, K, p, e or µ.
To measure the d yields, we apply a weight to each candidate to correct for detector and selection acceptance, and we then extract the yields from the d candidate energy loss distributions of the DCH and SVT using a weighted fit. Global trigger and event selection efficiencies are determined from simulated e + e − → Y → 2(N N )X events, where X corresponds to zero or more additional final-state particles, and Y = Υ (2S), π + π − Υ (1S), or qq in which the nucleons are produced promptly in fragmentation. Those efficiencies are assumed to be the same for Υ (2S), Υ (3S), and Υ (4S), so only the first is explicitly calculated and used also for the other resonances. Corrections for the kinematic selections are computed as a function of p CM using the fraction of d, in bins of CM momentum, which would pass the selection in the LAB frame. This fraction depends on the angular distribution of d's in the CM frame with respect to the beam axis, which we determine from MC generator coalescence studies with coalescence momentum p 0 = 160 MeV/c following a similar approach to [1]. Decays of Υ → ggg produce d's and d's isotropically, while in e + e − → qq there is a dependence on the CM polar angle.
The d reconstruction efficiency is determined in bins of p LAB and cos θ LAB from the "signal" MC samples. We compute an additional correction to this efficiency to account for the differing interactions of d and d in the material of the BABAR detector. Starting from the differing reconstruction probabilities of protons and antiprotons, as determined by Geant4 simulation of Υ (nS) decays, we estimate the effect of material interaction on d and d by rescaling for the larger d absorption cross sections, determined in Ref. [14]. As a cross-check on the result, the values obtained at cos θ LAB = 0 are found to be consistent with a prediction obtained from d inelastic cross sections from Ref. [15] and the known distribution of material in the BABAR detector. The final weight applied to each track is the inverse of the product of the trigger and event selection efficiencies, the kinematic acceptance fraction, and the d reconstruction efficiency.
Yields are computed using a fit to the distribution of the normalized residual of ionization energy loss. Specific ionization measurements from the DCH and SVT, their uncertainties, and their expected values from the Bethe-Bloch formula are calibrated using high-statistics control samples of particles of other species. The independent dE/dx measurements, in arbitrary units, from the DCH and SVT are averaged according to their respective uncertainties, after rescaling the SVT measurement and its uncertainty such that the expected value matches that of the DCH. The normalized residual is computed as the difference between the averaged dE/dx measurement and the expected value, divided by the uncertainty of the former. Ideally, this residual (shown in Fig. 1 for Υ (2S) data) has a Gaussian distribution centered at zero for d's and d's, while the value for other species is far from zero (the value is positive or negative for particles with higher or lower mass, respectively).
The probability density function (p.d.f.) for d's and d's is estimated using signal MC events, while the distribution for all other particles is taken from the distribution for negatively-charged particles for the simulated generic decays of Υ (nS), which contain no true d. The signal distribution is found to be well-modeled by a piecewise combination of a Gaussian function with an exponential tail toward lower values. After requiring that the piecewise function and its first derivative be continuuous, the addition of the tail adds only a single parameter as compared to a pure Gaussian. The residual distribution for other particles ("background") is more complicated, and extends into the signal region. This background distribution, obtained from simulated Υ (2S) events, is described well by the sum of a Gaussian function and an exponential function. Only the functional forms of the shapes used are extracted and validated using the MC samples, while almost all the p.d.f. parameters are estimated in the fit to the data. Very few candidates with large weights are removed since they could have an undue influence on the weighted fit [16]. An example of the fit to the residuals distribution for the d in the Υ (2S) data is shown in find a contribution from triton produced in material interactions. The tritons' distribution is similar to that of the d's and d's, so it is modeled using the same distribution with parameters allowed to float separately.
Candidates are divided into categories according to their charge, their CM momentum, and the type of dataset (Onpeak or Offpeak), and a weighted unbinned maximum likelihood fit is performed to all categories simultaneously. The d, d, triton and background distri-bution functions are the same for all categories, but the yields are floated separately. We divide the CM momentum range [0.35,2.25] GeV/c into nine bins containing approximately equal numbers of candidates, with no bin narrower than 100 MeV/c. To improve the quality of the fit the width of the Gaussian function for the background is also allowed to float separately for different bins of CM momentum, to account for significant difference in the distribution for different energy ranges. To achieve a more stable fit, if the fit results for a split parameter (i.e., one allowed to take different values in different sub-samples) are statistically compatible between two or more sub-samples, the parameter is forced to have the same floating value among those sub-samples. We perform a simultaneous fit to Onpeak and Offpeak datasets for Υ (2S) and Υ (3S), obtaining the number of d in each CM momentum bin. Yield values and their uncertainties in each bin are reported in Ref. [17]. We determine the final number of d by subtracting the yields for d in the Offpeak dataset from the Onpeak, after rescaling for the luminosity and the lowest-order 1/s correction of the e + e − → qq cross section.
The Υ (4S) decays almost exclusively to BB final states, and d production in B decays is kinematically disfavored, so the production from Υ (4S) decays is well below our sensitivity. Therefore, we proceed by combining the yields from Onpeak and Offpeak datasets to obtain the production rate for e + e − → qq. As a crosscheck, when subtracting the rescaled yields in the Offpeak dataset from the yields in the Onpeak one, the number of d's are compatible with zero in all bins.
Corrections due to the requirement on the number of DCH dE/dx samplings can be computed by comparing the distribution for this variable in data and signal simulation for d's in a narrow signal window that provides a clean and high-statistics control sample. We correct the fitted yields by −7% to account for this effect, which is not dependent on the CM momentum or polar angle.
To validate the fit procedure and check for possible biases in the d yields, we perform a series of fits to pseudodatasets generated according to the fit p.d.f.'s and we assign a systematic uncertainty based on any bias found (i.e., the deviation from zero of the mean of the normalized residuals distribution). The uncertainty from our choice of the background model distribution function is estimated by comparing to results obtained using a background model consisting of two Gaussian functions fixed to a common mean. Additional systematic uncertainties result from the weights used to correct for reconstruction efficiency in the detector and kinematic acceptance. We evaluate the effect of finite Monte Carlo statistics by allowing the weights to vary according to a Gaussian distribution centered at the nominal value with the width given by the uncertainties in the weights. We assign uncertainties in each bin from the distribution of fit results. The uncertainty in the correction for d material interaction is computed similarly, with the width of the Gaussian distribution set to the statistical uncertainty in the calculation plus the uncertainty assuming the correction distributed uniformly.
We estimate the contribution from secondary d's, described above, by fitting the DOCA distributions for data and simulated events in and outside the selected region. The simulation describes the data except for a slowly falling exponential component, which we ascribe to secondary d's. The fits indicate a small contribution to the selected events, which we take as a one-sided systematic uncertainty. We compare DOCA distributions of simulated antiprotons with a control sample of wellidentified antiprotons and we assign a +5.8% one-sided uncertainty. The uncertainty in the event selection efficiency is estimated by comparing the nominal efficiency with the efficiency for events for which at least one of the generated prompt (anti)nucleons is a (anti)proton within the detector angular acceptance. Finally, uncertainties in the cross sections and number of Υ (nS) mesons are propagated along with the other systematics. The contributions to the systematic uncertainties are summed in quadrature separately for the positive and negative sides, and their values for each contribution and the totals are summarized in Table II. Systematic values in each bin are reported in Ref. [17].
The numbers of d's extracted from the fit are corrected for the differences between data and MC samples mentioned above, and then are converted into branching fractions for Υ (2S) and Υ (3S) using the total number of Υ decays. Using the total luminosity from Onpeak and Offpeak Υ (4S) datasets, we compute the cross section for d production from e + e − → qq at a CM energy of ≈ 10.58 GeV.
The number of Υ (1S) decays is computed as where N fit Υ (1S) is the number of Υ (2S) → ππΥ (1S) events reconstructed inclusively, obtained by a fit to the invariant mass recoiling against a reconstructed π + π − system in Υ (2S) data. The quantities f sig and f sb are respectively the fraction of the fitted Υ (1S) recoil mass distribution in the signal and sideband regions used to subtract the contribution from background events, and ε filter is the average efficiency of the trigger and the event filter to accept Υ (2S) → ππΥ (1S) decays obtained from the Υ (2S) MC sample. After applying these factors, we find N Υ (1S) = (9.670 ± 0.023) × 10 6 . The final values for the differential branching fractions and e + e − → qq cross section are shown in Fig. 2. The total rates, presented in Table III, are obtained from the measured differential spectra by fits to the "fireball" model distribution [18] where E is the d CM energy and α and β are free parameters determined by the fit. The fits are shown in   Table III are the integral of these distributions from 0 GeV/c to 4 GeV/c, and the associated uncertainties are those in the integral taking into account the full covariance of α and β. We find that values for the parameter β in Υ decays are mutually compatible within 1σ, and the average value is β = (4.71 ± 0.19) GeV −1 . Fitting to the e + e − → qq spectrum yields a lower value of (3.92±0.22) GeV −1 , corresponding to a somewhat harder spectrum. As an additional cross-check on the cross section for e + e − → qq production, we fit separately the differential spectra from the Υ (4S) Offpeak dataset only, and obtain a cross section of 12.2 ± 1.8 fb, where the error is statistical only. Values for both this cross section and other parameters of the fit are in good agreement with the Onpeak plus Offpeak result.
In summary, we have performed measurements of inclusive d production in Υ (1, 2, 3S) decays and in e + e − → qq. These are the first measurements of d production in e + e − → qq at a CM energy of ≈ 10.58 GeV and in Υ (3S) decay, and the most precise measurement in Υ (2S) decay. Our total and differential rates for inclusived production in Υ (1S) and Υ (2S) decay are in good agreement with previous measurements [15], and with the expected spectral shapes from the coalescence model [20,21]. We additionally note an order of magnitude suppression of d production in quark-dominated e + e − → qq relative to the gluon-dominated Υ decays.
We are grateful for the excellent luminosity and machine conditions provided by our PEP-II colleagues, and for the substantial dedicated effort from the computing organizations that support BABAR. The collaborating institutions wish to thank SLAC for its support and kind hospitality. This work is supported by DOE