Observation of the $\psi(3686)$ decays into $\Sigma^{+}\bar{\Sigma}^{-}\omega$ and $\Sigma^{+}\bar{\Sigma}^{-}{\mathcal{\phi}}$

Based on $(27.08\pm 0.14)\times10^{8}$ $\psi(3686)$ events collected with the BESIII detector operating at the BEPCII collider, the $\psi(3686)\to\Sigma^{+}\bar{\Sigma}^{-}\omega$ and $\Sigma^{+}\bar{\Sigma}^{-}\phi$ decays are observed for the first time with statistical significances of 13.8$\sigma$ and 7.6$\sigma$, respectively. The corresponding branching fractions are measured to be $\mathcal{B}(\psi(3686)\to\Sigma^{+}\bar{\Sigma}^{-}\omega)=(1.90 \pm 0.18 \pm 0.21) \times 10^{-5}$ and $\mathcal{B}(\psi(3686)\to\Sigma^{+}\bar{\Sigma}^{-}\phi)=(2.96 \pm 0.54 \pm 0.41) \times 10^{-6}$, where the first uncertainties are statistical and the second systematic.


I. INTRODUCTION
Quantum Chromodynamics (QCD) is the theory of strong interactions and it has been rigorously tested at high energies [1]; in the low-energy regime, the non-Abelian nature of the theory requires a non-perturbative approach which must rely either on lattice QCD (LQCD) or on QCD-inspired models.Therefore, experimental measurements in the low-energy regime are crucial to verify the pertinent models, calibrate their parameters, and stimulate the development of new theoretical computations [2].As a result, many interesting properties associated with the strong decays of / and (3686) mesons have been investigated, with the aim to advance our knowledge about the QCD in the interplay between perturbative and non-perturbative strong interaction regime [3].BESIII has collected large data samples of vector charmonia from  +  − annihilations, such as / and (3686), which provide us the opportunity to conduct detailed experimental studies on rare decays of these states and further explore the intermediate structures inherent in these processes [4].

II. BESIII DETECTOR AND MONTE CARLO SIMULATION
The BESIII detector [15] records  +  − collisions provided by the BEPCII storage ring [16] in the center-ofmass energy range from 2.0 to 4.95 GeV, with a peak lu-minosity of 1×10 33 cm −2 s −1 achieved at √  = 3.77 GeV.BESIII has collected large data samples in this energy region [17][18][19].The cylindrical core of the BESIII detector covers 93% of the full solid angle and consists of a heliumbased 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 magnet is supported by an octagonal flux-return yoke with modules of resistive plate muon counters (MUC) interleaved with steel.The chargedparticle momentum resolution at 1 GeV/c is 0.5%, and the d/d resolution is 6% for the electrons from Bhabha scattering at 1 GeV.The EMC measures photon energy with a resolution of 2.5% (5%) at 1 GeV in the barrel (end-cap) region.The time resolution of the TOF barrel part is 68 ps, while that of the end-cap part 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, which benefits ∼ 83% of the data used in this analysis [20][21][22].
Monte Carlo (MC) simulated data samples produced with a geant4-based [23] software package, which includes the geometric description [24] of the BESIII detector and the detector response, are used to optimize the event selection criteria, estimate the signal efficiency and the level of background.The simulation models the beam energy spread and initial-state radiation in the  +  − annihilation using the generator kkmc [25,26].The inclusive MC sample includes the production of the (3686) resonance, the initial-state radiation production of the / meson, and the continuum processes incorporated in kkmc.Particle decays are generated by evtgen [27,28] for the known decay modes with BFs taken from the Particle Data Group (PDG) [29] and lundcharm [30,31] for the unknown ones.Final-state radiation from charged final-state particles is included using the photos package [32].
To determine the detection efficiency for each signal process, signal MC samples are generated with a modified data-driven generator BODY3 [27,28], to simulate contributions from different intermediate states in data for a given three-body final state, as discussed in Sec.IV.

III. EVENT SELECTION
In the analysis of the -mode, the Σ + ( Σ− ) and  particles are reconstructed via the Σ + ( Σ− ) →  0 (p 0 ) and  →  +  −  0 decays.Within these processes, the  0 is reconstructed through the  0 →  decay.For the -mode, to improve the detection efficiency, a partial reconstruction of (3686) → Σ + Σ−  →  0 p 0  +  − is performed, where only one  0 is required to be reconstructed while the other  0 is treated as a missing particle.
All charged tracks are required to satisfy |cos | < 0.93, where  is the polar angle defined with respect to the  axis, which is the symmetry axis of the MDC.The charged tracks not originating from Σ + ( Σ− ) decays are required to have their point of closest approach to the interaction point (IP) within 10 cm along the -axis, and within 1 cm in the transverse plane.For the charged tracks from the Σ + ( Σ− ) decays, the distance of closest approach to the IP must be less than 2 cm in the transverse plane, due to the hyperons' lifetime.The measurements of the flight time in the TOF and / in the MDC for each charged track are combined to compute particle identification (PID) confidence levels for the pion, kaon, and proton hypotheses.The tracks are assigned to the particle type with the highest confidence level.
Photon candidates are identified using showers in the EMC.The deposited energy of each shower must be greater than 25 MeV in the barrel region (| cos | < 0.80) or greater than 50 MeV in the end cap region (0.86 < | cos | < 0.92).To suppress electronic noise and energy depositions not associated with the event, the EMC cluster timing from the reconstructed event start time is further required to satisfy 0 ≤  ≤700 ns.The invariant mass of the  0 candidates reconstructed from a  pair is constrained to the  0 known mass [29] by a kinematic fit, and the  2 1C is further required to be less than 20.For the -mode, in order to suppress the potential backgrounds and improve the mass resolution, a seven-constraint (7C) kinematic fit is performed for the (3686) →  0 p 0  +  −  0 hypothesis by enforcing energy-momentum conservation and constraining the invariant mass of each of the three photon pairs to the nominal  0 mass.If there are more than one combination in the event, the one with the smallest  2 7C is chosen.Furthermore, the selection  2 7C <45 is applied, by optimizing the figure-of-merit (FOM) defined as  √ + , where  denotes the number of signal events obtained from the MC simulation, while  is the number of background events obtained from the inclusive MC sample.For the -mode, since we reconstruct only one  0 meson to enhance the reconstruction efficiency, a two-constraint (2C) kinematic fit is performed on the p +  −  0 combinations, requiring that the missing mass corresponds to the nominal  0 mass, as well as the invariant mass of the photon pair.If there is more than one combination in the event, the one with the smallest The /-related backgrounds are vetoed by requiring the  +  − ( 0  0 ) recoil mass outside the / mass window, and the -related backgrounds are suppressed by requiring the invariant mass of  0  0  0 outside the  mass window.All the mass windows are determined according to their mass resolutions and are listed in Table I, where  and  denote the invariant mass and the known mass [29], respectively, and  de-notes the recoil mass.For the -mode, to suppress the background channels with a number of photons different from six, the  2 4 (6p +  − )< 2 4 (5p +  − ) and  2 4 (6p +  − )< 2 4 (7p +  − ) requirements are also added.-mode After applying all the selection criteria, the twodimensional (2D) distributions of the p 0 invariant mass versus the  0 invariant mass in data are shown in Fig. 1, where clear Σ + Σ− signal peaks are seen in the two decay modes.The one-dimensional (1D) Σ + ( Σ− ) signal and sideband regions are defined as (1176, 1197) and (1144, 1165) or (1208, 1229) MeV/ 2 for the -mode, (1176, 1200) and (1140, 1164) or (1212, 1236) MeV/ 2 for the -mode, respectively.The 2D Σ + Σ− signal region is defined as the square region with both  0 and p 0 combinations lying in the 1D Σ + ( Σ− ) signal regions.The Σ + Σ− sideband I regions are defined as the square regions with either one of the  0 or p 0 combinations located in the 1D Σ + ( Σ− ) sideband regions and the other in the 1D signal region.The sideband II regions are defined as the square regions with both  0 and p 0 combinations located in the 1D Σ + ( Σ− ) sideband regions.The potential remaining backgrounds are investigated with the inclusive (3686) MC samples, using the eventtype analysis tool TopoAna [33].According to our strategy for extracting signal yields, we consider decay pro-cesses involving (for -mode) or (for -mode) as socalled peaking backgrounds.Furthermore, peaking backgrounds were found only for the -mode, mainly from the (3686) →  1,2 ,  1,2 → Λ Λ processes, and contaminations from these processes were found to be negligible.For the non-peaking backgrounds, we consider them as a smooth distribution in the invariant mass distribution.The possible non-Σ + ( Σ− ) and non-Σ + Σ− peaking backgrounds from (3686) decays are studied with sideband events, and will be considered when extracting the signal yields later.We determine the corresponding peaking contributions, N SD I = 63.0 ± 14.5 and N SD II = 61.8± 12.7 for the -mode, N SD I = 29.5 ± 8.4 and N SD II = 26.2± 8.5 for the -mode, by fitting the data in the 2D sideband I and sideband II regions.In the fitting process, the description of the peaking contribution and the smooth background is consistent with the strategy of fitting the Σ + Σ− signal region events in Section IV.
To investigate the possible quantum electrodynamics (QED) background, the same selection criteria are applied to data samples collected at center-of-mass energies of 3.65 GeV and 3.773 GeV, corresponding to values of integrated luminosity of 454 pb −1 and 2931.8 pb −1 [34], respectively.Only a peaking contribution of 1.6 is found for the -mode at 3.773 GeV, and no peaking background contribution is seen for the -mode.We also consider the system uncertainty caused by the interference effect between the resonance decay and the continuum processes, but the corresponding effect is negligible.

IV. SIGNAL YIELDS
The signal yields are determined by extended unbinned maximum likelihood fits on the  +  −  0 (-mode) or  +  − (-mode) invariant mass distributions in the Σ + Σ− signal regions.The total probability density function consists of a signal component and various background contributions.The signal component is modelled with the MC-simulated shape convolved with a Gaussian function to account for the possible difference in the mass resolution between data and MC simulation.
For the background of the 2D sideband I region and sideband II region, the shape is described using the MCsimulated shape of the signal, while the number of events is fixed to a normalization value, i.e. 1  2 According to the calculation in section III, the background contributions from all 2D sideband regions are fixed as 16.1 events for the -mode, and 8.2 events for the -mode.
The remaining smooth backgrounds for the two modes are modelled by a second-order polynomial function and an Argus function, respectively.The numbers of fitted signal events are 198.7 ± 18.9 and 55.2 ± 10.0 for -mode and -mode, respectively.The corresponding statistical significances are determined to be 13.8 and 7.6 for -mode and -mode, respectively.The statistical significance is estimated by the likelihood difference between the fits with and without the signal component, taking the change in the degrees of freedom into account.The potential intermediate states are investigated via Dalitz plots and their 1D projections.Figure 4 shows the invariant mass distributions of the different two-body combinations for the two decay modes, where the background contributions have been subtracted by using the normalized sideband events.No obvious structure can be found in any of the mass distributions.Nevertheless, the experimental distributions are not consistent with the signal MC sample generated according to the phase space (PHSP) distribution.To improve the accuracy of the detection efficiency, the PHSP model is replaced by the modified data-driven generator BODY3, where the MC-simulated events are sampled according to the Dalitz   distribution in the data.The comparisons are shown in Fig. 4. The MC sample of BODY3 describes the data better.

V. RESULTS FOR BRANCHING FRACTIONS
The BF for each signal decay is calculated by: Here,   is the number of signal events determined by the fit,  (3686) is the total number of (3686) events, V denotes the vector meson (), ℬ  is the BF of the -th intermediate state taken from the PDG [29],  is the detection efficiency, which is determined by the MC simulation based on the BODY3 model.The corresponding numerical values are listed in Table II.

VI. SYSTEMATIC UNCERTAINTIES
As sources of systematic uncertainties in the BF measurements we consider the tracking and PID efficiencies, the photon detection efficiency, the  0 reconstruction efficiency, the kinematic fit, the procedure for signal extraction and background subtraction, the MC simulation modeling, and so on.They are described as follows.
(i) Tracking and PID: The uncertainties due to the tracking efficiency are estimated with the control samples (3686) →  +  − /, / →  0   ±  ∓ , and / → p +  − , and are determined to be 1% [35][36][37] for each charged track.The PID uncertainties are determined to be 1% for each charged track as well, based on the same samples used to estimate the tracking uncertainties.
(ii) Photon detection: The difference in the photon detection efficiencies between data and MC simulation is studied using a control sample of  +  − →  +  − and found to be less than 0.5% for each photon.
(iii)  0 reconstruction: Based on the control samples (3686) →  0  0 /, / →  +  − and  +  − →  0 at √  = 3.773 GeV, the relative difference of the  0 reconstruction efficiencies between data and MC simulation has been estimated, and the results on the two datasets are consistent with each other.We have studied the relative difference as a function of the polar angle and of the total  0 momentum , and it decreases linearly as a function of the momentum following (0.06 − 2.41 × )%.The resulting systematic uncertainties of the  0 reconstruction efficiency are determined to be 2.9% and 0.8% for the -mode and the -mode, respectively.
(iv) Kinematic fit: The uncertainty associated to the kinematic fit is estimated by comparing the efficiencies with and without the helix parameter correction [38].
(v) Signal yields • Mass window: The systematic uncertainties related to each individual mass window requirement for background rejection are estimated by varying the size of each mass window by one standard deviation of the corresponding mass resolution.The larger variation in the BF for each mass requirement is considered as the related systematic uncertainty.For the 2D mass window of Σ + Σ− signal selection, the uncertainty is estimated by the control sample / → Σ + Σ− .The differences in the selection efficiencies between data and MC simulation from the control sample are taken as the corresponding systematic uncertainties.• Fitting range: The uncertainty due to the fitting range is estimated by randomly changing the range of the fits (for the -mode, the lower and the upper side in the range [630, 670] MeV/ 2 and [890, 930] MeV/ 2 , respectively; for the -mode, the upper side in the range [1080, 1120] MeV/ 2 ) and performing the fit 800 times, using a recalculated signal efficiency, as described in Ref. [39].The variance of the re-measured BFs is taken as the corresponding systematic uncertainty, respectively.• Signal shape: The uncertainty due to the signal shape is estimated by replacing the MCsimulated shape convolved with a Gaussian function with only the MC-simulated shape.
The difference in the measured BF is taken as the systematic uncertainty.• Non-peaking background: To estimate the uncertainty of the non-peaking background in the fit, we perform alternative fits by replacing the second-order Chebychev polynomial with a third-order polynomial function or the Argus function with the third-order Chebychev polynomial for the data.The difference in the measured BF is taken as the systematic uncertainty.• Δ-related background: The uncertainty due to the possible peaking background (3686) → Δ + Δ− () is estimated by including in the fit this channel, assuming ℬ(Δ + →  0 ) = 66.6% and ℬ((3686) → Δ + Δ − ()) equal to ℬ((3686) → Σ + Σ − ()), and the largest differences between these results and the nominal ones are taken as the systematic uncertainty, separately for the -mode and the -mode.
• Sideband background: The systematic uncertainty associated to the 2D sideband background is estimated by increasing and reducing the sideband regions by one standard deviation of the mass resolution, or change the normalized value of the background contribution by one standard deviation of the statistical uncertainty.The difference in the measured BF is taken as the systematic uncertainty.
• Interference: The effect due to the interference between the (3686) decay and the continuum production is estimated with the method described in Ref. [40].The cross section for a certain exclusive final state  can be written as where √  is the center-of-mass energy,    and    are the amplitudes of the continuum contribution and of the (3686) → Σ + Σ−  decay, respectively, and Φ is the relative phase between the two amplitudes [41].The systematic uncertainty due to the interference effect is assigned as the difference of    between the constructive (Φ = /2) and destructive (Φ = 3/2) interference, where    represents the ratio of cross section from the interference term with respect to the cross section of (3686) → Σ + Σ− .Since the prior distribution of the interference angle is assumed to be uniform, this difference divided by √ 12, namely 3.1%, is assigned as the corresponding systematic uncertainty for the -mode.As previously mentioned, this uncertainty for the -mode is found negligible.
(vi) MC simulation: The systematic uncertainty associated to the MC simulation is related to the use of the BODY3 generator and to the differences in the polar angle distributions of each final state particle between data and MC simulation.The first contribution is estimated by varying the bin size of the input Dalitz plot by ±25%, shifting the 1D sideband regions to the left or to the right by one standard deviation of the mass resolution, and changing in the BODY3 generator the background level in the input Dalitz plot by ±1, where  denotes the statistical uncertainty of the background level determined from the fit result.Combining the results from the three sources, the largest change to the nominal detection efficiency is taken as the systematic uncertainty.The second contribution is estimated by rescaling the corresponding distribution from the MC simulation to match the data.The difference in efficiency before and after the rescaling is taken as the systematic uncertainty.
(vii) Quoted BFs: The uncertainties for the quoted BFs are taken from PDG [29].

FIG. 1 .
FIG. 1. Distributions of  p 0 versus   0 of the accepted candidates for (A) (3686) → Σ + Σ−  and (B) (3686) → Σ + Σ− , where the red solid rectangle denotes the signal region, the pink dashed rectangles the sideband I region, and the green long-dashed rectangles the sideband II region.
Figures 2 and 3 show the fit results.

FIG. 2 .
FIG.2.The  ( +  −  0 ) distribution of the accepted events in the 2D signal region.The black points with uncertainties are data, the gray solid curve is the fit result, the orange shaded histogram represents the signal, the green histogram denotes the scaled 2D sideband contribution and the blue dotted curve represents the remaining smooth background.

FIG. 3 .
FIG.3.The  ( +  − ) distribution of the accepted events in the 2D signal region.The black points with uncertainties are data, the gray solid curve is the fit result, the orange shaded histogram represents the signal, the green histogram denotes the scaled 2D sideband contribution and the blue dotted curve represents the remaining smooth background.

TABLE I .
Mass selection windows for each mode.