Search for the reaction channel $e^+e^- \rightarrow \eta_c\eta \pi^+\pi^-$ at center-of-mass energies from 4.23 to 4.60 GeV

Using data collected with the BESIII detector operating at the Beijing Electron Positron Collider, we search for the process $e^+e^-\rightarrow \eta_c\eta \pi^+\pi^-$. The search is performed using five large data sets recorded at center-of-mass energies of 4.23, 4.26, 4.36, 4.42, and 4.60 GeV. The $\eta_c$ meson is reconstructed in 16 exclusive decay modes. No signal is observed in the $\eta_c$ mass region at any center-of-mass energy. The upper limits on the reaction cross sections are determined to be 6.2, 10.8, 27.6, 22.6 and 23.7 pb at the 90% confidence level at the center-of-mass energies listed above.

The BESIII collaboration has reported a charmonium like state Z c (3900) ± decaying into J/ψπ ± in the reaction e + e − → J/ψπ + π − at √ s = 4.26 GeV [1].This resonance was also observed by the Belle experiment [2] and confirmed using CLEO-c data [3].In the CLEO-c dataset, evidence was found for the neutral state Z c (3900) 0 in e + e − → J/ψπ 0 π 0 [3], which was later also observed by BESIII [4].The Z c (3900) can not be explained as a conventional meson, because of its decay to charmonia and the existence of its charged state.
An enhancement in the D D * system at a mass of 3890.3MeV/c 2 has been observed in the reaction channel e + e − → π + (D D * ) − at a center-of-mass energy of 4.26 GeV, which may be identified as the Z c (3900) ± [5].Because of the observation of this second decay, interpretations favor the resonance to be either a tetra-quark state or a D-meson molecule [6].In addition, a further charged resonance Z c (4020) ± was found in the subsystem h c π ± of the reaction e + e − → h c π + π − , also at a center-of-mass energy of 4.26 GeV, closely followed by the discovery of its isospin partner the Z c (4020) 0 [7,8].Furthermore, structures whose poles are compatible with the Z c (4020) have been observed by the BESIII collaboration in the reactions e + e − → π + (D * D * ) − and e + e − → π 0 (D * D * ) 0 [9,10].
The observations of the isospin triplets Z c (3900) decaying to J/ψπ and Z c (4020) decaying to h c π suggest the possibility of an unobserved triplet of Z ±,0 c states decaying to η c π ±,0 and an isospin-singlet state decaying to η c η. Ref. [11] predicts a tetra-quark state in the mass region of 3.7 to 3.9 GeV with J P C = 0 ++ that would satisfy this latter hypothesis.The observation of this resonance would, therefore, add important information to the puz-zle of new states and would improve the understanding of their internal structure.
We search for the reaction e + e − → η c ηπ + π − , as any events observed for this process will allow for studies of possible resonant structure in the η c η subsystem.The ππ system must have a relative angular momentum of L=1 to conserve C-parity.It is expected that this pion decay proceeds mainly via the ρ resonance (vector dominance model).This leads to suppression of the decay channel due to isospin conservation and, in addition, a limited phase space below center-of-mass energies of 4.3 GeV.Any observed events will also allow for studies of possible resonant substructures in the η c π ± subsystem.

A. Detector and Monte Carlo Simulation
The BESIII detector [12] is located at the BEPCII double-ring e + e − collider.The detector consists of a helium-based multi-layer drift chamber (MDC) with a momentum resolution of 0.5% for charged particles with a transverse momentum of 1 GeV/c, a plastic scintillator based time-of-flight (TOF) system with a time resolution of 68 ps in the barrel and 110 ps in the end caps, a CsI(Tl) electromagnetic calorimeter (EMC) with energy resolutions of 2.5% and 5.0% for 1 GeV photons in the barrel and end caps respectively, and a multilayer resistive-plate chamber muon-detection system.The BESIII detector operates in a 1 T magnetic field provided by a superconducting solenoid and has a geometrical acceptance of 93%.
To optimize selection criteria, estimate detector resolution and reconstruction efficiency, Monte Carlo (MC) simulations are used.The simulation of the BESIII de-tector is based on GEANT4 [13] which models the interaction of particles with the detector material.The initial interaction of the e + e − system is modeled with KKMC [14] generator which also handles initial state radiation.Subsequent particle decays are generated with EVTGEN [15].The generation of final state radiation is handled by PHOTOS [16].In the simulations the signal reaction channel e + e − → η c ηπ + π − is generated according to a phase-space distribution.The η c is reconstructed in the following 16 decay modes: pp, ppπ 0 and ppπ + π − .The branching ratios for these decays are taken from Ref. [17].The branching ratios of the decays π 0 → γγ, η → γγ, and K 0 S → π + π − are taken from the Particle Data Group (PDG) [18].For the optimization of the suppression of background reactions various simulated data sets are used, e.g.samples containing light quark and open charm and charmonium final states as well as e + e − or µ + µ − MC samples.

B. Data Analysis
The search for the reaction is performed at five different center-of-mass energies.The integrated luminosity of these data sets is given in Table I.During event reconstruction, the charged tracks are required to have a point of closest approach to the interaction point within a cylinder with a radius of V xy = 1 cm in the x-y plane and a length of V z = ±10 cm along the beam axis.In addition, the polar angle with respect to the beam axis has to be in the acceptance of the MDC, corresponding to |cos(θ)| < 0.93.For tracks originating from the decay of long lived particles like the K 0 S meson, the V xy requirement is omitted while V z is increased to ±20 cm.For each event, the net charge of all reconstructed tracks has to be zero.For particle identification, the joint probability from the energy loss in the MDC and the time-of-flight information of the TOF system is calculated for each particle species (π, K, p) and compared for the selection.Photons are reconstructed from clusters in the EMC.To suppress noise in the EMC, the reconstructed photon energy has to be larger than 25 MeV for |cos(θ)| < 0.80 and 50 MeV for 0.84 < |cos(θ)| < 0.92.Furthermore, the time difference between the event-start time and the EMC timing information has to be 0 ns ≤ t EMC ≤ 700 ns.To suppress clusters formed by split-off photons from charged particle tracks, the angle between a cluster in the EMC and the closest charged track has to be at least 10°.The number of reconstructed photons has to be at least equal to the number of photons expected for the final state in question.
π 0 and η mesons are reconstructed from combinations of photon pairs.To select π 0 candidates, the invariant mass of two photons must satisfy 110 MeV/c 2 ≤ m γγ ≤ 150 MeV/c 2 while for η candidates the invariant mass has to be in range 500 MeV/c 2 ≤ m γγ ≤ 570 MeV/c 2 .Candidate K 0 S mesons are reconstructed by applying a vertex fit to all pairs of oppositely charged particles assuming a pion hypothesis, but requiring no particle identification criteria.For these pairs the decay length L and its uncertainty σ L are calculated from the decay vertex and the primary vertex position.The pair is kept as a of the fit is smaller than 100.
In addition, the decay length of the K 0 S candidate has to satisfy L/σ L > 2. Finally, the invariant mass of the pion pair must lie within 15 MeV/c 2 of the nominal K 0 S mass.A vertex fit is applied to events passing these criteria, excluding the tracks originating from K 0 S candidates.In the cases that the vertex fit converges, a kinematic fit is performed to improve the momentum resolution.The fit is constrained by the initial four momentum of the e + e − pair and a mass constraint on the η mass.If the final state of the η c contains additional π 0 , η or K 0 S mesons, mass constraints are applied on the invariant masses of their daughter particles as well.The selection criteria on the χ 2 value from the kinematic fit is used to suppress poorly reconstructed events and is chosen for each final state to retain 90% of the signal events.The kinematic fit is not able to discriminate between pions from the initial reaction and pions from the subsequent η c decay as the total four momentum is identical for these two hypotheses.This can lead to multiple candidates per event for the whole reconstruction, with each candidate having the same χ 2 .In these cases all candidates are kept for further analysis.It was checked with signal MC data sets that these candidates which have the wrong assignment do not contribute to the signal yield as the form a smooth background distribution.Also all other background distributions show a smooth behaviour at the signal region.

C. Cross Section and Upper-Limit Determination
To determine the total event yield of the reaction channel e + e − → η c ηπ + π − a simultaneous extended maximum-likelihood fit is performed on the invariantmass distribution m of the η c candidates for all η c decay modes.The fit function f fit is given by for each η c decay mode separately.The signal shape for each decay mode is described by a Breit-Wigner function B, whose parameters are fixed to the nominal mass µ PDG ηc and width Γ PDG ηc of the η c meson from the PDG [18].This function is convolved with a Gaussian function G with mean µ and standard deviation σ to account for the detector resolution, which is extracted from signal MC simulation.The number of background events (combinatorial and physical), n b , for each η c decay mode is determined simultaneously in the fit.For the majority of η c decay modes the background is described by a n th -order Chebychev polynomial function where the single terms T k are weighted by the coefficients a k .For certain decay modes (η c → pp, K + K − η, K + K − π 0 and K s Kπ) it is found that an exponential background function provides a better description of the background distribution.The number of signal events n s in each η c decay mode is related to the cross section σ via the relation where L is the integrated luminosity, B(η c → X) the branching ratio of η c decaying to X, and ε tot,X the corresponding reconstruction efficiency, which is obtained by fitting the reconstructed η c invariant-mass distribution from signal MC simulation.Table I shows the sum over all η c final states of the products of the efficiency and branching ratios.The invariant mass distribution summed over all decay modes is shown in Fig. 1 together with the sum of the fitting curves.The resulting values for the observed cross section can be found in Table II.
The uncertainties are purely statistical and obtained by a likelihood scan using the MINOS tool [20].
As no significant signal is observed, an upper limit on the cross section is calculated.For this calculation a Bayesian approach is used.For the prior distribution π(σ), we assume that it is zero for negative values of the cross section and follows a flat distribution otherwise.With this assumption the upper limit is given by where L is the likelihood function of the simultaneous fit as depicted in Fig. 2. The derived upper limits at 90% confidence level can also be found in Table II.

D. Systematic Uncertainties
There are several sources of possible systematic bias in the analysis, for which uncertainties are assigned.These originate in discrepancies in the detector description between MC simulation and data, the knowledge of the η c branching ratios, the knowledge of the resonance parameters of the η c , the kinematic fit, and the background model and fit range used in the simultaneous fit.The systematic uncertainties are summarized in Table III.
The uncertainty associated with the understanding of the track reconstruction in the MDC is studied with the decays J/ψ → ppπ + π − and J/ψ → ρπ and is found to be 1% per charged track [21].An additional 1% per track is applied to account for the knowledge of the particle identification performance [22].The systematic uncertainty on the photon detection is estimated using the control samples ψ(3680) → π + π − J/ψ with J/ψ → ρ 0 π 0 and is determined to be smaller than 1% for each photon [23].The systematic uncertainty on the K 0 S reconstruction is estimated to be 1.2% using the control samples J/ψ → K * (892) ± K ∓ with K * (892) ± → K s π ± and J/ψ → φK s K ± π ∓ [24].The systematic uncertainty associated with π 0 and η reconstruction is estimated to be 1% per π 0 /η, following the studies reported in Ref. [25] using the control samples J/ψ → π + π − π 0 and J/ψ → ηpp.The influence of these uncertainties on the cross section extracted from the simultaneous fit is estimated by multiplying the reconstruction efficiency of each η c final state X with a correction factor α X , which is given by Each κ Y (with Y = T for tracks, γ for photons, π 0 /η and K 0 S for the reconstructed mesons) follows a Gaussian distribution centered at one and a width set to the corresponding uncertainty, while n Y is the number of tracks, photons etc. in each final state X.The simultaneous fit is performed 1000 times while changing the values for κ Y for each fit.The width of the resulting distribution normalized to the extracted cross-section is taken as the systematic uncertainty of the reconstruction efficiency.
The η c branching ratios entering the simultaneous fit are derived from the BESIII measurement in Ref. [17], using the following relation Here the branching ratio B(ψ(3680) → π 0 h c ; h c → γη c ) is obtained by combining two measurements performed by BESIII [26] and CLEO [27].To estimate the systematic uncertainty of the branching ratios for the η c final states a random number is drawn from a Gaussian distribution whose width is set to the total uncertainty of the combined measurement of the common denominator, and one for each of the 16 modes in the numerator separately.The simultaneous fit is performed again with the updated branching ratios.This is repeated 1000 times and the width of the obtained cross-section distribution normalized to the extracted cross-section is taken as systematic uncertainty.
During the simultaneous fit the mean and width of the signal Breit-Wigner distributions are fixed to the values given by the PDG [18].To account for the uncertainties on these values, 1000 fits are performed in which new values for the mean and width of the η c are randomly generated from two independent Gaussian probability distributions, with the parameters of these distributions set according to the central values and uncertainties of the PDG.The standard deviations of the resulting cross-section distributions are assigned as a systematic uncertainty in the measurement.
To improve the agreement of the χ 2 distribution of the kinematic fit between signal MC and data, the helix parameters of the charged tracks are smeared using the decay J/ψ → φf 0 (980).The systematic uncertainty associated with the kinematic fit is estimated by switching off this correction, repeating the simultaneous fit and assigning the difference in the cross sections as the uncertainty.
The influence of the mass range over which the fit is performed is studied by narrowing and increasing the range of the fit by 5 MeV/c 2 .The systematic uncertainty is calculated by taking the maximum difference of the nominal fit value and the values obtained by varying the fit range.
The systematic uncertainty associated with the description of the background shape is estimated by increasing the order of the Chebychev polynomials by one.For those cases where the baseline description of the background is an exponential function, this is replaced by a second-order Chebychev polynomial.The difference between the cross section obtained with the new fit and that with the nominal background model is taken as the systematic uncertainty.
The integrated luminosity is determined by using Bhabha events.The systematic uncertainty of the luminosity measurement has been studied in Ref. [28] and a relative uncertainty of 1% is assigned for each centerof-mass energy.
To include the systematic uncertainties into the calculation of the upper limits, the likelihood is folded with a Gaussian distribution with a width set to the size of the systematic uncertainties The likelihood graphs from this procedure are shown in Fig. 2. The upper limits for the observed cross sections including the systematic uncertainty are listed in Table II.
To obtain upper limits for the Born cross sections σ Born , the observed cross sections have to be corrected for initial-state radiation (ISR) and vacuum polarization.The equation for the number of signal events in each η c decay mode then reads where ε X ,tot is the total efficiency for final-state X, L is the integrated luminosity, δ ISR and δ vp are the correction factors for initial-state radiation and vacuum polarization, respectively.The ISR correction factor is given by Here x is the fraction of the beam energy carried away by the ISR photon, ε(x) the corresponding reconstruction efficiency, σ(x) is the line shape of a single resonance, which is assumed to have Breit-Wigner shape in the calculations, and σ 0 and ε 0 are their counterparts in the absence of initial-state radiation.W (x) is the so-called radiator function [29].The value of δ ISR has a strong dependence on the parameters of the Breit-Wigner line shape, with the correction being largest for narrow resonances.As no resonances are observed, we make the conservative assumption that any resonance present has a width of 10 MeV/c 2 which is well below the measured parameters of, for example, the Y(4260).For the determination of the upper limit the most conservative approach is taken by assuming this small resonance is located such that the correction factor is largest.The value of δ ISR is estimated to be 0.64, independent of the collision energy and the η c final state.This is shown in Table IV, together with the values of δ vp , which are energy dependent and calculated with alphaQED [30].The upper limits for the Born cross section are given in the right column of Table II   GeV (e).The interval correponding to the upper limit at 90% confidence level is indicated as gray area.

FIG. 1 .
FIG. 1. Result of the simultaneous fit to 16 ηc decay modes.Shown is the sum of all modes at √ s of 4.23 GeV (a), 4.26 GeV (b), 4.36 GeV (c), 4.42 GeV (d), and 4.60 GeV (e).Black points are data, blue line is the sum of the fitting functions.

FIG. 2 .
FIG. 2. Likelihood curves convoluted with the Gaussian function representing the systematic uncertainties as a function of the cross section at √ s of 4.23 GeV (a), 4.26 GeV (b), 4.36 GeV (c), 4.42 GeV (d), and 4.60GeV (e).The interval correponding to the upper limit at 90% confidence level is indicated as gray area.

TABLE I .
[19]grated luminosity of the used data samples and sum over all ηc final states of the products of the efficiencies and branching ratios at the different center-of-mass energies.The center-of-mass energies are taken from Ref.[19].

TABLE II .
Observed cross section σ and upper limits (UL) for the reaction e + e − → ηcηπ + π − at the five center-of-mass energies.3rd column: UL without systematic uncertainties, 4th column: UL with systematic uncertainties, 5th column: UL with systematic uncertainties plus ISR and vacuum polarization correction.

TABLE IV .
[7,31]es for the vacuum polarization and ISR corrections for the different data sets.Calculations of the vacuum polarization correction are based on alphaQED[30].We perform a search for the process e + e − → η c ηπ + π − at √ s = 4.23, 4.26, 4.36, 4.42, and 4.60 GeV with data collected by the BESIII detector.The cross section at each center-of-mass energy is determined by a simultaneous fit to the invariant mass of the η c meson for 16 decay modes.The observed cross sections are determined to be σ 4.23 GeV = −5.39+3.15−2.83±2.05 pb σ 4.26 GeV = −0.98 +4.11 −3.53 ± 1.67 pb σ 4.36 GeV = 8.59 +6.72 −6.03 ± 2.04 pb σ 4.42 GeV = 3.07 +5.36 −5.12 ± 5.83 pb σ 4.60 GeV = 3.16 +6.91 −6.51 ± 4.19 pb.These upper limits are of the same order of magnitude as the measured cross sections of the processes e + e − → J/ψπ + π − and e + e − → h c π + π −[7,31].As no significant e + e − → η c ηπ + π − signal is seen in the current data set it is not yet possible to conclude about possible resonant structures in the final-state subsystems.