Search for baryon and lepton number violation decay $D^{\pm}\to n(\bar{n})e^{\pm}$

Using a data set of electron-positron collisions corresponding to an integrated luminosity of ${\rm 2.93~fb^{-1}}$ taken with the BESIII detector at a center-of-mass energy of 3.773 GeV, a search for the baryon ($B$) and lepton ($L$) number violating decays $D^{\pm}\to n(\bar{n})e^{\pm}$ is performed. No signal is observed and the upper limits on the branching fractions at the $90\%$ confidence level are set to be $1.43\times10^{-5}$ for the decays $D^{+(-)}\to \bar{n}(n)e^{+(-)}$ with $\Delta|B-L|=0$, and $2.91\times10^{-5}$ for the decays $D^{+(-)}\to n(\bar{n})e^{+(-)}$ with $\Delta|B-L|=2$ , where $\Delta|B-L|$ denotes the change in the difference between baryon and lepton numbers.

Using a data set of electron-positron collisions corresponding to an integrated luminosity of 2.93 fb −1 taken with the BESIII detector at a center-of-mass energy of 3.773 GeV, a search for the baryon (B) and lepton (L) number violating decays D ± → n(n)e ± is performed. No signal is observed and the upper limits on the branching fractions at the 90% confidence level are set to be 1.43 × 10 −5 for the decays D +(−) →n(n)e +(−) with ∆|B − L| = 0, and 2.91 × 10 −5 for the decays D +(−) → n(n)e +(−) with ∆|B − L| = 2 , where ∆|B − L| denotes the change in the difference between baryon and lepton numbers.

I. INTRODUCTION
The search for new physics beyond the Standard Model (SM) is one of the major goals of particle physics. The matter-antimatter asymmetry of the universe is one prominent observations that cannot be explained within the SM, and as such is a serious challenge to our understanding of nature. This asymmetry suggests the existence of baryon number violation (BNV) [1]. While proton decay has been searched for decades but not yet observed, the search for decays of heavy mesons and baryons that are forbidden in the SM can provide an alternative probe to search for BNV. In most grand unified theories (GUTs) [2][3][4][5][6] and some SM extension models [7,8], baryon-number and lepton-number violation (LNV) is allowed, but the difference of baryon and lepton numbers is conserved (∆|B − L| = 0). Dimension-six operators allow processes with ∆|B − L| = 0 to proceed, mediated by heavy gauge bosons X with charge 4 3 or Y with charge 1 3 , as shown in Feynman diagrams in Fig. 1 (a) and (b) for D meson decays. Furthermore, there is another BNV process possible under dimensionseven operators, mediated by an elementary scalar field φ, as shown in Fig. 1 (c). In this process, the differ-ence of baryon and lepton number is changed by 2 units (∆|B − L| = 2). Reference [9] argues that the decay amplitudes of these two kinds of BNV processes are expected to be of comparable strength. Thus, experimental searches for these BNV decays probe new physics effects and test different models beyond the SM. The CLEO, BABAR, and CLAS experiments searched for BNV processes in D, B meson and hyperon decays [10][11][12], respectively, without finding evidence of a signal. Upper limits (ULs) were set on the decay branching fractions in the range of 10 −5 ∼ 10 −8 at the 90% confidence level (CL). Recently, the BESIII experiment searched for D + meson decays to a hyperon and an electron, i.e. D + →Λ(Σ 0 )e + with ∆(B − L) = 0 and D + → Λ(Σ 0 )e + with ∆(B − L) = 2. No signal was found and ULs of around 10 −6 were set on the decay branching fractions at the 90% CL [13]. It is natural to extend the search to D + meson decays to a (anti-)neutron and electron pair. A higher generation SUSY model [14] predicts the branching fraction of D 0 →pℓ + (ℓ + = e + , µ + ) to be less than 4.0 × 10 −39 , thus the decay D + →nℓ + is also expected to be of a comparable magnitude because it differs only by the change of a spectator quark.
In this paper, we report the first search for the BNV process D +(−) →n(n)e +(−) with ∆|B − L| = 0, and D +(−) → n(n)e +(−) with ∆|B−L| = 2 by using 2.93 fb −1 of electron-positron collision data taken at a center-ofmass energy of √ s = 3.773 GeV. Throughout this paper, the presence of charge-conjugated processes are implied unless explicitly stated otherwise.

II. BESIII DETECTOR AND MONTE CARLO SIMULATION
The BESIII detector is a magnetic spectrometer [15] located at the Beijing Electron Positron Collider (BEPCII) [16]. The cylindrical core of the BE-SIII detector consists of a helium-based multilayer drift chamber (MDC), a plastic scintillator time-of-flight system (TOF), and a CsI(Tl) electromagnetic calorimeter (EMC), which are all enclosed in a superconducting solenoidal magnet providing a 1.0 T magnetic field. The solenoid is supported by an octagonal flux-return yoke with resistive plate counter muon identifier modules interleaved with steel. The acceptance of charged particles and photons is 93% over 4π solid angle. The chargedparticle momentum resolution at 1 GeV/c is 0.5%, and the specific energy loss (dE/dx) resolution is 6% for the electrons from Bhabha scattering. The EMC measures photon energies with a resolution of 2.5% (5%) at 1 GeV in the barrel (end-cap) region. The TOF measures flight time of charged particle with a resolution 68 ps in the barrel region, and 110 ps in the end-cap region when the data sets in this analysis were collected.
Monte Carlo (MC) simulation samples, generated with a geant4-based [17] package [18], including the geometric and material description of the BESIII detector, are used to determine the detection efficiency, optimize the selection criteria and estimate the backgrounds. The analysis is performed in the framework of the BESIII Offline Software System [19], which takes care of the detector calibration, event reconstruction and date storage. The simulation includes the beam energy spread and initial state radiation (ISR) in the e + e − annihilations modeled with the generator kkmc [20]. The inclusive MC samples contain the production of DD pairs, the non-DD decays of the ψ(3770), the ISR production of the J/ψ and ψ(3686) states, and continuum processes, in which the known decay modes are modeled with evtgen [21] using branching fractions taken from the Particle Data Group (PDG) [22], and the remaining unknown decays from charmonium states are modeled with lundcharm [23,24]. The final-state radiation from charged particles is incorporated using photos [25]. In the signal MC sample, D + D − pairs are generated by the VSS model from evtgen [21] and the signal process is generated with a uniform momentum distribution in the phase space (PHSP) according to the conservation of angular momentum.

A. Analysis Method
At √ s = 3.773 GeV, D ± mesons are produced in pairs without the presence of any additional fragmentation particles. This property provides an ideal environment for investigating D ± meson decays with a doubletag (DT) method [26]. In this approach, the singletag (ST) D + meson is reconstructed in six hadronic decay modes D + → K − π + π + , D + → K − π + π + π 0 , D + → K 0 S π + , D + → K 0 S π + π 0 , D + → K 0 S π + π + π − and D + → K + K − π + , all of which have relatively large branching fractions and low background contamination. The DTs are then formed by reconstructing the other charm meson in the event in its decay to the signal mode. The decay branching fractions of the four signal modes (D +(−) →n(n)e +(−) and D +(−) → n(n)e +(−) ) are determined independently by where N tot ST and N DT are respectively the yields of ST and DT events in data summed over all ST decay modes. The effective signal detection efficiency in the presence of the ST D ± meson is calculated by are the corresponding detection efficiencies of the ST and DT method for the i th tag mode, respectively, and i sums over all ST decay modes.
To avoid possible bias, a blind analysis technique is performed in which the data in the interesting phase space region are viewed only after the analysis strategy is validated with MC simulation or data in the control region and then fixed.

B. Event Selection
Charged particles, including kaon, pion and electron/positron candidates, are reconstructed from the hit information in the MDC. The charged tracks, apart from in the case of pions from the decays of candidate K 0 S mesons, are required to have a distance of closest approach to the interaction point (IP) within ±10 cm in the beam direction and 1 cm in the plane perpendicular to the beam. The polar angle θ of charged tracks with respect to the z-axis of the MDC must satisfy |cosθ| < 0.93. Particle identification (PID), based on the information from the dE/dx and TOF measurements, is applied to the charged tracks and the CLs for the kaon and pion hypotheses (CL K,π ) are calculated. A kaon is identified by requiring CL K > CL π , and a pion by requiring CL π > CL K . To identify electrons, the deposited energy in the EMC is utilized, in addition to the dE/dx and TOF information, and the CLs are calculated for the electron, pion, kaon, and proton hypotheses (CL e,π,K,p ), individually. An electron is identified by requiring CL e /(CL e + CL π + CL K + CL p ) > 0.8.
Photon candidates are reconstructed from EMC showers and are required to have energy greater than 25 MeV in the barrel region (|cosθ| < 0.8), and 50 MeV in the end-cap region (0.86 < |cosθ| < 0.92). To suppress showers from electronic noise and those unrelated to the event under analysis, the EMC shower time is required to be within 700 ns of the start time of the event. The minimum opening angle between the photon candidate and all charged tracks is required to be greater than 10 • to avoid contamination from charged tracks showering in the EMC detector. The π 0 candidates are reconstructed from photon pairs by requiring the invariant mass (M γγ ) to be in the range (0.115,0.150) GeV/c 2 . To improve the kinematic resolution, a kinematic fit constraining M γγ to the nominal π 0 mass [22] is applied to the π 0 candidate. The kinematic variables after the kinematic fit are used in the subsequent analysis.
The K 0 S candidates are reconstructed from two charged tracks with opposite charges, polar angle in the range |cosθ| < 0.93 and points of closest approach to the IP within ±20 cm along the beam direction. No requirement on the distance of closest approach in the plane perpendicular to the z direction is applied. A vertex fit is performed on the two tracks on the assumption that they are pions from a common decay point. A further secondary vertex fit, which constrains the K 0 S to come from the beamspot, is applied to suppress background with the requirement L/σ L > 2, where L is the decay length, defined as the distance between the primary and secondary vertexes, and σ L is the corresponding resolution. The invariant mass (M ππ ) must satisfy [22].

C. Single-Tag Events
The ST D ± mesons are reconstructed in the six ST hadronic-decay modes and separated from background using two variables: the energy difference ∆E = E D − E beam and the beam constrained mass M BC = E beam − p 2 D , where E beam is beam energy, E D and p D are the energy and momentum of the ST D ± meson candidates in the rest frame of the e + e − system, respectively. When multiple candidates for a specific ST mode are found, the one with minimum |∆E| is retained. The ST candidate events are further required to have ∆E within (−55, +40) MeV for ST modes including a π 0 , and within (−25, +25) MeV otherwise.
To determine the yields of D ± meson for each ST decay modes, binned maximum-likelihood fits are performed to the M BC distributions in the range from 1.8365 to 1.8865 GeV/c 2 , as illustrated in Fig. 2 for the D + meson. In the fit, the signal is modeled by the MCsimulated shape convolved with a double-Gaussian function to take the resolution difference between data and MC simulation into account. The means and widths of the double-Gaussian function are free independent parameters in the fit. The combinatorial background is described by an ARGUS function [27]. Candidate events within M BC ∈ (1.863, 1.877) GeV/c 2 are kept for further analysis. The corresponding yields of the ST D ± mesons are determined by integrating the fitted signal lineshape in the same M BC range, as summarized in Table I The red and green dashed lines are signal and background, respectively. The blue solid lines are the sum of signal and background. The black points are data. DT and ǫ i,c2 DT are DT efficiencies for signal channels with ∆|B − L| = 0 and ∆|B − L| = 2, respectively.

D. Double-Tag Events
DT signal candidates are selected from the sample of ST D ± events by requiring an electron candidate and no additional charged tracks in the event. A two-constraint (2C) kinematical fit is performed by imposing energy and momentum conservation, and constraining the invariant mass of ST D ± candidates as well as the mass of the electron-(anti-)neutron system to be the known mass of D ∓ meson [22], in which the (anti-)neutron is regarded as a missing particle with unknown mass. The fit is required to converge, but no further selection on the χ 2 of the fit is applied. The momentum and invariant mass of the (anti-)neutron obtained from the kinematic fit (p n/n and M n/n ) are recorded for the subsequent analysis. To suppress backgrounds, candidate events are required to possess a shower in the EMC around the fitted direction of the (anti-)neutron within an opening angle of 30 • . If there are several showers in this region, the one with the largest energy is selected.
MC studies with a generic event-type analysis tool, TopoAna [28], indicate that the backgrounds in the selected samples are dominated by semi-leptonic decays of D ± meson with K 0 L and π 0 mesons in the final state. Taking into account the result of the MC simulation, the selected showers in the EMC are further required to lie no more than 10 • (15 • ) from the direction of the neutron (anti-neutron) candidate. In addition, a Multivariate Data Analysis (MVA) based on the shower shape in the EMC is performed based on a Gradient Boosted Decision Trees (GBDT) algorithm. The utilized variables include the total deposited energy E tot , the number of hit crystals in the EMC N hit , and the A20 and A42 Zernike moments as defined in Ref. [29]. In order to train and test the MVA, high-purity training and testing samples, including the (anti-)neutron signal, as well as K 0 L and photon backgrounds, are obtained from data using a selection that is independent of any EMC information. The (anti-)neutron sample comes from the decay process J/ψ →pnπ + + c.c., the K 0 L sample from J/ψ → KπK 0 L and the photon sample from J/ψ → ρπ 0 with π 0 → γγ. Studies show that the distributions of the shower-shape variables have a significant dependence on the momentum of the (anti-)neutron. Therefore, the MVA is performed in separate (anti-)neutron momentum bins of width 100 MeV/c. For a specific (anti-)neutron momentum bin, the training and testing background samples are reweighted according to their expected momentum lineshapes, which are obtained from the inclusive MC samples. The distributions of GBDT values in the different momentum bins for the anti-neutron and neutron as well as for the backgrounds are shown in Fig. 3. The selection on the GBDT values, optimized by maximizing the quantity ǫ/(1.5 + √ B) [30], where ǫ is the relative efficiency in the MVA signal sample and B is the number of background events normalized to match the luminosity of data in the inclusive MC sample, are applied and shown in Fig. 3.
The detection efficiencies for finding a matching shower, and for the selection on the opening angle and GBDT value are evaluated from data using the large and high purity control sample of J/ψ → pnπ decays. The efficiencies of finding a matching shower and for the requirement on the opening angle are studied as a function of two variables: the (anti-)neutron momentum and cosθ, following the procedure described in Ref. [31]. The above efficiencies are directly applied to the signal MC sample with a sampling approach.

E. Signal Extraction and Fitting
The mass distributions of (anti-)neutron M n/n from the kinematic fit, after all selection cuts, are shown in Fig. 4 for the four decay processes, where the upper two plots are for the processes D + →ne + and D − → ne − with ∆|B − L| = 0, and the lower two plots are for the processes D − →ne − and D + → ne + with ∆|B − L| = 2. No obvious signal is observed. The DT detection efficiencies for the different ST modes are determined to be 18.65% for D + →ne + , 19.92% for D − → ne − , 19.68% for D + → ne + and 18.85% for D − →ne − . An unbinned maximum-likelihood fit is performed to each M n/n distribution as shown in Fig. 4, individually. In the fit, signal and background are modeled by the MC-simulated shapes obtained from signal and inclusive MC samples, respectively. The yields of signal and background are left free in the fit and the returned values are shown in the plots. Since no significant signal is observed, conservative ULs are set by combining the two processes with ∆|B − L| = 0 and those with ∆|B − L| = 2, respectively, as described below.

F. Assignment of Systematic Uncertainties
The systematic uncertainties related to the efficiency for reconstructing the tag side cancel due to the DT method. The sources of possible systematic bias that remain include those associated with the DT-side selection efficiency, and the ST and DT yields extraction.
The uncertainties associated with the DT-side eventselection efficiency include those from the electron tracking and PID efficiencies, from the 2C kinematic fit, from the efficiency of finding a matched shower, and from requirements on the angle and the GBDT values.
The uncertainties on the tracking and PID efficiencies for electrons and positron are studied in control samples, as described in Ref. [32]. These efficiencies are studied in two-dimensional bins of momentum versus cosθ for data and MC simulation, individually. The average relative differences on the efficiencies between data and MC simulation, which is calculated by weighting the corresponding values according to the distribution of the electron/positron signal, are assigned as the systematic uncertainties.
The uncertainty associated with the 2C kinematic fit is studied with a high-purity control sample D → K 0 S eν e → ππeν e decays by mimicking the K 0 S ν e as a missing system. The same kinematic fit is performed on this sample, and the efficiency of decays surviving this fit is measured and compared with the corresponding efficiency in MC simulation. This comparison is made as a function of the invariant mass of the K 0 S ν e system, and the relative difference around the neutron mass is taken as the systematic uncertainty.
The detection efficiencies for finding a matching shower and for the angle requirement in the EMC for (anti-)neutron are estimated from a control sample of J/ψ → pnπ decays in two-dimensional bins of momentum versus cosθ, and then applied directly to the signal MC samples with a sampling approach. The dominant source of potential bias from this approach arises from the different physics environment in the EMC between the signal and control sample, as well as the statistical uncertainty associated with the size of the control sample. To estimate the size of this potential bias, we compare the efficiencies between the signal MC sample and the MC-simulated control sample J/ψ → pnπ, and assign the small differences observed as the systematic uncertainties. The uncertainty associated with the sample size is determined by standard error propagation.
The efficiency of the requirement on the GBDT value is determined from the control sample of J/ψ → pnπ decays in different (anti-)neutron momentum bins. Two sources of potential bias are considered: background in the control sample and the choice of momentum binning. The amount of possible contamination is estimated by fitting the missing mass in the control sample in the different momentum bins, and its full effect is determined on the efficiency measurement and taken as the corresponding systematic efficiency. The possible bias associated with the momentum binning is evaluated by varying the bin size and taking the maximum change in the measured efficiency as the systematic uncertainty.
The uncertainties associated with DT selection efficiency are summarized in Table II, and the total uncertainties are the quadratic sum of the individual values. The uncertainties on the ST yields are estimated to be 0.5% by studying the variation in results when using a different fit range of (1.8415,1.8865) GeV/c 2 , describing the combinatorial background with a 3-order polynomial rather than an ARGUS function, and by imposing a different endpoint of 1.8863 or 1.8867 GeV/c 2 for the ARGUS function [33].
The uncertainty associated with the fit of the M n/n distribution arises from the imperfect knowledge of the background shape and the choice of fit range, which will be considered in the determination of the upper limits.

G. Determination of the Upper Limits
Since no signal is observed, an UL is set on ∆|B − L| = 0 processes by performing a fit to the M n/n distributions of D + →ne + and D − → ne − , similar to that of Sec. III E, but with the signal yields fixed. A UL is also set on ∆|B−L| = 2 processes from a fit to the D − →ne − and D + → ne + distributions.
The fixed signal yields in the fits correspond to different branching fraction assumptions, according to the effective detection efficiencies, ST yields and the uncertainties. Likelihood values are obtained as a function of the branching fraction, where the effects of systematic uncertainties associated with DT efficiencies and ST yields are included by convolving the likelihood distribution with Gaussian functions of mean zero and width equal to their absolute uncertainties, as described in Refs. [34]. The ULs on the branching fraction at the 90% CL are calculated by integrating the likelihood distribution from zero to 90% of the total curve. To take into account the effects on the imperfect knowledge of background shape and the fit range on the M n/n distributions, alternative fits are performed with different assumptions for background lineshape and fit range. The most conservative ULs obtained at the 90% CL are shown in Fig. 5, which are taken as the final results. These ULs are B D + →ne + < 1.43×10 −5 and B D + →ne + < 2.91×10 −5 . Studies of ensembles of simulated experiments ('toy MC') that contain no signal give results that are consistent with these measurements within 2σ.

IV. SUMMARY
In summary, by analyzing e + e − collision data with an integrated luminosity of 2.93 fb −1 at √ s = 3.773 GeV taken with the BESIII detector, we search for the BNV decays D + →ne + with ∆|B − L| = 0 and D + → ne + with ∆|B − L| = 2 for the first time. No signal is found and the ULs on branching fraction are determined to be B D + →ne + < 1.43 × 10 −5 and B D + →ne + < 2.91 × 10 −5 for the processes with ∆|B − L| = 0 and ∆|B − L| = 2 , respectively. More data at this collision energy is being collected, up to an integrated luminosity of around 20 fb −1 [35]. With this larger sample, and assuming no signal, it will be possible to improve the ULs by around a factor of three.

ACKNOWLEDGMENTS
The BESIII collaboration thanks the staff of BEPCII and the IHEP computing center for their strong support. This work is supported in part by National Key R&D Program of China under Contracts Nos.