Feasibility study for the measurement of $\pi N$ TDAs at PANDA in $\bar{p}p\to J/\psi\pi^0$

The exclusive charmonium production process in $\bar{p}p$ annihilation with an associated $\pi^0$ meson $\bar{p}p\to J/\psi\pi^0$ is studied in the framework of QCD collinear factorization. The feasibility of measuring this reaction through the $J/\psi\to e^+e^-$ decay channel with the PANDA (AntiProton ANnihilation at DArmstadt) experiment is investigated. Simulations on signal reconstruction efficiency as well as the background rejection from various sources including the $\bar{p}p\to\pi^+\pi^-\pi^0$ and $\bar{p}p\to J/\psi\pi^0\pi^0$ reactions are performed with PandaRoot, the simulation and analysis software framework of the PANDA experiment. It is shown that the measurement can be done at PANDA with significant constraining power under the assumption of an integrated luminosity attainable in four to five months of data taking at the maximum design luminosity.

The exclusive charmonium production process inpp annihilation with an associated π 0 mesonpp → J/ψπ 0 is studied in the framework of QCD collinear factorization. The feasibility of measuring this reaction through the J/ψ → e + e − decay channel with the PANDA (AntiProton ANnihilation at DArmstadt) experiment is investigated. Simulations on signal reconstruction efficiency as well as the background rejection from various sources including thepp → π + π − π 0 andpp → J/ψπ 0 π 0 reactions are performed with PandaRoot, the simulation and analysis software framework of the PANDA experiment. It is shown that the measurement can be done at PANDA with significant constraining power under the assumption of an integrated luminosity attainable in four to five months of data taking at the maximum design luminosity.

I. INTRODUCTION
Understanding of the hadronic structure in terms of the fundamental degrees of freedom of QCD is one of the fascinating questions of the present day physics. Lepton beam initiated reactions, allowing to resolve individual quarks and gluons inside hadrons, proved to be a handy tool for this issue. The factorization property established for several classes of hard (semi-)inclusive and exclusive processes allows to separate the short distance dominated stage of interaction and the universal non-perturbative hadronic matrix elements. Some of the matrix elements which have been the subject of significant interest include the Parton Distribution Functions (PDFs) [1], Generalized Parton Distributions (GPDs) [2,3], Transverse Momentum Dependent Parton Distribution Functions (TMD PDFs) [4], (Generalized) Distribution Amplitudes ((G)DAs) [5] and Transition Distribution Amplitudes (TDAs) [6,7] encoding valuable information on the hadron constituents.
Alongside with the study of lepton beam induced reactions, one can get access to the same non-perturbative functions in a complementary way by considering the cross conjugated channels of the corresponding reactions. For example, protonantiproton annihilation into a lepton pair and a photon (or a meson) can be seen as the cross conjugated counterpart of the leptoproduction of photons (or mesons) off protons, and provides access to nucleon GPDs and/or nucleon-to-photon (nucleon-to-meson) TDAs.
Such investigations have been hindered up to now by the limitations of antiproton beam luminosities. However, very significant results on the electromagnetic form factors in the time-like region using the pp → e + e − reaction were obtained by the E835 experiment at FNAL (Fermilab National Accelerator Laboratory) [8]. But inclusive lepton pair production and hard exclusive channels still remain unexplored. This situation will be largely improved in the next decade with the availability of the high intensity antiproton beam at FAIR (Facility for Antiproton and Ion Research) with momentum up to 15 GeV/c. The PANDA experiment [9] will dedicate an important part of its physics program to the investigation of the nucleon structure in antiproton-proton annihilation reactions. It includes the detailed study of the time-like electromagnetic nucleon form factors employing both the e + e − and µ + µ − production channels in a broad kinematic range. It is also planned to access PDFs through the Drell-Yan mechanism, by measuring inclusive e + e − production and GPDs considering γ * γ and γ * π 0 exclusive channels at large production angles. Finally, the reactionspp → γ * M → e + e − M and pp → J/ψM → e + e − M, where M stands for a light meson M = {π 0 , η, ρ 0 , ω, . . .}, are proposed to study nucleon-tomeson TDAs.
Nucleon-to-meson (and particularly nucleon-to-pion) TDAs were introduced as a further generalization of the concepts of both GPDs and nucleon light-cone wave functions (DAs). They describe partonic correlations inside nucleons and allow to access the non-minimal Fock components of the nucleon light-cone wave function with additional quark-antiquark pair seen as a light meson. Therefore, in particular, πN TDAs provide information on the nucleon's pion cloud.
Nucleon-to-pion TDAs arise in the collinear factorized description of several hard exclusive reactions such as backward electroproduction of pions off nucleons [10,11], which can be studied at JLab [12] and COMPASS in the space-like regime, while PANDA will provide access to the same nonperturbative functions in the time-like regime [13,14].
On the theory side, the possibility to study nucleon-tomeson TDAs is provided by the collinear factorization theorem similar to the well known collinear factorization theorem for hard meson electroproduction [15], giving rise to the description in terms of GPDs. However, the collinear factorization theorem for the TDA case has never been proven explicitly. Therefore, one of the important experimental tasks is to look for experimental evidence of the validity of the factorized description of the corresponding reactions in terms of nucleon-to-meson TDAs. This can be done either by verifying the appropriate scaling behavior or by checking the angular dependence of the produced lepton pair specific for the dominant reaction mechanism. Bringing trustworthy evidence for the validity of the factorized description of a new class of hard exclusive reaction will, by itself, represent a major experimental achievement of PANDA.
Recently, a detailed study of the access to πN TDAs in the reactionpp → γ * π 0 → e + e − π 0 following the cross section estimates of Ref. [13] with PANDA has been presented in Ref. [16]. The investigation of the reactionpp → J/ψπ 0 → e + e − π 0 constitutes a natural complement to the latter study. The resonant case presents the noticeable advantage of a larger cross section, and a cleaner signal selection due to the resonant e + e − production. The simultaneous measurement of both resonant and non-resonant channels provides constraints to the πN TDAs in different kinematic ranges, and allows to test the universality of πN TDAs. While the non-resonant pp → γ * π 0 → e + e − π 0 has never been measured, some scarce data exist forpp → J/ψπ 0 → e + e − π 0 [17][18][19], which can be used to constrain the predictions. Production of a J/ψ with an associated π 0 inpp collisions has indeed been investigated in the past by the E760 experiment at FNAL, since it constitutes a background in the search for charmonium states via their decay into J/ψπ 0 . An important part of the PANDA program will also focus on such studies, as described in Ref. [20]. This brings additional motivation for the detailed measurements of thepp → J/ψπ 0 → e + e − π 0 reaction. From the theory point of view, access to the πN TDA in the J/ψ production channel is also more favorable, since one can take advantage of the known J/ψ →pp decay width [21] in order to reduce ambiguities related to the choice of the phenomenological parametrization for the relevant nucleon DA [22].
The aim of the present study is therefore to explore the feasibility of the measurement of the reactionpp → J/ψπ 0 → e + e − π 0 with PANDA at different incident momenta of the antiproton beam, based on the cross section estimates of Ref. [22]. The paper is organized as follows: Section II outlines the design of the PANDA experimental setup with a focus on the most relevant components to the analysis. Section III covers the properties of thepp → J/ψπ 0 → e + e − π 0 reaction which constitutes the signal. In Section IV, the different background contributions are discussed. Section V is devoted to the description of the simulation and analysis procedure. In Section VI, the expected precision on differential cross section measurements is presented.

A. FAIR Accelerator Complex
The FAIR accelerator complex which is under construction to extend the existing GSI (Gesellschaft für Schwerionenforschung) facilities in Darmstadt, Germany, will provide beam for four experimental pillars, one of which is the PANDA experiment dedicated to hadronic physics. FAIR will use the existing SIS18 synchrotron as an injection ring into a new larger synchrotron SIS100. The SIS100 ring will generate an intense pulsed beam of protons with energies reaching up to 29 GeV that can be directed at an antiproton production target. Time averaged production rates in the range of 5.6 × 10 6 to 10 7p s −1 are expected. Antiprotons are collected and phase space cooled in the CR (Collector Ring), then transferred to the RESR (Recycled Experimental Storage Ring) accumulator, and then injected into the HESR (High Energy Storage Ring), equipped with stochastic and electron cooling, where they will be used by the PANDA experiment in a fixed target setup. This full setup is designed to provide beams with up to 10 11 antiprotons, and peak instantaneous luminosities reaching up to 2 × 10 32 cm −2 s −1 . Such a scenario will allow the accumulation of an integrated luminosity of 2 fb −1 in about five months, which will be used as the basis for all results shown in this analysis. However, the more likely scenario currently is a staged construction with a reduced setup at the start of operations without the RESR until the full design can be realized. In this case, the HESR will be used as an accumulator in addition to its original task of cooling the beam and storing it for experiments with internal targets. This will result in a luminosity that is about a factor ten lower than the full design goal during the initial phases of operating FAIR.

B. The PANDA Detector
The proposed PANDA detector is depicted in Fig. 1. The discussion here will focus on the subsystems that are particularly relevant for the presented analysis. The PANDA detector consists of the target spectrometer surrounding the target area and the forward spectrometer designed to detect particles in the forward rapidity region. The target spectrometer is divided into a barrel region with polar angle reach from 22 • to 145 • , and an endcap region that covers polar angles below 22 • , down to 10 • in the horizontal plane and 5 • in the vertical plane. Particles with polar angles below the endcap coverage are detected by the forward spectrometer. In addition, PANDA will be equipped with a Luminosity Monitor Detector (LMD) at very forward angles, built for precise determination of both absolute and relative time integrated luminosities.
Two technologies are being developed to provide a hydrogen target with sufficient density that allows to reach the design luminosity within the restricted space available [9]. A target thickness of 4 × 10 15 hydrogen atoms per cm 2 is required to achieve a peak luminosity of 2 × 10 32 cm −2 s −1 assuming 10 11 stored antiprotons in the HESR. The Cluster-Jet Target system operates by pumping pressurized cold hydrogen gas into vacuum through a Laval-type nozzle, leading to a condensation of hydrogen molecules into a narrow jet of hydrogen clusters with each cluster containing 10 3 -10 6 hydrogen molecules. The main advantages of this setup are the homogeneous density profile and the ability to focus the antiproton beam at the highest areal density point. The Pellet Target system creates a stream of frozen hydrogen micro-spheres (pellets) of diameter 25 -40 µm, that cross the antiproton beam perpendicularly. The Pellet Target will be equipped with an optical tracking system that can determine the vertex position of individual events with high precision.
The innermost part of the barrel region is occupied by charged particle tracking detectors, which in turn are surrounded by particle identification detectors, followed by a solenoid magnet that generates a nearly uniform 2 T field pointing in the direction of the beam. The innermost layers of tracking are provided by the Micro Vertex Detector (MVD) [23], based on silicon pixel detectors for the innermost two layers, and a double sided strip detectors for the remaining two layers. The Straw Tube Tracker (STT) [24] is constructed from aluminized mylar tubes with gold-plated tungsten anode wires running along the axis. The tubes operate with an active gas mixture composed of argon and CO 2 held at a pressure of 2 bar allowing them to be mechanically self-supporting. The STT adds only about X/X 0 ≈ 1.2% to the total radiation length of the tracking system on top of the ≈ 10% expected from the MVD. The MVD and the STT also measure ionization energy loss by charged particles in their layers. A truncated mean of the specific energy loss of the tracks measured by the MVD and STT layers is used to estimate the energy loss dE/dx for each track.
The barrel region tracking subsystems are immediately surrounded by various dedicated particle identification (PID) detectors. The innermost PID detector that is used in this analysis is the barrel DIRC [25] (Detection of Internally Reflected Cherenkov light), where particles are identified by the size of the Cherenkov opening angle. The DIRC is followed by the barrel Electromagnetic Calorimeter (EMC) [26], constructed from lead tungstate (PbWO 4 ) doped crystals, operated at a temperature of -25 • C to optimize the light yield. Photons from each crystal in the barrel EMC are detected by a pair of APDs (Avalanche Photodiodes). The EMC constitutes the most powerful detector for the identification of electrons through the momentum-energy correlation, particularly at momenta higher than ≈ 1 GeV/c. The DIRC, together with the MVD and STT dE/dx measurements, ensure coverage at lower momenta where the EMC electron identification (EID) capacity is weaker.
The design of the barrel spectrometer of PANDA also includes a Muon Range System (MRS) surrounding the Solenoid for the identification of muons as well as Time of Flight (TOF) detectors between the tracking layers and the DIRC for general PID. The MRS and TOF are not used in the analysis presented here.
For the endcap section of the target spectrometer, tracking points are provided by the STT as well as a set of four disc shaped MVD layers, and three chambers of Gas Electron Multiplier (GEM) trackers. PID is performed using information from the endcap EMC and the endcap DIRC [27] (also called Disc DIRC in reference to its geometrical shape). Apart from its location, and angular coverage, the Disc DIRC operates using the same basic principle as the barrel DIRC. The endcap EMC is instrumented using the same PbWO 4 crystals as the barrel EMC, however photon detection is performed using APDs only for the outer lying crystals. The crystals close to the beam axis are readout by VPTs (Vacuum Phototriodes) because of the stringent requirements on radiation hardness there.
The design of PANDA also provides for coverage at angles below those of the target spectrometer (< 5 • ), through the Forward Tracking System (FTS), a Ring Imaging Cerenkov (RICH) detector system and a Shashlyk calorimeter. Charged particles traversing the forward tracking system are subject to a field integral of 2 Tm generated by a dipole magnet, allowing for momentum determination. For this analysis, the forward Shashlyk calorimeter was not included in the simulations. As a result, our efficiency prediction is underestimated for events whose kinematics leads to charged particles requiring energy measurement for identification in the extremely forward direction. With the full PANDA setup, the performance for such events will be better than that reported in this paper.
Precise determination of integrated luminosity is a critically important ingredient for the whole PANDA physics program. The LMD is a detector that has been designed to provide both absolute and relative time integrated luminosity measurements with 5% and 1% systematic uncertainty, respectively [28]. The LMD will rely on the determination of the differential cross section of antiproton -proton elastic scattering into a polar angle (laboratory reference frame) range of 3.5 -8 mrad and full azimuth to achieve this goal. The LMD tracks antiprotons using four planes of HV-MAPS (High Voltage Monolithic Active Pixel Sensors) tracking stations, a setup chosen to fulfil the constraints of high spatial resolution and low material budget. The LMD will be placed 10.5 m downstream from the interaction point within a vacuum sealed enclosure to reduce systematic uncertainties from multiple scattering.

C. Simulation and Analysis Software Environment
The simulation and analysis software framework called PandaRoot [29,30] is used for the feasibility study described here. PandaRoot is a collection of tools used for the simulation of the transport of particles through a GEANT4 [31] implementation of the PANDA detector geometry, as well as detailed response simulation and digitization of hits in the various detector elements that takes into account electronic noise. Software for the reconstruction of tracks based on the simulated tracking detector hit points is implemented in PandaRoot, as well as the association of reconstructed tracks to signal in outer PID detectors. A PID probability is assigned to each track based on the response in all the outer detectors, using a simulation of five possible particle species: e ± , µ ± , π ± , K ± , p ± . Clusters in the electromagnetic calorimeter that are not associated to a reconstructed track are designated as neutral candidates, and used for the reconstruction of photons. The simulation studies are done in the e + e − decay channel for the J/ψ due to much higher EID efficiency as compared to muon identification efficiency, at a fixed pion rejection probability. This is particularly pertinent forpp → J/ψπ 0 , due to the very high pionic background event rates, as will be discussed in Section IV.
The tracking points from the tracking detectors are used for pattern recognition to find charged particle tracks. Points that are found to belong to the same track are in a first step fitted to a simplified helix for an initial estimate of the momentum, which is then used as a starting point for an iterative Kalman Filter procedure relying on the GEANE [32] track follower. The output from the Kalman Filter is a more refined estimate of the momentum of tracks that takes into account multiple scattering as well as changes in curvature due to energy loss in the detector material.
Since the Kalman Filter does not take into account non-Gaussian alterations of track parameters, it can not correctly handle changes of track momentum through Bremsstrahlung energy loss. This is particularly pernicious for electrons that can lose on average up to 10% of their total energy through photon emission. To correct this effect event by event, a procedure was developed [33]. For each track this algorithm looks for potential Bremsstrahlung photon candidates in the EMC, and adds that energy to the track. This method was demonstrated to work over a wide range of momenta and angles including those relevant for this analysis.

III. THEORETICAL DESCRIPTION OF THE SIGNAL CHANNEL
The feasibility study presented here is carried out at three values of the square of the center-of-mass (c.m.) energy s: 12.3 GeV 2 , 16.9 GeV 2 and 24.3 GeV 2 . The first value is chosen to coincide with existing data from E835 forpp → J/ψπ 0 [17][18][19]. The remaining two values are chosen at the incidentp momenta of 8 GeV/c and 12 GeV/c, respectively, to explore the kinematic zone between the first point and the maximum availablep momentum at FAIR of 15 GeV/c.

A. Kinematics
In order to present the cross section estimates within the description based on the πN TDAs employed for our feasibility study, we would like to review briefly the kinematics of the signal reaction: making special emphasis on the kinematic quantities employed in the collinear factorization approach. The natural hard scale for the reaction (1) is introduced by the c.m. energy squared s = (p N + pN) 2 and the charmonium mass squared M 2 ψ . The collinear factorized description is supposed to be valid in the two distinct kinematic regimes, corresponding to the generalized Bjorken limit (large s and Q 2 ≡ M 2 ψ for a given s/Q 2 ratio and small cross channel momentum transfer squared): • the near-forward kinematics |t| ≡ |(p π − pN) 2 | s, M 2 ψ ; it corresponds to the pion moving almost in the direction of the initial antinucleon in the NN center-ofmass system (CMS); • the near-backward kinematics |u| ≡ |(p π − p N ) 2 | s, M 2 ψ corresponding to the pion moving almost in the direction of the initial nucleon in the NN CMS.
Due to the charge-conjugation invariance of the strong interaction there exists a perfect symmetry between the nearforward and near-backward kinematic regimes of the reaction (1). These two regimes can be considered in exactly the same way, and the amplitude of the reaction within the uchannel factorization regime can be obtained from that within the t-channel factorization regime, with the help of the obvious change of the kinematic variables (see Eq. (8) below). In the NN CMS these two regions look perfectly symmetric. However, we note that PANDA operates with the antibaryon at beam momentum and the baryon at rest in the lab frame. Consequently, the symmetry between the near-forward and near-backward kinematics is not seen immediately in the PANDA detector. Moreover, this introduces acceptance differences between the two regimes which will be explored in Section 6 as a function of the incident pN momentum.
For definiteness below, we consider the near-forward (tchannel) kinematic regime. The detailed account of the relevant kinematic quantities is presented in Appendix A of Ref. [14]. It is convenient to choose the z-axis along the nucleon-antinucleon colliding path, selecting the direction of the antinucleon as the positive direction. Introducing the lightcone vectors p t and n 2 (2p t · n t = 1), one can perform the Sudakov decomposition of the particle momenta. Neglecting the small mass corrections (assuming m π = 0 and m N √ s) we obtain: where ξ t stands for the t-channel skewness variable, which characterizes the t-channel longitudinal momentum transfer: We also introduce the transverse (with respect to the selected z-axis) t-channel momentum transfer squared (∆ t T ) 2 . This quantity can be expressed in terms of the skewness variable ξ from Eq. (3) and the t-channel momentum transfer squared (∆ t ) 2 ≡ t as: For (∆ t T ) 2 = 0 the momentum transfer is purely longitudinal and hence the pion is produced exactly in the forward direction. This corresponds to the maximal possible value of the momentum transfer squared: The t -channel skewness variable can be expressed through the reaction invariants as: In the present study, following Ref. [22] we neglect all t/s and m 2 N /s corrections, and employ the simple expression for the skewness variable: In order to apply the same formalism for the u-channel (near-backward) kinematic regime, it suffices to perform the following variable transformations in the relevant formula: Therefore, in what follows we omit the superscript referring to the kinematic regime for the kinematic variables.
The generalized Bjorken limit, in which the validity of the collinear factorized description of the reaction (1) is assumed, is defined by the requirement ∆ 2 s, Q 2 ≡ M 2 ψ . There is no explicit theoretical means to specify quantitatively the condition |∆ 2 | s, Q 2 . However, the common practice coming from studies of the similar reactions suggests |∆ 2 | < 1 GeV 2 can be taken as a reasonable estimate. For the fixed value of the skewness parameter this condition can be translated into the corresponding kinematic cut for ∆ 2 T or, equivalently, to cuts in the pion scattering angle in the CMS and (after the appropriate boost transformation) in the lab frame. This allows to specify the span of the "forward" and "backward" cones in which the collinear factorized description is supposed to be valid for the reaction (1).
However, once such a kinematic cut has been implemented, one has to be prudent. Indeed, the kinematic formulas derived in the Appendix A of Ref. [14] represent the approximation, which is valid in the generalized Bjorken limit s, Q 2 → ∞ while its validity for a given kinematic setup is not necessarily ensured. Certainly it is useless to employ the approximate kinematic formulas in the region where the approximation does not provide a satisfactory description of the kinematic quantities.
Therefore, it is instructive to compare the approximate result for (∆ t T ) 2 given by Eq. (4) with ξ t given by Eq. (6) and by the less accurate expression from Eq. (7) with the general result obtained with the help of the exact kinematic relation for the CMS pion scattering angle θ * π : Here, denotes the CMS momentum of the produced pion, where: is the usual Mandelstam function, and the s-channel CMS scattering angle is expressed as: .
(10) Figure 2 shows the comparison of (∆ t T ) 2 computed using the approximate formula from Eq. (4) with the exact result from Eq. (9) for the three selected values of s. We plot (∆ t T ) 2 as a function of ∆ 2 for ∆ 2 ≤ ∆ 2 max , where ∆ 2 max is the maximal kinematically accessible value of ∆ 2 , corresponding to ∆ T = 0 (i.e., the pion produced exactly in the forward direction). The validity limit of our kinematic approximation is shown in Fig. 2 by the solid vertical lines. It is calculated by imposing the maximum allowable deviation of 20% on the value of (∆ t T ) 2 from the exact result from Eq. (9). The final validity range ∆ 2 min is determined by picking the more conservative limit of the kinematic approximation and the standard constraint |∆ 2 | ≤ 1 GeV 2 , ensuring the smallness of the |∆ 2 | comparing to s and Q 2 = M 2 ψ in the generalized Bjorken limit. Table I summarizes the valid ranges of ∆ 2 in which we are going to apply the factorized description of the reaction (1) for the three selected energies considered here. The lower limit comes from the applicability of collinear factorization (Bjorken limit) for the lowest beam momentum (5.5 GeV/c) and from the shortcoming of the approximation employed for the kinematic quantities (kinematic constraints) for the higher two beam momenta (8 and 12 GeV/c). The last two columns show how these limits translate to the polar angle of the π 0 in the lab frame, θ π 0 lab , namely the maximum (minimum) valid θ π 0 lab in the near-forward (near-backward) validity range. At the other end, the minimum (maximum) valid polar angle in the near-forward (near-backward) validity range is 0 • (180 • ) for all energies. TABLE I. The kinematic range in which we assume the validity of the factorized description of the signal channel in terms of πN TDAs and nucleon DAs for the three values of the incidentp momentum employed in the present study. The last two columns give the limits in terms of the polar angle of the π 0 in the lab frame for the nearforward and near-backward regimes.

B. Cross Section Estimates within the Collinear Factorization Approach
The calculation of the N +N → J/ψ + π cross section within the collinear factorization approach follows the same main steps as those in the calculations of the J/ψ →pp decay width in the perturbative QCD (pQCD) approach [34][35][36]. The small and large distance dynamics is factorized, and the corresponding amplitude is presented as the convolution of the hard part, computed in the pQCD, with the hadronic matrix elements of the QCD light-cone operators (πN TDAs and nucleon DAs) encoding the long distance dynamics (see Fig. 3). The hard scale, which justifies the validity of the perturbative description of the hard subprocess, is provided by the mass of Collinear factorization of the annihilation process N(pN )N(p N ) → J/ψ(p ψ )π(p π ). Top panel: near-backward kinematics (u ∼ 0). Bottom panel: near-forward kinematics (t ∼ 0).N(N) DA stands for the distribution amplitude of antinucleon (nucleon); πN(πN) TDA stands for the transition distribution amplitude from a nucleon (antinucleon) to a pion. Figures reproduced from Ref. [22].
Within the approximation to order leading twist-three, only the transverse polarization of the J/ψ is relevant, and the cross section with the suggested reaction mechanism for either the near-forward or the near-backward kinematic regimes reads [22]: where the squared matrix element |M T | 2 is expressed as: Here the factor C reads: where f π = 93 MeV is the pion weak decay constant; f ψ = (413 ± 8) MeV is the normalization constant of the non-relativistic light-cone wave function of heavy quarkonium; α s stands for the strong coupling. The two functions I (ξ , ∆ 2 ) and I (ξ , ∆ 2 ) denote the convolutions of the hard scattering kernels with the πN TDAs and (anti)nucleon DAs. In our studies we employ the estimate of the cross section from Eq. (11) within the simple cross channel nucleon exchange model for the πN TDAs [37]. In this model the convolution integrals I , I read as: where f N = (5.0 ± 0.5) GeV 2 is the nucleon wave function normalization constant; g πNN 13 is the phenomenological pion-nucleon coupling; M 0 is the standard convolution integral of nucleon DAs occurring in the expression for the J/ψ →pp decay width within the pQCD approach of Ref. [36]: In order to compute the value of the cross section given by Eq. (11), one has to employ the phenomenological solutions for the leading twist nucleon DAs to compute the convolution integral M 0 and to specify the appropriate value of the strong coupling α s for the characteristic virtuality of the process. Unfortunately, some controversy on this issue exists in recent literature (see e.g. the discussion in Chapter 4 of Ref. [38]). There are several classes of phenomenological solutions for the leading twist nucleon DAs: • one class (usually referred as the Chernyak-Zhitnitskytype solutions) contains nucleon DAs which differ considerably from the asymptotic form of the leading twist nucleon DAs. Such DAs require α s ∼ 0.25 to describe the experimental charmonium decay width Γ(J/ψ → pp) from Eq. (15); • another class of solutions (e.g. the Bolz-Kroll solution [39] and Braun-Lenz-Wittmann NLO model [40]) contains nucleon DAs which instead are rather close to the asymptotic form and require α s ∼ 0.4 to reproduce the experimental Γ(J/ψ →pp).
For our rough cross section estimates, intended for the feasibility studies, we follow the prescription from Ref. [22], and employ the results for thepp → J/ψπ 0 cross section with the value of α s fixed by the requirement that the given phenomenological solution for the nucleon DAs reproduces the experimental J/ψ →pp decay width. In this case the simple nucleon pole model for πN TDAs results in the samē pp → J/ψπ 0 cross section predictions for any input phenomenological nucleon DA solution.  Figure 4 shows the prediction of the s dependence of the differential cross section ofpp → J/ψπ 0 at |∆ 2 T | = 0 which decreases only by a factor of about 4 between the lowest to highest CMS collision energies that will be available at PANDA. Figure 5 shows the ∆ 2 T dependence of the cross section for the three collision energies which are used for the studies addressed here, each of them limited to its validity range. It is possible to compare this prediction to the E835 measurement of thepp → J/ψπ 0 → e + e − π 0 cross section [17][18][19], taken at an incidentp momentum of 5.5 GeV/c that corresponds to s = 12.3 GeV 2 or √ s = 3.5 GeV. This value of √ s is 25 MeV below the threshold for the production of the h c resonance. The reported values lie in the range from 90 pb to 230 pb. Integrating the differential cross section from the TDA model at s = 12.3 GeV 2 over its validity range (−0.092 < t [GeV 2 ] < 0.59, which is smaller than the full kinematically accessible range), we find 206.8 pb, combining near-forward and near-backward kinematics. This value is on the upper end of the measured range by E835. However, given that the majority of the production rate from the TDA model prediction lies within the validity range, this result gives a degree of confidence that the model can be used as a basis for a feasibility study within its validity range.
Let us stress that the validity of the factorized description of the reaction (1) in terms of πN TDAs and nucleon DAs has only been conjectured. The corresponding collinear factorization theorem has never been proven explicitly. Therefore, one of the important experimental challenges is to establish evidence for the validity of this description. In general, there are several essential marking signs for the onset of the collinear factorization regime for a given hard exclusive reaction. The most obvious one is the characteristic scaling behavior of the cross section with the relevant virtuality 1/Q 2 . However, for the reaction (1) this feature is of little use, since the virtuality is fixed by the mass of the heavy quarkonium. Another opportunity is to look for the specific polarization dependence. For the case of the nucleon-antinucleon annihilation into J/ψ in association with a forward (or backward) neutral pion it is the transverse polarization of the J/ψ that is dominant within the collinear factorized description in terms of πN TDAs. This dominating contribution manifests through the characteristic (1 + cos 2 θ * ) distribution of the decay leptons in the lepton pair CMS scattering angle θ * . The dominance of the corresponding polarization has to be verified by means of a dedicated harmonic analysis.

C. Event Generator
A Monte Carlo (MC) event generator for the signal reaction was implemented by relying on Eq. (13). The angular distribution of π 0 s in the lab frame from this generator is shown in Fig. 6. The entire near-forward kinematic validity range of the reaction is concentrated in a polar angle window below ≈ 30 • (and even smaller window at higher collision energy), whereas the near-backward kinematic validity range occupies a window above ≈ 45 • (larger at higher collision energy) extending all the way to 180 • . This has implications for the signal reconstruction efficiency as will be shown in Section V I.

IV. BACKGROUND PROPERTIES
In the context of the PANDA detector setup, the signal reactionpp → J/ψπ 0 → e + e − π 0 has multiple background sources with varying degree of importance. These background sources are the subject of discussion in this section. For each source, rate estimates are given and in cases where a full MC is warranted, the details about the event generators that were used to simulate events are presented.
A. Three Pion Productionpp → π + π − π 0 Thepp → π + π − π 0 reaction is an important background source since the cross section is orders of magnitude larger than that of the signal and the possibility to misidentify the charged pion pair as an electron-positron pair. This reaction has been studied in the past at various incidentp momenta. Despite the limited statistics that were collected, the results from these early measurements tabulated in Ref. [41] provide a valuable benchmark for this study. As shown in Fig. 7, these results point to a steep decline of the total cross section as a function of incidentp momentum. For our background simulations, we used total cross sections slightly larger than the interpolation between the nearest existing data points. Very few experiments collected enough statistics to give high quality spectra for this reaction. As a result, the feasibility study relies on the hadronic event generator DPM (Dual Parton Model) [42] to simulate the shape of the spectra, since the model is constrained by taking into account the sparse experimental differential distributions. In particular, it includes the production of ρ and f 2 resonances in agreement with experimental observation [43]. Table II shows the J/ψπ 0 and π + π − π 0 cross sections in a mass window of 2.8 -3.3 GeV/c 2 for the charged pion pair, within the validity ranges of the TDA model (cf. Table I). The validity range includes both the near-forward and nearbackward kinematic approximation zones. The last column of the table gives the approximate signal to background ratio of produced event rates taking into account the J/ψ → e + e − branching fraction of 5.69%. The ratio of signal to background event production rate is of the order 10 −7 to 10 −6 . The fact that the signal e + e − invariant mass distribution is peaked allows to gain a factor of about 10 in signal to background ratio (S/B) before any PID information has been used. The rejection of the rest of the background has to come mostly from PID. The resulting strong requirement on the electron-pion discrimination power will make PID cuts a crucial component of the analysis. Multi-pion (N π ≥ 4) final states with at least one π + π − pair are also potential sources of background to thepp → J/ψπ 0 → e + e − π 0 channel, if one charged pion pair is misidentified as an e + e − pair. They have cross sections that are up to a factor 15 higher than three pion production, but larger rejection factors can be achieved due to the different kinematics. To confirm this, we performed detailed simulation studies for the π + π − π 0 π 0 and π + π − π + π − π 0 channels, and used the results to conclude on the rejection capability for channels with higher number of pions. The cross sections TABLE II. Signal (pp → J/ψπ 0 → e + e − π 0 ) cross section and production yields, background (pp → π + π − π 0 ) cross sections and signal over background ratio. The cross sections have been integrated over the validity range of the model (both near-forward and nearbackward). The branching ratio for J/ψ → e + e − and the mass cut on the π + π − pair in the range 2.8 -3.3 GeV/c 2 have been taken into account. The production yields for the signal are integrated for a luminosity of 2 fb −1 , corresponding to about five months at full luminosity.
pp ∆σ for the simulation of multi-pion final states is estimated by scaling the π + π − π 0 cross sections discussed in the previous section, with the scaling factor derived from DPM simulations. Table III gives the ratio of cross sections between those multi-pion final states (π + π − π 0 π 0 and π + π − π + π − π 0 ), and the three-pion final state (π + π − π 0 ).
For the cross section of thepp → J/ψπ 0 π 0 → e + e − π 0 π 0 channel, which is not predicted by the DPM model, we assumed a scaling to the J/ψπ 0 channel according to: where the π + π − π 0 π 0 to π + π − π 0 ratios were calculated based on the DPM event generator output. The results forpp → J/ψπ 0 π 0 → e + e − π 0 π 0 are 35.3 pb, 52.7 pb and 40.7 pb for beam momenta of 5.5 GeV/c, 8.0 GeV/c and 12.0 GeV/c, respectively. Although there is no existing measurement of thepp → J/ψπ 0 π 0 cross section to confirm these assumptions, we provide arguments below, which suggest that they are reasonably conservative. The E760 collaboration reported in Ref. [17] the nonobservation of a signal forpp → J/ψπ 0 π 0 → e + e − π 0 π 0 at c.m. energies close to the h c (1P) mass of 3.5 GeV/c 2 . This determines an upper limit for the cross section of 3 pb, which is about a factor 10 below our assumption at this energy. Calculations of non-resonant channels for thepp → J/ψπ + π − reaction at the X(3872) energy have been performed with an hadronic model [44], yielding a cross section of about 60 pb forpp → J/ψπ + π − → e + e − π + π − at √ s around 3.872 GeV. With the assumption of a factor two smaller cross section for pp → J/ψπ 0 π 0 → e + e − π 0 π 0 based on isospin coefficients, this calculation is consistent with our assumption. These calculations, outlined in Ref. [44], are however likely to have been overestimated due to the absence of vertex cut-off form factors, as in Ref. [45] for the case ofpp → J/ψπ 0 . Finally, the cross section for the production of thepp → J/ψπ 0 π 0 → e + e − π 0 π 0 channel via feed-down from a resonance might in some cases be expected to be lower than, or of the same magnitude, as our assumptions. This is obviously the case for resonances with charge conjugation C = 1, e.g. the X(3872), which do not decay into J/ψ plus any number of π 0 s. But charmonium states with C = 1 might also contribute to the J/ψπ 0 π 0 channel with cross sections that are lower or of the same order of magnitude. For example, a prediction of 30 pb for the reactionpp → Y (4260) → J/ψπ 0 π 0 → e + e − π 0 π 0 at √ s = 4.260 GeV was provided in Ref. [9], with the assumption that, having the same quantum numbers as the ψ(2S), the Y (4260) decays into J/ψπ 0 π 0 with the same branching ratio as what was reported in Ref. [46]. The Y (4260) resonance will probably have a very low branching fraction for decays into the J/ψπ 0 and π + π − π 0 channels, and will therefore not contribute significantly to the signal and to the other backgrounds.
D. Di-electron Continuum:pp → γ * π 0 → e + e − π 0 The production of a γ * with an associated π 0 is another process which can be used to demonstrate the universality of the TDAs [16]. When the invariant mass of the electronpositron pair is near the J/ψ mass, this channel represents a background source to thepp → J/ψπ 0 → e + e − π 0 channel, which can not be rejected in the analysis procedure since the particles in final state are identical to the signal. A similar TDA formalism to the one used here for the prediction of the J/ψπ 0 channel can also be used to estimate the differential cross section forpp → γ * π 0 → e + e − π 0 , as demonstrated in Ref. [16]. The same predictions can be used to integrate the cross section numerically over the corresponding validity domains at each collision energy (cf. Table I).
After setting the range 8.4 < Q 2 [GeV 2 ] < 10.1 to match the window around the J/ψ mass, which will be used for the analysis, cross sections of 13.6 fb, 21.6 fb and 24.8 fb are obtained at s = 12.3 GeV 2 , s = 16.9 GeV 2 and s = 24.3 GeV 2 , respectively. Taking into account the branching ratio of J/ψ → e + e − (5.94%), this results in contamination on the 10 −3 level and therefore has not been considered for further simulations.

E. Hadronic Decays of J/ψ
The reactionpp → J/ψπ 0 , with the J/ψ decaying into π + π − , where the π + π − pair is subsequently misidentified as a e + e − pair, is another potential source of background. It is highly suppressed by the branching fraction of J/ψ into π + π − (≈ 10 −4 ), and the low probability of misidentifying the pions as electrons (cf. Section V B).
Similarly,pp → J/ψπ 0 events with a hadronic J/ψ decay that can mimic the signal's final state (for example J/ψ → π + π − π 0 or J/ψ → γπ + π − ) is heavily suppressed by the probability to identify the π + π − pair. With the same production cross sections as the signal, and with very low probability of misidentifying the π + π − pair as an e + e − pair, such final states have negligible detection rates, and therefore will not be fully simulated.

V. SIMULATIONS AND ANALYSIS
For the present feasibility study, full MC simulations were performed on events from both the signal and background event generators described in the previous sections. In this section details of the analysis procedure will be provided. After a brief overview of the analysis methodology, a discussion of the PID and selection procedure, the reconstruction of π 0 and e + e − pairs, as well as the use of kinematic fits to gain further rejection for one class of background reaction where PID alone is not sufficient is provided. The section concludes with global signal to background ratios and signal purity for each background type included in the simulation study.

A. Brief Description of the Method
The analysis starts with the generation of signal and background events as described in Section III, followed by the transport of tracks in GEANT4, where the geometry of the PANDA detector has been implemented. The simulation of the detector response and digitization of the signals follows this step, at which point tracks and neutral particles can be reconstructed.
Signal events are then passed to the full event selection chain, including PID and analysis cuts. This ensures the most realistic description possible of the reconstructed signal, including statistical fluctuations. The number of signal events to simulate was picked to correspond exactly to the expected signal counts shown in Table II within the validity domain. By directly plotting spectra with tracks that pass all identification cuts, the plots will reflect the expected statistical accuracy.
For the background events, a different approach was followed. To check the feasibility of the measurement, the residual background contamination needs to be understood with a good precision in each bin used for the extraction of the physics observable (e.g. differential cross section distribution in t). This would not be possible if the cuts are directly applied due to the small number of background events that would pass all cuts. Therefore, a weight proportional to the product of single charged track misidentification probabilities is applied to each event instead of direct application of cuts as in the case of signal event simulation. The charged pion misidentification probability is parameterized as a function of momentum of the track. The photon selection is done by direct application of the cuts, just as in the case of the signal event analysis. Figure 8 shows the polar angle versus momentum distributions of reconstructed final state tracks in the background (top row) and signal (bottom row) simulations at the three different p incident momenta used in this study.

B. PID Efficiency
For the investigation of PID efficiency, single particle simulations of electrons, positrons, positive and negative pions were performed in the three kinematic zones depicted by red, black and blue boxes in Fig. 8. This is to focus the simulation to values of p and θ that matter most for this study. In total, 5×10 6 single e + and e − events each were used for the estimation of the electron efficiency. To obtain a good precision on the very low pion misidentification probability, a larger statistics of 25 × 10 6 each for single π + and π − events were used. The following section discusses in more detail how these efficiencies were determined. For the sake of simplicity in the remainder of the discussion, the word electron is used to refer to both electrons and positrons, except when the distinction is necessary.
As mentioned above, one of the critical aspects of this analysis is the reduction of the hadronic background using PID cuts, since the π + π − π 0 final state is kinematically very similar to the J/ψπ 0 → e + e − π 0 final state when the π + π − invariant mass is close to the J/ψ invariant mass. To this end, the PID should be used in the most effective way possible to maximize pion rejection until the cost to electron efficiency is prohibitive. This section describes the selection cuts of the analysis.

Combined PID Probability
The starting point for the application of global PID is the detector-by-detector PID probability information. To simplify the process of combining information from various detectors, the relevant PID variables from a given detector subsystem i ∈ {subsys} = {EMC, DIRC, DISC, STT, MVD} is used to determine an estimate of the probability P ID i ( j) that a given track is one of the five charged particles species that can be identified by PANDA, where j ∈ {e ± , µ ± , π ± , K ± , p ± }. The combined probability that a given track is of type j is then given by: where: Equation (17) ensures the proper normalization of the probabilities assigned for each track: ∑ j P ID comb ( j) = 1. With the combined probability in hand, a cut is applied at an appropriate value for the identification of any given species. In the present case, a sufficiently high threshold on P ID comb (e ± ) is imposed to select electrons and reject all other species. Here we show in Fig. 9 the global PID efficiency ε e ± eid (θ MC , p MC ) as a function of the true MC polar angle θ MC and momentum p MC with a cut of P ID comb (e ± ) > 0.9, based on a simulation of 10 7 electrons and positrons. Figure 10 shows the pion misidentification probability ε π ± eid (p MC ) as a function of p MC for the same cut based on a simulation of 5 × 10 7 charged pions.

Optimization of the Cut on
To determine an optimal cut on P ID comb (e ± ) for EID, the relationship between average efficiency of EID ε(e ± ) versus average misidentification probability for charged pions ε(π ± ) was studied as a function of the cut on P ID comb (e ± ), and plotted in Fig. 11 as a ROC (Receiver Operating Characteristics) curve, which shows the performance of the classifier as the discrimination threshold is varied. ε(e ± ) and ε(π ± ) are determined by taking the bin-by-bin weighted average of the electron efficiency shown in Fig. 9 and pion misidentification probability shown in Fig. 10. The weight for each bin is set to the content of the corresponding bin in the kinematic distributions shown in Fig. 8. One can observe that there is a significant gain in charged pion rejection with relatively small loss in EID efficiency up to a cut of P ID comb (e ± ) > 90%, beyond which the rejection gain no longer justifies the associated loss in efficiency. The cut of P ID comb (e ± ) > 90% is therefore chosen for the remainder of this analysis. It should be noted that the difference of the ROC curves between different incidentp momenta comes from the differences in the momentum and angular distributions of single tracks in the signal and background events that were used as a weight to average the  The application of the cut on P ID comb (e ± ) is straightforward for the signal simulations. The cut is applied on a track-bytrack basis to every reconstructed electron, and the spectra are constructed from those tracks that pass the cut. This approach is however not realistic for the background simulations, since the rejection is very high (≈ 10 −4 per charged pion). An unrealistically large number of background events need to be simulated to produce sufficient statistics in each bin after all cuts are applied. For this reason a different approach is adopted. The single pion misidentification probability shown in Fig. 10 is parameterized by a function f ε π ± (p) to smooth out the statistical fluctuation, and subsequently used as a weight for each background event based on the product of the values of f ε π ± (p) at the respective true MC momenta of the two identified pions, p π + and p π + : where N MC evt is the number of fullpp → π + π − π 0 reaction events that were simulated, and N BG is the number of background events that are expected from 2 fb −1 integrated luminosity based on the cross sections given in column 4 of Table II. The constant factor N BG /N MC evt ensures the proper normalization of the background spectra.

C. π 0 Reconstruction
Neutral pions are reconstructed through their two-photon decay channel, in the invariant mass spectrum formed by combining all photons within an event into γγ pairs. A cluster from the EMC with a minimum reconstructed energy of 3 MeV is considered to originate from a photon if there is no charged track candidate whose extrapolation to the EMC falls within a 20 cm radius from the EMC cluster. The invariant mass spectra show a contribution from combinatorial γγ pairs, which can be reduced by relying on the kinematic correlation of π 0 decay photons that the combinatorial γγ pairs do not display. Figure 12 shows the correlation of the reconstructed average photon energy to the opening angle in the lab frame between two photons, from all γγ pairs including those from combinatorial pairs on the left, and for γγ pairs that decay from a π 0 on the right, inpp → J/ψπ 0 → e + e − π 0 events simulated within the validity domains of the TDA model.
The following cuts are applied to the data: with: where OA is the opening angle between the two photons, E γ1 and E γ2 are the energies of the two photons, and a L i and a U i , with i = 0, 1, 2 are coefficients of the parametrization determined independently for each collision energy. Thus, it is possible to reduce the combinatorial background to a few percent while keeping an efficiency larger than 90% for pairs where both photons stem from π 0 decays. The effect of this cut is shown in Fig. 13 in a simulation ofpp → J/ψπ 0 → e + e − π 0 events, for all photon pairs on the left and for photons originating from π 0 decays on the right. In addition to this, an invariant mass cut of 110 < m γγ < 160 MeV/c 2 is applied on the two photon system.

D. J/ψ Reconstruction
The reconstruction of J/ψ candidates is accomplished by pairing all positive charged candidates passing EID cuts with all negative charged candidates that also pass the EID cuts. Figure 14 shows the invariant mass spectrum of e + e − pairs after application of the EID cuts, while requiring the presence of at least one reconstructed π 0 in the event. The mass distribution can be described satisfactorily by a Crystal Ball function [47]. A fit to the mass distribution has a peak at (3.088 ± 0.001) GeV/c 2 and a width of (51.3 ± 1.0) MeV/c 2 . The solid vertical lines show the mass window 2.8 < M e + e − < 3.3 GeV/c 2 which is used as a selection to reconstruct e + e − pairs from J/ψ.

E. J/ψπ 0 Signal Reconstruction
Finally, the full event is reconstructed by pairwise combining all reconstructed π 0 s with all J/ψ candidate e + e − pairs in the same event. Due to the presence of combinatorial background in both the π 0 reconstruction (random γγ pairs) and the J/ψ reconstruction (random e + e − pairs), there could be more than one candidate J/ψπ 0 pair per event. In such events, the angle between the π 0 and J/ψ in the CMS is calculated for each pair and the combination closest to 180 • (the most back-to-back pair) is selected. Figure 15 depicts the invariant mass distributions of the signal reaction as well as all the simulated background sources after the selection of the π 0 and e + e − pair. The sum of contributions (S+B) from signal (S) and all background sources (B) is shown in the same figure as the black histogram. Table IV shows the signal to background ratios of the different background sources simulated at this stage of the analysis, and the first column of Table V displays the efficiency of the signal and the rejection powers of different sources of background. The channels with a charged pion pair are sup-pressed with rejection powers of the order 10 7 . The main effect (roughly 10 6 ) comes from PID, while the remaining factor of 10 comes from the cut on the charged pair invariant mass. As a result, background events can be rejected to levels where they can be subtracted when needed in the e + e − invariant mass spectra by a sideband analysis, in which invariant mass regions to the right and left side of the J/ψ peak are used to estimate the background contribution under the peak. This is discussed in more detail in Section V H. As expected, the J/ψπ 0 π 0 channel is selected with an efficiency similar to the J/ψπ 0 channel. It is therefore now the dominant background source, roughly a factor four larger than the signal. This ratio is a direct result of the conservatively high cross section assumption discussed in Section IV C. A dedicated analysis of the J/ψπ 0 π 0 channel will be possible with PANDA, allowing for measurement of the cross section at the same c.m energy as the J/ψπ 0 signal. In the mean time, we stick to the cross section inputs as described in Section IV C, and propose further analysis cuts described in the following section to reduce this background to the percent level, keeping in mind that that they will later be adjusted to match realistic values of the J/ψπ 0 π 0 cross section.

F. Kinematic Fit
The J/ψπ 0 π 0 background will be reduced by exploiting the kinematic differences to the J/ψπ 0 channel. Figure 16 shows the χ 2 distributions of a kinematic fit with four-momentum conservation enforced as a constraint by assuming an exclusive γγe + e − event (denoted as χ 2 2γe + e − ). The plots are shown for the three beam momenta of this study. The insets show the same plot on a linear scale with a restricted range on the χ 2 axis, demonstrating that χ 2 2γe + e − is peaked at values compatible with the four degrees of freedom of the kinematic fit to thepp → J/ψπ 0 hypothesis. In contrast, the J/ψπ 0 π 0 events show a significantly flatter χ 2 2γe + e − distribution extending to very large values, providing a powerful tool for further rejection. As shown in the second column of Table V, by applying a maximum cut-off on χ 2 2γe + e − of 20, 50 and 100 at pp = 5.5, 8.0 and 12.0 GeV/c, respectively, it is possible to reduce the J/ψπ 0 π 0 contamination to less than 8%, while keeping the corresponding loss in signal efficiency to ≈ 15 -30%, depending on pp.
Despite the significant improvement of the S/B ratio, additional cuts are needed to bring the contamination by J/ψπ 0 π 0 to under a few percent at all energies. A comparison of the χ 2 value of the kinematic fit for the 2γe + e − hypothesis (χ 2 2γe + e − ) with the χ 2 value for the 4γe + e − hypothesis (χ 2 4γe + e − ) can provide additional discrimination power. The correlation between χ 2 2γe + e − and χ 2 4γe + e − is shown in Fig. 17. By requiring χ 2 4γe + e − > χ 2 2γe + e − for those reconstructed J/ψπ 0 events with an additional γγ pair in the event lowers the J/ψπ 0 π 0 contamination down to less than 5%, depending on pp, with no significant loss in signal efficiency.
Finally, by requiring that the number of photons in the event with energy above 20 MeV (N γ(>20MeV ) , shown in Fig. 18 for J/ψπ 0 and J/ψπ 0 π 0 events) to be less than or equal to three, it is possible to achieve a J/ψπ 0 π 0 contamination of less than 2% with an additional loss of efficiency of only about 10 -15%. Table V summarizes key quantities, signal efficiency, contamination from J/ψπ 0 π 0 and π + π − π 0 and purity from combinatorial background as well as S/B ratio including all background sources after each step in the analysis procedure described above, starting from the e + e − π 0 selection.

G. Signal to Background Ratio
The expected background contamination from all sources considered in this study is plotted in Fig. 19 as the shaded histogram for each validity range at the three beam momenta. In the worst scenario (at the lowest energy studied here at pp = 5.5 GeV/c), the S/B ratio will be at least of the order of 15 at all values of t. At higher energies, the background contamination frompp → π + π − π 0 is about or less than a percent, even with the conservative estimates of the background cross section used. As already mentioned, PANDA will provide dedicated measurements of these background channels, hence allowing for a subtraction of the corresponding contributions. In addition, background channels with no J/ψ in the final state can be suppressed using a sideband analysis of the Individual contributions from the various sources of background discussed in the text are also plotted with the line style depicted in the legend. π + π − π 0 , π + π − π 0 π 0 and π + π − π + π − π 0 are generated using DPM, whereas J/ψπ 0 π 0 is generated using the phase-space (PHSP) model. The combined contribution (S+B) of the signal (S) and all background (B) reactions is shown by the solid line. Events of the signal channel J/ψπ 0 are generated using the TDA model-based generator. invariant mass distribution of the charged pair, as will be discussed in the next section. This gives us confidence that the measurement should be readily feasible from the point of view of background rejection. V. The efficiency for signal events and rejection factor (R) of the two significant background contributions, together with background contamination (C ), signal purity from combinatorial background (P comb ) and total signal to background ratio. The quantities are tabulated with successive application of kinematic cuts starting from the selection of the π 0 and e + e − pair (column e + e − π 0 sel.). The signal and background count rates are determined within a M e + e − window of 2.8 to 3.3 GeV/c 2 in the invariant mass of the charged pair. The beam momentum for each table is given in the top left cell.

H. Sideband Background Subtraction
The number of reconstructed signal events is extracted by subtracting the background contribution from the total number of entries within the range 2.8 < M e + e − < 3.3 GeV/c 2 . The background contribution is obtained by integrating the polynomial component of the Crystal Ball (see Section V D) plus third order polynomial function fitted to the sum of signal and background e + e − invariant mass histograms. The fitted range was set to 2.5 < M e + e − < 3.8 GeV/c 2 . The parameter for the position of the maximum and the width of the Crystal Ball component were fixed to values extracted from the fit to the integrated J/ψ yield shown in Fig. 14, whereas the remaining parameters of the function were allowed to vary freely during fitting. The loss of signal events that fall outside of the counting window due to the Bremsstrahlung tail of the J/ψ peak is accounted for as an analysis cut efficiency loss in the correction procedure that will be outlined in Section V I.
The result of the subtraction procedure is summarized in Fig. 19 as a function of squared momentum transfer for the near-forward and near-backward approximations. For comparison, the count rates of e + e − pairs within the same J/ψ 2γe + e − and χ 2 4γe + e − for J/ψπ 0 (left) and J/ψπ 0 π 0 (right) events at beam momentum of pp = 5.5 GeV/c. mass window from the signal simulation (with no background added) are shown as open markers. The sideband subtraction procedure overestimates the signal count rate by roughly the amount of contamination that comes from J/ψπ 0 π 0 events, since it can only take into account background sources such as misidentified π + π − π 0 events that vary smoothly.

I. Efficiency Correction
This section describes a procedure for efficiency correction of the reconstructed signal count rate to obtain differential cross sections, and compare the result to the TDA model that was used to generate the signal events. The result will also serve as a guide to the statistical uncertainties expected for this measurement. For illustration, the t variable in the nearforward kinematic approximation is used, but the same arguments hold for the near-backward kinematic approximation if t is replaced by u. The fully corrected data points at a given value of t are calculated using the equation: where L int is the integrated luminosity (2 fb −1 ), and BR(J/ψ → e + e − ) indicates the branching ratio of the J/ψ → e + e − decay channel (5.94%). The quantity N J/ψπ 0 is the number of reconstructed J/ψπ 0 events in a particular bin in t, and ∆t is the width of the bin. The location of the point on the t axis is determined by the mean t value of all the entries in the given t-bin at the generator level. For the efficiency calculation, a separate simulation data set of the signal channel was used. The efficiency correction ε(t) in a given bin [t min , t max ] is defined as the ratio of the number of reconstructed events to the number of the generated events within that window: .
We note that in the numerator, the number of signal events is counted in a bin determined by the reconstructed value of t, denoted by t REC , whereas in the denominator the generated (true MC) value of t, denoted by t MC , is used. The signal reconstruction efficiencies calculated in this manner for the three incidentp momenta are shown in Fig. 20 for the full kinematic range accessible at each energy regardless of the validity range of the TDA model. The efficiency as a function of u is just the mirror image of this distribution, where the point of reflection sits midway on the full domain covered at the particular energy. Therefore, the efficiency for the backward π 0 emission can be deduced from this plot by looking at the most negative values of t. For clarity, the validity ranges of the TDA model on the t axis for both the near-forward and near-backward kinematic approximation regimes are shown as arrows at the bottom of the plot with a color that corresponds to the efficiency histogram. As far as count rates are concerned, this efficiency plot illustrates that our study, which is based on the TDA model, can be extrapolated to any other model. In particular, the reactionpp → J/ψπ 0 was recently studied in a hadronic model including intermediate baryonic resonances [45], where detailed predictions have been made for the center of mass energy dependence of the cross section and the π 0 angular distribution, which can be compared to the future PANDA data. While the efficiency in the nearforward and near-backward regimes for the lowest beam momentum (5.5 GeV/c) are more or less comparable, this is not true for the higher beam momenta, where the efficiency in the backward regime tails off to lower values. This is due to the increasing probability for one of the leptons from the J/ψ decay to fall in the forward spectrometer region outside the acceptance of the detectors included in this analysis.

VI. SENSITIVITY FOR TESTING TDA MODEL
A. Differential Cross Sections Figure 21 shows a comparison between the expected precision of the measured differential cross section for an integrated luminosity of 2 fb −1 to the prediction of Ref. [22] that was used as the basis for the signal event generator. The measurements have a satisfactory precision of about 8 -10% relative uncertainty. This level of precision will allow a quantitative test of the prediction of TDA models. The angular distribution of the e + and e − in the J/ψ reference frame constitutes a key observable that could allow to test the validity of the factorization approach. The leading twist description of the TDA model predicts a specific form procedure using the sideband method, as well as from uncertainties related to the estimation of residual background from J/ψπ 0 π 0 events, which will amount to about 1% in the pessimistic scenario assumed in this work. Finally, there will be a contribution to the systematic uncertainty coming from the calculation of the efficiency correction, but this can only be determined by how well the simulation describes the performance of the fully constructed detector. For the current study, we simulated 6 × 10 6 events to calculate the efficiency at each beam momentum; as a result, the statistical uncertainty on the efficiency is negligible.

VII. CONCLUSIONS
In this study, the feasibility of measuring J/ψπ 0 production inpp annihilation reactions with PANDA at the future FAIR facility was investigated. The study is based on the TDA model for the description of the signal. A combination of available data and the established hadronic event generator DPM was used to describe the background with only pions in the final state. Non-resonant J/ψπ 0 π 0 was also considered, for which a PHSP event generator was used for the event distributions and conservative estimates were made for the cross sections given the lack of constraints from existing data. Generated events were simulated and reconstructed using the Pan-daRoot software based on the proposed PANDA setup.
Using events selected to include a π 0 and an e + e − pair with invariant mass inside a broad window around the J/ψ mass, background reactions with no J/ψ in the final state can be reduced to the level of a few percent. The residual back-ground can be subtracted using the procedure outlined. The signal efficiency decreases from 18% for beam momentum of 5.5 GeV/c to 9% for 12.0 GeV/c. To reject the J/ψπ 0 π 0 background reaction to the 2% level, additional selection based on kinematic fits is needed, which further reduces the signal efficiency by about an additional 30%. The severeness of these cuts will however be adjusted depending on the measured yield for this background channel.
The cross section plots were produced with the assumption of a 2 fb −1 integrated luminosity which could be accumulated in approximately five months of data taking at each beam momentum setting with the full design luminosity of the FAIR accelerator complex. The resulting uncertainties to measure the differential cross section as a function of the four momentum transfer in the two validity regions is expected to be of the order of 5 -10%. In addition, the angular distribution of the leptons in the J/ψ center of mass frame provide important information to test the leading-twist approximation used in the TDA model. The distribution of these observables can also serve as a test of other models, as the recent hadronic model of Ref. [45], which will be particularly interesting outside the validity range of the TDA model at large emission angles.
We conclude that the measurement proposed here can be performed with PANDA, even at c.m. energies corresponding to resonances with a forbidden or weak decay to J/ψπ 0 . Together withpp → γ * π 0 → e + e − π 0 , the measurement of thē pp → J/ψπ 0 reaction will enable to test TDA models and opens exciting perspectives for the study of nucleon and antinucleon structure.