c : Observation of

: Using data taken at 29 center-of-mass energies between 4.16 and 4.70 GeV with the BESIII detector at the Beijing Electron Positron Collider corresponding to a total integrated luminosity of approximately 18.8 , the process is observed for the first time with a statistical significance of . The average Born cross sections in the energy ranges of (4.160, 4.380) GeV, (4.400, 4.600) GeV and (4.610, 4.700) GeV are measured to be fb, fb and fb, respectively, where the first uncertainties are statistical and the second are systematic. The line shapes of the and invariant mass spectra are consistent with phase space distributions, indicating that no hexaquark or di-baryon state is observed.


INTRODUCTION
One of the most fundamental questions in hadron physics is related to the mechanism of color confinement in Quantum Chromodynamics (QCD). Colorsinglet hadronic configurations of quarks and gluons can form bound states or resonances. Besides the well-known combinations of qq for mesons and qqq for baryons, other combinations, such as gqq for hybrid states [1], multi-gluons for glueball states [2], qqqq for tetraquark states [3], qqqqq for pentaquark states [4] and hexaquark states (qqqqqq), are also allowed by QCD. Di-baryon and hexaquark states have been searched for in a range of nucleon-nucleon scattering reactions. Recently, an isoscalar resonant structure was observed in the isoscalar two-pion fusion process pn → dπ 0 π 0 [5] by the WASA Collaboration and was later confirmed in the other twopion fusion processes pn → dπ + π − [6] and pp → dπ + π 0 [7], and the two-pion non-fusion process pn → ppπ 0 π − [8] and pn → pnπ 0 π 0 [9]. This state was denoted by d * (2380) following the convention used for nucleon excitations. These observations indicate the possibility of the existence of hexaquark and di-baryon configurations. In 2021, the BESIII Collaboration reported the search for hexaquark and di-baryon states in examining the invariant mass spectra of two baryons in the process e + e − → 2(pp) [10], and no significant signal was observed.
Analyzing data sets corresponding to a total integrated luminosity of approximately 18.8 fb −1 taken at center-of-mass energies √ s between 4.16 and 4.70 GeV with the BESIII detector, we present in this paper the first measurement of the cross section of the process e + e − → pppnπ − + c.c.. We search for the d * (2380) and other possible hexaquark or di-baryon states with the data samples with energies above 4.60 GeV, where thē pn system with a mass around 2.4 GeV for d * (2380) is kinematically accessible. The mass range of thepn system around 2.4 GeV/c 2 , in which the d * (2380) might contribute, is covered by the data samples with energies above 4.60 GeV. Throughout this paper, charge conjugation is always implied unless explicitly stated, and in discussing systematic uncertainties.

THE BESIII DETECTOR AND DATA SAMPLES
The BESIII detector [11] records symmetric e + e − collisions provided by the BEPCII storage ring [12], which operates in the center-of-mass energy range from 2.0 to 4.95 GeV. BESIII has collected large data samples in this energy region [13]. The cylindrical core of the BESIII detector covers 93% of the full solid angle and 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 identification modules interleaved with steel. The charged-particle momentum resolution at 1 GeV/c is 0.5%, and the specific energy loss (dE/dx) resolution is 6% for 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 time resolution in the TOF barrel region is 68 ps, while that in the end cap region is 110 ps. The end cap TOF system was upgraded in 2015 using multi-gap resistive plate chamber technology, providing a time resolution of 60 ps [14].
The data sets were collected at 29 center-of-mass energies between 4.16 and 4.70 GeV. The nominal energies of the data sets from 4.16 to 4.60 GeV are measured by the di-muon process e + e − → (γ ISR/FSR )µ + µ − [15,16], where the subscript ISR/FSR stands for the initial-state or final-state radiation process, respectively. The data sets from 4.61 to 4.70 GeV are calibrated by the process e + e − → Λ + cΛ − c [17]. The integrated luminosity L int is determined using large-angle Bhabha scattering events [17,18]. The total integrated luminosity of all data sets is approximately 18.8 fb −1 .
The response of the BESIII detector is modeled with Monte Carlo (MC) simulations using the software framework boost [19] based on geant4 [20], which includes the geometry and material description of the BESIII detector, the detector response and digitization models, as well as a database that keeps track of the running conditions and the detector performance. Large MC simulated event samples are used to optimize the selection criteria, evaluate the signal efficiency, and estimate background contributions.
Inclusive MC simulation samples are generated at different center-of-mass energies to study potential background reactions. These samples consist of open charm processes, the ISR production of vector charmonium and charmonium-like states, and the continuum processes incorporated in kkmc [21]. The known decay modes are modeled with evtgen [22] using branching fractions taken from the Particle Data Group (PDG) [23], and the remaining unknown decays of the charmonium states are simulated with lundcharm [24]. Final-state radiation from charged final-state particles is incorporated with photos [25]. The signal MC simulation sample of e + e − → pppnπ − at each energy point is generated with the events being uniformly distributed in phase space.

DATA ANALYSIS
Events with two positive and two negative charged tracks are selected. For each charged track, the polar angle in the MDC with respect to the z direction must satisfy |cosθ| < 0.93. All charged tracks are required to originate from the interaction region, defined as R xy < 1 cm and |V z | < 10 cm, where R xy and |V z | are the distances from the point of closest approach of the tracks to the interaction point in the x − y plane and in the z direction, respectively. The combined dE/dx and TOF information are used to calculate particle identification (PID) confidence levels for the pion, kaon, and proton hypotheses. Each track is assigned as the particle hypothesis with the highest confidence level. The final state in the e + e − → pppnπ − process is reconstructed with three (anti-)protons and one π − .
Since the neutron can not be well reconstructed with the BESIII detector, the signal process is determined via the recoiling mass of the reconstructed charged particles, defined as where E e + e − and P e + e − are the center-of-mass energy and the momentum of the e + e − system, respectively; E pppπ − and P pppπ − are the total reconstructed energy and total momentum of the pppπ − system, respectively. Events with M rec greater than 0.8 GeV/c 2 are kept for further analysis.
Studies based on the inclusive MC simulation samples [26] show that no peaking background events survive the selection criteria. To further suppress background events, two additional selection criteria are imposed on the accepted candidate events. First, the invariant mass M pπ − of the reconstructed pπ − system is required to be outside the Λ signal region, i.e. |M pπ − − 1.115| > 0.010 GeV/c 2 , to remove the possible background associated with Λ decays. Here, 1.115 GeV/c 2 is the known Λ mass [29], and 0.010 GeV/c 2 corresponds to about three times the mass resolution. Second, the invariant mass of ppp (M ppp ) must be less than 3.6 GeV/c 2 due to the remaining neutron and pion in the event.
The M rec distribution of the accepted candidates after the above selection criteria from the combined data sets is displayed in Fig. 1, where a significant neutron signal is observed. The signal yield is determined by a maximum likelihood fit to this distribution. In the fit, the signal is represented by the luminosity weighted MC-simulated shape convolved with a Gaussian function and the remaining background is described by a linear function. From the fit, the signal yield is determined to be 123 ± 14. The statistical significance of the signal is determined to be 11.5σ, which is evaluated as −2 ln(L 0 /L max ), where L max is the maximum likelihood of the nominal fit and L 0 is the likelihood of the fit without involving the signal component. The change of the degree of freedom is 1. The neutron signal region is de-fined as M rec ∈ (0.925, 0.968) GeV/c 2 and the corresponding sideband regions are defined as M rec ∈ (0.857,0.900) ∪ (0.990,1.033) GeV/c 2 . Figure 2 shows the comparisons of the momentum and polar angle distributions of the neutron of the accepted candidate events between data and signal MC simulation samples, where the data distribution is from the combined data sets and the MC simulation distribution has been weighted by the signal yields in data. The invariant mass of any two or three particles, the momentum and cosθ distributions of the other final state particles have also been examined. The agreement between data and MC simulation allows to determine the detection efficiency with the signal MC simulation events generated uniformly distributed in the five-body phase space.
To search for hexaquark and di-baryon states, thepn invariant mass spectrum is examined. Figure 3 shows the ppπ − andpn invariant mass spectra of the candidate events for the reaction e + e − → pppnπ − . In the fit to Mpn, the signal is represented by the luminosity weighted phase space MC simulation shape and the remaining combinatorial background is described by a linear function. The goodness-of-fit is χ 2 /ndf = 2.10/2. Here, ndf is the number of degrees of freedom. Compared to the phase space hypothesis, no obvious structure is observed.

AVERAGE CROSS SECTIONS
In each data set, only a few events have been observed in the neutron signal region, with a statistical significance of less than 3σ. To obtain significant neutron signals the data sets are combined into three sub-samples in the energy ranges of (4.160, 4.380), (4.400, 4.600) and (4.610, 4.700) GeV for further analysis.
The average observed cross section for e + e − → pppnπ − is calculated by where N sig j is the number of signal events from the j-th combined data set, L i and ǫ i are the integrated luminosity and efficiency of data set i, , respectively, i represents the i-th energy point in j-th sub-data set. The detection efficiency is corrected by the PID and tracking efficiencies correction factors, f PID and f trk , which are determined to be 0.92 and 0.98 by weighting the differences between data and MC simulation efficiencies in different momentum ranges, respectively. Inserting the numbers which are listed in Table 1 into Eq. 2 yields the average observed cross sections (19.4±5.1±1.0) fb, (42.8±9.8±2.3) fb and (54.2 ± 8.6 ± 2.9) fb for the three sub-data sets, respectively, where the first uncertainties are statistical and the second are systematic.
To measure the average Born cross section of e + e − → pppnπ − , a similar lineshape as that of e + e − → 2(pp) [10] is assumed to determine the ISR and vacuum polariza-5 No. X Observation of e + e − → pppnπ − + c.c. 6 The obtained Born cross sections are then used as input in the generator and the cross section measurements are iterated with the updated detection efficiencies. This process is repeated until the (1+δ γ )·ǫ values become stable at all energies, i.e. the difference of (1+δ γ )·ǫ between the last two iterations is less than 4%. Figure 4 shows the obtained average Born cross sections in the defined subsamples. The average Born observed cross sections are calculated with Eq. 3, and the results are (21.5±5.7±1.2) fb, (46.3±10.6±2.5) fb and (59.0±9.4±3.2) fb for the three sub-data sets, respectively, where the first uncertainties are statistical and the second are systematic. Two different functions are used to compare the trend of the average Born cross section to a reaction where a similar behaviour is expected. The first one is a simple five-body energydependent phase space lineshape [10,27] and the second one is an exponential function [10,28], which are shown in Figure 4. The exponential function is constructed as where p 0 and p 1 are free parameters, M th = (3m p + m n + m π − ), m p , m n , and m π − are the known masses of p, n, and π − taken from the PDG [29]. This is similar to the one used for the cross section lineshape of e + e − → 2(pp) in Ref. [10], as they are similar reactions where one of thep has been exchanged bynπ − . However, it should be noted that the two functions in Figure 4 are not fit results, but drawn with arbitrary scale factors for comparison since a qualitative fit is not possible due to the limited statistics. The systematic uncertainties in the cross section measurements will be discussed in the next section.

SYSTEMATIC UNCERTAINTY
In the cross section measurements, the systematic uncertainties mainly comes from the integrated luminosity, tracking efficiency, PID efficiency, ISR correction, M rec fit, and veto of background events associated with Λ decays.
The integrated luminosity of the data set is measured by large-angle Bhabha scattering events, and the uncertainty in the measurement is 1.0% [18], which is dominated by the precision of the MC generator used for efficiency correction. The tracking and PID efficiencies have been studied with high purity control samples of J/ψ → ppπ + π − and ψ(3686) → π + π − J/ψ → π + π − pp decays [30,31]. The differences of the tracking and PID efficiencies between data and MC simulation in different transverse momentum and total momentum ranges are obtained separately. The averaged differences for the tracking (PID) efficiencies are corrected by the factors f trk (f PID ) as mentioned in Sec. 4. The uncertainties of the tracking and PID efficiencies are reweighted by the p/p and π + /π − momenta of the signal MC simulation events. The reweighted uncertainties for tracking (PID) efficiencies, 0.1% (0.3%) per p, 0.1% (0.4%) perp, 1.0% (0.5%) per π + and 0.8% (0.4%) per π − , are assigned as the systematic uncertainties. Adding them linearly gives the total systematic uncertainties due to the tracking and PID efficiencies to be 1.1% and 1.6% for the process e + e − → pppnπ − , and 1.3% and 1.9% for the process e + e − → pppnπ + , respectively. The input Born cross sections in the generator are iterated until the (1 + δ γ ) · ǫ values converge. The largest difference of (1 + δ γ ) · ǫ between the last two iterations at all energy points, 3.2%, is taken as the corresponding systematic uncertainty.
Three different tests were performed to estimate the uncertainty associated with the M rec fit. The fit range is increased or decreased by 5 MeV/c 2 . The background shape is replaced with a second-order Chebychev polynomial function, and the signal shape is replaced with an MC simulation-derived shape convolved with a Gaussian function. The quadrature sum of these changes, 3.6%, is taken as the relevant uncertainty.
The systematic uncertainty due to the veto of Λ background events is estimated by changing the Λ veto mass window from ±3σ to ±5σ, where σ is the invariant mass resolution and the value is 3 MeV/c 2 . The change of the measured cross section, 0.03%, is assigned as the uncertainty.
Adding the above systematic uncertainties summarized in Table 2 in quadrature yields the total systematic uncertainties of 5.3% and 5.4%, for the processes e + e − → pppnπ − and e + e − → pppnπ + , respectively. The average systematic uncertainty, 5.35%, is taken as the total systematic uncertainty in the cross section measurement for the process e + e − → pppnπ − + c.c..

SUMMARY
By using the data sets taken at the center-of-mass energies between 4.16 and 4.70 GeV, the process e + e − → pppnπ − + c.c. has been observed for the first time with a statistical significance of 11.5σ. The average Born cross sections in the three energy ranges of (4.160, 4.380), (4.400, 4.600) and (4.610, 4.700) GeV are measured to be (21.5±5.7±1.2) fb, (46.3±10.6±2.5) fb and (59.0±9.4±3.2) fb, respectively, where the first uncertainties are statistical and the second systematic. The Born cross section close to threshold is larger than would be expected from five-body phase space. The lineshape of the average Born cross sections for the process e + e − → pppnπ − +c.c. shows similar behaviour to that of the process e + e − → 2(pp). The shape of the invariant-mass spectra ofpn and ppπ − are in good agreement with the phase-space distributions, thereby indicating no hexaquark or di-baryon state observed with the current data sample size.
The BESIII collaboration thanks the staff of BEPCII and the IHEP computing center for their strong support.