Measurement of $e^{+}e^{-}\rightarrow \pi^{+}\pi^{-}\psi(3686)$ from 4.008 to 4.600$\sim$GeV and observation of a charged structure in the $\pi^{\pm}\psi(3686)$ mass spectrum

TL;DR

This paper measures the cross section of the process $e^{+}e^{-} ightarrow \pi^{+}\pi^{-}\psi(3686)$ across energies from 4.008 to 4.600 GeV, discovering two resonances and a charged structure, with improved precision over previous results.

Contribution

First observation of a lower-mass resonance in $e^{+}e^{-} ightarrow \\pi^{+}\pi^{-}\\psi(3686)$ with high significance, and identification of a charged charmonium-like structure in the mass spectrum.

Findings

Two resonances at 4209.5 MeV and 4383.8 MeV with widths around 80 MeV.

First observation of the lower-mass resonance with 5.8 sigma significance.

Charged structure observed in the \\pi^{\\pm}\\psi(3686) spectrum at 4.416 GeV.

Abstract

We study the process $e^{+}e^{-} \rightarrow \pi^{+}\pi^{-}\psi(3686)$ using 5.1 $fb^{-1}$ of data collected at 16 center-of-mass energy ($\sqrt{s}$) points from 4.008 to 4.600$\sim$GeV by the BESIII detector operating at the BEPCII collider. The measured Born cross sections for $e^{+}e^{-}\rightarrow \pi^{+}\pi^{-}\psi(3686)$ are consistent with previous results, but with much improved precision. A fit to the cross section shows contributions from two structures: the first has $M=4209.5\pm7.4\pm1.4$$\sim$MeV/$c^{2}$ and $\Gamma=80.1\pm24.6\pm2.9$$\sim$MeV, and the second has $M=4383.8\pm4.2\pm0.8$$\sim$MeV/c$^{2}$ and $\Gamma=84.2\pm12.5\pm2.1$$\sim$MeV, where the first errors are statistical and the second systematic. The lower-mass resonance is observed in the process $e^{+}e^{-}\rightarrow \pi^{+}\pi^{-}\psi(3686)$ for the first time with a statistical significance of $5.8\sigma$. A charged charmonium-like structure is observed in the $\pi^{\pm}\psi(3686)$ invariant mass spectrum for data at $\sqrt{s} = 4.416~$GeV. A fit with an $S$-wave Breit-Wigner function yields a mass $M=4032.1\pm2.4$~$MeV$^{2}$, where the errors are statistical only. However, there are still unresolved discrepancies between the fit model and data. The width of the intermediate state varies in a wide range for different kinematic regions within the data set. Therefore no simple interpretation of the data has been found, and a future data sample with larger statistics and more theoretical input will be required to better understand this issue.

Tables

Table 1. Table 1: Summary of the measurement of the Born cross section σ B superscript 𝜎 𝐵 \sigma^{B} at individual c. m. energies. The subscript 1 1 1 or 2 2 2 denotes mode I or II, respectively. The first uncertainties are statistical, and the second systematic. An upper limit at the 90% confidence level (C.L.) is determined by a profile likelihood method  uplimit for data samples with low signal significance.

$\sqrt{s}$ (GeV)	$L$ (pb-1)	$N_1^{\text{obs}}$	${ϵ}_1$ (%)	$N_2^{\text{obs}}$	${ϵ}_2(\%)$	$\left(1+{\delta }^r\right)$	$\left(1+{\delta }^v\right)$	${\sigma }_1^B$ (pb)	${\sigma }_2^B$ (pb)	${\sigma }^B$ (pb)
4.008	482	$0.0\pm 0.6$	22.6	$0.2\pm 2.3$	4.7	0.70	1.056
4.085	52.6	$4.0\pm 2.0$	36.1	$1.0\pm 1.0$	20.9	0.75	1.056	$6.5\pm 3.2\pm 0.9$	$3.9\pm 3.9\pm 0.3$	$5.4\pm 2.5\pm 0.6$
4.189	43.1	$3.8\pm 2.0$	39.2	$2.8\pm 2.2$	27.7	0.76	1.056	$6.8\pm 3.6\pm 0.7$	$9.9\pm 7.8\pm 1.5$	$7.3\pm 3.3\pm 0.7$
4.208	54.6	$8.9\pm 3.3$	40.7	$7.0\pm 3.0$	27.5	0.76	1.057	$12.2\pm 4.5\pm 1.7$	$20.0\pm 8.4\pm 1.9$	$14.0\pm 4.0\pm 1.5$
4.217	54.1	$13.0\pm 3.6$	40.9	$0.0\pm 0.7$	27.3	0.76	1.057	$17.8\pm 4.9\pm 1.5$		$17.8\pm 4.9\pm 1.5$
4.226	1092	$315\pm 18$	39.2	$141\pm 14$	28.1	0.76	1.056	$22.3\pm 1.3\pm 1.7$	$19.4\pm 1.9\pm 2.0$	$21.3\pm 1.1\pm 1.6$
4.242	55.6	$11.0\pm 3.3$	41.4	$7.9\pm 3.1$	28.0	0.76	1.053	$14.6\pm 4.4\pm 1.3$	$21.5\pm 8.4\pm 1.8$	$16.0\pm 3.9\pm 1.2$
4.258	826	$241\pm 16$	40.3	$84\pm 11$	23.5	0.76	1.054	$22.0\pm 1.4\pm 1.7$	$18.3\pm 2.5\pm 1.8$	$20.9\pm 1.2\pm 1.5$
4.308	44.9	$17.0\pm 4.2$	41.6	$15.0\pm 4.1$	27.3	0.74	1.053	$28.2\pm 6.9\pm 2.6$	$53.2\pm 14.5\pm 7.4$	$32.1\pm 6.2\pm 2.8$
4.358	540	$439\pm 21$	41.2	$275\pm 19$	29.8	0.79	1.051	$57.8\pm 2.8\pm 4.4$	$69.8\pm 4.8\pm 5.2$	$61.0\pm 2.4\pm 4.3$
4.387	55.2	$56.6\pm 7.6$	39.4	$25.7\pm 6.1$	29.4	0.86	1.051	$70.1\pm 9.4\pm 6.7$	$59.4\pm 14.1\pm 5.6$	$66.4\pm 7.8\pm 5.5$
4.416	1074	$693\pm 27$	37.8	$415\pm 24$	27.4	0.96	1.053	$41.0\pm 1.6\pm 3.2$	$47.3\pm 2.7\pm 3.7$	$42.8\pm 1.4\pm 3.0$
4.467	110	$15.1\pm 4.2$	32.6	$8.3\pm 4.2$	23.9	1.10	1.055	$8.8\pm 2.5\pm 1.0$	$9.2\pm 4.7\pm 1.9$	$8.9\pm 2.2\pm 0.9$
4.527	110	$13.4\pm 4.0$	29.1	$7.0\pm 3.6$	20.8	1.25	1.055	$7.7\pm 2.3\pm 0.9$	$7.8\pm 4.0\pm 1.2$	$7.7\pm 2.0\pm 0.8$
4.575	47.7	$4.5\pm 2.3$	28.3	$5.7\pm 3.2$	20.2	1.23	1.055	$6.2\pm 3.2\pm 0.8$	$15.4\pm 8.7\pm 1.8$	$7.3\pm 3.0\pm 0.8$
4.600	567	$106\pm 11$	31.8	$71\pm 10$	21.5	1.08	1.055	$12.6\pm 1.3\pm 1.2$	$17.2\pm 2.4\pm 1.6$	$14.6\pm 1.1\pm 1.1$

Table 2. Table 2: Results of the fit to the e + ​ e − → π + ​ π − ​ ψ ​ ( 3686 ) → superscript 𝑒 superscript 𝑒 superscript 𝜋 superscript 𝜋 𝜓 3686 e^{+}e^{-}\to\pi^{+}\pi^{-}\psi(3686) cross section. The error is statistical only.

Parameters	Solution I	Solution II
$M\left(Y4220\right)$   (MeV/$c^2$)	$4209.5\pm 7.4$
$\Gamma \left(Y\left(4220\right)\right)$   (MeV)	$80.1\pm 24.6$
$ℬ{\Gamma }^{e^{+}{e}^{-}}\left(Y\left(4220\right)\right)$   (eV)	$1.61\pm 1.27$	$1.80\pm 1.41$
$M\left(Y4390\right)$   (MeV/$c^2$)	$4383.8\pm 4.2$
$\Gamma \left(Y\left(4390\right)\right)$   (MeV)	$84.2\pm 12.5$
$ℬ{\Gamma }^{e^{+}{e}^{-}}\left(Y\left(4390\right)\right)$   (eV)	$7.25\pm 2.8$	$10.96\pm 3.8$
${\varphi }_1$   (rad)	$3.3\pm 1.0$	$2.8\pm 0.4$
${\varphi }_2$   (rad)	$0.8\pm 0.9$	$4.7\pm 0.1$

Table 3. Table 3: Systematic uncertainty in the measurement of resonance parameters.

Sources	$Y\left(4220\right)$		$Y\left(4390\right)$
	$M$ (MeV/$c^2)$	$\Gamma$ (MeV)	$M$ (MeV/$c^2)$	$\Gamma$ (MeV)
Energy scale	0.1	0.0	0.1	0.0
Energy spread	1.3	0.1	0.7	1.6
Fit range	0.2	1.8	0.2	0.5
$Y\left(4660\right)$	0.6	2.3	0.3	1.2
Total	1.4	2.9	0.8	2.1

Table 4. Table 4: Masses and widths of the vector charmonium states observed from different processes at BESIII in the mass region between 4.2 and 4.4 GeV/ c 2 superscript 𝑐 2 c^{2} . The subscript 1 or 2 denotes the lower mass state or higher mass state.

Process	$M_1$ (MeV/$c^2$)	${\Gamma }_1$ (MeV)	$M_2$ (MeV/$c^2$)	${\Gamma }_2$ (MeV)
$e^{+}{e}^{-}\to \omega {\chi }_{c0}$	$4230\pm 8\pm 6$	$38\pm 12\pm 2$ omegachicj
$e^{+}{e}^{-}\to {\pi }^{+}{\pi }^{-}J/\psi$	$4220.0\pm 3.1\pm 1.4$	$44.1\pm 4.3\pm 2.0$	$4320.0\pm 10.4\pm 7.0$	${101.4}_{-19.7}^{+25.3}\pm 10.2$pipijpsi
$e^{+}{e}^{-}\to {\pi }^{+}{\pi }^{-}{h}_c$	${4218.4}_{-4.5}^{+5.5}\pm 0.9$	${66.0}_{-8.3}^{+12.3}\pm 0.4$	${4391.5}_{-6.8}^{+6.3}\pm 1.0$	${139.5}_{-20.6}^{+16.2}\pm 0.6$ pipihc
$e^{+}{e}^{-}\to {\pi }^{+}{D}^0{D}^{\ast -}+c.c$	$4224.8\pm 5.6\pm 4.0$	$72.3\pm 9.1\pm 0.9$	$4400.1\pm 9.3\pm 2.1$	$181.7\pm 16.9\pm 7.4$ piDD
$e^{+}{e}^{-}\to {\pi }^{+}{\pi }^{-}\psi \left(3686\right)$	$4209.5\pm 7.4\pm 1.4$	$80.1\pm 24.6\pm 2.9$	$4383.8\pm 4.2\pm 0.8$	$84.2\pm 12.5\pm 2.1$

Table 5. Table 5: Parameters of the resolution for M 2 ​ ( π ± ​ ψ ​ ( 3686 ) ) superscript 𝑀 2 superscript 𝜋 plus-or-minus 𝜓 3686 M^{2}(\pi^{\pm}\psi(3686)) , M 𝑀 M and σ 𝜎 \sigma are the mean and resolution of Gaussian function, respectively. The subscript 1 1 1 or 2 2 2 denote Gaussian function one or two. f ​ r ​ a 𝑓 𝑟 𝑎 fra is the fraction of Gaussian function one.

$\sqrt{s}$ (GeV)	$M_1$ (GeV/c2)2	${\sigma }_1$ (GeV/c2)2	$M_2$ (GeV/c2)2	${\sigma }_2$ (GeV/c2)2	$fra$
4.416 (Mode I)	$0.00088\pm 0.00017$	$0.01611\pm 0.00016$	$0.0377\pm 0.0014$	$0.0543\pm 0.0010$	$0.7862\pm 0.0064$
4.416 (Mode II)	$0.00014\pm 0.00047$	$0.01918\pm 0.00031$	$0.0581\pm 0.0046$	$0.0311\pm 0.0023$	$0.841\pm 0.019$
4.358 (Mode I)	$0.00114\pm 0.00015$	$0.01585\pm 0.00014$	$0.0376\pm 0.0013$	$0.0567\pm 0.0010$	$0.7809\pm 0.0057$
4.358 (Mode II)	$0.00011\pm 0.00037$	$0.01839\pm 0.00029$	$0.0508\pm 0.0049$	$0.0335\pm 0.0021$	$0.839\pm 0.020$
4.258 (Mode I)	$0.00037\pm 0.00012$	$0.01562\pm 0.00011$	$0.0326\pm 0.0010$	$0.05346\pm 0.00079$	$0.7974\pm 0.0048$
4.258 (Mode II)	$-0.001167\pm 0.00026$	$0.01724\pm 0.00023$	$0.0432\pm 0.0032$	$0.0344\pm 0.0013$	$0.841\pm 0.014$

Equations

$${\sigma }_i^{\mathrm{B}}=\frac{N_i^{\mathrm{obs}}}{{\mathcal{L}}_{\mathrm{int}}(1+{\delta }^r)(1+{\delta }^v)\mathcal{B}{r}_i{ϵ}_i}(i=1,2),$$

$$\mathcal{A}=f_1{e}^{i{\varphi }_1}+f_2+f_3{e}^{i{\varphi }_2},$$

$$f_i=\frac{M_i}{M}\frac{\sqrt{12\pi {\mathcal{B}}_i{\Gamma }_i^{e^{+}{e}^{-}}{\Gamma }_i}}{M^2-M_i^2+i{M}_i{\Gamma }_i}\sqrt{\frac{\Phi (M)}{\Phi (M_i)}},$$

$$\frac{p\cdot q/c^2}{(M_R^2-x{)}^2+M_R^2\cdot {\Gamma }^2/c^4}+\frac{p\cdot q/c^2}{(M_R^2-y{)}^2+M_R^2\cdot {\Gamma }^2/c^4},$$

$$ϵ(x,y)\cdot (\frac{p\cdot q}{(M_R^2-x{)}^2+M_R^2\cdot {\Gamma }^2}+\frac{p\cdot q}{(M_R^2-y{)}^2+M_R^2\cdot {\Gamma }^2})\otimes \sigma (x,y),$$

Full text

I INTRODUCTION

In the past decade, a series of new vector ($J^{PC}=1^{--}$) charmoniumlike states, e.g., the $Y(4260)$, $Y(4360)$, and $Y(4660)$, have been observed in $e^{+}e^{-}$ annihilation through dipion hadronic transitions to low mass charmonium states Y4260BaBar ; Y4260Cleo ; Y4260Belle ; Y4360Belle ; Y4360BaBar . Many theoretical interpretations have been proposed to understand the underlying structure of the $Y$ family of states, such as hybrid charmonium Y4360hybrid , tetraquark Y4360tetra , or hadronic molecule Y4360molecule etc. The $Y(4360)$ was first observed in $e^{+}e^{-}\to\gamma_{\rm ISR}Y(4360)\to\gamma_{\rm ISR}\pi^{+}\pi^{-}\psi(3686)$ by BaBar Y4360BaBar and subsequently confirmed by Belle Y4360Belle . Recently, two resonant structures were observed in the processes $e^{+}e^{-}\to\pi^{+}\pi^{-}J/\psi$ at BESIII pipijpsi , indicating the $Y(4260)$ resonance reported by previous experiments Y4260BaBar ; Y4260Cleo ; Y4260Belle actually consists of two structures. There are two resonances observed in the process $e^{+}e^{-}\to\pi^{+}\pi^{-}h_{c}$ pipihc in the mass region between 4.2 and 4.4 GeV/$c^{2}$. The lower mass state, $Y(4220)$, has a mass consistent with the lower mass state observed in $e^{+}e^{-}\to\pi^{+}\pi^{-}J/\psi$ pipijpsi , but has a somewhat larger width. The higher mass state, $Y(4390)$, does not match any known vector charmonium or charmonium-like states. Therefore, a more precise measurement of the cross section of $e^{+}e^{-}\to\pi^{+}\pi^{-}\psi(3686)$ at BESIII will help to clarify the spectrum of vector particles.

In recent years, another mysterious and interesting new pattern of charmoniumlike states, the $Z_{c}^{\pm}$’s, have been observed both in final states containing a charged pion and a low mass charmonium state Y4260Belle ; ZC3900BES ; ZC3900CLEO ; ZC4020BES , and in pairs of charmed mesons $(D^{(\ast)}\bar{D}^{\ast})^{\pm}$ ZC3885BES ; ZC3885BES2 ; ZC4025BES . A similar pattern of states is also observed in the bottomonium system zb . All of these states contain a minimum of four quarks as they are electrically charged and their constituents include a $c\bar{c}$ pair. They are good candidates for meson-meson molecules or tightly bound tetraquark states. To clarify their true nature, it is essential to search for this kind of exotic state in other final states. Belle’s recently updated result of $e^{+}e^{-}\to\gamma_{\rm ISR}\pi^{+}\pi^{-}\psi(3686)$, using their full data samples, shows evidence for a new $Z_{c}^{\pm}$ state in the mass spectrum of $\pi^{\pm}\psi(3686)$ with a mass around $4.05$ GeV/$c^{2}$ and a statistical significance of $3.5\sigma$ Y4360Belle . It is important to confirm the existence of this new $Z_{c}^{\pm}$ candidate in other experiments.

In this paper, we present a study of $e^{+}e^{-}\to\pi^{+}\pi^{-}\psi(3686)$, with two decay modes $\psi(3686)\to\pi^{+}\pi^{-}J/\psi$ (mode I) and $\psi(3686)\to\text{neutrals}+J/\psi$ (mode II), where ‘neutrals’ refers to $\pi^{0}\pi^{0}$, $\pi^{0}$, $\eta$ and $\gamma\gamma$, at center-of-mass (c.m.) energies $\sqrt{s}$ from 4.008 to 4.600 GeV. The $J/\psi$ is reconstructed in its prominent decay mode $J/\psi\to\ell^{+}\ell^{-}$ ($\ell=e/\mu$). The actual c.m. energies are measured by studying the process $e^{+}e^{-}\to\mu^{+}\mu^{-}$ cme . The data samples used in this analysis were collected with the BESIII detector at the BEPCII collider at 16 different c.m. energy points with a total integrated luminosity of 5.1 fb*-1* luminosity .

II THE BESIII EXPERIMENT AND THE DATA SETS

BEPCII is a double-ring $e^{+}e^{-}$ collider running at c.m. energies between 2.0 and 4.6 GeV, reaching a peak luminosity of $1.0\times 10^{33}$cm*-2s-1* at a c.m. energy of 3770 MeV. The cylindrical BESIII detector has an effective geometrical acceptance of $93\%$ of 4$\pi$ and is divided into a barrel section and two endcaps. It contains a small cell, helium-based (60$\%$ He, 40$\%$ C3H8) main drift chamber (MDC) which provides a momentum measurement of charged particles with a resolution of $0.5\%$ at a momentum of 1 GeV/c in a magnetic field of 1 Tesla. The energy loss measurement ($dE/dx$) provided by the MDC has a resolution better than $6\%$. A time-of-flight system (TOF) consisting of 5-cm-thick plastic scintillators can measure the flight time of charged particles with a time resolution of 80 ps in the barrel and 110 ps in the endcaps. An electromagnetic calorimeter (EMC) consisting of 6240 CsI (Tl) in a cylindrical structure and two endcaps is used to measure the energies of photons and electrons. The energy resolution of the EMC is $2.5\%$ in the barrel and $5.0\%$ in the endcaps for photon/electrons with an energy of 1 GeV. The position resolution of the EMC is 6 mm in the barrel and 9 mm in the endcaps. A muon system (MUC) consisting of 1000 m2 of Resistive Plate Chambers (RPC) is used to identify muons and provides a spatial resolution better than 2 cm. A detailed description of the BESIII detector can be found in Ref. besint .

The Geant4-based geant4 Monte Carlo (MC) simulation software package Boost Deng is used to generate the signal and background MC samples. The $e^{+}e^{-}$ collision is simulated with the KKMC kkmc generator taking into account the beam energy spread and Initial State Radiation (ISR), where the line shape of the $e^{+}e^{-}\to\pi^{+}\pi^{-}\psi(3686)$ cross section is taken from the latest Belle results Y4360Belle . The QED Final State Radiative (FSR) correction for $J/\psi\to\ell^{+}\ell^{-}$ is incorporated with Photos photos . Large signal MC samples $e^{+}e^{-}\to\pi^{+}\pi^{-}\psi(3686)$ as well as the dominant backgrounds are generated exclusively at each c.m. energy, where the decay $\psi(3686)\to\pi^{+(0)}\pi^{-(0)}J/\psi$ is simulated with the Jpipi model jpipi , and the decay $e^{+}e^{-}\to\pi^{+}\pi^{-}\psi(3686)$ is generated with both the Jpipi model and a model distributed uniformly in phase space (PHSP). The Jpipi model is constructed with an effective Lagrangian for the hadronic decays and the coupling constants in the model are obtained by a fit to the invariant $\pi^{+}\pi^{-}$ mass distributions in $\psi(3686)\to\pi^{+}\pi^{-}J/\psi$. The generic inclusive MC samples at $\sqrt{s}=4.258$ and 4.358 GeV, with integrated luminosity equivalent to that of data, are generated to study the potential background contributions, where the known decay modes are generated with Evtgen evtgen with branching fractions set to values in the Particle Data Group (PDG) pdg , and the remaining unknown ones are generated with Lundcharm lundcharm .

III Event selection

The final states of the decay $e^{+}e^{-}\to\pi^{+}\pi^{-}\psi(3686)$ have topologies $\pi^{+}\pi^{-}\pi^{+}\pi^{-}\ell^{+}\ell^{-}$ and \pi^{+}\pi^{-}$$\ell^{+}\ell^{-}+N\gamma ($N\geq 2)$ for mode I and II, respectively. A good charged track must have a polar angle $|\cos\theta|<0.93$ in the MDC, and have the point of closest approach to the interaction point within 1 cm in the plane perpendicular to the beam and within $\pm 10$ cm along the beam. Candidate events are required to have five (with $\pm 1$ net charge) or six (with zero net charge) charged tracks for mode I and four (with zero net charge) charged tracks for mode II. Photon candidates are reconstructed by clustering EMC energy deposits. They must have a minimum energy of 25 MeV for $|\cos\theta|<$ 0.8 or 50 MeV for 0.86 $<|\cos\theta|<$ 0.92, and be isolated from all charged tracks by an angle $>10^{\circ}$. EMC timing is used to suppress electronic noise and energy deposition unrelated to the event. At least two photon candidates are required for mode II. Since pions and leptons are kinematically well separated for the signal process, charged tracks with momentum $p>1.0$ GeV/$c$ are assumed to be leptons, while those with p<0.65\;\mbox{GeV/c} are pions. Electron and muon separation is performed with the deposited energy $E$ in the electromagnetic calorimeter. Electrons must have $E/pc>0.7$, while muons have $E<0.45\;\mbox{GeV}$. Candidate events are required to have one lepton pair $\ell^{+}\ell^{-}$, and the other tracks must be pions. The number of hits in the muon counter is further required to be not less than 5 for $J/\psi\to\mu^{+}\mu^{-}$ candidates.

In mode I, a four-constraint (4C) kinematic fit imposing energy-momentum conservation under the $e^{+}e^{-}\to\pi^{+}\pi^{-}\pi^{+}\pi^{-}\ell^{+}\ell^{-}$ hypothesis is applied for the event candidates with six charged tracks, and a one-constraint (1C) kinematic fit with one missing pion is applied for the candidates with five charged tracks. The corresponding $\chi^{2}$ of the fit must be $\chi^{2}_{4C}<60$ or $\chi^{2}_{1C}<15$, respectively. The $J/\psi$ signal appears with low background in the spectrum of the $\ell^{+}\ell^{-}$ invariant mass, $M(\ell^{+}\ell^{-})$, calculated with the corrected momentum after the kinematic fit, and is required to be within 3.05<M(\ell^{+}\ell^{-})<3.15\;\mbox{GeV/c^{2}}. To improve the mass resolution of the $\psi(3686)$, the selected candidates with six or five good charged tracks are further fed into a five-constraint (5C) or two-constraint (2C) kinematic fit with an additional constraint on $M(\ell^{+}\ell^{-})$ to the nominal $J/\psi$ mass ($M(J/\psi)$) pdg , respectively. The $\psi(3686)$ signal is reconstructed with the $\pi^{+}\pi^{-}J/\psi$ system (four combinations per event) whose invariant mass $M(\pi^{+}\pi^{-}J/\psi)$, calculated with the corrected momentum after the 5C or 2C fit, is closest to the nominal $\psi(3686)$ mass ($M(\psi(3686))$) pdg .

In mode II, no kinematic fit is applied. The $J/\psi$ signal is extracted by requiring 3.05<M(\ell^{+}\ell^{-})<3.15\;\mbox{GeV/c^{2}}. The $\psi(3686)$ signal is reconstructed using the mass recoiling against the $\pi^{+}\pi^{-}$ system, $M^{\rm recoil}(\pi^{+}\pi^{-})$. The requirement |M^{\rm corr}_{\psi(3686)}-M(\psi(3686))|>8\;\mbox{MeV/c^{2}} is imposed to veto the background $e^{+}e^{-}\to\text{neutrals}+\psi(3686)$ ($\psi(3686)\to\pi^{+}\pi^{-}J/\psi$), where $M^{\rm corr}_{\psi(3686)}$ = $M(\pi^{+}\pi^{-}\ell^{+}\ell^{-})-M(\ell^{+}\ell^{-})+M(J/\psi)$. A requirement |M(\gamma\gamma\pi^{+}\pi^{-})-M(\eta)|>50\;\mbox{MeV/c^{2}} is used to eliminate the background $e^{+}e^{-}\to\eta J/\psi$ ($\eta\to\pi^{+}\pi^{-}\pi^{0}$), where $\gamma\gamma$ are the two photons with largest energy, and $M(\eta)$ is the nominal $\eta$ mass pdg . A requirement $\cos\theta_{\pi^{+}\pi^{-}}<0.9$ is also applied to remove the radiative Bhabha and dimuon background in which a photon converts into a $e^{+}e^{-}$ pair and is misidentified as a $\pi^{+}\pi^{-}$ pair.

IV Extraction of the Born cross section

With the above selection criteria, a prominent $\psi(3686)$ signal over a small background is observed in the $M(\pi^{+}\pi^{-}J/\psi)$ and $M^{\rm recoil}(\pi^{+}\pi^{-})$ spectra for Mode I and Mode II, respectively, as shown in Fig. 1 for data at $\sqrt{s}=4.416\;\mbox{GeV}$ as an example. Potential backgrounds without a $J/\psi$ in the final state are explored with events within the $J/\psi$ sideband regions. Backgrounds including a $J/\psi$ in the final state are studied with both inclusive and exclusive MC samples. From these studies, no peaking background is expected in the range of interest in the two $\psi(3686)$ decay modes. To determine the signal yields, unbinned maximum likelihood fits are performed to the $M(\pi^{+}\pi^{-}J/\psi)$ and $M^{\rm recoil}(\pi^{+}\pi^{-})$ spectra for the two modes, respectively. Here, the signal probability density function (PDF) is described with the MC simulated shape convoluted with a Gaussian function to account for the mass resolution difference between data and MC simulation, and the PDF for background is presented as a linear function. Simultaneous fits to data at 16 different c.m. energies are performed. The c.m. energies share a common background shape and a common data/MC resolution difference.

The Born cross section is determined from

[TABLE]

where $N^{\rm obs}_{i}$ is the signal yield extracted from the fit, $\mathcal{L}_{\rm int}$ is the integrated luminosity, $1+\delta^{r}$ is the ISR correction factor, $1+\delta^{v}$ is the vacuum polarization factor, $\mathcal{B}r$ is the product of the decay branching fractions of the intermediate states in the cascade decay (4.11% for mode I and 2.95% for mode II), and $\epsilon$ is the detection efficiency. The various numbers used in the Born cross section calculation are summarized in Table 1.

Several sources of systematic uncertainty are considered in the measurement of the Born cross section. The uncertainty in the tracking efficiency for leptons and pions are 1.0% per track trkuncertainty . The uncertainty in the photon efficiency is 1.0% per photon phouncertainty . The uncertainty related to the requirement on the number of hits in the muon counter is 4.2%, determined by studying a control sample of $e^{+}e^{-}\to\pi^{+}\pi^{-}J/\psi$. The uncertainty associated with the kinematic fit is estimated with the same technique as in Ref. kinematic . The uncertainty associated with the $J/\psi$ mass requirement is estimated by smearing the $M(\ell^{+}\ell^{-})$ distribution of MC samples according to the resolution of $M(\ell^{+}\ell^{-})$ in data, and the difference in MC efficiencies is taken as the uncertainty. In mode II, the uncertainties associated with the requirements on $M^{corr}_{\psi(3686)}$, $M(\gamma\gamma\pi^{+}\pi^{-})$ and $\cos\theta_{\pi^{+}\pi^{-}}$ are evaluated by changing the requirements. The uncertainty associated with the fit procedure is investigated by varying the fit range, replacing the linear function by a second-order polynomial function for the background, and varying the width of the Gaussian function for signal. For data sets with large luminosity, the detection efficiency is estimated with a MC sample weighted according to the observed Dalitz plot $M^{2}(\pi^{+}\psi(3686))$ versus $M^{2}(\pi^{+}\pi^{-})$. The corresponding uncertainty is estimated by varying the weight factors according to the statistical uncertainty of data in each bin on the Dalitz plane. For data sets with low luminosity, the detection efficiency is estimated using MC samples with the Jpipi model; the corresponding uncertainty is evaluated with large luminosity data sets, by taking the largest difference of efficiencies between the Jpipi model MC sample and their nominal values with weighted MC. The uncertainty related with the ISR correction factor is studied by replacing the input cross section line shape with the latest results from BaBar Y4360BaBar in KKMC, and the change on $(1+\delta^{v})~{}\epsilon$ is taken as an uncertainty. The uncertainty in the vacuum polarization factor is 0.5% taken from a QED calculation vacuum . The uncertainty on the integrated luminosity is 1%, determined with large angle Bhabha events luminosity . The uncertainties in the decay branching fractions of intermediate states are quoted from the PDG pdg . The uncertainty from others sources, such as lepton separation, trigger efficiency, and FSR are negligible, and are conservatively taken to be 1.0%. Assuming all sources of systematic are independent, the total uncertainties are obtained by adding the individual values in quadrature, and are in a range between $7.7\%$ to $14.1\%$ and $7.4\%$ to $20.1\%$, depending on c.m. energy, for mode I and II, respectively.

The measured Born cross sections of e^{+}e^{-}\to\pi^{+}\pi^{-}$$\psi(3686) at individual c.m. energies for the two $\psi(3686)$ decay modes are consistent with each other within their uncertainties. The measurements are therefore combined by considering the correlated and un-correlated uncertainties between the two modes, according to Refs. combine1 ; combine2 . The comparison of the combined Born cross section of $e^{+}e^{-}\to\pi^{+}\pi^{-}\psi(3686)$ with those from previous experimental results is shown in Fig. 2. The results are consistent with former experiments, and have much improved precision.

V Fit to the cross section

To study the possible resonant structures in $e^{+}e^{-}\to\pi^{+}\pi^{-}\psi(3686)$, a binned $\chi^{2}$ fit is applied to describe the cross section obtained in this analysis in a energy region from 4.085 GeV to 4.600 GeV. Assuming that three resonances exist, the PDF can be parameterized as

[TABLE]

where $f_{1}$ is for the $Y(4220)$, $f_{2}$ is for the $Y(4390)$, $f_{3}$ is for the $Y(4660)$, $\phi_{1}$ is the phase angle between $Y(4390)$ and $Y(4220)$, and $\phi_{2}$ is the phase angle between $Y(4390)$ and $Y(4660)$.

The amplitude $f_{i}~{}(i=1,2,3)$ for each resonance is a $P$-wave Breit-Wigner function, defined as

[TABLE]

where $M_{i}$ is the mass of the resonance, $\Gamma^{e^{+}e^{-}}_{i}$ is the partial width to $e^{+}e^{-}$, $\Gamma_{i}$ is the total width, $\mathcal{B}_{i}$ is the branching fraction of the resonance decaying to $\pi^{+}\pi^{-}\psi(3686)$, and $\Phi(M_{i})$ is the three-body PHSP factor. In the fit, $M_{i}$, $\Gamma_{i}$, the product $\mathcal{B}_{i}\Gamma^{e^{+}e^{-}}_{i}$, and the relative phase angle $\phi_{i}$ $(i=1,2)$ are free parameters. The parameters of the $Y(4660)$ are fixed to Belle’s latest results Y4360Belle .

The fit result is shown in Fig. 2 as the solid curve and the resulting parameters of the resonances are summarized in Table 2. There are two solutions with equally good fit quality, $\chi^{2}/\text{ndf}=9.98/7$, where ndf is the number of degrees of freedom. An alternative fit without $Y(4220)$ is performed as shown by the dashed curve in Fig. 2, yielding $\chi^{2}/ndf=54.54/11$. The significance of the $Y(4220)$ is determined to be $5.8\sigma$.

The systematic uncertainties accociated with the resonance parameters include the following. The uncertainty from the c.m. energy measurement is studied by taking the uncertainty of the c.m. energy $\Delta(\sqrt{s})$=0.8 GeV cme into consideration. The uncertainty due to the c.m. energy spread is estimated by convoluting the fit PDF with a Gaussian function with a width of 1.6 MeV, which is the energy spread determined by the Beam Energy Measurement System espread . The uncertainty from the fit range is investigated by including the first energy point $\sqrt{s}=4.008$ GeV or excluding the energy point $\sqrt{s}=4.085$ GeV. The uncertainty from the $Y(4660)$ resonance is studied by varying the parameters of the $Y(4660)$ within its uncertainties. Table 3 summarizes the systematic uncertainties for the mass and total width of the resonances $Y(4220)$ and $Y(4390)$.

A comparison of the masses and widths of the $Y(4220)$ and $Y(4390)$ determined in this work to those of other vector states observed at BESIII in the mass region between 4.2 and 4.4 GeV/$c^{2}$ are shown in Table 4.

VI Study of intermediate states

Intermediate states in the process $e^{+}e^{-}\to\pi^{+}\pi^{-}\psi(3686)$ are investigated in the data samples that have large integrated luminosity. A requirement of 3.68<M(\pi^{+}\pi^{-}J/\psi)(M^{\rm recoil}(\pi^{+}\pi^{-}))<3.70\;\mbox{GeV/c^{2}} is applied to extract the $\psi(3686)$ signal, and the sideband regions, 3.63<M(\pi^{+}\pi^{-}J/\psi)(M^{\rm recoil}(\pi^{+}\pi^{-}))<3.65\;\mbox{GeV/c^{2}} or 3.73<M(\pi^{+}\pi^{-}J/\psi)<3.75\;\mbox{GeV/c^{2}}, are used to explore the potential non-$\psi(3686)$ backgrounds, where only the left side sideband region is used in mode II since a long tail appears on the right side of the $\psi(3686)$ signal due to ISR. The non-$\psi(3686)$ backgrounds are found to be small, and do not produce peaks in the various distributions.

With the above selection criteria, the Dalitz plots of $M^{2}(\pi^{+}\pi^{-})$ versus $M^{2}(\pi^{\pm}\psi(3686))$ and the corresponding one-dimensional projections are shown in Fig. 3 for data samples at $\sqrt{s}=4.416,4.358,4.258$ and $4.226\;\mbox{GeV}$, individually, where the plots include the candidates of the two $\psi(3686)$ decay modes. It should be noted that the two higher energy points are from the $Y(4390)$ peak region, while the other two are from the $Y(4220)$ peak region. These four points are marked with vertical arrows in Fig. 2.

For data at $\sqrt{s}=4.416\;\mbox{GeV}$, a prominent narrow structure is observed around 4030\;\mbox{MeV/c^{2}} in the $M(\pi^{\pm}\psi(3686))$ spectrum. The structure is also evident in the corresponding Dalitz plot, but it appears to be more complex when looking at high and low $M(\pi^{+}\pi^{-})$ ranges. In the low $M(\pi^{+}\pi^{-})$ region, there seem to be two separate $\pi^{\pm}\psi(3686)$ structures, presumably corresponding to a physical structure and its kinematic reflection. But in the high $M(\pi^{+}\pi^{-})$ region, only one broad $\pi^{\pm}\psi(3686)$ structure is seen. For data at $\sqrt{s}=4.358\;\mbox{GeV}$, there is no obvious structure observed in the $M(\pi^{\pm}\psi(3686))$ spectrum, but a cluster of events appear in the low $M(\pi^{+}\pi^{-})$ region on the corresponding Dalitz plot. It is worth noting that, at this c.m. energy, a physical structure with a mass of 4030\;\mbox{MeV/c^{2}} in the $M(\pi^{\pm}\psi(3686))$ spectrum has a reflection at the same mass position. For data at $\sqrt{s}=4.258\;\mbox{GeV}$, two bumps around 3900 and 4030 MeV/$c^{2}$ are visible in both the Dalitz plot and in the $M(\pi^{\pm}\psi(3686))$ spectrum. For data at $\sqrt{s}=4.258\;\mbox{GeV}$, the possible structures with masses of 3900 and 4030 MeV/$c^{2}$ in the $M(\pi^{\pm}\psi(3686))$ spectrum have kinematic reflections at each other’s mass positions. For data at $\sqrt{s}=4.226\;\mbox{GeV}$, no structure is clearly seen, which is very different from the behavior at the energy point close by, $\sqrt{s}=4.258~{}\mbox{GeV}$. A further striking feature for data at $\sqrt{s}=4.226\;\mbox{GeV}$ is a very different $M(\pi^{+}\pi^{-})$ distribution from those at the other three energy points.

To characterize the structure observed on the $M(\pi^{\pm}\psi(3686))$ spectrum for data at $\sqrt{s}=4.416\;\mbox{GeV}$, an unbinned maximum likelihood fit is carried out on the Dalitz plot of $M^{2}(\pi^{+}\psi(3686))$ versus $M^{2}(\pi^{-}\psi(3686))$ (denoted as $x$ and $y$ in formula 4). Assuming an intermediate state with spin parity $1^{+}$, the Dalitz plot is parameterized by the coherent sum of the process with an intermediate state and the direct process $e^{+}e^{-}\to\pi^{+}\pi^{-}\psi(3686)$. The PDF of the intermediate state is described with an $S$-wave Breit-Wigner function without considering interference among the charged conjugate modes,

[TABLE]

where $p$ ($q$) is the $\psi(3686)$ (intermediate state) momentum in the $\pi^{\pm}\psi(3686)$ (initial $e^{+}e^{-}$) rest frame, and $M_{R}$ and $\Gamma$ are the mass and width of the intermediate state. The 2-dimensional mass resolution and the detection efficiency, determined from MC simulations, are incorporated in the PDF for the intermediate states in the fit. The PDF of the direct process $e^{+}e^{-}\to\pi^{+}\pi^{-}\psi(3686)$ is taken from a MC-simulated shape using the Jpipi model, and that of the non-$\psi(3686)$ background is described with the distribution of events in the $\psi(3686)$ sideband region. A simultaneous fit constraining the mass and width of the intermediate state is carried out by maximize the product of the likelihood values of the two $\psi(3686)$ decay modes. The fit process is validated using MC samples. The distributions of mass resolution and detection efficiency are provided in the appendix.

The fit yields a mass of M=4032.1\pm 2.4\;\mbox{MeV/c^{2}} and a width of $\Gamma$ = $26.1\pm 5.3\;\mbox{MeV}$ for the intermediate state with a significance of 9.2$\sigma$, evaluated by comparing the likelihood values with or without the intermediate states included. The fit projections on $M^{2}(\pi^{\pm}\psi(3686))$ and $M^{2}(\pi^{+}\pi^{-})$ for data at $\sqrt{s}=4.416\;\mbox{GeV}$ are shown in Fig. 3. It can be seen that the overall fit curve does not match the peaking structure on the $M(\pi^{\pm}\psi(3686))$ spectrum in data, and the corresponding confidence level (C.L.) of the fit is only 8%, as estimated by toy-MC tests. Alternative fits with different assumptions of the spin parity of the intermediate state, including the interference among the charge conjugated modes, and including the contribution of $Z_{c}(3900)^{\pm}$ are explored. In these fits the parameters of the intermediate state are close to the norminal fit result and the fit qualities are not improved significantly. As shown in the Dalitz plot, the behavior of the structure is very different between the high $M(\pi^{+}\pi^{-})$ region and the low $M(\pi^{+}\pi^{-})$ region. A similar fit to data with the additional requirement M^{2}(\pi^{+}\pi^{-})>0.3\;(\mbox{GeV/c^{2}})^{2} is performed, which yields a mass of M=4030.3\pm 0.1\;\mbox{MeV/c^{2}} and a width of $\Gamma=5.1\pm 0.2\;\mbox{MeV}$. The corresponding projection of the fit and data on the $M^{2}(\pi^{\pm}\psi(3686))$ distribution is shown in Fig. 4, and the fit C.L. is 50%.

Similar fits are carried out to data at $\sqrt{s}=4.358$ and $4.258\;\mbox{GeV}$, where the parameters of the intermediate state are fixed to those from the fit to data at $\sqrt{s}=4.416\;\mbox{GeV}$. The projections on $M^{2}(\pi^{\pm}\psi(3686))$ and $M^{2}(\pi^{+}\pi^{-})$ are also shown in Fig. 3. The statistical significance of the intermediate state is $3.6\sigma$ and $9.6\sigma$ for data at $\sqrt{s}=4.358$ and $4.258\;\mbox{GeV}$, respectively. For data at $\sqrt{s}=4.358$ GeV, as shown in the Dalitz plot, a cluster of events appears in the $M^{2}(\pi^{\pm}\psi(3686))$ spectrum at low $M^{2}(\pi^{+}\pi^{-})$, which also shows a hint of an intermediate structure. As mentioned previously, for data at $\sqrt{s}=4.258\;\mbox{GeV}$, the structures with masses of 3900 and 4030 MeV/$c^{2}$ on the $M(\pi^{\pm}\psi(3686))$ spectrum are kinematic reflections of each other, so the fit results are strongly dependent on whether the $Z_{c}(3900)^{\pm}$ is included in the fit or not. For data at $\sqrt{s}=4.226\;\mbox{GeV}$, an intermediate state with a mass of 4030 MeV and its reflection are very close to the kinematic boundary of the three-body decay, so no obvious peak is observed in the $M(\pi^{\pm}\psi(3686))$ spectrum. The anomalous distribution on the $M(\pi^{+}\pi^{-})$ spectrum is, however, discussed in Ref. zoubs . For the other energy points with high statistics, such as $\sqrt{s}=4.387$ and 4.600 GeV, the Dalitz plots and the distribution of $M^{2}(\pi^{\pm}\psi(3686))$ and $M^{2}(\pi^{+}\pi^{-})$ are shown in Fig. 5.

VII SUMMARY

In summary, based on 5.1 fb*-1* of $e^{+}e^{-}$ collision data with c.m. energies between 4.008 and 4.600 GeV, collected with the BESIII detector, the Born cross sections of $e^{+}e^{-}\to\pi^{+}\pi^{-}\psi(3686)$ are measured. The measured cross sections are in good agreement with those from Babar and Belle Y4360BaBar ; Y4360Belle , but with significantly improved precision. The cross section is fit with three coherent Breit-Wigner functions. The parameters of two resonances, the $Y(4220)$ and $Y(4390)$, are determined from the fit to be $M=4209.5\pm 7.4\pm 1.4$ MeV/$c^{2}$ and $\Gamma=80.1\pm 24.6\pm 2.9$ MeV for the $Y(4220)$, and $M=4383.8\pm 4.2\pm 0.8$ MeV/c2 and $\Gamma=84.2\pm 12.5\pm 2.1$ MeV for the $Y(4390)$, where the first errors are statistical and the second systematic. The resonance $Y(4220)$ is observed in the process $e^{+}e^{-}\to\pi^{+}\pi^{-}\psi(3686)$ for the first time with a significance of $5.8\sigma$. Both of the resonances $Y(4220)$ and $Y(4390)$ are consistent with the vector resonances observed in $e^{+}e^{-}\to\pi^{+}\pi^{-}h_{c}$ pipihc . The $Y(4220)$ state is possibly the same state as observed in $e^{+}e^{-}\to\pi^{+}\pi^{-}J/\psi$ pipijpsi and $e^{+}e^{-}\to\omega\chi_{c0}$ omegachicj , and an attempt was made with a combined fit to the $e^{+}e^{-}\to\omega\chi_{c0}$, $\pi^{+}\pi^{-}J/\psi$, $\pi^{+}\pi^{-}h_{c}$, and $\pi^{+}D^{0}D^{*-}+c.c$ cross sections measured by BESIII to determine the resonant parameters of the $Y(4220)$ combinefit . Peaking at a similar mass was also observed in the process $e^{+}e^{-}\to\eta J/\psi$ etajpsi at BESIII.

A charged charmonium-like structure is observed in the $M(\pi^{\pm}\psi(3686))$ spectrum for data at $\sqrt{s}=4.416\;\mbox{GeV}$. A fit to the structure, assuming the existence of a spin-parity $1^{+}$ charmonium-like state, yields a mass M=4032.1\pm 2.4\;\mbox{MeV/c^{2}}, where the errors are statistical only. However, the fit cannot describe data well. The width of the intermediate state varies over a wide range for different kinematic regions of data at $\sqrt{s}=4.416$ GeV. Similar fits are carried out to data at $\sqrt{s}=4.258$ and $4.358\;\mbox{GeV}$, where there is also evidence for the new structure at 4030 MeV/$c^{2}$ in the $M(\pi^{\pm}\psi(3686))$ spectra. For data at $\sqrt{s}=4.226\;\mbox{GeV}$, no obvious structure in the $M(\pi^{\pm}\psi(3686))$ spectrum is observed, but the $M(\pi^{+}\pi^{-})$ spectrum shows an anomalous distribution. The measured mass of the intermediate state deviates from that of the structure observed by Belle with their full data sample Y4360Belle by over $3\sigma$. Larger data sets, as well as additional theoretical input, are necessary for further understanding of the intermediate structures in $e^{+}e^{-}\to\pi^{+}\pi^{-}\psi(3686)$.

VIII ACKNOWLEDGMENTS

The BESIII collaboration thanks the staff of BEPCII, the IHEP computing center and the supercomputing center of USTC for their strong support. This work is supported in part by National Key Basic Research Program of China under Contract No. 2015CB856700; National Natural Science Foundation of China (NSFC) under Contracts Nos. 11235011, 11322544, 11335008, 11375170, 11275189, 11425524, 11475164, 11475169, 11625523, 11605196, 11605198, 11635010; the Chinese Academy of Sciences (CAS) Large-Scale Scientific Facility Program; the CAS Center for Excellence in Particle Physics (CCEPP); Joint Large-Scale Scientific Facility Funds of the NSFC and CAS under Contracts Nos. U1332201, U1532257, U1532258, U1532102; CAS under Contracts Nos. KJCX2-YW-N29, KJCX2-YW-N45, QYZDJ-SSW-SLH003; 100 Talents Program of CAS; National 1000 Talents Program of China; INPAC and Shanghai Key Laboratory for Particle Physics and Cosmology; German Research Foundation DFG under Contracts Nos. Collaborative Research Center CRC 1044, FOR 2359; Istituto Nazionale di Fisica Nucleare, Italy; Joint Large-Scale Scientific Facility Funds of the NSFC and CAS; Koninklijke Nederlandse Akademie van Wetenschappen (KNAW) under Contract No. 530-4CDP03; Ministry of Development of Turkey under Contract No. DPT2006K-120470; National Natural Science Foundation of China (NSFC); National Science and Technology fund; The Swedish Resarch Council; U. S. Department of Energy under Contracts Nos. DE-FG02-05ER41374, DE-SC-0010118, DE-SC-0010504, DE-SC-0012069; University of Groningen (RuG) and the Helmholtzzentrum fuer Schwerionenforschung GmbH (GSI), Darmstadt; WCU Program of National Research Foundation of Korea under Contract No. R32-2008-000-10155-0.

Appendix A Distribution of $\epsilon(x,y)$ and $\sigma(x,y)$

In Eq. 5, the fit procedure can be carried out by the Breit-Wigner function for $M^{2}(\pi^{+}\psi(3686))$ and $M^{2}(\pi^{-}\psi(3686))$ (denoted as $x$ and $y$) convoluted with a resolution function $\sigma(x)$, $\sigma(y)$, respectively, and then each of them multiplies by a flat polynomial on $y$ and $x$ to form a 2-dimensional function. The signal shape can then be obtained by multiplying a 2-dimensional resolution function $\epsilon(x,y)$ after adding the above two 2-dimensional functions together.

[TABLE]

The resolution for distributions $M^{2}(\pi^{+}\psi(3686))$ and $M^{2}(\pi^{-}\psi(3686))$ are shown in Fig. 6 for c.m. energies $\sqrt{s}=$4.416, 4.358, 4.258 and 4.226 GeV. It can be parameterised by a double-Gaussian with the parameters listed in Table 5. The distributions of 2-dimensional efficiency function are shown in Fig. 7.