Observation of $\psi(3686)\to p\bar{p}\phi$

Using a data sample of $4.48\times10^{8}$ $\psip$ events collected with the BESIII detector, we present a first observation of $\psi(3686)\to p\bar{p}\phi$, and we measure its branching fraction to be $[6.06\pm0.38 ($stat.$) \pm 0.48 ($syst.$)]\times10^{-6}$. In contrast to the earlier discovery of a threshold enhancement in the $p\bar{p}$-mass spectrum of the channel $J/\psi\to\gamma p\bar p$, denoted as $X(p\bar{p})$, we do not find a similar enhancement in $\psi(3686)\to p\bar{p}\phi$. An upper limit of $1.82\times10^{-7}$ at the $90\%$ confidence level on the branching fraction of $\psi(3686)\to X(p\bar{p})\phi\to p\bar{p}\phi$ is obtained.

Using a data sample of 4.48 × 10 8 ψ(3686) events collected with the BESIII detector, we present a first observation of ψ(3686) → ppφ, and we measure its branching fraction to be [6.06 ± 0.38(stat.) ± 0.48(syst.)] × 10 −6 . In contrast to the earlier discovery of a threshold enhancement in the pp-mass spectrum of the channel J/ψ → γpp, denoted as X(pp), we do not find a similar enhancement in ψ(3686) → ppφ. An upper limit of 1.82 × 10 −7 at the 90% confidence level on the branching fraction of ψ(3686) → X(pp)φ → ppφ is obtained.

I. INTRODUCTION
An intriguing enhancement near the pp-mass threshold, referred to as the X(pp), was discovered by BES in the channel J/ψ → γpp [1] and subsequently con-firmed by CLEO [2] and BESIII [3]. A more recent partial-wave amplitude analysis of J/ψ → γpp [4] supports the existence of the structure and concludes to a spin-parity assignment of J P C = 0 −+ . There is no experimental evidence of such an enhancement in radia-tive Υ(1S) → γpp [5] decay nor in the J/ψ → ωpp decay [6]. It is tempting to associate this enhancement with the X(1835), a resonance that was recently confirmed by BESIII [7] after it was first observed in J/ψ → γπ + π − η ′ decay [8]. Whether or not the pp-mass threshold enhancement and the X(1835) are related to the same source still needs further study. As a result, lots of theoretical speculations have been proposed to interpret the nature of this structure, including the quasibound nuclear baryonium [9,10], a multiquark resonance [11] or an effect caused by final-state interaction (FSI) [12,13] near the proton-antiproton production threshold.
Most recently, BESIII reported the study of J/ψ → ppφ [14] and no evidence of a near-threshold enhancement in the pp-mass spectrum was found. Moreover, no significant signatures of resonances in the pφ orpφ mass spectra were observed. For the decay of ψ(3686) → ppφ, BES reported an upper limit on the branching fraction B(ψ(3686) → ppφ) of 2.6 × 10 −5 at the 90% confidence level (C.L.) [15]. The latest measurement came from CLEO [16], who reported an upper limit on the branching fraction B(ψ(3686) → ppφ) of 2.4 × 10 −5 at the 90% C.L.. These experimental observations, together with similar results found in different decays, give rise to a discussion on the nature of the threshold effect and stimulate theoretical developments.
In this work, we report on the data analysis of the charmonium decay ψ(3686) → ppφ. The data have been obtained with the BESIII detector at the BEPCII storage ring at which a total of (4.481 ± 0.029) × 10 8 ψ(3686) events [17] were produced in electron-positron annihilations. The aim of this work is to search for a near-threshold enhancement in the pp-mass spectrum and to search for pφ(pφ) resonances that might hint to the existence of pentaquarks with hidden strangeness. Moreover, we measured the branching fraction of the process ψ(3686) → ppφ which allows us to inspect the '12% rule' proposed in 1975 [18]. The rule is based on perturbative quantum chromodynamics (QCD) calculations, in which the ratio of the branching fractions of ψ(3686) and J/ψ into the same final hadronic state is given by

II. DETECTOR AND MONTE CARLO SIMULATION
The BESIII detector is a magnetic spectrometer [19] located at the Beijing Electron Positron Collider (BEPCII) [20]. The cylindrical core of the BESIII 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 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 time resolution of the TOF barrel part is 68 ps, while that of the end cap part is 110 ps.
Simulated samples produced with the geant4based [21] Monte Carlo (MC) package which includes the geometric description of the BESIII detector and the detector response, are used to determine the detection efficiency and to estimate the backgrounds. The simulation includes the beam energy spread and the initialstate radiation (ISR) in the e + e − annihilations modelled with the generator kkmc [22]. The inclusive MC sample consists of the production of the ψ(3686) resonance, and the continuum processes incorporated in kkmc. The known decay modes are modelled with evtgen [23] using branching fractions taken from the Particle Data Group (PDG) [24], and the remaining unknown decays from the charmonium states with lundcharm [25]. The final-state radiation (FSR) from charged final-state particles is incorporated with the photos package [26]. The background is studied using a sample of 5.06 × 10 8 inclusive ψ(3686) MC events. The analysis is performed in the framework of the BESIII offline software system (BOSS) [27] incorporating the detector calibration, event reconstruction and data storage.

A. Event selection and background analysis
The ψ(3686) → ppφ reaction is identified with the φ subsequently decaying into K + K − resulting in a final state of four charged tracks, namely ppK + K − . The charged tracks must have been detected in the active region of the MDC, corresponding to | cos θ| < 0.93, where θ is the polar angle of the charged track with respect to the beam direction. Moreover, the tracks are required to pass within ±10 cm of the interaction point in the beam direction and within ±1 cm in the plane perpendicular to the beam. Two of the charged tracks are identified as a proton and an antiproton by using combined TOF and dE/dx information. To improve the detection efficiency, the events with at least one K + (K − ) are selected for further analysis. Thus, the candidate events are required to have three or four charged tracks. A one-constraint (1C) kinematic fit is subsequently performed under the hypothesis of ψ(3686) → ppK + K − , where K + or K − is treated as a missing particle with the nominal mass of a kaon. For the events with both kaons detected, two 1C kinematic fits are performed assuming a missing K + or K − . The one with the least χ 2 1C is retained. To suppress background events, the χ 2 1C is required to be less than 10.
The potential backgrounds are investigated using the inclusive ψ(3686) MC sample. Besides the irreducible backgrounds from the non-resonant decay ψ(3686) → ppK + K − , the reducible backgrounds are dominated by the processes involving Λ(Λ) intermediate states. To suppress the above backgrounds, all other charged tracks except for the selected proton, antiproton and kaon candidates are assumed to be pions, and events are excluded if any combination of pπ − orpπ + has an invariant mass lying in the range There are also some background events found originating from the process ψ(3686) →pK + Λ(1520) + c.c. with Λ(1520) → pK. A MC sample is generated to describe its shape, and the number of background events of ψ(3686) →pK + Λ(1520) is expected to be 40 ± 21, which is estimated by a fit to the measured pK − invariant-mass spectrum. The signal shape of the Λ(1520) → pK is modeled with a Breit Wigner (BW) function, and the background is described with a second-order Chebychev polynomial function. Only the background from the continuum process e + e − → ppφ was found to have a peaking structure underneath the φ-signal region. This contribution from this background is studied using the offresonance samples taken at √ s = 3.773 GeV, and its absolute magnitude is determined according to the formula is the number of events which remained in the off-resonance samples after applying the same event selections that are used to identify ψ(3686) → ppφ. L, σ and ε refer to the integrated luminosities (L ψ(3686) = 668.55 pb −1 [17], L 3773 = 2931.8 pb −1 [28], the cross sections and the detection efficiencies of the data samples taken at the two corresponding center-of-mass energies, respectively. Figure 1 shows the K + K − invariant-mass spectrum after applying all the selection criteria mentioned above. Note that a clear signal corresponding to the decay φ → K + K − is visible in the spectrum. Figure 2 shows the Dalitz plot of ψ(3686) → ppφ for the events with a K + K − invariant-mass that falls within the φ-mass region (1.005 GeV/c 2 < M K + K − < 1.035 GeV/c 2 ). The data show no evident resonance structures. Figure 3 shows its projections on the pφ andpφ invariant-mass distributions. These distributions show that the data are well described by a phase-space distribution of the signal channel together with the continuum background and non-peaking background.  ψ(3686) → ppφ, is generated according to a phase-space assumption. The background contribution from the continuum process e + e − → ppφ is obtained as discussed above and its shape and yield are fixed in the fit. The other background events are parameterized by a modified ARGUS function [29]. The parameters of the Gaussian function and the ARGUS function are left free in the fit. The fit, shown in Fig. 1, yields N obs = 753 ± 47 signal events. The statistical significance is found to be 21 σ, which is determined from the change in −2 ln L in the fits of mass spectrum with and without assuming the presence of a signal while considering the change in degrees of freedom of the fits.  with, where N obs is the number of the observed signal events which comes from the fit. N ψ(3686) is the total number of ψ(3686) events. The branching fraction of φ → K + K − , B(φ → K + K − ) = (49.2 ± 0.5)%, is taken from the PDG [24]. ε is the detection efficiency. To obtain a reliable detection efficiency, the MC sample of ψ(3686) → ppφ, distributed according to a phase-space assumption, is weighted to match the distribution of the backgroundsubtracted data with the mass distribution of pp, and the average detection efficiency is determined to be 56.4%. The branching fraction, B(ψ(3686) → ppφ), is measured to be (6.06 ± 0.38 ± 0.48) × 10 −6 , where the uncertainties are the statistical and systematic uncertainty, respectively. The systematic uncertainties will be discussed in detail in the following section.

C. Systematic uncertainties
The systematic uncertainties that affect the branchingfraction measurement can be divided into two categories. The first category is given by the uncertainties in the track reconstruction, the particle identification (PID), 1C kinematic fit, and Λ/Λ veto efficiency. The other category comprises the uncertainties which originate from the fit of the mass spectrum, the weighting procedure, the cited branching fraction of the decay of the intermediate state, and the total number of ψ(3686) events.
The difference in the efficiencies of the track reconstruction for p/p between MC and data is studied using a clean sample of J/ψ → ppπ + π − and found to be less than 1.0% per track. For the K ± , the systematic uncertainty is studied using a clean control sample of J/ψ → K 0 S K ± π ∓ . 1.0% per tracking is taken as the systematic uncertainty for the tracking efficiency [30].
The PID efficiency of p/p is also studied from the same data sample of J/ψ → ppπ + π − . The results indicate that the p/p PID efficiency for data agrees with the MC simulation within 1%. The PID efficiency for the kaon is measured in the clean channel J/ψ → K + K − η. It is found that the difference between the PID efficiency of data and MC is less than 1% for each kaon. In this analysis, three charged tracks are required to be identified as a proton, an anti-proton and a kaon. Hence, 3% is taken as the systematic uncertainty associated with the PID.
With a clean control sample of ψ(3686) → pK −Λ + c.c., the systematic uncertainty of the 1C kinematic fit is estimated to be 3.4% by calculating the difference of ratio of signal yields with χ 2 1C cut and without 1C kinematic fit between MC simulation and data.
The φ-signal yields are obtained by fitting the K + K − invariant-mass spectrum. Systematic uncertainties related to the fit have been estimated by using different signal and background shapes, alternative fit ranges, and by taking into consideration an additional resonant structure. To estimate the uncertainty from the modeling of the φ-signal shape, an alternative fit with an acceptancecorrected BW function to describe the φ-signal has been performed. To estimate the uncertainty due to the back- In the K + K − invariant-mass distribution, we observed a small bump around 1 GeV/c 2 . Although this structure might be due to statistical fluctuations, we considered the possibility of an additional resonance. We, therefore, fitted the distribution with an extra BW function convolved with a Gaussian function. The change of signal yield in the different fit is taken as the corresponding systematic uncertainty. The quadratic sum of the four individual uncertainties is taken as the systematic uncertainty related with the mass spectrum fit, and it is found to be 5.5%.
To obtain a reliable detection efficiency, the MC sample modeled using a phase-space distribution is weighted to match the distribution of the background-subtracted data. To consider the effect on the statistical fluctuations of the signal yield in the data, a set of toy-MC samples are used to estimate the detection efficiencies. With the reweighting, a maximum deviation in detection efficiencies of 1.0% is found and quoted as the corresponding systematic uncertainty.
The branching fraction uncertainty of the intermediate decay φ → K + K − , 1.0%, is taken from the PDG and the uncertainty of the number of ψ(3686) events is 0.6% [17].
In Table I, a summary is shown of all contributions to the systematic uncertainties on the branching fraction measurements. The total systematic uncertainty is given by the quadratic sum of the individual contributions, assuming all sources to be independent.   Figure 4 depicts the pp invariant-mass distribution for the events with a K + K − invariant mass that falls within the φ mass region (1.005 GeV/c 2 < M K + K − < 1.035 GeV/c 2 ), where no evident enhancement near the pp-mass threshold is visible. It is found that the events from the phase-space process together with other background components provide a good description of the data, which is shown in Fig. 4. Therefore, an upper limit for the X(pp) production rate can be measured. For J/ψ → ppφ, we divide pp invariant-mass spectrum into 9 bins in the region of [1.876, 2.056] GeV/c 2 . With the same procedure as described above, the number of the φ events in each bin can be obtained by fitting to the corresponding K + K − -mass spectrum. Subsequently, the nonφ-background-subtracted M pp distribution is obtained as shown in Fig. 5, where the errors are statistical only, and m p is the nominal mass of proton [24]. The spin (J) and parity (P ) of X(pp) have been determined by an amplitude analysis of J/ψ → γpp decay and resulted in J P C = 0 −+ [4]. In our analysis, we parametrize the X(pp) signal by an efficiency-weighted S-wave BW function, where M is the pp invariant mass, the parameter f F SI accounts for the effect of the FSI, q is the momentum of the proton in the pp rest frame, κ is the momentum of φ in the ψ(3686) rest frame, L = 0 is the relative orbital angular momentum of the pp system, M 0 and Γ 0 are the mass and width of X(pp), respectively, which are fixed to those in Ref. [4]. ε rec (M ) is the mass-dependent detection efficiency which is obtained from MC simulations of ψ(3686) → X(pp)φ → ppφ. We ignore possible interference effects of the X(pp) resonance with non-resonant background contributions.
To determine the upper limit on the size of the pp enhancement, a series of binned least-χ 2 fits are performed to the background-subtracted pp-mass spectrum with the expected signal. Fit-related uncertainties are included by considering the following three aspects: (a) the X(pp) signal is described by excluding the FSI factor with f F SI = 1 or taking into account the Jülich FSI value as described in Ref. [13]; (b) the non-resonant background is represented by the shape obtained from the ψ(3686) → ppφ MC simulation or parameterized by a function of f (δ) = N (δ 1/2 + a 1 δ 3/2 + a 2 δ 5/2 ) (δ = M pp − 2m p , a 1 and a 2 are free parameters); and (c) the fit is performed in the range of  limit on the branching fraction is determined by, where N UL is the maximum number of X(pp) events. To be conservative, the multiplicative uncertainties listed in Table I are considered by convoluting the normalised χ 2 distribution with a Gaussian function. The detection efficiency, ε, is obtained from MC simulations, and is determined to be 58.9%. The upper limit on the branching fraction of ψ(3686) → X(pp)φ → ppφ at the 90% C.L. is calculated to be 1.82 × 10 −7 .