Cross section measurement of $e^+e^- \to p\bar{p}\eta$ and $e^+e^- \to p\bar{p}\omega$ at center-of-mass energies between 3.773 GeV and 4.6 GeV

Based on $14.7~\textrm{fb}^{-1}$ of $e^+e^-$ annihilation data collected with the BESIII detector at the BEPCII collider at 17 different center-of-mass energies between $3.7730~\textrm{GeV}$ and $4.5995~\textrm{GeV}$, Born cross sections of the two processes $e^+e^- \to p\bar{p}\eta$ and $e^+e^- \to p\bar{p}\omega$ are measured for the first time. No indication of resonant production through a vector state $V$ is observed, and upper limits on the Born cross sections of $e^+e^- \to V \to p\bar{p}\eta$ and $e^+e^- \to V \to p\bar{p}\omega$ at the $90\%$ confidence level are calculated for a large parameter space in resonance masses and widths. For the current world average parameters of the $\psi(4230)$ of $m=4.2187~\textrm{GeV}/c^{2}$ and $\Gamma=44~\textrm{MeV}$, we find upper limits on resonant production of the $p\bar{p}\eta$ and $p\bar{p}\omega$ final states of $7.5~\textrm{pb}$ and $10.4~\textrm{pb}$ at the $90\%$ CL, respectively.

indication of resonant production through a vector state V is observed, and upper limits on the Born cross sections of e + e − → V → ppη and e + e − → V → ppω at the 90 % confidence level are calculated for a large parameter space in resonance masses and widths. For the current world average parameters of the ψ(4230) of m = 4.2187 GeV/c 2 and Γ = 44 MeV, we find upper limits on resonant production of the ppη and ppω final states of 7.5 pb and 10.4 pb at the 90% CL, respectively. INTRODUCTION In recent years, several unexpected states have been discovered in the charmonium sector. Notable examples are the χ c1 (3872) discovered by Belle [1], the charged charmonium-like Z c (3900) discovered by BESIII [2] and the vector states ψ(4220) and ψ(4390), originally discovered by BaBar [3] as a single broad resonance named Y (4260) in the e + e − → γ ISR π + π − J/ψ process, and later found to be two distinct states by BESIII [4]. Charmonium-like vector states similar to the ψ(4220) and ψ(4390) were observed by BESIII in the processes e + e − → π + π − h c [5], π + π − ψ(3686) [6], D 0 D * − π + [7], ωχ c0 [8] and ηJ/ψ [9][10][11]. However, so far no observations have been made of decays into light mesons or baryons for either the ψ(4220) or the ψ(4390). Cross sections have been determined for the processes e + e − → ppπ 0 [12], φφφ, φφω [13], pK 0 Sn K − [14], K 0 S K ± π ∓ [15], K 0 S K ± π ∓ π 0 and K 0 S K ± π ∓ η [16] without significant detection of resonant production. Channels that include a proton anti-proton pair are especially interesting, since the partial width of decays of the type V → pph, where V is a vector-state in the charmonium region and h is a light meson, can be related to the production cross section pp → V h [17]. In light of the upcoming PANDA experiment at the FAIR facility [18], it is important to obtain the production cross sections of potentially exotic resonances in the charmonium sector.
Multiple theoretical explanations have been offered concerning the nature of the ψ(4220) and ψ(4390) states. Possible interpretations include D 1D molecules, hybrid charmonia or baryonium states, and their compatibility with experimental data has recently been discussed in detail in Ref. [19]. Additional information is needed from experiments in order to discriminate between the different hypotheses. Thus, the search for new decay modes of the ψ(4220) and ψ(4390) is important.
In this work, we report measurements of the Born cross sections of the processes e + e − → ppη and e + e − → ppω for data collected at 17 different center-of-mass energies between √ s = 3.7730 GeV and 4.5995 GeV with the BESIII detector. The center-of-mass energy dependence of both cross sections is investigated for possible resonant contributions V → ppη and V → ppω.

II. BESIII DETECTOR AND MONTE CARLO SIMULATIONS
The BESIII detector is a magnetic spectrometer [20] located at the Beijing Electron Positron Col-lider (BEPCII) [21]. 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 the 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. The end cap TOF system was upgraded in 2015 with multi-gap resistive plate chamber technology, providing a time resolution of 60 ps [22]. This improvement affects data at 11 of the 17 center-of-mass energy points.
A Monte Carlo (MC) simulation of the full BESIII detector, based on geant4 [23], is used to optimize selection requirements, determine the product of detector acceptance and reconstruction efficiency and study and estimate possible background contributions. These simulations also account for the observed beam energy spread.
In addition, an inclusive MC sample at a center-ofmass energy of √ s = 4.1784 GeV, corresponding to the dataset with the largest integrated luminosity (see Table I), is used to study potential background contributions. This sample includes open charm processes, ISR production of vector charmonium(-like) states and continuum qq (where q is a u, d, s quark) processes. Known decay modes are modeled with evtgen [26] using branching fractions taken from the PDG [25], whereas unknown processes are modeled by the lundcharm model [27]. Final state radiation (FSR) from charged final state particles is incorporated with the photos package [28]. The inclusive MC sample at √ s = 4.1784 GeV corresponds to 40 times the luminosity for the data at this center-of-mass energy.

III. EVENT SELECTION
Two different final states are studied in this work, namely ppγγ (for e + e − → ppη with η → γγ) and ppπ + π − γγ (for e + e − → ppη with η → π + π − π 0 and e + e − → ppω with ω → π + π − π 0 , both with a subsequent π 0 → γγ decay). In the analysis, the invariant mass distributions m(γγ) and m(π + π − π 0 ) will be used to identify and quantify η and ω contributions. The polar angle θ of each charged track detected in the MDC has to satisfy | cos θ| < 0.93, and its point of closest approach to the interaction point must be within ±10 cm in the beam direction and within 1 cm in the plane perpendicular to the beam direction. For particle identification (PID), the TOF information and the specific energy deposit dE/dx in the MDC are combined to calculate a likelihood P (h) for the particle hypotheses h = π, K, p. The particle type with the largest likelihood is assigned to each track. In addition, we require a minimum likelihood of P (h) > 10 −5 to suppress background.
For photons, a minimum energy deposit in the calorimeter of 25 MeV in the barrel region (| cos θ| < 0.80) or of 50 MeV in the end-cap regions (0.86 < | cos θ| < 0.92) is required. In addition, the time information from the shower in the calorimeter relative to the event start time has to be less than 700 ns. Showers within 10 • of the impact point of any charged track are discarded.
Only events containing exactly one good proton and one good anti-proton candidate (and exactly one good π + and one good π − candidate in case of the ppπ + π − γγ final state) and at least two good photon candidates are retained. A four-(five-)constraint kinematic fit is performed to the ppγγ (ppπ + π − π 0 with π 0 → γγ) hypothesis requiring four-momentum conservation between initial and final states and, if relevant, an additional mass constraint for the π 0 → γγ decay. If more than one γγ combination in an event satisfies the above requirements, only the combination with the lowest kinematic fit χ 2 is kept for further analysis. The resulting invariant mass spectra for the decays η → γγ, η → π + π − π 0 and ω → π + π − π 0 are displayed for the data at √ s = 4.1784 GeV in Fig. 1. According to the inclusive MC sample, backgrounds containing tracks misidentified as (anti-)proton candidates are very rare and can be neglected. In the case of the ppγγ final state, the main background channels are the processes e + e − → pp , ppπ 0 , ppω with a subsequent ω → π 0 γ decay and e + e − → γ ISR J/ψ with a subsequent J/ψ → pp decay. For the ppπ + π − γγ final state, the main background channels are e + e − → ∆∆π, e + e − → p∆ρ, e + e − → p∆ππ and e + e − → ppρπ (the various possible charges of the ∆,∆, ρ and π and overall charge conjugation are taken into account), which can all lead to the signal final state. No peaking backgrounds in the invariant mass distributions m(γγ) and m(π + π − π 0 ) are found.
The selection criteria with respect to the kinematic fit are optimized according to s √ s+b , where s and b are the number of signal and background events in the inclusive MC sample, that has been scaled to data, after requiring χ 2 < χ 2 sel . In case of the ppγγ final state, χ 2 < 30 is found as an optimal selection condition, whereas for the other two reactions the background is dominated by other processes leading to the ppπ + π − π 0 system so that no requirement is made on the χ 2 of the kinematic fit. This procedure is repeated with s and b being signal and background contributions determined directly from data using the η (ω) signal-and sidebands. The resulting selection condition agrees with the one determined from the inclusive MC sample within statistical uncertainties. As both the data and the inclusive MC samples at the other center-of-mass energies are significantly smaller than those at √ s = 4.1784 GeV, we do not repeat this optimization, and the requirement χ 2 < 30 for the ppγγ final state is applied for all center-of-mass energies.
The number of signal events is determined from a fit to the invariant mass spectra (see Fig. 1). In the fit, the signal is described by a shape determined from signal MC simulations convolved with a Gaussian function in order to account for a possible underestimation of the mass resolution in MC simulation. The background is described by first-and second-order polynomial functions in the η → γγ and η(ω) → π + π − π 0 channels, respectively. In the first step, we perform a global binned maximum likelihood fit to the sum of the data at all center-of-mass energies, determining a width of the Gaussian smearing functions as well as a background shape for each of the three channels. A binned maximum likelihood fit is then performed to each dataset and channel individually, with only the signal and background yields as free parameters. The signal and background shape are fixed to the result of the global fit in order to obtain reliable results from all sample sizes. The number of signal events in channel i at center-of-mass energy √ s is then defined as N i sig (s) = n sr −n bg,sr , where n sr is the number of events in the signal region (defined as the symmetric region around the nominal η (ω) meson mass containing 95% of all signal events according to the signal shape) and n bg,sr is the number of background events in the signal region according to the polynomial background function (see Fig. 1).

IV. EFFICIENCY DETERMINATION
We define the efficiency, i (s), according to where N i acc (s) is the number of reconstructed signal events and N i gen (s) is the total size of the signal MC sample in channel i for a center-of-mass energy √ s. If the efficiency is not constant over the full n-particle phase-space, Eq. 1 only holds if the signal MC sample properly reflects data in all relevant coordinates x = {p p , θ p , φ p , pp, θp, φp, ...}. Since the data distribution is a priori unknown, we will perform a partial wave analysis of the data in order to re-weight our MC sample.
The isobar model [29] is used in the partial wave analysis by decomposing the full e + e − → γ * → ppη and e + e − → γ * → ppω processes into a sequence of two-body decays. Each two-body decay is described in the helicity formalism [30]. The η meson is treated as a stable particle in the amplitude analysis, whereas the three-body decay of the ω meson is described by a three-body am-plitude according to Ref. [31]. Blatt-Weisskopf barrier factors [30] are used for both the production γ * → ad and the two-body decay a → bc according to Ref. [31] . The dynamical part of the amplitude is described by relativistic Breit-Wigner amplitudes. Only the line-shape of the η and ω mesons is not described in this way. Here we employ the signal MC simulations that are used for normalization in the partial wave analysis for the lineshapes. For the two e + e − → ppη channels, we include as intermediate pp resonances the J/ψ (using a Voigt distribution) and possible contributions of a J P C = 1 −− and a J P C = 3 −− resonance. We also include one J P = 1 2 −(+) and one J P = 3 2 +(−) resonance contribution for the ηp(p) system. For of the e + e − → ppω channel, only intermediate states that decay to the pp system are included. For each of the three possible quantum numbers J P C = 0 −+ , 0 ++ , 2 ++ we include a phase-space contribution that is constant as a function of m(pp) as well as two resonant contributions. All masses and widths of intermediate states are free parameters in the fit, apart from the J/ψ contribution, where mass and width are fixed to the values given in the PDG [25]. Note that the aim of this partial wave analysis is only to describe the data accurately enough to enable an accurate determination of the efficiency.
The partial wave analysis is performed as an unbinned maximum likelihood fit using the software package PAW-IAN [32]. Details on likelihood construction in PAW-IAN can be found in Refs. [31][32][33]. The remaining background events underneath the η (ω) peaks are accounted for in the partial wave analysis by adding the sidebands defined in Fig. 1 to the likelihood with negative weights in such a way that the combined sideband weight is equal to the integral of the background function in the signal region.
The data for the two decay channels η → γγ and η → π + π − π 0 are fitted simultaneously with all amplitudes fully constrained between the two channels apart from an overall scaling factor. The results of the partial wave analysis for the different channels are displayed in Fig. 2 for the high statistics data at a center-of-mass energy of √ s = 4.1784 GeV. In general, the fits have a similar quality for all center-of-mass energies. For each energy point, we obtain event weights, w( x), from the partial wave analysis as a function of the coordinates in the n-particle phase-space, and they are used to determine the efficiency i (s) as The efficiencies obtained in this way are summarized in Tables I and II.  → γγ (a -b) and η → π + π − π 0 decays (c -d) and the e + e − → ppω process (e -f) for the data at a center-of-mass energy of √ s = 4.1784 GeV. The left column shows the invariant mass of the pp system, the right column the invariant mass of the pη (pω) system. Black points correspond to data, full (red) lines show the result of the amplitude analysis.

V. DETERMINATION OF BORN CROSS SECTIONS
The Born cross section of the process e + e − → ppη (e + e − → ppω) is given by where N (s) is the number of signal events observed in the data sample at center-of-mass energy √ s, L(s) is the corresponding integrated luminosity determined us-ing Bhabha scattering [34], δ r (s) and 1 |1−Π| 2 are corrections accounting for initial state radiation and vacuum polarization, (s) is the efficiency and B is the product of branching ratios involved in the decay. The correction 1 |1−Π| 2 is calculated with the alphaQED software package [35] with an accuracy of 0.5%. As initial state radiation depends on the shape of the cross section and can, in general, have an effect on the efficiency, it is treated in an iterative procedure starting from a flat energy dependence of the Born cross section σ B (s). We consider two successive iterations converged if κ i /κ i−1 = 1 within statistical uncertainties, where κ(s) = (s) · (1 + δ r (s)) is the product of the ef-ficiency and a radiative correction factor 1 + δ r (s) obtained from the ConExc MC generator. The product of branching fractions is given by B = Br(η → γγ) for the e + e − → ppη(→ γγ) channel, B = Br(η → π + π − π 0 ) · Br(π 0 → γγ) for the e + e − → ppη(→ π + π − π 0 ) channel and B = Br(ω → π + π − π 0 ) · Br(π 0 → γγ) for the e + e − → ppω(→ π + π − π 0 ) channel. A combined Born cross section σ B is determined for the two different η decay modes by using a weighted least squares method [36]. The resulting cross sections are displayed in Fig. 3 and all necessary values for their calculation are summarized in Tables I and II. Table I. Summary of the Born cross sections σB of the process e + e − → ppη for the datasets at different center-of-mass energies √ s, integrated luminosity L, radiative corrections 1 + δr, vacuum polarization correction 1 |1−Π| 2 , number of observed events Ni, efficiency i, and cross section σi in the two channels (1) η → γγ and (2)

VI. SYSTEMATIC UNCERTAINTIES
Various sources of systematic uncertainties contributing to the measurement of the e + e − → ppη and e + e − → ppω Born cross sections have been considered.
The uncertainty of the integrated luminosity determined using Bhabha scattering is 1% [34]. The systematic uncertainty of the tracking efficiency has been determined using a J/ψ → ppπ + π − control sample in Ref. [37] as 1% per track. Similarly, systematic uncertainties of photon detection efficiencies have been studied using a J/ψ → ρπ control sample [38] and were found to be 1% per photon. For PID efficiency, a systematic uncertainty of 1% per proton and 1% per pion are taken from Ref. [12,39]. For multiple particles, each of the track-finding, PID, and photon efficiency uncertainties are added linearly [12,[37][38][39]. Uncertainties on the branching fractions are taken from the PDG [25]. With regard to the kinematic fit, where a selection condition of χ 2 < 30 is applied in case of the e + e − → ppη(→ γγ) mode, the selection condition is varied between χ 2 < 5 and χ 2 < 55 in steps of δχ 2 = 5 and the resulting Born cross section is determined and compared with the nominal value R = σstep σnom . We take the standard deviation of a weighted sample of the ratio R as the systematic uncertainty due to potential differences in the χ 2 distributions between data and MC simulation. Here, 1/δR is used as the weight, where δR is the uncertainty taking into account the sizable correlation between the event samples. The nominal symmetric signal region containing 95% of the total signal is altered to a set of both smaller and larger signal regions and we determine the resulting Born cross sections. As outlined above, we take the standard deviation of a sample of ratios R weighted by the inverse of the statistical uncertainty as the systematic uncertainty resulting from the choice of signal region. For the background description, the polynomial shapes were increased by one order from the nominal first order polynomial used for the η → γγ invariant mass spectrum, and the second order polynomial in the case of the η → π + π − π 0 and ω → π + π − π 0 invariant mass spectra. The fits are then repeated and the difference to the nominal results is taken as a systematic uncertainty. For the radiative correction factor, we performed five additional iterations and found no difference beyond the statistical uncertainty. This contribution to the systematic uncertainty is therefore neglected.
The systematic uncertainties are summarized in Table III for the data at a center-of-mass energy of √ s = 4.1784 GeV. The total systematic uncertainty is obtained by adding each contribution in quadrature. Correlated systematic uncertainties in the two e + e − → ppη channels are accounted for in the calculation of the combined Born cross section following Ref. [36].

VII. SEARCH FOR RESONANT CONTRIBUTIONS
The final Born cross sections for the e + e − → ppη and e + e − → ppω processes are displayed in Fig. 3. In order to search for possible e + e − → V → ppη (e + e − → V → ppω) resonant contributions, we perform two different fits. In the first fit, only a non-resonant contribution of the type defined in Ref. [12] is used. The second fit includes a single Breit-Wigner amplitude of the form that is coherently added to the non-resonant term. Unbinned maximum likelihood fits are performed where the likelihood L(x; Θ) given the data x and the fit parameters Θ is defined as the product L(x; Θ) = i,j L ij (Θ), where L ij is a set of likelihood functions, one each for each dataset i and decay mode j. These likelihood functions are transformed such that they only depend on the expected number of signal events N ≡ N ij (Θ) which can be calculated for each dataset according to Eq. 3. The likelihood L ij (N ) is then obtained from data via a likelihood scan of the number of signal events in the invariant mass distributions of the meson decay systems. These likelihood scans are parameterized by asymmetric Gaussian distributions. Incorporating the systematic uncertainties of dataset i and channel j, the likelihood is In the fit, all systematic uncertainties apart from the one on the branching ratio of the meson decays are considered uncorrelated between the different c.m. energies. While a correlation of a systematic uncertainty between two c.m. energies can not in general be ruled out, our assumption of a vanishing correlation leads to the most conservative upper limit estimation. We find no evidence for a resonant contribution from the fits and set upper limits at the 90% confidence level. As the resonant contribution is added coherently, the fit finds two ambiguous solutions for constructive and destructive interference by construction. The upper limits are obtained by integrating L(x; Θ) = i,j L ij (Θ) according to where the prior π(Θ) is given by The procedure outlined above is repeated with a step size of 1 MeV for different masses m in the range 4 GeV/c 2 < m < 4.4 GeV/c 2 and widths Γ in the range 40 MeV < Γ < 300 MeV for a potential resonant contribution. The results are shown in Fig. 4. The most stringent upper limits are found for resonant contributions with mass m = 4.389 GeV/c 2 and width Γ = 40 MeV (Γ = 296 MeV) in the ppη (ppω) channel with values of 4.03 pb and 7.88 pb at the 90% CL, respectively. The upper limits for a resonant contribution of the ψ(4230), using current world average values for mass (m = 4.2187 GeV/c 2 ) and width (Γ = 44 MeV) [25] are 7.5 pb and 10.4 pb at the 90% CL, respectively.

VIII. SUMMARY
The processes e + e − → ppη and e + e − → ppω have been studied using 14.7 fb −1 of electron-positron annihilation data at 17 different center-of-mass energies between 3.7730 GeV and 4.5995 GeV. Both processes are clearly identified at all center-of-mass energies and Born cross sections are determined. We find no evidence for a resonant contribution from a fit to the e + e − → ppη and e + e − → ppω Born cross sections, and set upper limits at the 90% confidence level for a wide range of resonance parameters m and Γ. Using the approach outlined in Ref. [17], these upper limits will serve as valuable input for model calculations of the processes pp → V η and pp → V ω for the upcoming PANDA experiment.