Cross section measurements of $e^{+}e^{-}\rightarrow p \bar{p} \pi^{0}$ at center-of-mass energies between 4.008 and 4.600 GeV

Based on $e^+e^-$ annihilation data samples collected with the BESIII detector at the BEPCII collider at 13 center of mass energies from 4.008 to 4.600 GeV, measurements of the Born cross section of $e^{+}e^{-}\rightarrow p \bar{p} \pi^{0}$ are performed. No significant resonant structure is observed in the measured energy-dependent cross section. The upper limit on the Born cross section of $e^{+}e^{-}\rightarrow Y(4260) \rightarrow p \bar{p} \pi^{0}$ at the 90\% C.L. is determined to be 0.01 pb.

e Also at the Novosibirsk State University, Novosibirsk, 630090, Russia

Introduction
The Born cross section of e + e − → ppπ 0 in the vicinity of the ψ(3770) has been measured recently by BESIII [1]. Information on the cross section of e + e − → ppπ 0 at higher energies is however still lacking. The experimental data on the cross section of e + e − → hadrons can be used as an input to calculate the hadronic vacuum polarization via dispersion integrals [2,3,4,5].
The charmonium-like state Y (4260) was first observed in its decay to π + π − J/ψ [6]. So far, there is no evidence of the Y (4260) in the measured open charm decay channels [7,8] and R value scans [9,10,11,12,13,14,15]. Many theoretical models have been proposed to interpret the nature of Y (4260), e.g. as a tetraquark state [16], a D 1 D or D 0 D * hadronic molecule [17], a hybrid charmonium [18,19], or a baryonium state [20]. Searches for new decay modes of the Y (4260) may provide information that can shed light on the nature of Y (4260). In particular, the hybrid model [18] predicts a sizable coupling between the Y (4260) and charmless decays.
In this analysis, we report measurements of the cross section of e + e − → ppπ 0 based on the e + e − annihilation samples collected with the BESIII detector at 13 center of mass energies in the range √ s = 4.008−4.600 GeV as shown in Table 1. Results of the measurements can be used to estimate the cross section of pp → X cc π 0 , which is of high importance for the planned PANDA experiment [21] at FAIR in Darmstadt, Germany.

BESIII detector and Monte-Carlo simulation
The BESIII detector [22] is a magnetic spectrometer operating at BEPCII, a double-ring e + e − collider with center-of-mass energies between 2.0 and 4.6 GeV and a peak luminosity of 10 33 cm −2 s −1 near the ψ(3770) mass. The cylindrical core of the BESIII detector consists of a helium-based main drift chamber (MDC), a plastic scintillator time-of-flight system (TOF), and a CsI(Tl) electromagnetic calorimeter (EMC) that 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 for charged particles and photons is 93% of the 4π solid angle, and the charged-particle momentum resolution is 0.5% for transverse momenta of 1 GeV/c. The energy resolution for showers in the EMC is 2.5 (5%) for 1 GeV photons in the barrel (endcaps) region.
A GEANT4-based [23] Monte Carlo (MC) simulation software package is used to optimize the event selection criteria, estimate backgrounds and determine the detection efficiency. For each energy point, we generate 200,000 signal MC events of e + e − → ppπ 0 uniformly in phase space. Effects of initial state radiation (ISR) are simulated with KKMC [24], where the line shape of the production cross section of e + e − → ppπ 0 is taken from results of the measured cross section iteratively. Effects of final state radiation off charged particles are simulated with PHOTOS [25].
To study possible backgrounds, a MC sample of inclusive Y (4260) decays, equivalent to an integrated luminosity of 825.6 pb −1 , is also generated at √ s = 4.26 GeV. In these simulations, the Y (4260) is allowed to decay generically, with the main known decay channels being generated using EVTGEN [26] with branching fractions set to world average values [27]. The remaining events associated with charmonium decays are generated with LUNDCHARM [28], while the continuum hadronic events are generated with PYTHIA [29]. QED events (e + e − → e + e − , µ + µ − , and γγ) are generated with KKMC [24]. The sources of backgrounds at other energy points are assumed to be similar.

Event selection
The final state in this decay is characterized by two charged tracks and two photons. Two charged tracks with opposite charge are required. Each track is required to have its point of closest approach to the beam axis within 10 cm of the interaction point in the beam direction and within 1 cm of the beam axis in the plane perpendicular to the beam. The polar angle of the track is required to be within the region of | cos θ| < 0.93.
The time-of-flight and the specific energy loss dE/dx of a particle measured in the MDC are combined to calculate particle identification probabilities for pion, kaon, and proton hypotheses. For each track, the particle type yielding the largest probability is assigned. In this analysis, one charged track is required to be identified as a proton and the other one as an anti-proton.
Photon candidates are reconstructed using clusters of energy deposited in the EMC. The energy deposited in nearby TOF counters is included in EMC measurements to improve the reconstruction efficiency and the energy resolution. Photon candidates are selected by requiring a minimum energy deposition of 25 MeV in the barrel EMC (| cos θ| < 0.8) or 50 MeV in the end cap EMC (0.86 < | cos θ| < 0.92). To reject photons radiated from charged particles, the angle between the photon candidate and the proton is required to be greater than 10 degrees. A more stringent cut of 30 degrees between the photon candidate and the anti-proton is applied to exclude the large number of photons from antiproton annihilation.
For events with one proton, one anti-proton, and at least two photons, a kinematic fit (4C) with the total four momenta of all particles constrained to the energy and three momentum-components of the initial e + e − system is applied. When more than two photons are found in an event, all possible ppγγ combinations are considered and the one yielding the smallest χ 2 4C is retained for further analysis. The χ 2 4C is required to be less than 30. After selecting the ppγγ candidate, the π 0 candidates are selected by requiring |M (γγ) − m π 0 | < 15 MeV/c 2 , where m π 0 is the nominal π 0 mass [27].
The Dalitz plot for the events passing the above selection criteria for data at √ s = 4.258 GeV is shown in Fig. 1(a). The corresponding invariant mass spectra of pp, pπ 0 andpπ 0 are shown in Fig. 1(b), (c) and (d), respectively.
The potential backgrounds for e + e − → ppπ 0 are studied using the inclusive MC sample at √ s = 4.26 GeV. After imposing all event selection requirements, the remaining background events are found to have the final state topologies e + e − → γpp, γγpp and γγγpp.
No other background survives. The non-π 0 background events can be evaluated from events in the π 0 sidebands. The π 0 sideband regions are defined as 0.07 < M (γγ) < 0.10 GeV/c 2 and 0.17 < M (γγ) < 0.20 GeV/c 2 . The background contamination estimated using π 0 sidebands at √ s = 4.258 GeV is 0.3%. The background contributions are neglected in the subsequent analysis.

Study of intermediate structures by Partial Wave Analysis
As shown in Fig. 1, a prominent structure near the threshold in the pp mass spectrum is visible. Structures are also seen in the pπ 0 andpπ 0 mass spectra. To evaluate the detection efficiencies of the decay e + e − → ppπ 0 properly, a partial wave analysis (PWA) is performed with the e + e − → ppπ 0 candidates to study the intermediate states present.
For the process e + e − → ppπ 0 , the isospin of the ppπ 0 system can be I = 0 or I = 1. The quasitwo-body decay amplitudes in the sequential decay pro- [30,31]. All 1 −− and 3 −− states above pp threshold, N * and ∆ * states with spin up to 5/2, listed in the summary tables of the PDG [27], are considered in this analysis. According to the framework of soft π meson theory [32], the off-shell decay process should be included. Thus, N (940) with a mass of 940 MeV/c 2 and zero width representing a virtual proton which could emit a π 0 is considered as a possible component. No isoscalar vector meson is considered, since there is no candidate above the pp threshold in the summary tables of the PDG. The ρ * states are parameterized by a constant-width relativistic Breit-Wigner (BW) propagator with barrier factors included. The N * and ∆ * states are parameterized by a BW propagator as described in Ref. [30]. The resonance parameters are fixed according to previous measurements [27] due to limited statistics. The complex coefficients of the amplitudes are determined by an unbinned maximum likelihood fit. The details of the likelihood function construction can be found in Ref. [33].
For ρ * states with J = 1, the pp final state interaction (FSI) effect using the Jülich model [34] is taken into consideration by factorizing the partial wave amplitude into the amplitude without the FSI effect and the S wave pp scattering amplitude in the scattering length approximation given in Ref. [34]. The direct process of e + e − → ppπ 0 can be modeled by 1 −− or 3 −− phase space of the pp system (1 −− or 3 −− PHSP). All combinations of the components in Ref. [35] are evaluated. The changes in the negative log-likelihood (NLL) and the number of free parameters in the fit with and without a resonance are used to evaluate its statistical significance. Resonances with significance greater than 5σ are retained in the PWA solution. The selection of PWA components is performed at the energy points with the high statistics, i.e. at √ s = 4.008, 4.226, 4.258 and 4.416 GeV, as shown in Table 1. The selected components are used to describe the data at other nearby energy points. The data at √ s = 4.189 − 4.600 GeV can be described by the N (1440), ρ(2150), ρ 3 (1990) and 1 −− PHSP amplitudes. The data at √ s = 4.008 − 4.085 GeV can be described by the N (1520), N (2570), ρ(2150), ρ 3 (1990) and 1 −− PHSP amplitudes. The N (940) is not included in the fits since its significance is less than 5σ. If we perform an alternative PWA fit with N (1440), ρ(2150), ρ 3 (1990) and 1 −− PHSP at √ s = 4.008 GeV, the NLL worsens by 37.8. The change of efficiency determined with the alternative fit with respect to the nominal value is considered as a source of systematic uncertainty. Comparisons of the data and the fit projection (weighted by MC efficiencies) in terms of the invariant mass spectra of pp, pπ 0 andpπ 0 at √ s = 4.258 GeV are shown in Fig. 1(b), (c) and (d), respectively. The χ 2 over the number of bins is displayed in those figures.

Cross section for e + e − → ppπ 0
The Born cross section for e + e − → ppπ 0 is determined as where N obs is the number of observed events; L is the integrated luminosity; ǫ is the detection efficiency derived from MC events generated according to the results of the PWA fit; (1 + δ r ) is the radiative correction factor, which is taken from a QED calculation taking the line shape of the cross section e + e − → ppπ 0 of data as input in an iterative procedure; (1 + δ v ) is the vacuum polarization factor, including leptonic and hadronic contributions, taken from a QED calculation with an accuracy of 0.5% [36]; and B π 0 is the branching fraction of π 0 decaying to γγ according to the PDG [27]. The measured Born cross section of e + e − → ppπ 0 at each energy point is listed in Table 1. Uncorrelated systematic uncertainties in the Born cross section measurements mainly originate from the π 0 mass window requirement, kinematic fit and the intermediate states in PWA. The systematic uncertainty from the requirement on the π 0 signal region is estimated by smearing the invariant mass of the γγ pair in the signal MC with a Gaussian function to compensate for the resolution difference between data and MC. The parameters for smearing are determined by fitting the π 0 distribution of data with the MC shape convoluted with a Gaussian function. The difference in the detection efficiency between signal MC samples with and without the extra smearing is taken as the systematic uncertainty. The systematic uncertainty due to the kinematic fit is estimated by correcting the helix parameters of charged tracks for the signal MC sample according the method described in Ref. [37]. The difference in the detection efficiency between the MC samples with and without this correction is taken as the systematic uncertainty. The systematic uncertainty from the intermediate Table 1: The results on e + e − → ppπ 0 . Shown in the table are the integrated luminosity L, the radiative correction factor (1 + δ r ), the vacuum polarization factor (1 + δ v ), the number of observed events N obs , the detection efficiency ǫ and the Born cross section σ B (e + e − → ppπ 0 ) at each energy point. The errors of ǫ are from the PWA fit. The first errors of σ B are statistical, and the second ones are systematic. states in PWA includes those from the BW parametrization, resonance parameters and extra resonances. Uncertainties from the BW parametrization of intermediate states are estimated by replacing the BW formula of N (1440) and N (1520) as used in Ref. [30] with a constant BW formula and replacing those of ρ(2150) and ρ 3 (1990) with the BW formula with the Gounaris-Sakurai (GS) model [38]. In the PWA fit, the resonance parameters are fixed according to the previous measurements [39,40]. Alternative fits are performed in which the resonance parameters are set as free parameters and the changes in the results are taken as systematic uncertainties. Uncertainties from additional resonances are estimated by adding the most significant additional resonance among each J P assignment in Ref. [35] into the PWA solution individually, and their influences on the cross section measurements are taken as the systematic uncertainties.
Correlated systematic uncertainties among the different energy points include those from luminosity measurement (1.0%) [41], MDC tracking (2% for two charged tracks) [42], particle identification (2% in total for proton and anti-proton) [43], photon detection efficiency (2%) [44] and radiative correction. The difference in ǫ(1+δ r ) between the third and fourth iteration is taken as the systematic uncertainty due to the radiative correction, as the radiative correction dependent quantity ǫ(1 + δ r ) converges after three iterations.
The total systematic uncertainty of the different energy points is calculated by adding the individual uncertainties in quadrature as shown in Table 2.
(2) to the calculated cross sections. In Eq. (2), σ con and σ Y represent the continuum cross section and resonant cross section, respectively, and σ con can be described by a function of s, σ con = C/s λ , where the exponent λ is a priori unknown. The parameter φ describes the phase between resonant and continuum production amplitudes. The mass m and width Γ of the Y (4260) are fixed to the PDG values [27]. The values of C, λ, σ Y , and the interference phase φ are free in the fit. The uncorrelated systematic uncertainties in the Born cross section measurements are directly considered in the fit and the effect of the correlated systematic uncertainties on the final results is estimated by the method in Ref. [45], in which the error propagation is determined from shifting the data by the aforementioned correlated uncertainties and adding the deviations in quadrature. In addition, the uncertainties for the beam energy measurements of all the data points taken from Ref. [46] are considered in the fit. The best fit function is shown in Fig. 2 as the solid line. The dashed line represents the fit with σ Y = 0. The optimal value of σ Y is (1.6±5.9)×10 −3 pb with a statistical significance of 0.5σ. The significance is calculated based on the changes in the χ 2 value and the number of free parameters in the fit with and without the assumption of existence of the Y (4260) resonance. The result for the phase between resonant and continuum production amplitudes is φ = 3.4±1.0. The parameters describing the slope of the continuum cross section are C = (5.4 ± 5.3) · 10 5 GeV 2λ pb and λ = 4.2 ± 0.4. The upper limit on σ Y at the 90% C.L., σ up Y , is determined by σ up Y 0 G(σ Y , σ σY )dx/ ∞ 0 G(σ Y , σ σY )dx = 0.9, where G(σ Y , σ σY ) is a Gaussian function with mean value σ Y = 1.6 × 10 −3 pb and standard deviation σ σY = 5.9 × 10 −3 pb. The uncertainties from mass and width of the Y (4260) are considered by varying them by one standard deviation according to the PDG values [27] and the most conservative σ up Y is taken as the final result. The obtained upper limit is 0.01 pb.

Summary
Based on 13 data samples between √ s = 4.008 and 4.600 GeV collected with the BESIII detector, the process e + e − → ppπ 0 is studied. The Born cross section of e + e − → ppπ 0 is measured. No resonant structure is observed in the shape of the cross section. The upper limit on the Born cross section of e + e − → Y (4260) → ppπ 0 at the 90% C.L is estimated to be 0.01 pb.