Observation of $D^+_s\rightarrow p\bar{n}$ and confirmation of its large branching fraction

The baryonic decay $D^+_s\rightarrow p\bar{n}$ is observed, and the corresponding branching fraction is measured to be $(1.21\pm0.10\pm0.05)\times10^{-3}$, where the first uncertainty is statistical and second systematic. The data sample used in this analysis was collected with the BESIII detector operating at the BEPCII $e^+e^-$ double-ring collider with a center-of-mass energy of 4.178~GeV and an integrated luminosity of 3.19~fb$^{-1}$. The result confirms the previous measurement by the CLEO Collaboration and is of greatly improved precision, which may deepen our understanding of the dynamical enhancement of the W-annihilation topology in the charmed meson decays.

The baryonic decay D + s → pn is observed , and the corresponding branching fraction is measured to be (1.21 ± 0.10 ± 0.05) × 10 −3 , where the first uncertainty is statistical and second systematic. The data sample used in this analysis was collected with the BESIII detector operating at the BEPCII e + e − double-ring collider with a center-of-mass energy of 4.178 GeV and an integrated luminosity of 3.19 fb −1 . The result confirms the previous measurement by the CLEO Collaboration and is of greatly improved precision. This result will improve our understanding of the dynamical enhancement of the W-annihilation topology in the charmed meson decays. The decay D + s → pn is the only kinematically allowed baryonic decay of the three ground-state charmed mesons D 0 , D + , and D + s . It provides a unique probe of hadronic dynamics and is of great importance to the study of weak annihilation decays of charmed mesons [1][2][3][4][5]. At the short-distance level, under the vacuum-insertion approximation, its branching fraction (BF) is predicted to be very small (of the order 10 −6 ) owing to chiral suppression by the factor (m π /m Ds ) 4 which follows from the partially conserved axial current [4]. This physically corresponds to the mechanism of helicity suppression.
The CLEO Collaboration reported evidence for the decay D + s → pn with 13.0 ± 3.6 signal events, resulting in an anomalously large BF of B(D + s → pn) = (1.30 ± 0.36 +0. 12 −0.16 ) × 10 −3 [6]. This unexpectedly large BF stimulates the interest of theorists. Many phenomenological possibilities have been proposed to explain the apparent discrepancy between theoretical predictions and the experimental measurement, e.g., the not well justified factorization ansatz due to the light mass of charm quark and the complicated final state interaction at the threshold of pn production [2][3][4], a contribution of additional decay mechanisms such as final state scattering [4], or the effect of the time-like baryonic form factors from the axial vector currents [7]. Experimentally, the confirmation of the observation of the decay D + s → pn by different experiments is highly desirable, and a much improved precision on its decay BF is necessary to distinguish between different phenomenological models and understand the decay dynamics of charmed mesons. The e + e − annihilation sample collected at √ s = 4.178 GeV with the Beijing Spectrometer (BESIII) in 2016, which corresponds to an integrated luminosity of 3.19 fb −1 and is 5 times larger in statistics compared to the CLEO data, provides a good opportunity for this measurement.
BESIII is a general-purpose detector with 93% coverage of the full solid angle. Details of the detector can be found in Ref. [8]. In 2015, BESIII was upgraded by replacing the two endcap time-of-flight (TOF) systems with a new detectors that use multi-gap resistive plate chambers (MRPC), which achieve a time resolution of 60 ps [9].
A geant4-based [10] Monte Carlo (MC) simulation software package, which includes the description of the BESIII detector geometry and its response, is used to generate MC simulated event samples. The simulation includes the beam energy spread and initial state radiation (ISR) in the e + e − annihilations modeled with the generator ConExc [11]. The final state radiation (FSR) from charged tracks is incorporated with the photos package [12]. The generic MC samples, consisting of the production of open charm processes, the ISR return to low-mass charmonium (ψ) states, and continuum processes (quantum electrodynamics processes and continuum production of light quarks qq, q = u, d, s), have a size corresponding to an integrated luminosity 35 times larger than that of the data. The known particle decays are generated using evtgen [13] with the BFs taken from the Particle Data Group (PDG) [14], and the remaining unknown decays of low mass ψ states are generated with lundcharm [15]. We also generate a signal MC sample of 4 × 10 6 e + e − → D * + s D − s → D + s γD − s events, in which D + s decays to pn, while the D − s is set to decay generically, while the charge conjugated modes is also included; this sample is used to obtain the shapes of kinematic variables in signal decays and to estimate systematic uncertainties. Throughout the Letter, charge conjugated modes are implicitly implied, unless otherwise noted.
In this analysis, the D + s sample is predominantly produced in the reaction e + e − → D * + s D − s → D + s γD − s . We fully reconstruct a D − s meson, named "single tag (ST)", in eleven decay modes that correspond to 25% of the total decay width [14]: s sample, we further require an isolated photon consistent with D * + s decay and reconstruct the D + s → pn signal in the side recoiling against the D − s candidate, referred to as the "double tag (DT)". Both the D + s directly produced in the e + e − annihilation and the one from D * + s decay are considered. Thus, the numbers of ST (N i ST ) and DT (N i DT ) candidates for a specific tag mode i are where N tot is the total number of e + e − → D * + s D − s + c.c. events in the data, B i , B D * + s →γD + s , and B D + s →pn are the BFs for D − s tag mode i, D * + s → γD + s , and D + s → pn, respectively, and ǫ i ST(DT) is the ST (DT) detection efficiency. The factor 2 indicates that the signal D + s is either directly produced in the e + e − annihilation or from D * + s decay. Based on Eqs. (1) and (2), combining the eleven ST modes leads to the expression where N tot DT is the total number of DT signal events reconstructed from all ST modes.
All charged tracks are reconstructed from hits in the main drift chamber (MDC) with a polar angle θ (with respect to the beam direction) within | cos θ| < 0.93. Charged tracks, except for those from K 0 S decays, are required to have a point of closest approach to the interaction point (IP) within ±10 cm along the beam direction and within 1 cm in the plane perpendicular to the beam axis. Particle identification (PID) is performed by combining the specific energy loss dE/dx measured in the MDC and the TOF information. A charged π(K) candidate is identified by requiring the PID likelihood value Photon candidates are reconstructed with energy deposits in the electromagnetic calorimeter (EMC) that are not associated with reconstructed charged tracks. The photon is required to have an energy larger than 25 MeV in the barrel region (| cos θ| < 0.8), or 50 MeV in the endcap region (0.86 < | cos θ| < 0.92). To suppress electronic noise and energy deposits unrelated to the events, the shower time in the EMC must be within 700 ns of the event start time [16]. The π 0 and η candidates are reconstructed from γγ pairs with an invariant mass M γγ within (0.115, 0.150) GeV/c 2 and (0.50, 0.57) GeV/c 2 , respectively. Candidates with both photons in the endcap regions are rejected due to the bad energy resolution. To improve the momentum resolution, a 1C kinematic fit is performed, constraining M γγ to the nominal π 0 or η mass [14] and requiring χ 2 < 30. The updated momentum of each photon from the kinematic fit is used in the further analysis.
The K 0 S candidates are reconstructed via the decay K 0 S → π + π − by performing a vertex-constrained fit to all oppositely charged track pairs without PID requirements applied. The charged tracks must be within | cos θ| < 0.93, and have a point of closest approach to the IP within ±20 cm along the beam direction; no requirement is placed on the point of closest approach in the plane perpendicular to the beam. The χ 2 of the vertex fit must be less than 100. To suppress the combinatorial background, a secondary vertex fit is performed, constraining the direction of the K 0 S momentum to point back to the IP, and requiring χ 2 < 20. The flight length L, defined as the distance between the common vertex of the π + π − pair and the IP, is obtained in the secondary vertex fit and required to satisfy L > 2σ L for accepted K 0 S candidates, where σ L is the uncertainty of L. The four-momenta after the secondary vertex fit are used in the subsequent analysis. The K 0 S candidate is required to have a mass within the range (0.487, 0.511) GeV/c 2 , corresponding to three standard deviations on the mass distribution.
In the ST mode D − s → K − π + π − , the π + π − invariant mass is required to be outside the range (0.480, 0.515) GeV/c 2 to avoid double counting with the ST mode D − s → K 0 S K − . For a given ST mode, the D − s candidates are reconstructed by all possible combinations of selected K ± , π ± , K 0 S , π 0 , η and η ′ candidates in an event, and are identified with the corresponding invariant mass M tag . To suppress the background from the non-strangeness excited D * decay D * → πD, the π ±(0) candidates from D + s decays must have a momentum larger than 100 MeV/c. To further suppress the non D * + s D − s backgrounds, a variable that represents the invariant mass of the system recoiling against the selected D − s candidate is defined as where E cm is the center-of-mass energy, p Ds is the momentum of the selected D − s candidate in the center-ofmass system, and M Ds is the nominal D − s mass [14]. In the process e + e − → D * + s D − s → D + s γD − s , the selected D − s candidates are produced either directly in the e + e − annihilation or from the decay D * − s → γD − s . The corresponding M rec distribution for the former case peaks at the nominal D * + s mass M D * + s [14] smeared by the mass resolution, and that for the latter case has a relatively flat distribution between 2.05 and 2.18 GeV/c 2 . The D − s candidates are accepted by requiring 2.05 < M rec < 2.18 GeV/c 2 . The e + e − → D + s D − s process is highly suppressed by this requirement. For an event with multi-D − s candidates for a specific tag mode per charge, only the one with minimum |M rec − M D * + s | is kept. The M tag distributions of the events passing the above selection criteria are shown in Fig. 1 for all ST modes. The ST yields are determined by performing a binned maximum likelihood fit. In the fit, the D − s signal is described by the MC-simulated line-shape convolved with a Gaussian function representing the resolution difference between data and MC simulation, where the parameters of the Gaussian functions are free parameters the fit. The background is described by Chebychev polynomial functions of the first kind of first or second order. The fit results are superimposed on the data in Fig. 1. For further study, we require that M tag is within 2.5 times the resolution around the D − s peak. The requirements on M tag , the ST yields, and the corresponding ST detection efficiencies obtained with the generic MC samples are summarized in Table I for each individual ST mode. The signal D + s → pn and the isolated photon from the D * s decay are reconstructed from the remaining tracks and photons that are not used in the ST D − s reconstruction. Exactly one remaining charged track with opposite charge to the ST D − s meson and at least one remaining good photon are required. The charged track is identified as a proton by requiring L(p) ≥ L(K), L(p) ≥ L(π) and L(p) ≥ 0.001. The angle between this isolated photon and the nearest charged track is required to be larger than 10 • .
To improve the resolution and the likelihood of associating the correct photon candidate from the D * s decay, we perform a kinematic fit with constraints on the masses of the ST D − s , signal D + s , intermediate state D * ± s , and the initial four-momentum. The two hypotheses, i.e., e + e − → D * + s (γ + pn)D − s (ST ) or e + e − → D + s (pn)D * − s (γ + ST ), are tested, and the one with the smaller fit χ 2 is chosen. In the fit, the antineutron is treated as a missing particle with unknown mass, thus there are a total of 7 − 4 = 3 constraints. The χ 2 of the kinematic fit is required to be less than 200. This requirement retains most of the signal events, but removes 50% of background. For an event with more than one remaining photon, we try all photon candidates in the kinematic fit, and the one with the smallest χ 2 is selected. The updated momenta after the kinematic fit are used in the subsequent analysis. The resulting mass of the missing particle M miss , using all ST modes, is shown in Fig. 2. A prominent antineutron signal is visible. The uncertainties are statistical only. The BFs of π 0 /η → γγ, K 0 S → π + π − , η ′ → π + π − η and η ′ → γπ + π − are not included in efficiencies.

ST mode
Mtag The potential backgrounds are classified into (a) non-D − s background and (b) real-D − s background. The background (a) is dominated by continuum processes with proton and antineutron in the final state and can be estimated with the events in the M tag sideband region (3.5∼5.0σ away from the D s peak). The corresponding M miss distribution of background (a) is shown as the shaded histogram in Fig. 2. No obvious peak is observed in the vicinity of the anti-neutron signal. Since D + s → pn is the only baryonic decay mode for the D + s meson, no peaking background is expected for background (b). The properties of the backgrounds are validated by studying the generic MC samples. The total DT signal yield is determined by performing an unbinned maximum likelihood fit to the M miss distribution in Fig. 2, where the signal is described by an MC-simulated line-shape convolved with a Gaussian function representing the resolution difference between data and MC simulation; the background is modeled by an ARGUS function [17]. The fit shown in Fig. 2 returns 193 ± 17 D + s → pn signal events. The DT efficiencies for the individual ST mode are estimated by performing the same procedure on the generic MC samples, and are summarized in Table I. Based on Eq. (3), inserting all the numbers reported above and incorporating the world-average value for B(D * + s → γD + s ) [14], we obtain B(D + s → pn) = (1.21 ± 0.10) × 10 −3 , where the uncertainty is statistical only.
With a DT technique, the systematic uncertainties on detecting the ST D − s meson largely cancel. For the reconstruction of the isolated photon and the signal D + s → pn, the following sources of systematic uncertainties are studied, resulting in a total systematic uncertainty of 4.4% when the individual contributions are summed in quadrature.
The efficiencies for proton tracking and PID are studied as function of cos θ and momentum using the control sample e + e − → π + π − pp. The results are then weighted by the cos θ and momentum distributions of the proton in the signal MC. The average efficiency difference between data and MC simulation combined for tracking and PID is 3.2%, which is taken as the systematic uncertainty.
We study the uncertainties associated with the photon detection and the kinematic fit simultaneously with a control sample of D + s → K 0 S K + decays produced in the process e + e − → D * + s D − s → D + s γD − s . The resultant difference on the efficiencies between data and MC simulation is 2.4%, which is assigned as the systematic uncertainty from this source.
The proton and antineutron may produce additional showers in the EMC that might then affect the efficiency of detecting D + s → pn decays. To estimate this effect, we examine the detection efficiencies determined with two different signal MC samples that are produced with and without the neutron interaction effect in the EMC, respectively. Conservatively, we assign half of the difference between the two efficiencies, 0.9%, as the uncertainty.
The uncertainty sources associated with the fit to the M miss distribution include the background parameterization and the fit range. The corresponding uncertainties are estimated by performing fits with alternative background shape obtained with the events in the ST M tag sideband region and various fit ranges. The resultant changes on the signal yields are regarded as the corresponding uncertainties. The sum of the three uncertainties above in quadrature is 0.7%, which is taken as the associated systematic uncertainty.
For the ST D − s yields, there is a contribution from the process e + e − → γ ISR D + s D − s , which causes a tail falling into the M rec windows. We estimate this background contributes to our ST yields by at most 0.3% based on the MC simulation. We take this upper limit as the systematic uncertainty from this source.
According to Eq. (3), the uncertainty related to the ST efficiency is expected to be canceled. However, due to the different multiplicities, the ST efficiencies estimated with the generic and the signal MC samples are expected to differ slightly. Thus, the uncertainty associated with the ST efficiency is not canceled fully, which results in a so called "tag bias" uncertainty. We study the tracking/PID efficiencies in different multiplicities, and take the combined differences between data and MC simulation, 0.6%, as the corresponding uncertainty.
The uncertainties associated with the quoted BF of D * + s → γD + s and the limited MC statistics are also considered, which lead to 0.8% and 1.1%, respectively.
In summary, using an e + e − collision data sample corresponding to 3.19 fb −1 collected at √ s = 4.178 GeV with the BESIII detector, we report the observation of D + s → pn and measure the absolute BF to be (1.21 ± 0.10 ± 0.05) × 10 −3 , where the first uncertainty is statistical and second systematic. The decay D + s → pn is confirmed and the precision of the BF measurement is much better than that of the previous measurement [6]. The anomalously large BF for D + s → pn explicitly shows that the weak annihilation process featured as a shortdistance dynamics is not the driving mechanism for this transition, while the hadronization process driven by non-perturbative dynamics determines the underlying physics. The measurement is important since similar annihilation effect is also present in other hadronic decays of charmed mesons. Relating this baryonic decay rate to the leptonic rate should provide important clues on how baryons are produced in hadronic interactions. The improved measurement also sets up the non-perturbative