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

TL;DR

This study reports evidence for the decay of the exotic $Z_c^{\u00b1}$ states into $ ho^{b1}\u03b5_c$, measuring the cross section at specific energies and setting upper limits at others, aiding understanding of their nature.

Contribution

First evidence of $Z_c^{b1} o ho^{b1}\u03b5_c$ decay, with measured cross section and ratios to other decay modes, providing new insights into the structure of these exotic states.

Findings

Evidence for $Z_c^{b1} o ho^{b1}\u03b5_c$ at 4.226 GeV with 3.9$\sigma$ significance.

Measured cross section times branching fraction: $(48 \u00b1 11 \u00b1 11)$ pb.

Upper limits set at other energies; ratios of decay modes used to test theoretical models.

Abstract

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$.

Tables

Table 1. Table 1: The Born cross section ( σ B superscript 𝜎 B \sigma^{\rm B} ) for the e + ​ e − → π + ​ π − ​ π 0 ​ η c → superscript 𝑒 superscript 𝑒 superscript 𝜋 superscript 𝜋 superscript 𝜋 0 subscript 𝜂 𝑐 e^{+}e^{-}\to\pi^{+}\pi^{-}\pi^{0}\eta_{c} process and the numbers that enter the calculation (see Eq. ( 1 )). Here, s 𝑠 \sqrt{s} is in GeV, ℒ ℒ \mathcal{L} is in pb − 1 superscript pb 1 \rm pb^{-1} , ∑ ε ​ ℬ 𝜀 ℬ \sum\varepsilon{\cal B} is in % and σ B superscript 𝜎 B \sigma^{\rm B} is in pb.

$\sqrt{s}$	$ℒ$	$N_{\mathrm{sig}}$	$\left(1+\delta \right)$	$\frac{1}{{\left|1-\Pi \right|}^2}$	$\sum \epsilon ℬ$	${\sigma }^B$(${\sigma }_{U.L.}^B$)
4.226	1091.7	${324}_{-80}^{+83}$	0.74	1.056	0.82	${46}_{-11}^{+12}\pm 10$
4.258	825.7	${157}_{-68}^{+73}$	0.76	1.054	0.80	${30}_{-13}^{+14}\pm 9$ ( )
4.358	539.8	${32}_{-24}^{+62}$	1.03	1.051	0.62	$9_{-7}^{+17}\pm 2$ ( )
4.416	1073.6	${19}_{-18}^{+82}$	1.15	1.053	0.49	$3_{-3}^{+13}\pm 1$ ( )
4.600	566.9	$0_{-0}^{+28}$	1.32	1.055	0.31	$0_{-0}^{+12}\pm 13$ ( )

Table 2. Table 2: Born cross sections of e + ​ e − → π ∓ ​ Z c ​ ( 3900 ) ± → π ∓ ​ ρ ± ​ η c → superscript 𝑒 superscript 𝑒 superscript 𝜋 minus-or-plus subscript 𝑍 𝑐 superscript 3900 plus-or-minus → superscript 𝜋 minus-or-plus superscript 𝜌 plus-or-minus subscript 𝜂 𝑐 e^{+}e^{-}\to\pi^{\mp}Z_{c}(3900)^{\pm}\to\pi^{\mp}\rho^{\pm}\eta_{c} and e + ​ e − → π ∓ ​ Z c ​ ( 4020 ) ± → π ∓ ​ ρ ± ​ η c → superscript 𝑒 superscript 𝑒 superscript 𝜋 minus-or-plus subscript 𝑍 𝑐 superscript 4020 plus-or-minus → superscript 𝜋 minus-or-plus superscript 𝜌 plus-or-minus subscript 𝜂 𝑐 e^{+}e^{-}\to\pi^{\mp}Z_{c}(4020)^{\pm}\to\pi^{\mp}\rho^{\pm}\eta_{c} . 𝒮 𝒮 \mathcal{S} is the statistical significance of the signal. Other parameters are defined in the same way as those in Table  1 . Here, Z c ​ ( 3900 ) subscript 𝑍 𝑐 3900 Z_{c}(3900) is labeled as Z c subscript 𝑍 𝑐 Z_{c} and Z c ​ ( 4020 ) subscript 𝑍 𝑐 4020 Z_{c}(4020) is labeled as Z c ′ superscript subscript 𝑍 𝑐 ′ Z_{c}^{\prime} .

$\sqrt{s}$ (GeV)	$N_{\mathrm{obs}}^{Z_c}$	$N_{\mathrm{obs}}^{Z_c^{\prime }}$	$\left(1+\delta \right)$	$\frac{1}{{\left|1-\Pi \right|}^2}$	$\sum {\epsilon }^{Z_c}ℬ$ (%)	$\sum {\epsilon }^{Z_c^{\prime }}ℬ$ (%)	${\sigma }^{{\mathrm{BZ}}_c}$ (pb)	${\sigma }_{U.L.}^{{\mathrm{BZ}}_c}$ (pb)	${\sigma }_{U.L.}^{{\mathrm{BZ}}_c^{\prime }}$ (pb)	${𝒮}^{Z_c}$ ($\sigma$)	${𝒮}^{Z_c^{\prime }}$ ($\sigma$)
4.226	${240}_{-54}^{+56}$	${21}_{-11}^{+15}$	0.74	1.056	0.59	0.52	${48}_{-11}^{+11}\pm 11$	…		4.3	1.0
4.258	${92}_{-43}^{+48}$	$0_{-0}^{+11}$	0.76	1.054	0.50	0.56	${28}_{-13}^{+15}\pm 8$			2.0	…
4.358	${12}_{-8}^{+40}$	$0_{-0}^{+15}$	1.03	1.051	0.44	0.42	$5_{-3}^{+16}\pm 2$			0.3	…
4.416	${101}_{-44}^{+48}$	$6_{-4}^{+17}$	1.15	1.053	0.35	0.34	${22}_{-10}^{+10}\pm 5$			2.2	…
4.600	$0_{-0}^{+11}$	$0_{-0}^{+10}$	1.32	1.055	0.20	0.21	$0_{-0}^{+7}\pm 1$			…	…

Table 3. Table 3: Comparison of the measured R Z c ​ ( 3900 ) subscript 𝑅 subscript 𝑍 𝑐 3900 R_{Z_{c}(3900)} and R Z c ​ ( 4020 ) subscript 𝑅 subscript 𝑍 𝑐 4020 R_{Z_{c}(4020)} with the theoretical predictions.

Ratio	Measurement	Tetraquark	Molecule
$R_{Z_c\left(3900\right)}$	$2.3\pm 0.8$  Collaboration:2017njt	${230}_{-140}^{+330}$ Esposito:2014hsa	${0.046}_{-0.017}^{+0.025}$ Esposito:2014hsa
		${0.27}_{-0.17}^{+0.40}$ Esposito:2014hsa	$1.78\pm 0.41$ Goerke:2016hxf
		0.66 Faccini:2013lda	$6.84\times {10}^{-3}$ Patel:2014zja
		$0.56\pm 0.24$ Agaev:2016dev	$0.12$ Ke:2013gia
		$0.95\pm 0.40$ Dias:2013xfa
		$1.08\pm 0.88$ Wang:2017lot
		$1.28\pm 0.37$ Goerke:2016hxf
		$1.86\pm 0.41$ Goerke:2016hxf
$R_{Z_c\left(4020\right)}$	Ablikim:2013wzq	${6.6}_{-5.8}^{+56.8}$ Esposito:2014hsa	${0.010}_{-0.004}^{+0.006}$ Esposito:2014hsa

Equations

$${\sigma }^{\mathrm{B}}(e^{+}{e}^{-}\to {\pi }^{+}{\pi }^{-}{\pi }^0{\eta }_c)=\frac{N_{\mathrm{sig}}}{\mathcal{L}(1+\delta )\frac{1}{|1-\Pi {|}^2}{\sum }_i{\epsilon }_i{\mathcal{B}}_i},$$

Full text

The charged charmonium-like states $Z_{c}(3900)^{\pm}$ Ablikim:2013mio ; Liu:2013dau ; Ablikim:2013xfr and $Z_{c}(4020)^{\pm}$ Ablikim:2013wzq ; Ablikim:2013emm were first observed in 2013. Although their observed properties indicate they are not conventional mesons consisting of a quark-antiquark pair, their exact quark configurations are still unknown. Several models have been developed to describe their inner structure reviews , including loosely bound hadronic molecules of two charmed mesons Voloshin:1976ap , compact tetraquarks Maiani:2004vq ; Wang:2013vex , and hadro-quarkonium Voloshin:2007dx ; Dubynskiy:2008mq .

It has recently been suggested that the relative decay rate of $Z_{c}$ states, such as $Z_{c}(3900)\to\rho\eta_{c}$ to $\pi J/\psi$ (or $Z_{c}(4020)\to\rho\eta_{c}$ to $\pi h_{c}$), can be used to discriminate between the tetraquark and meson molecule scenarios Esposito:2014hsa . In Ref. Esposito:2014hsa , the predicted ratio $R_{Z_{c}(3900)}={\cal B}(Z_{c}(3900)\to\rho\eta_{c})/{\cal B}(Z_{c}(3900)\to\pi J/\psi)$ is $230_{-140}^{+330}$ or $0.27_{-0.17}^{+0.40}$ based on the diquark-antidiquark tetraquark model, depending on how the spin-spin interaction outside the diquarks is treated. On the other hand, using Non-Relativistic Effective Field Theory techniques, this ratio is only $0.046_{-0.017}^{+0.025}$ if we assume the $Z_{c}(3900)$ is a meson molecule state. Similarly, the predicted ratio of $R_{Z_{c}(4020)}={\cal B}(Z_{c}(4020)\to\rho\eta_{c})/{\cal B}(Z_{c}(4020)\to\pi h_{c})$ is $6.6_{-5.8}^{+56.8}$ in the tetraquark model, but only $0.010_{-0.004}^{+0.006}$ in the meson molecule model Esposito:2014hsa . However, the well-separated predictions for $R_{Z(3900)}$ and $R_{Z(4020)}$, shown above, could move closer or even overlap according to different theoretical approaches. Within QCD sum rule approaches Faccini:2013lda ; Agaev:2016dev ; Dias:2013xfa ; Wang:2017lot and convaraint quark model approaches Goerke:2016hxf to the tetraquark scenario, the predicted value of $R_{Z_{c}(3900)}$ can vary from 0.66 to 1.86. Furthermore, different approaches to the meson molecule model Patel:2014zja ; Ke:2013gia ; Goerke:2016hxf can lead to predictions for $R_{Z_{c}(3900)}$ from $6.8\times 10^{-3}$ to $1.8$. Consequently, the capability to separate the molecular and tetraquark models is currently model-dependent. In the hadron-charmonium model, the $Z_{c}(3900)$ is treated as a $J/\psi$ embedded in an S-wave spinless excitation of light-quark matter and consequently the transition $Z_{c}(3900)\to\rho\eta_{c}$ is expected to be suppressed compared to $Z_{c}(3900)\to\pi J/\psi$. A search for the decays of $Z_{c}(3900)$ or $Z_{c}(4020)$ to $\rho\eta_{c}$ thus offers an important opportunity to discrimate among the wide range of theoretical predictions.

In this Letter, we first report a search for the process $e^{+}e^{-}\to\pi^{+}\pi^{-}\pi^{0}\eta_{c}$. Then, based on the first step, we study the subprocesses $e^{+}e^{-}\to\pi Z_{c}(3900)^{\pm}$; $Z_{c}(3900)^{\pm}\to\rho^{\pm}\eta_{c}$ and $e^{+}e^{-}\to\pi Z_{c}(4020)^{\pm}$; $Z_{c}(4020)^{\pm}\to\rho^{\pm}\eta_{c}$. We use data samples collected with the BESIII detector bes3-detector at center-of-mass (c.m.) energies above 4 GeV, as listed in Table 1. The c.m. energies are measured using the $e^{+}e^{-}\to\mu^{+}\mu^{-}$ process with an uncertainty of $\pm 0.8$ MeV Ablikim:2015zaa . The beam spread is measured to be 1.6 MeV.

The design and performance of the BESIII detector are given in Ref bes3-detector . A geant4-based geant4 Monte Carlo (MC) simulation software package 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 evtgen . Initial state radiation (ISR) is simulated with kkmc kkmc , and final state radiation is handled with photos photos .

Charged tracks, photons and $K_{S}^{0}$ candidates are reconstructed using the standard criteria of the BESIII experiment Ablikim:2017pg . Candidate $\pi^{0}$ and $\eta$ decays to $\gamma\gamma$ are reconstructed from pairs of photons with invariant mass in the range [0.120, 0.145] GeV/$c^{2}$ for the $\pi^{0}$ and [0.50, 0.57] GeV/$c^{2}$ for the $\eta$. 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 $\pi^{0}$ or $\eta$ mass Tanabashi:2018oca .

The $\eta_{c}$ candidates are reconstructed using nine hadronic decays: $p\bar{p}$, $2(K^{+}K^{-})$, $K^{+}K^{-}\pi^{+}\pi^{-}$, $K^{+}K^{-}\pi^{0}$, $p\bar{p}\pi^{0}$, $K_{S}^{0}K^{\pm}\pi^{\mp}$, $\pi^{+}\pi^{-}\eta$, $K^{+}K^{-}\eta$, and $\pi^{+}\pi^{-}\pi^{0}\pi^{0}$. All combinations with invariant mass in the range [2.7, 3.2] GeV/$c^{2}$ are kept within each event. The signal region for the $\eta_{c}$ candidates is defined as [2.95, 3.02] GeV/$c^{2}$ and the sidebands as [2.78, 2.92] and [3.05, 3.19] GeV/$c^{2}$.

After the above selection, a four-constraint (4C) kinematic fit is performed for each event, and the $\chi^{2}$ of the fit ($\chi^{2}_{\rm 4C}$) is required to be less than 40 to suppress backgrounds. In each event, the mass of each track (excluding $K_{S}^{0}$ daughters) is taken to be that of the kaon, pion or proton, depending on the decay mode under study. Finally, only the combination of mass assignments with the minimum $\chi^{2}_{min}\equiv\chi^{2}_{\rm 4C}+\chi^{2}_{\rm 1C}+\chi^{2}_{\rm PID}+\chi^{2}_{\rm vertex}$ is kept. Here, $\chi^{2}_{\rm 1C}$ is the $\chi^{2}$ of the 1C fit for $\pi^{0}$ ($\eta$), $\chi^{2}_{\rm PID}$ is the sum of the $\chi^{2}$ for the PID of all charged tracks, and $\chi^{2}_{\rm vertex}$ is the $\chi^{2}$ of the $K_{S}^{0}$ 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^{-}\to\pi^{+}\pi^{-}\pi^{0}\eta_{c}$ are classified into two categories. They are events from (1) charmonium(-like) states decays (most of which include open-charm decays, e.g. $\psi\to D^{(*)}\bar{D}^{(*)}$); and (2) the continuum process, $e^{+}e^{-}\to q\bar{q}$, with $q=u$, $d$, and $s$.

By analyzing 600,000 $e^{+}e^{-}\to\pi^{+}\pi^{-}h_{c}$ MC simulation events with $h_{c}$ decaying inclusively, a small enhancement in the $\eta_{c}$ signal region is found. Using the measured cross section given in Ref Ablikim:2013wzq and the luminosity of data, its contribution, $N_{\rm bkg}^{\rm peaking}$, is estimated to be $8.7\pm 2.0$ at $\sqrt{s}=4.226$ GeV. The contributions at other energies are estimated in a similar way.

To suppress background events with charmed mesons, events are rejected if a $D$ meson candidate is reconstructed in one of its five decay modes: $D^{0}\to K^{\pm}\pi^{\mp}$, $D^{0}\to K^{\pm}\pi^{\mp}\pi^{0}$, $D^{\pm}\to K^{\pm}\pi^{\mp}\pi^{\pm}$, $D^{\pm}\to K_{S}^{0}\pi^{\pm}$, and $D^{\pm}\to K_{S}^{0}\pi^{\pm}\pi^{0}$. To accomplish this, we require the invariant mass of $D^{0}$($D^{\pm}$) candidates to be outside the region $m(D^{0})\pm 24$ MeV ($m(D^{\pm})\pm 10$ MeV). To reduce the continuum background, events with a $K^{*}(892)\to K\pi$, an $\omega\to\pi^{+}\pi^{-}\pi^{0}$, or an $\eta\to\pi^{+}\pi^{-}\pi^{0}$ candidate are removed by requiring $|M(K\pi)-m(K^{*})|>32$ MeV, $|M(\pi^{+}\pi^{-}\pi^{0})-m(\omega)|>26$ MeV, and $|M(\pi^{+}\pi^{-}\pi^{0})-m(\eta)|>10$ MeV, respectively. Here, $m(D^{0})$, $m(D^{\pm})$, $m(K^{*})$, $m(\omega)$ and $m(\eta)$ are the nominal masses of the corresponding states.

The mass windows for the background veto mentioned above and the $\chi^{2}$ requirement of the 4C kinematic fit are determined by optimizing the figure-of-merit (FOM), which is defined as ${\rm FOM}=S/\sqrt{S+B}$. Here, $S$ is the number of signal events from the MC simulation assuming $\sigma(e^{+}e^{-}\to\pi^{+}\pi^{-}\pi^{0}\eta_{c})=50~{}{\rm pb}$, which is evaluated from a measurement with unoptimized selection criteria. $B$ is the number of background events obtained from the $\eta_{c}$ sidebands in the data and extrapolated to the signal region linearly. The optimization is performed through iterations until all the selection criteria become stable.

To obtain the $\pi^{+}\pi^{-}\pi^{0}\eta_{c}$ yield, the invariant mass distributions of the $\eta_{c}$ candidates in the nine decay modes are fitted simultaneously using an unbinned maximum likelihood method. In the fit, the $\eta_{c}$ signal shape is determined from MC simulation and is described with a constant-width Breit-Wigner function (mass and width are fixed to the world average values Tanabashi:2018oca ) convolved with a Crystal Ball function, which represents instrumental resolution. The background is described with a second order Chebyshev Polynomial (CP). 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 Ablikim:2017pg . The total signal yield of the nine channels is labeled $N_{\rm obs}$, which is shared for all the channels and required to be positive. The free parameters in the fit include $N_{\rm obs}$ and the background yield and shape parameters for each decay mode. Figure 1(left) shows the fit results at $\sqrt{s}$ = 4.226 GeV projected onto the sum of events from all nine $\eta_{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\sigma$, which is obtained by comparing the change of the log-likelihood value $\Delta(-\ln L)=9.0$ with and without the $\pi^{+}\pi^{-}\pi^{0}\eta_{c}$ signal in the fit with one degree of freedom. The same selection criteria are applied to the other data sets, but no significant signals are observed.

The Born cross section of the $e^{+}e^{-}\to\pi^{+}\pi^{-}\pi^{0}\eta_{c}$ reaction is calculated using

[TABLE]

where $N_{\rm sig}=N_{\rm obs}-N_{\rm bkg}^{\rm peaking}$ is the number of signal events after the peaking background subtraction; $\mathcal{L}$ is the integrated luminosity; $(1+\delta)$ is the ISR correction factor, assuming the $\pi^{+}\pi^{-}\pi^{0}\eta_{c}$ signal is from $Y(4260)$ decays Tanabashi:2018oca ; and $\frac{1}{|1-\Pi|^{2}}$ is the vacuum-polarization factor Jegerlehner:1985gq . The cross sections and the numbers used for their calculation are listed in Table 1 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 $\sigma^{\rm B}$. The systematic uncertainties are incorporated into the upper limit by smearing the probability density function of the cross section Ablikim:2017pg . The corresponding results for $\sigma^{\rm B}_{\rm U.L.}$ are also listed in Table 1.

The $Z_{c}(3900)^{\pm}$ and $Z_{c}(4020)^{\pm}$ signals are examined after requiring that the invariant mass of an $\eta_{c}$ candidate is within the $\eta_{c}$ signal region [2.95, 3.02] GeV/$c^{2}$ and the invariant mass of $\pi^{\pm}\pi^{0}$ is within the $\rho$ signal region [0.675, 0.875] GeV/$c^{2}$. Here, we don’t distinguish the pions from $\eta_{c}$ decay or from collision and $\rho$ decay, therefore all possible combinations in one event are kept to avoid bias. To suppress the combinatorial background, the momenta of the pions from the $\rho$ decays are required to be less than 0.8 GeV/$c$. The events in the $\eta_{c}$ sidebands and $\rho$ sideband, which is defined as [0.475, 0.675] GeV/$c^{2}$, are investigated and no peaking structure is found. In addition, the simulated background events are studied (Fig. 2 left) and show good agreement with data both in the $\eta_{c}$ signal (Fig. 3 top) and sideband regions (Fig. 2 right). In the data sample, the $Z_{c}(3900)^{\pm}$ signal is apparent, but there is no statistically significant $Z_{c}(4020)^{\pm}$ signal.

To obtain the yields of $e^{+}e^{-}\to\pi^{\mp}Z_{c}(3900)^{\pm}\to\pi^{\mp}\rho^{\pm}\eta_{c}$ and $e^{+}e^{-}\to\pi^{\mp}Z_{c}(4020)^{\pm}\to\pi^{\mp}\rho^{\pm}\eta_{c}$, the invariant mass of $\rho^{\pm}\eta_{c}$ candidates in the nine $\eta_{c}$ decay channels are fitted simultaneously using the same method as for $e^{+}e^{-}\to\pi^{+}\pi^{-}\pi^{0}\eta_{c}$. In the fit, a possible interference between the signal and the background is neglected. The mass and width of the $Z_{c}(3900)^{\pm}$ are fixed to the values from the latest measurement Collaboration:2017njt and those of the $Z_{c}(4020)^{\pm}$ are fixed to world average values Tanabashi:2018oca . The mass resolution is obtained from MC simulation and parameterized as a Crystal Ball function Oreglia:1980cs . The background is described with a second order CP function. To validate the fit model, we perform a fit with the same model on the simulated background as shown in Fig. 2 (left). The signal yields of $Z_{c}(3900)^{\pm}$ and $Z_{c}(4020)^{\pm}$ are $48\pm 46$ and $0\pm 4$, respectively, and the statistical significance of the $Z_{c}(3900)^{\pm}$ is $0.6\sigma$. We also fit the sideband events both from data and MC with the second order CP function and the function can describe the sidebands well as shown in Fig. 2 (right). After the validation, we apply the fit model to data. Figure 3 shows the fit to the data set taken at $\sqrt{s}=4.226$ GeV. The total $Z_{c}(3900)^{\pm}$ signal yield is $240^{+56}_{-54}$ events with a statistical significance of $4.3\sigma$, and that of the $Z_{c}(4020)^{\pm}$ is $21^{+15}_{-11}$ events with a statistical significance of $1.0\sigma$. The signals at the other c.m. energies are not statistically significant.

The Born cross section for $e^{+}e^{-}\to\pi^{\mp}Z_{c}^{\pm}$ with $Z_{c}^{\pm}\to\rho^{\pm}\eta_{c}$ is calculated using the same equation as shown in Eq. (1). The numbers used in the calculation and the results are listed in Table 2.

The systematic uncertainties in the $\sigma^{\rm B}(e^{+}e^{-}\to\pi^{+}\pi^{-}\pi^{0}\eta_{c})$ measurement originate from the uncertainty of each factor in Eq. (1). The integrated luminosity has an uncertainty of $1.0\%$ Ablikim:2015nan . The uncertainty due to the subtraction of the $e^{+}e^{-}\to\pi^{+}\pi^{-}h_{c}$ peaking background events includes both the uncertainty due to the cross section and the statistical error of the MC sample. To estimate the uncertainty due to ISR correction, the c.m. energy dependent cross section of $e^{+}e^{-}\to\pi^{+}\pi^{-}J/\psi$ measured by the BESIII experiment Ablikim:2016qzw is used instead of Y(4260). The uncertainty from the signal shape consists of the mass resolution discrepancy between data and MC simulation and the uncertainty of the $\eta_{c}$ resonant parameters. The former is studied using an $e^{+}e^{-}\to\gamma_{ISR}J/\psi$ Ablikim:2017ove sample and the latter is estimated by varying the $\eta_{c}$ mass and width by $\pm 1\sigma$ around the world average values Tanabashi:2018oca . The uncertainty for the background shape is estimated by changing the order of the CP function and adjusting the fit boundaries. The methods for estimating the uncertainties due to the vacuum polarization and $\sum_{i}\varepsilon_{i}{\cal B}_{i}$ are the same as those described in Ref. Ablikim:2017pg . Furthermore, the uncertainty due to the $e^{+}e^{-}\to\pi^{+}\pi^{-}\pi^{0}\eta_{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 systematics uncertainty is that of $\sum_{i}\varepsilon_{i}{\cal B}_{i}$. The total systematic uncertainties are listed in Table 1.

For the $\sigma^{\rm B}(e^{+}e^{-}\to\pi^{\mp}Z_{c}(3900)^{\pm}(Z_{c}(4020)^{\pm})\to\pi^{\mp}\rho^{\pm}\eta_{c})$ measurement, the uncertainties on $\mathcal{L}$, ISR factors, $\sum_{i}\varepsilon_{i}{\cal B}_{i}$ and the vacuum polarization factor are studied following the methods described in the measurement of $\sigma^{\rm B}(e^{+}e^{-}\to\pi^{+}\pi^{-}\pi^{0}\eta_{c})$. Moreover, additional systematic uncertainties arise from the $\rho$ and $\eta_{c}$ selections, and the fit of the invariant mass spectrum of $\rho^{\pm}\eta_{c}$. The uncertainty due to the $M(\pi^{\pm}\pi^{0})$ mass window is estimated by comparing the invariant mass of $M(\omega\to\pi^{+}\pi^{-}\pi^{0})$ in data and MC assuming the mass resolution of $M(\pi^{+}\pi^{-}\pi^{0})$ is larger than $M(\pi^{\pm}\pi^{0})$. The discrepancy is found to be negligible. The uncertainty of the $\eta_{c}$ line shape is estimated by the varying the mass and width of the $\eta_{c}$ within the errors given by world average values Tanabashi:2018oca . The uncertainties affecting the fit to the $Z_{c}(3900)^{\pm}$ ($Z_{c}(4020)^{\pm}$) are estimated with the same methods as in the $\pi^{+}\pi^{-}\pi^{0}\eta_{c}$ case. All these sources and those in the $\sigma^{\rm B}(e^{+}e^{-}\to\pi^{+}\pi^{-}\pi^{0}\eta_{c})$ measurement are assumed to be independent and added in quadrature. The uncertainties related to the fit of invariant mass of $\eta_{c}\to$ hadrons are excluded because they don’t affect the $e^{+}e^{-}\to\pi Z_{c}$ measurement. The largest systematics uncertainty comes from $\sum_{i}\varepsilon_{i}{\cal B}_{i}$. The total systematic uncertainties are listed in Table 2.

To evaluate the effect of the systematic uncertainty on the signal significance at $\sqrt{s}=4.226$ GeV, we vary the signal shape, background parametrization, and fit range, or free the $Z_{c}$ mass, then repeat the fit. We find that the statistical significance of the $Z_{c}(3900)$ is always larger than $3.9~{}\sigma$.

In summary, using the $e^{+}e^{-}$ annihilation data at $\sqrt{s}=4.226$, $4.258$, $4.358$, $4.416$, and $4.600$ GeV, we study the $e^{+}e^{-}\to\pi^{+}\pi^{-}\pi^{0}\eta_{c}$ process for the first time. Evidence of this process is observed at $\sqrt{s}=4.226$ GeV with a significance of $4.2\sigma$ and the Born cross section $\sigma^{\rm B}(e^{+}e^{-}\to\pi^{+}\pi^{-}\pi^{0}\eta_{c})$ is measured to be $(46^{+12}_{-11}\pm 10)\rm\,pb$, excluding the processes $e^{+}e^{-}\to\omega\eta_{c}$ and $\eta\eta_{c}$. Evidence for the $\rho^{\pm}\eta_{c}$ decay mode of the charged charmonium-like state $Z_{c}(3900)^{\pm}$ is found in the process $e^{+}e^{-}\to\pi^{\mp}Z_{c}(3900)^{\pm}$ with $Z_{c}(3900)^{\pm}\to\rho^{\pm}\eta_{c}$ from the same data set. The measured cross section times branching ratio $\sigma^{\rm B}(e^{+}e^{-}\to\pi^{\mp}Z_{c}(3900)^{\pm})\times{\cal B}(Z_{c}(3900)^{\pm}\to\rho^{\pm}\eta_{c})$ is $(48\pm 11\pm 11)\rm{\,pb}$. This result indicates that the $e^{+}e^{-}\to\pi^{+}\pi^{-}\pi^{0}\eta_{c}$ process is dominated by the subprocess $e^{+}e^{-}\to\pi^{\mp}Z_{c}(3900)^{\pm}\to\pi^{\mp}\rho^{\pm}\eta_{c}$ (and implicitly $e^{+}e^{-}\to\pi^{0}Z_{c}(3900)^{0}\to\pi^{0}\rho^{0}\eta_{c}$). The significance of $Z_{c}(3900)^{\pm}\to\rho^{\pm}\eta_{c}$ is $3.9\sigma$ including the systematical uncertainty. No significant signal of $e^{+}e^{-}\to\pi^{+}\pi^{-}\pi^{0}\eta_{c}$ is observed at $\sqrt{s}=$4.258, 4.358, 4.416, and 4.600 GeV and no significant signal of $e^{+}e^{-}\to\pi^{\mp}Z_{c}(4020)^{\pm}$ with $Z_{c}(4020)^{\pm}\to\rho^{\pm}\eta_{c}$ is found in any of the data sets. Upper limits are deterimened at 90% C.L.

Using the results from Refs. Ablikim:2013wzq and Collaboration:2017njt , we calculate the ratios $R_{Z_{c}(3900)}={\cal B}(Z_{c}(3900)^{\pm}\to\rho^{\pm}\eta_{c})/{\cal B}(Z_{c}(3900)^{\pm}\to\pi^{\pm}J/\psi)$ and $R_{Z_{c}(4020)}={\cal B}(Z_{c}(4020)^{\pm}\to\rho^{\pm}\eta_{c})/{\cal B}(Z_{c}(4020)^{\pm}\to\pi^{\pm}h_{c})$. The results obtained from the measurements at $\sqrt{s}=$4.226, 4.258, and 4.358 GeV are listed in Table 3, together with the theoretical predictions for comparison.

The measured $R_{Z_{c}(3900)}$ is closer to the calculation of the tetraquark model than to that of the meson molecule model in Ref. Esposito:2014hsa . The measurement is also consistent with several other independent calculations based on the tetraquark scenario Faccini:2013lda ; Agaev:2016dev ; Dias:2013xfa ; Wang:2017lot ; Goerke:2016hxf . For the molecule model, as we mentioned before, the calculated $R_{Z_{c}(3900)}$ is highly model dependent  Patel:2014zja ; Ke:2013gia ; Goerke:2016hxf . Therefore, it is necessary to narrow down the theoretical uncertainty in the molecular framework to have a better comparison with the measurement. In the hadron-charmonium model, the ${\cal B}(Z_{c}(3900)\to\rho\eta_{c})$ is suppressed compared with ${\cal B}(Z_{c}(3900)\to\pi J/\psi)$ and therefore inconsistent with the measurement Voloshin:2013dpa . Furthermore, this model predicts a new resonance $W_{c}(3785)$, which can be produced via $e^{+}e^{-}\to\rho W_{c}\to\rho\pi\eta_{c}$, the same final state we analyzed here. As we found that the $e^{+}e^{-}\to\pi^{+}\pi^{-}\pi^{0}\eta_{c}$ process is saturated by $e^{+}e^{-}\to\pi Z_{c}(3900)\to\rho\pi\eta_{c}$, we can conclude that the production of the $W_{c}$, if present, is small compared to $e^{+}e^{-}\to\pi Z_{c}(3900)$.

For $R_{Z_{c}(4020)}$, we can only report upper limits, but they are smaller than the value calculated based on the tetraquark model. On the other hand, the upper limits are not in contradiction with the molecule model calculation, which is about two orders of magnitude smaller than the current upper limits Esposito:2014hsa .