Evidence for $Z_{c}^{\pm}$ decays into the $\rho^{\pm} \eta_{c}$ final state

We study $e^{+}e^{-}$ collisions with a $\pi^{+}\pi^{-}\pi^{0}\eta_{c}$ final state using data samples collected with the BESIII detector at center-of-mass energies $\sqrt{s}=4.226$, $4.258$, $4.358$, $4.416$, and $4.600$~GeV. Evidence for the decay $\zcpm\to\rhopm\etac$ is reported with a statistical significance of $3.9\sigma$ with various systematic uncertainties taken into account at $\sqrt{s} = 4.226$~GeV, and the Born cross section times branching fraction $\sigma^{B}(\EE\to \pimp\zcpm)\times \BR(\zcpm\to\rhopm\etac)$ is measured to be $(48 \pm 11 \pm 11)\,\rm{pb}$. The $\zcpm\to \rhopm\etac$ signal is not significant at the other center-of-mass energies and the corresponding upper limits are determined. In addition, no significant signal is observed in a search for $\zcppm\to \rho^{\pm}\etac$ with the same data samples. The ratios $R_{\zc}=\BR(\zcpm\to \rho^{\pm} \etac)/\BR(\zcpm\to \pi^{\pm} \jpsi)$ and $R_{\zcp}=\BR(\zcppm\to \rho^{\pm} \etac)/\BR(\zcppm\to \pi^{\pm} \hc)$ are obtained and used to discriminate between different theoretical interpretations of the $\zcpm$ and $\zcppm$.

It has recently been suggested that the relative decay rate of Z c s with the same charge, such as Z ±,0 c (3900) → ρ ±,0 η c to π ±,0 J/ψ (or Z ±,0 c (4020) → ρ ±,0 η c to π ±,0 h c ), can be used to discriminate between the molecule and tetraquark scenarios [11]. In the tetraquark scenario, the predicted ratio of B(Z c (3900) → ρη c )/B(Z c (3900) → πJ/ψ) is 230 or 0.27, depending on whether or not the spin-spin interaction outside the diquarks is kept. In the meson molecule framework, on the other hand, this ratio is only 0.046. Similarly, the predicted ratio of B(Z c (4020) → ρη c )/B(Z c (4020) → πh c ) is 6.6 in the tetraquark model, but only 0.010 in the meson molecule model [11]. Therefore, a search for the decays of Z c (3900) or Z c (4020) to ρη c offers an important opportunity to understand their internal structure.
The design and performance of the BESIII detector are given in Ref [12]. The GEANT4-based [15] Monte Carlo (MC) simulation software package, which includes the geometric description of the BESIII detector and the detector response, is used to optimize event selection criteria, determine the detection efficiencies, and estimate the backgrounds. At each energy, the signal events are generated according to phase space using EVTGEN [16]. The ISR is simulated with KKMC [17], and final state radiation is handled with PHOTOS [18].
The charged tracks, photons and K 0 S candidates are reconstructed with nominal criteria used at the BESIII experiment and can be found in Ref. [19]. To reconstruct π 0 and η, the invariant mass of a photon pair is required to be in the range [0.120, 0.145] GeV/c 2 for π 0 and [0.50, 0.57] GeV/c 2 for η, which are approximately equivalent to a ±3σ window around the nominal mass of π 0 or η. To improve the resolution, a one-constraint (1C) kinematic fit is imposed on the selected photon pairs to constrain their invariant mass to the nominal π 0 or η mass [20]. The entire decay sequence is fully reconstructed for each η c decay channel. After the above selection, a four-constraint (4C) kinematic fit is performed for each event, and the χ 2 of the fit (χ 2 4C ) is required to be less than 40 to suppress backgrounds. In the kinematic fit, all charged tracks are assigned to the assumed particle species if there is only one type of particle in the final state. Otherwise, each track is assigned to all possible hypotheses. If there is more than one combination in an event, the one with minimum χ 2 min ≡ χ 2 4C + χ 2 1C + χ 2 PID + χ 2 vertex is kept for further analysis. Here, χ 2 1C is the χ 2 of the 1C fit for π 0 (η), χ 2 PID is the sum of the χ 2 for the PID of all charged tracks, and χ 2 vertex is the χ 2 of the K 0 S secondary vertex fit. Inclusive MC samples with the same statistics as the data are studied to understand the backgrounds. The major backgrounds to e + e − → π + π − π 0 η c are classified into two categories. They are events from (1) charmonium(-like) states decays (most of which include open-charm decays, e.g. ψ → D ( * )D( * ) ); and (2) the continuum process, e + e − → qq, with q = u, d, and s. It is found that about two thirds of the background events originate from the continuum process, and the rest are from resonance decays.
Background events with charmed mesons are rejected if a D meson candidate is reconstructed in one of its five major decay modes: To accomplish this, we require the invariant mass of D 0 (D ± ) candidates to be outside the region m(D 0 ) ± 24 MeV (m(D ± ) ± 10 MeV). Events with a K * (892) → Kπ, an ω → π + π − π 0 , or an η → π + π − π 0 candidate are removed to reduce the continuum background by requiring , m(ω) and m(η) are the nominal masses of the corresponding states.
By analyzing 600,000 e + e − → π + π − h c MC simulation events with h c decaying inclusively, a small enhancement in the η c signal region is found. Its contribution, N peaking bkg , is estimated to be 8.7 ± 2.0 at √ s = 4.226 GeV using results from Ref [4]. The contributions at other energies are estimated in a similar way. Except this mode, no peaking background is found.
The mass windows for the background veto and the χ 2 requirement of the 4C kinematic fit are determined by optimizing the figure-of-merit (FOM), which is defined as FOM = S/ √ S + B. Here, S is the number of signal events from the MC simulation (assuming σ(e + e − → π + π − π 0 η c ) = 50 pb, which is evaluated from a measurement using unoptimized selection criteria), and B is the number of background events obtained from the η c sidebands in the data. The optimization is performed through iterations until all the selection criteria become stable.
To obtain the π + π − π 0 η c yield, the invariant mass distributions of the η c candidates in the nine decay modes are fitted simultaneously using an unbinned maximum likelihood method. In the fit, both the signal and background shapes are channel dependent but the relative signal yields among all the channels are constrained by branching fractions and efficiencies. The η c signal is described with a constant width Breit-Wigner function (mass and width are fixed to the world average values [20]) convolved with instrumental resolution determined from the MC simulation. The background is parameterized with a second order Chebyshev Polynomial (CP) function. The total signal yield of the nine channels is N obs , which is shared for all the channels and required to be positive. N obs × f i is the signal yield of the ith channel. Here, the weight factor is f i ≡ ε i B i / i ε i B i , where B i denotes the branching fraction of η c decays to the ith mode [21] and ε i represents the corresponding efficiency. The free parameters in the fit include N obs and the background parameters for each decay mode. Figure 1(left) shows the fit results at √ s = 4.226 GeV projected onto the sum of events from all nine η c decay modes. Figure 1(right) shows the background subtracted distribution. The total signal yield is 333 +83 −80 with a statistical significance of 4.2σ, which is obtained by comparing the change of the log-likelihood value ∆(− ln L) = 9.0 and degrees of freedom ∆dof = 1 with and without the π + π − π 0 η c signal in the fit. The same selection criteria are applied to the data sets at √ s =4.258, 4.358, 4.416, and 4.600 GeV, but no significant signals are observed.
The Born cross section of the e + e − → π + π − π 0 η c reaction is calculated using where N sig = N obs − N peaking bkg is the number of signal events after the peaking background subtraction; L is the integrated luminosity; (1 + δ) is the ISR correction factor, assuming the π + π − π 0 η c signal is from Y (4260) decays [20]; and 1 |1−Π| 2 is the vacuum-polarization factor [22]. The cross sections and the numbers used in its calculation are listed in Table I for all energy points. The upper limits of the cross sections at 90% confidence level (C.L.) are determined using a Bayesian method, assuming a flat prior in σ B . The systematic uncertainties are incorporated into the upper limit by smearing the probability density function of the cross section [19]. The corresponding results for σ B U.L. are also listed in Table I. The Z c (3900) ± (Z c (4020) ± ) signals are examined after requiring that the invariant mass of an η c candidate is within the η c signal region [2.95, 3.02] GeV/c 2 and the invariant mass of π ± π 0 is within the ρ signal region [0.675, 0.875] GeV/c 2 . Here, all possible combinations in one event are kept to avoid bias. To suppress the combinatorial background, the momenta of the pions from ρ decays are required to be less than 0.8 GeV/c. The events in the η c and ρ sidebands are used to estimate the background. The recoil mass spectrum of the remaining π ∓ is shown in Fig. 2 for the data at √ s = 4.226 GeV, together with the contribution from the η c and ρ sideband events. Here, the ρ sideband is defined as [0.475, 0.675] GeV/c 2 . In Fig. 2, the Z c (3900) ± signal is apparent, but there is no statistically significant Z c (4020) ± signal.
To obtain the yields of e + e − → π ∓ Z c (3900) ± → π ∓ ρ ± η c and e + e − → π ∓ Z c (4020) ± → π ∓ ρ ± η c , the invariant masses of ρ ± η c candidates in the nine η c decay channels are fitted simultaneously using the same method as for e + e − → π + π − π 0 η c . In the fit, a possible interference between the signal and the background is neglected. The mass and width of the Z c (3900) ± are fixed to the values from the latest measurement [23] and those of the Z c (4020) ± are fixed to world average values [20]. The mass resolution is obtained from MC simulation and parameterized as a Crystal Ball function [24]. The background is described with a second order CP function. Figure 2(left) shows the fit to the invariant mass of ρ ± η c summed over the nine η c decays at √ s = 4.226 GeV. Figure 2(right) is the background subtracted distribution. The total Z c (3900) ± signal yield is 240 +56 −54 events with a statistical significance of 4.3σ, and that of the Z c (4020) ± is 21 +15 −11 events with a statistical significance of 1.0σ. The signals at the other c.m. energies are not statistically significant.
The Born cross section for e + e − → π ∓ Z ± c with Z ± c → ρ ± η c is calculated using where N Zc obs is the total signal yield of the Z c and ε Zc i is the corresponding detection efficiency. The definitions of other factors are the same as in Eq. (1). The numbers used in the calculation and the results are listed in Table II. The systematic uncertainties in the σ B (e + e − → π + π − π 0 η c ) measurement originate from the uncertainty of each factor in Eq. (1). The integrated luminosity has an uncertainty of 1.0% [13]. The number of e + e − → π + π − h c , h c → γη c peaking background events is estimated with MC simulation of e + e − → π + π − h c , so both the cross section uncertainty and the statistical error of the peaking background from the MC sample are taken into account to estimate this uncertainty. For the ISR correction, we assume that the e + e − → π + π − π 0 η c events mainly come from the Y (4260). Alternatively, the c.m. energy dependent cross section of e + e − → π + π − J/ψ measured by the BESIII experiment [25] is used to estimate this uncertainty. The uncertainty from the signal shape consists of the mass resolution discrepancy between data and MC simulation and the uncertainty of the η c resonant parameters. The former is studied using an e + e − → γ ISR J/ψ [26] sample and the latter is estimated by varying the η c parameters by ±1σ around the world average values [20]. The uncertainty for the background shape is estimated by changing the order of the CP function. The uncertainty in the fitting range is obtained by adjusting the boundaries. The methods to estimate the uncertainties due to the vacuum polarization and i ε i B i are the same as described in Ref. [19]. Furthermore, when we estimate the uncertainty of i ε i B i , the uncertainty due to the e + e − → π + π − π 0 η c decay dynamics is obtained by comparing the simulations with and without the Z c resonance. All of the sources are assumed to be independent and added in quadrature and the largest contribution to the systematics is the uncertainty of i ε i B i . The total systematic uncertainties are listed in Table I. For the σ B (e + e − → π ∓ Z c (3900) ± (Z c (4020) ± ) → π ∓ ρ ± η c ) measurement, the uncertainties on L, ISR factors, εB and the vacuum polarization factor are studied following the methods described in the measurement of σ B (e + e − → π + π − π 0 η c ). Moreover, additional sources of systematic un-   ) for the e + e − → π + π − π 0 ηc process and the numbers that enter the calculation (see Eq. (1)). σ B U.L. is the upper limit of the cross section at the 90% C.L., and S is the statistical significance of the signal.  TABLE II. Born cross sections of e + e − → π ∓ Zc(3900) ± (Zc(4020) ± ) → π ∓ ρ ± ηc (numbers for Zc(4020) ± are in brackets). The parameters are defined in the same way as those in Table I. certainties come from the ρ and η c selections, and the fit to the Z c (3900) ± (Z c (4020) ± ) resonances. The uncertainty due to the M (π ± π 0 ) mass window is estimated by comparing the invariant mass of M (ω → π + π − π 0 ) in data and MC assuming the mass resolution of M (π + π − π 0 ) is larger than M (π ± π 0 ). The discrepancy is found to be negligible. The uncertainty of the η c line shape is estimated by the variation of the mass and width of the η c within the errors given by world average values [20]. The uncertainties due to the fit to the Z c (3900) ± (Z c (4020) ± ) are estimated with the same methods as in the π + π − π 0 η c case. All these sources and those in the σ B (e + e − → π + π − π 0 η c ) measurement (except the uncertainties of the fit) are assumed to be independent and added in quadrature, and the total systematic uncertainties are listed in Table II. To consider the effect of the systematic uncertainty on the signal significance at √ s = 4.226 GeV, we vary the signal shape, background parametrization, and fit range, or free the Z c mass in the fit. We find the statistical significance of the Z c (3900) is always larger than 3.9 σ in all the checks.
In summary, we report the first evidence for the ρη c decay mode of the charged charmonium-like state Z c (3900) ± in the e + e − annihilation data at √ s = 4.226 GeV, and measure the cross section times branching ratio σ B (e + e − → π ∓ Z c (3900) ± ) × B(Z c (3900) ± → ρ ± η c ) = (48 ± 11 ± 11) pb. This result is very close to the cross section of e + e − → π + π − π 0 η c , which is (46 +12 −11 ± 10) pb. This indicates that the e + e − → π + π − π 0 η c process is dominated by the process e + e − → π ∓ Z c (3900) ± → π ∓ ρ ± η c . The statistical significance of Z c (3900) ± → ρ ± η c is 4.3σ (3.9σ including the systematical uncertainty). No signal is observed at √ s =4.258, 4.358, 4.416, and 4.600 GeV and the upper limits of the production cross sections at 90% C.L. are determined. No significant signal of Z c (4020) ± → ρ ± η c is found in any of the data sets and the upper limits of the production cross sections are determined.
The measured R Zc(3900) is closer to the calculation of the tetraquark model than that of the molecule model [11]. The measurement is also consistent with several other independent calculations based on the tetraquark scenario [27][28][29][30]. As for the molecule model, we notice that the calculated R Zc(3900) is highly model dependent and varies from 6.7 × 10 −3 to 1.8 in different approaches [30][31][32][33][34]. Therefore, it is neccessary to narrow down the theoretical uncertainty in the molecular framework to have a better comparison with the measurement. In the hadron-charmonium model, the B(Z c (3900) → ρη c ) is suppressed compared with B(Z c (3900) → πJ/ψ) and therefore inconsistent with the measurement [35]. Furthermore, this model predicts a new resonance W c (3785), which can be produced via e + e − → ρW c → ρπη c , the same final state we analyzed here. As we found that the e + e − → π + π − π 0 η c process is saturated by the e + e − → πZ c (3900) → ρπη c , the production of the W c is minor compared with e + e − → πZ c (3900) if it exists.
For R Zc(4020) , we can only report upper limits, but it is smaller than the calculations based on the tetraquark model. Meanwhile, it is not in contradiction with the molecule model calculation, which is about two orders of magnitude smaller than the current upper limit [11]. TABLE III. Comparison of the measured R Zc(3900) and R Zc(4020) with the theoretical predictions from Ref. [11]. "Type-I" and "Type-II" represent the tetraquark model predictions with and without taking into account the spin-spin interaction outside the diquarks, and "Molecule" is the molecule model prediction.