Study of two-photon decays of pseudoscalar mesons via $J/\psi$ radiative decays

Using a sample of $4.48\times10^{8}$ $\psi(3686)$ events collected with the BESIII detector at the BEPCII collider, we study the two-photon decays of the pseudoscalar mesons $\pi^0$, $\eta$, $\eta^\prime$, $\eta(1405)$, $\eta(1475)$, $\eta(1760)$, and $X(1835)$ in $J/\psi$ radiative decays using $\psi(3686)\to\pi^{+}\pi^{-}J/\psi$ events. The $\pi^0$, $\eta$ and $\eta^\prime$ mesons are clearly observed in the two-photon mass spectra, and the branching fractions are determined to be $B(J/\psi\to\gamma\pi^{0}\to3\gamma)=(3.57\pm0.12\pm0.16)\times10^{-5}$, $B(J/\psi\to\gamma\eta\to3\gamma)=(4.42\pm0.04\pm0.18)\times10^{-4}$, and $B(J/\psi\to\gamma\eta'\to3\gamma)=(1.26\pm0.02\pm0.05)\times10^{-4}$, where the first errors are statistical and the second systematic. No clear signal for $\eta(1405)$, $\eta(1475)$, $\eta(1760)$ or $X(1835)$ is observed in the two-photon mass spectra, and upper limits at the $90\%$ confidence level on the product branching fractions are obtained.


I. INTRODUCTION
Within the framework of quantum chromodynamics (QCD), the two-photon decay width of a meson plays a crucial role in understanding the nature of the meson and helps to distinguish glueballs from conventional mesons since glueballs are believed to have a relatively small twophoton decay width [1]. Therefore, experimental studies of the two-photon decays of mesons are very important to help interpret the meson spectrum.
The η(1405)/η(1475) pseudoscalar meson was once regarded as a glueball candidate since it was copiously produced in J/ψ radiative decays [2] and was not observed in two-photon collisions [3]. However, the measured mass is much lower than the prediction of Lattice QCD for a pseudoscalar glueball, which lies above 2GeV/c 2 [4,5]. Later, the experiments found two different pseudoscalar states, η(1405) and η(1475), with the former mainly decaying to a 0 (980)π and KKπ and the latter mainly to K * (892)K [6]. At present, the one state assumption and the nature of η(1405)/η(1475) are still controversial. Another pseudoscalar meson, the η(1760), has been proposed as a mixture of a glueball with a conventional qq state [7], rather than a pure qq meson or a glueball, and this hypothesis is supported by the large production rate of the η(1760) in J/ψ → γωω decays [8,9]. The nature of the X(1835) is still an open question although a number of theoretical interpretations have been proposed, including an NN bound state [10], baryonium with sizable gluon content [11], a pseudoscalar glueball [12], a radial excitation of the η ′ [13], and an η c -glueball mixture [14]. None of these interpretations have been completely ruled out or confirmed.

II. DETECTOR AND MONTE CARLO SIMULATION
BEPCII is a double-ring e + e − collider running at center-of-mass energies from 2.0 to 4.6 GeV. The BE-SIII [18] detector at BEPCII, with a geometrical acceptance of 93% of 4π solid angle, operates in a 1.0 T magnetic field provided by a superconducting solenoid magnet. The detector is composed of a helium-based drift chamber (MDC), a plastic-scintillator time-of-flight (TOF) system, a CsI(Tl) electromagnetic calorimeter (EMC) and a resistive plate chamber (RPC)-based muon chamber (MUC) in the iron flux return yoke of the magnet. The spatial resolution of the MDC is better than 130 µm, the charged-particle momentum resolution is 0.5% at 1.0 GeV/c, and the specific energy loss (dE/dx) resolution is better than 6% for electrons from Bhabha events. The time resolution of the TOF is 80 ps in the barrel and 110 ps in the endcaps. The energy resolution of the EMC at 1.0 GeV/c is 2.5% (5%) in the barrel (endcaps), and the position resolution is better than 6 mm (9 mm) in the barrel (endcaps). The position resolution in the MUC is better than 2 cm.
Monte Carlo (MC) simulations are used to estimate backgrounds and determine the detection efficiencies. The GEANT4-based [19] simulation software BOOST [20] includes the geometric and material description of the BESIII detector, detector response, and digitization models, as well as the tracking of the detector running conditions and performance. Production of the charmonium state ψ(3686) is simulated with KKMC [21,22], while the decays are generated with EVTGEN [23,24] for known decay modes with branching fractions taken from the Particle Data Group (PDG) [6] and by LUNDCHARM [25] for the remaining unknown decays. We use a sample of 5.06 × 10 8 simulated ψ(3686) events, in which the ψ(3686) decays generically ('inclusive MC sample'), to study the backgrounds. The analysis is performed in the framework of the BESIII offline software system (BOSS) [26] which incorporates the detector calibration, event reconstruction, and data storage.

III. DATA ANALYSIS
In this paper, the two-photon decays of the pseudoscalar mesons are investigated with J/ψ radiative decays. Hence the candidate events for the reconstruction of ψ(3686) → π + π − J/ψ, J/ψ → 3γ are required to have two oppositely charged tracks and at least three photon candidates. Each charged track, reconstructed using hits in the MDC, is required to be in the polar angle range | cos θ| < 0.93 and pass the interaction point within ±10 cm along the beam direction, and within ±1 cm in the plane perpendicular to the beam. Both charged tracks are assumed to be pion candidates.
Photon candidates are reconstructed from clusters of energy deposited in the EMC, and the deposited energy of each is required to be larger than 25 MeV in the barrel region (| cos θ| < 0.80) or 50 MeV in the endcap region (0.86 < | cos θ| < 0.92). The opening angle between a shower and the nearest charged track must be greater than 15 • , and timing requirements in the EMC are used to suppress electronic noise and energy deposits unrelated with the collision event. Events that satisfy the above requirements are retained for further analysis.
A four-constraint (4C) kinematic fit imposing energy and momentum conservation is performed under the hypothesis of π + π − γγγ. If the number of photon candidates in an event is larger than three, the combination with the smallest χ 2 4C from the kinematic fit is selected, and χ 2 4C is further required to be less than 50. The distribution of the γγγ invariant mass, M γγγ , of selected candidate events is shown in Fig. 1, where a very clean J/ψ peak is seen with very low background. A mass window requirement |M γγγ −m J/ψ | < 0.08 GeV/c 2 is applied to select the J/ψ signal, where m J/ψ is the nominal mass of the J/ψ meson [6].
After the above requirements, the distribution of the two-photon invariant mass M γγ is shown in Fig. 2, where the photon momenta from the 4C kinematic fit are used to calculate M γγ and there are three entries per event.
The backgrounds without the J/ψ intermediate state (non-J/ψ background) can be estimated from the events within the J/ψ sideband regions, defined  as 3.072 GeV/c 2 < M γγγ < 3.080 GeV/c 2 and 3.114 GeV/c 2 < M γγγ < 3.122 GeV/c 2 , which are indicated in Fig. 1. The backgrounds from ψ(3686) → π + π − J/ψ with J/ψ decaying to neutral particle final states (J/ψ background) are investigated with the inclusive MC sample of 5.06×10 8 ψ(3686) events. One prominent background is ψ(3686) → π + π − J/ψ, with J/ψ → γπ 0 π 0 , which produces a peak around the π 0 mass region in the M γγ distribution. To estimate its contribution, a dedicated MC sample of ψ(3686) → π + π − J/ψ, J/ψ → γπ 0 π 0 is produced incorporating the amplitude analysis result of J/ψ → γπ 0 π 0 [27]. With the same selection criteria and taking into account the number of ψ(3686) events as well as the branching fractions of ψ(3686) → π + π − J/ψ [6] and J/ψ → γπ 0 π 0 [27], the corresponding distribution of M γγ is shown as the solid histogram in Fig. 2. The number of peaking background events in the π 0 signal region is expected to be 32 ± 2, which is estimated by a fit to the γγ invariant mass spectrum of the above MC sample, where the π 0 signal is modeled with the sum of a Crystal Ball (CB) [28] function and a Gaussian function, and the other J/ψ non-peaking background is described with a second order Chebychev polynomial function.  The signal yields of J/ψ → γ(π 0 , η, η ′ ) → 3γ are obtained from unbinned maximum likelihood fits to the γγ invariant mass spectra. In the fits, the signal shapes are modeled with the sum of a CB function and a Gaussian function. The non-J/ψ backgrounds are estimated with the events in the J/ψ sideband region, assuming the M γγγ distribution to be flat in the vicinity of the J/ψ. Their yields and shapes are fixed in the fit. The nonpeaking J/ψ background is parameterized with a secondorder Chebychev polynomial function. The fit results are shown in Fig. 3. The signal yields from the fit and the MC determined detection efficiencies are summarized in Table I, where the MC simulation is performed using an angular distribution of 1+cos 2 θ γ for the radiative photon in the J/ψ rest frame.
No obvious signals for the pseudoscalar mesons η(1405), η(1475), η(1760) or X(1835) are observed in the M γγ distributions. Upper limits on the signal yields are obtained by fits to the M γγ distributions in the vicinity of the corresponding signal region, as shown in Fig. 4. In the fits, the line shapes of the η(1405), η(1475), η(1760) and X(1835) signals are parameterized by Breit Wigner (BW) functions convolved with Gaussian functions to account for the mass resolution, where the mass and width of BW functions are fixed to the world average values taken from the PDG [6] and the mass resolutions are obtained from MC simulation. The background shapes are described by second-order Chebychev polynomial functions. We derive the upper limits from these fits using a Bayesian approach with a flat prior as input. The distribution of minimized likelihood values for a series of input signal event yields is taken as the probability den-sity function (PDF) for the expected number of events. The number of events at 90% of the integral of the PDF from 0 to the given number of events is defined as the upper limit at the 90% confidence level (C.L.). To take into account the systematic uncertainties related to the fits, alternative fits with different fit ranges and background shapes are also performed, and the maximum upper limit among these cases is selected.

IV. SYSTEMATIC UNCERTAINTIES
Systematic uncertainties in the branching fraction measurements mainly originate from efficiency differences between data and MC simulation in the MDC tracking, the photon detection, the kinematic fitting efficiency and the J/ψ mass window requirement. Additional uncertainties associated with the fit range, the background shape, the sideband regions, the MC statistics, the branching fraction of ψ(3686) → π + π − J/ψ, and the total number of ψ(3686) events are also considered.
The tracking efficiency of charged pions has been investigated using control samples of J/ψ → ppπ + π − [29]. The difference in tracking efficiency between data and MC simulation is found to be 1% per track, which is taken as the uncertainty from the tracking efficiency. The photon detection efficiency is studied with a clean sample of J/ψ → ρ 0 π 0 [30]. The result shows that the difference of detection efficiency between data and MC simulation is 1% per photon.
The systematic uncertainties associated with the 4C kinematic fit are studied with the track helix parameter correction method, as described in Ref. [31]. In this analysis, we take the efficiencies with correction as the nominal values, and the differences with respect to those without corrections are taken as the systematic uncertainties associated with the 4C kinematic fit.
Due to the difference in the mass resolution between data and MC, the uncertainty related to the J/ψ mass window requirement is investigated by smearing the MC simulation in accordance with the signal shape of data. The change of the detection efficiency is assigned as the systematic uncertainty for the J/ψ mass window requirement.
To study the uncertainty from the fit range, the fit is repeated with different fit ranges, and the resultant largest differences in the signal yields are taken as the systematic uncertainties.
To estimate the uncertainty associated with the background shape, alternative fits with first-order or thirdorder Chebychev polynomial functions for the background are performed, and the maximum differences in signal yields with respect to the nominal values are taken as the systematic uncertainties.
The uncertainties from the J/ψ sideband region is estimated by using alternative sideband regions. The maximum differences in signal yields are taken as the uncertainties.
The uncertainty from the decay branching fractions of ψ(3686) → π + π − J/ψ is taken from the PDG [6], and the systematic uncertainty due to the number of ψ(3686) events is determined to be 0.7% according to Ref. [17]. Table II summarizes the systematic uncertainties from all sources for each decay. The systematic uncertainties associated with the statistics of MC samples are also included. The total systematic uncertainty is obtained by adding all individual uncertainties in quadrature, assuming all sources to be independent.
To estimate the upper limits on product decay branching fractions for un-observed pseudoscalar mesons, the systematic uncertainties are taken into consideration by convolving the PDF of likelihood values in each decay with a Gaussian function G(µ, σ) = G(0, N σ sys ), where N is the signal yield and σ sys is the corresponding relative systematic uncertainty listed in Table II. The upper limits on the number of events and the branching frac-