Study of the resonance structures in $e^{+}e^{-}\rightarrow\pi^{+}\pi^{-}J/\psi$ process

A series of charmonium-like states, commonly referred to as $XYZ$ states, has been observed in the last two decades. The $Y(4260)$ was first observed by BaBar using an initial state radiation (ISR) technique in the process $e^{+}e^{-}\rightarrow\gamma_{\rm{ISR}}\pi^{+}\pi^{-}J/\psi$ [1] and was soon confirmed by Belle in the same process [2]. In 2017, the BESIII collaboration was able to resolve the $Y(4260)$ structure into a combination of two resonances, the $Y(4220)$ and the $Y(4320)$, using the world’s largest sample of $e^{+}e^{-}\rightarrow\pi^{+}\pi^{-}J/\psi$ events [3]. Moreover, similar structures have been observed in the processes $e^{+}e^{-}\rightarrow\pi^{+}\pi^{-}\psi(2S)$ [4], $\pi^{+}\pi^{-}h_{c}$ [5], $\pi^{+}D^{0}D^{\ast-}$ [6], $\eta J/\psi$ [7] and $\omega\chi_{cJ}$ $(J=0,\ 1,\ 2)$ [8, 9] with the BESIII data. The parameters of the two resonances are to a large extent consistent among reactions. However, to understand whether the structures observed in different final states are indeed the same, more investigations are needed.

Conventional charmonium states, such as the $\psi(4040)$, $\psi(4160)$, and $\psi(4415)$, mainly decay into open charm final states ($D^{(\ast)}\bar{D}^{(\ast)}$), while the $Y$ states show strong coupling to the hidden-charm final states [10]. The number of observed vector states in this energy region exceeds that of the predicted vector charmonium states [11]. These features suggest that some of these supernumerary vector states are candidates of an exotic nature, such as hybrid, tetraquark states, or mesonic molecules  [12, 13, 14, 15, 16, 17]. To clarify the nature of these states and to distinguish between the different theoretical models, precise measurements of the production cross section and of the resonance parameters are essential.

In this paper, an improved measurement of the energy-dependent $e^{+}e^{-}\rightarrow\pi^{+}\pi^{-}J/\psi$ cross section at center-of-mass energies ($\sqrt{s}$) between 3.7730 to 4.7008 GeV with a total integrated luminosity of about 23 fb^-1 is presented. The $XYZ$ data and the $R$-scan data (between $\sqrt{s}=4.410$ and $4.590$ GeV) analyzed in Ref. [3] are also analyzed here, and thus these measurements are correlated and can not be used in combination. Improvements in this paper include the use of additional data in the $Y(4220)$/$Y(4320)$ mass region, allowing us to study these two states in more detail, as well as the inclusion of more detailed background studies and systematic studies of the event selection.

The BESIII detector [18] records $e^{+}e^{-}$ collisions provided by the BEPCII storage ring [19] in a center-of-mass energy range from 2.0 to 4.9 GeV. The cylindrical core of the detector covers 93% of the full solid angle and consists of a helium-based multilayer drift chamber (MDC), a plastic scintillator time-of-flight system (TOF), and a CsI(Tl) electromagnetic calorimeter (EMC), which are all enclosed in a superconducting solenoidal magnet providing a 1.0 T magnetic field. The solenoid is supported by an octagonal flux-return yoke with resistive plate counter muon identification modules interleaved with steel. The charged-particle momentum resolution at $1~{}{\rm GeV}/c$ is $0.5\%$, and the $\mathrm{d}E/\mathrm{d}x$ resolution is $6\%$ for electrons from Bhabha scattering. The EMC measures photon energies with a resolution of $2.5\%$ ($5\%$) at $1$ GeV in the barrel (end cap) region. The time resolution in the TOF barrel region is 68 ps, while that in the end cap region is 110 ps. The end cap TOF system was upgraded in 2015 using multi-gap resistive plate chamber technology, providing a time resolution of 60 ps [20].

The data samples analyzed in this article are listed in appendix. They include 40 energy points (referred to as the $XYZ$ data sample) from 3.7730 to 4.7008 GeV with an integrated luminosity of more than 40 $\rm{pb^{-1}}$ at each $\sqrt{s}$, and 13 energy points (the $R$-scan data sample) from 4.410 to 4.590 GeV with a luminosity of $7-9$ $\rm{pb^{-1}}$ at each $\sqrt{s}$ to study possible structure around 4.5 GeV. The integrated luminosities and $\sqrt{s}$ are measured using Bhabha scattering, radiative di-muon events, and $\Lambda_{c}$ pairs [21, 23, 22].

The geant4-based [24] Monte Carlo (MC) simulation software packages boost [25] and evtgen [26] are used to determine detection efficiencies and to estimate the background contributions. The generator kkmc [27] is used to model the beam energy spread and the ISR emission in $e^{+}e^{-}$ annihilations. Final state radiation (FSR) from charged (final state) particles is incorporated using the photos package [28]. The signal MC samples of $e^{+}e^{-}\rightarrow\pi^{+}\pi^{-}J/\psi$, with $J/\psi\rightarrow\ell^{+}\ell^{-}~{}(\ell=e,~{}\mu)$, are generated using a phase-space (PHSP) model and weighted following the results of a partial wave analysis (PWA) of the data. The amplitude model and the PWA framework are the same as in Refs. [29, 30, 31]. Potential background contaminations are studied via inclusive MC samples described in Refs. [3, 32]. Bhabha, di-muon, and $e^{+}e^{-}\rightarrow 2(\pi^{+}\pi^{-})$ events are chosen as control samples to study the tracking efficiencies of $e^{\pm}$, $\mu^{\pm}$, and $\pi^{\pm}$, respectively. An exclusive MC sample for $e^{+}e^{-}\rightarrow e^{+}e^{-}\mu^{+}\mu^{-}$ (with $\gamma^{\ast}\gamma^{\ast}\rightarrow\mu^{+}\mu^{-}$, two-photon process) [33, 34] is generated for Boosted Decision Tree (BDT) training [35] used for background suppression.

Events with exactly four charged tracks and zero net charge are selected. In order to be considered, a charged track reconstructed in the MDC is required to satisfy $|\rm{cos}\theta|<0.93$, with $\theta$ being the polar angle, and the distance of closest approach to the $e^{+}e^{-}$ interaction point must be within $\pm 10$ cm in the beam direction and within 1 cm in the plane perpendicular to the beam direction. If there are only two tracks with momentum greater than 1.0 GeV$/c$, then the two tracks are judged as leptons; if there are three tracks with momentum greater than 1.0 GeV$/c$, then the two tracks with opposite charges and with an invariant mass closer to the $J/\psi$ mass [36] are regarded as leptons. The other two tracks are regarded as pions. The energy deposited in the EMC is used to distinguish between electrons and muons. For muon candidates, the deposited energy is required to be less than 0.4 GeV, while for electrons, it is required to be larger than 1.1 GeV. In order to suppress the background contribution and to improve the energy and momentum resolution, a four-constraint (4C) kinematic fit is applied to the event with the hypothesis $e^{+}e^{-}\rightarrow\pi^{+}\pi^{-}\ell^{+}\ell^{-}$, which constrains the total four-momentum of the final state particles to that of the initial colliding beams. The $\chi^{2}$ of the kinematic fit is required to be less than 60. The cosine of the opening angle of the pion-pair ($\rm cos(\pi^{+}\pi^{-})$) and, only in the $J/\psi\rightarrow e^{+}e^{-}$ decay channel, of the pion-electron pairs ($\rm cos(\pi^{\pm}e^{\mp})$) are required to be less than 0.98 to suppress gamma conversion background events of the radiative Bhabha and di-muon background contributions [3].

In order to reduce the background contributions from low momentum electrons, pions are identified using the $\mathrm{d}E/\mathrm{d}x$ information recorded by the MDC. A discriminator $\chi_{\pi^{+/-}}\equiv(\mu_{\rm m}-\mu_{\rm exp})/\sigma_{\rm m}$ is defined by combining the measured $\mathrm{d}E/\mathrm{d}x$ value ($\mu_{\rm m}$), the measurement uncertainty ($\sigma_{\rm m}$), and the expected value under a pion hypothesis ($\mu_{\rm exp}$). The conditions $\chi_{\pi^{+}}<2.8~{}(3.2)$ and $\chi_{\pi^{-}}<3.0~{}(4.10)$ for the $e^{+}e^{-}$ mode ($\mu^{+}\mu^{-}$ mode) are found to provide an optimal balance between signal efficiency loss and background rejection power. The polar angle distributions of $e^{+}$ and $e^{-}$ in the $J/\psi\rightarrow e^{+}e^{-}$ mode (Fig. 1) clearly show a significant contribution from the $e^{+}e^{-}\rightarrow e^{+}e^{-}\mu^{+}\mu^{-}$ process (a two-photon process corresponding to the yellow filled area in Fig. 1). A boosted decision tree method, implemented within the TMVA framework [35], is trained to efficiently suppress this background.

The energy deposited in the EMC, the time of flight from TOF, $\mathrm{d}E/\mathrm{d}x$, and the opening angle of the pion and electron candidates (Fig. 2) are used as input variables to the BDT, as they are mostly uncorrelated, as can be seen from Figs. 3(a) and (b). The BDT is trained using the MC simulation of the process $e^{+}e^{-}\rightarrow e^{+}e^{-}\mu^{+}\mu^{-}$. The signal MC and background MC simulation of multiple energy points are combined to train the BDT model. Figure 3(c) shows that the training samples and test samples in the BDT model are in good agreement, and the model can effectively distinguish between signal and background. The response of the BDT is required to maximize the $S/\sqrt{S+B}$ distribution (where the $S$ and $B$ indicate the number of signal and background events passing the condition) shown in Fig. 3(d).

After applying the event selection criteria mentioned above, a clear $J/\psi$ signal peak is observed in the invariant mass distribution of the lepton pairs ($M(\ell^{+}\ell^{-})$), as shown in Figs. 4(a)($e^{+}e^{-}$ mode) and (b)($\mu^{+}\mu^{-}$ mode), Figs. 4(c)(data) and (d)(MC) show the two-dimensional distribution of the invariant mass of the $\pi^{+}\pi^{-}$ and $\pi^{+}J/\psi$ pairs. To determine the signal yields, an unbinned maximum-likelihood fit to $M(\ell^{+}\ell^{-})$ is performed. The signal probability density function (PDF) is defined as the convolution between the MC signal and a Gaussian function, while the background contribution is parameterized with a linear term.

The dressed Born cross section of $e^{+}e^{-}\rightarrow\pi^{+}\pi^{-}J/\psi$ at a given $\sqrt{s}$ is calculated with

where $N^{\rm obs}$ is the number of signal events, $\mathcal{L}_{\rm int}$ is the integrated luminosity, $\epsilon$ is the selection efficiency, $\mathcal{B}$ is the branching fraction of $J/\psi\rightarrow\ell^{+}\ell^{-}$ [36] and ($1+\delta$) is the ISR correction factor. The final cross section is given by the average of the $e^{+}e^{-}$ and the $\mu^{+}\mu^{-}$ modes, weighted with the inverse of the respective statistical uncertainties. The parameters of Eq. 1 and the measured cross sections at the different $\sqrt{s}$ are summarized in the appendix.

The ISR correction factor is evaluated using an iterative procedure, in which the $\sigma(\sqrt{s})$ is initially assumed to be simply flat and iteratively recomputed until the difference between the final two iterations is less than 0.1%. The $(i+1)^{\rm th}$ iteration $(1+\delta)^{i+1}$ is given by [37]:

where $W_{j}^{i}$ is the weighting factor for the $j^{\rm th}$ event. The ISR correction factor $(1+\delta)^{0}$ corresponds to the $0^{\rm th}$ iteration, i.e. to an assumed flat shape of the cross section; $\sigma^{i}$ is the line-shape given by the $i^{\rm th}$ iteration and $\sigma^{0}$ is calculated with the correction factor $(1+\delta)^{0}$, which is calculated using the kkmc [27] program; $m_{j}$ is the invariant mass of the final state of the event $j$; and $N$ is the total number of generated MC events. Considering that the sizes of the $R$-scan data samples are small, only the $XYZ$ data samples are used in the iteration. The ISR correction factors for the $R$-scan data samples are calculated from the model describing the cross section in Sect. VI.

The main contributions to the systematic uncertainty of the cross section measurement are related to the measurement of the luminosity, the MC model, the tracking efficiency, the ISR correction, the branching fractions of $J/\psi$ decays, the kinematic fit, the fit to $M(\ell^{+}\ell^{-})$, and the BDT method. The integrated luminosities of all data sets are measured using large angle Bhabha scattering events with an uncertainty of 0.66% [21]. The branching fraction of $J/\psi\rightarrow\ell^{+}\ell^{-}$ is taken from PDG [36], and the uncertainties are 0.6% for each mode. The uncertainty related to the kinematic fit is estimated by adjusting the helix parameters of the charged tracks in MC simulation to match the $\chi^{2}$ distribution of data and MC [38]. The uncertainty is determined as the difference between the results before and after the helix parameters correction, resulting in 0.76% (1.34%) for the $\mu^{+}\mu^{-}$ ($e^{+}e^{-}$) mode. The uncertainty on the signal yield arising from the fitting of the $M(\ell^{+}\ell^{-})$ distribution is obtained by varying the fit range and by changing the background modeling from a first to a second order polynomial, which leads to a difference of 2.44% (1.16%) that is taken as the systematic uncertainty introduced by the fit method.

In the evaluation of the systematic uncertainties due to the tracking efficiency, the MC simulation model, the BDT method, and the ISR correction, 300 sets of Gaussian samplings according to the related central value ($\mu_{0}$) and statistical uncertainties, obtained from the control samples, the PWA fitted parameters, the data-MC simulation ratios of the BDT training variables, and the ISR correction factors, are generated. A Gaussian function ($\mu_{1},~{}\sigma_{1}$) is used to fit the resulting distributions from the 300 MC simulation samples. The final uncertainty of each source is given by $(|\mu_{1}-\mu_{0}|+\sigma_{1})/\mu_{0}\times 100\%$.

The uncertainty of the tracking efficiency comes from the uncertainties of the correction factors obtained from control samples, which are the Bhabha, di-muon and $e^{+}e^{-}\rightarrow 2(\pi^{+}\pi^{-})$ processes. The correction factors are given by the ratio of data to MC simulation in the two dimensions of transverse momentum and the cosine of the polar angle. Therefore, the Gaussian sampling is based on the correction value and its uncertainty, and obtains the efficiencies distributions with respect to the new correction factors. In the MC simulation model, the uncertainties come from the uncertainties of the amplitude parameters given by the PWA fit to the data. To estimate the uncertainties of the parameters, toy MC samples are generated according to the fitted parameters and the error matrix to obtain the efficiency distribution for each energy point. The difference in the efficiency determined with pseudoexperiments generated according to the amplitude model parameters with and without the three-body PHSP process is regarded as the uncertainty of the partial wave model.

The differences in line shape caused by statistical uncertainties lead to different ISR correction factors. Therefore, for each energy point, the cross section and error serve as the parameters of the Gaussian function of sampling, and after a new round of iteration, the cross section distribution is obtained for each point. Apart from the systematic uncertainty of the ISR correction, the systematic uncertainties of the other contributions are approximately the same for every energy point, and thus the average of each is taken, as listed in Table 1. The ISR correction uncertainty of $R$-scan data samples are estimated by the value of closest $XYZ$ data point, and the uncertainty of each point is listed in the appendix. Assuming all the sources to be independent, the total systematic uncertainties are obtained by adding them in quadrature.

To study the possible resonant structures in the $e^{+}e^{-}\rightarrow\pi^{+}\pi^{-}J/\psi$ process, a maximum likelihood fit is preformed to the measured cross section. The likelihood is constructed assuming the number of events satisfies a Gaussian distribution in $XYZ$ data and Poisson distribution in $R$-scan data. The cross section is parameterized with a coherent sum of Breit-Wigner (BW) functions. Due to the lack of data near the $\psi(3770)$ resonance, it is not feasible to determine the relative phase between the $\psi(3770)$ amplitude and other amplitudes. The cross section line shape is described by

where $R_{\psi(3770)}$ is used to describe the $\psi(3770)$ resonance and its mass and width are fixed to the world average values [36]. The $i$ is the imaginary unit. $R_{j}$ represents the amplitude to describe a given resonant structure and $\phi_{j}$ is the corresponding phase. The phase $\phi_{0}$ is set to zero and the other phases are given relative to the $R_{0}$. For the structure near 4.0 GeV, two different parameterization methods are applied, Model I: a BW function, and Model II: an exponential function (Exp) of the form $R_{0}(\sqrt{s})=PS(\sqrt{s})e^{-p_{0}(\sqrt{s}-M_{\rm{threshold}})}p_{1}$ [39], with $M_{\rm{threshold}}=m_{\pi^{+}}+m_{\pi^{-}}+m_{J/\psi}$, $PS(\sqrt{s})$ is the PHSP factor of the three-body decay $R_{j}\rightarrow\pi^{+}\pi^{-}J/\psi$ [36], and $p_{0}$ and $p_{1}$ are free parameters determined by the fit. The number of resonances is denoted by $n$, comprising the known $Y(4220)$ and $Y(4320)$ as well as further possible structures. The amplitude $R_{j}$ is defined as

where $M_{j}$, $\Gamma^{\rm tot}_{j}$ and $\Gamma^{\rm ee}_{j}$ are the mass, full width and electronic width of resonance $R_{j}$, respectively, and $\mathcal{B}(R_{j})$ is the branching fraction for $R_{j}\rightarrow\pi^{+}\pi^{-}J/\psi$.

In case of considering the states $Y(4220)$ and $Y(4320)$ ($n=2$), multiple sets of solutions are obtained given by the two models (Model I: BW, Model II: Exp) at 4.0 GeV. The fit results are shown in Fig. 5, and the fit parameters are summarized in Table 2. Sizable differences between the fit results of Model I and Model II appear mainly in the energy region between 3.7730 and 4.1574 GeV. The difference of $\chi^{2}/\rm{ndf}$ is 3.72, where $\rm{ndf}$ is the number of degrees of freedom. Therefore, Model I is chosen to be the default model for the final cross section fit result. The cross section fit shows larger fluctuations at $\sqrt{s}=$ 3.8713 GeV. These might be due to the influence of the $X(3872)$ [36] resonance which was not included in the model since the $X(3872)$ is very narrow and there are not sufficient data points around its nominal mass.

Considering the distribution of the pull ($\chi$) values, the above two models do not describe the interval from 4.4 to 4.6 GeV very well. Therefore, a third BW function ($n=3$) is added to study whether this deviation is caused by possible additional structures. When the (fit) parameters of the third BW function are floated, two possible solutions are obtained, one with a mass close to the $\psi(4415)$ and the other one close to 4.5 GeV. Compared with the mode of $n=2$, the significance of these two solutions are $4.0\sigma$ ($3.6\sigma$) and $2.1\sigma$ ($2.7\sigma$), respectively. The alternative fits using the parameters of the $\psi(4415)$ [36] and the newly discovered $Y(4500)$ structure [40] have also been attempted to describe the structure at 4.5 GeV, and led to results with significance of $2.6\sigma$ ($3.1\sigma$) and $3.3~{}\sigma$ ($3.3~{}\sigma$), respectively. The numbers in the brackets correspond to an alternative fit, in which the BW function(Model I) is replaced by an exponential function(Model II) to describe the structure near 4.0 GeV, as also shown in Table 2 and Table 3.

Due to the limited data samples around 4.4 GeV, the mass and width of the new additional BW function are fixed to $\psi(4415)$ or $Y(4500)$, the fit result in Table 3. With larger/new data samples in that region, the structure of the $\psi(4415)$ or $Y(4500)$ can be studied further. The fit results are shown in Fig. 6, and the parameters are summarized in Table 3. In conclusion, the parameters of the $Y(4220)$ state are stable for the different models, while the ones of $Y(4360)$ heavily depend on the presence of an additional structure close to 4.5 GeV.

Relevant systematic uncertainties of the fit parameters are caused by the fit model, the $\sqrt{s}$ energy measurement and its spread, and the PHSP factor. The systematic uncertainty related to the fit model is evaluated as the difference of the mass and width results of Model I and Model II. The $\sqrt{s}$ of all data sets have been measured with di-muon events with an uncertainty of 0.6 MeV that propagates directly to the uncertainty of the mass of the resonances. The uncertainties included by the $\sqrt{s}$ spread are obtained by convolving the resonant PDF with a Gaussian function whose width is taken to be 1.6 MeV, equal to the spread obtained from the Beam Energy Measurement System [41]. The uncertainty of the PHSP factor, due to the existence of intermediate states, is estimated by considering the PHSP of cascade two-body decays of $e^{+}e^{-}\rightarrow RJ/\psi$ (with $R=\sigma,~{}f_{0}(980),~{}f_{0}(1370)$) and $e^{+}e^{-}\rightarrow\pi^{\pm}Z_{c}(3900)^{\mp}$ , and the maximum value of the difference with respect to the result obtained when using the three-body PHSP factor is taken as the systematic uncertainty. The deviation of the resonant parameters introduced by the uncertainties of the $\psi(3770)$ resonance parameters are less than 0.1 MeV, and thus can be neglected. Assuming all of the systematic uncertainties are independent, adding them in quadrature delivers the total error as listed in Table 4.

In summary, a precision measurement of the energy-dependent cross section for the process $e^{+}e^{-}\rightarrow\pi^{+}\pi^{-}J/\psi$ from $\sqrt{s}=3.7730$ GeV to $\sqrt{s}=4.7008$ GeV at BESIII is performed. This measurement improves upon the precision of the previous results from the BESIII collaboration in the same channel [3] by about 60% at the points of same statistics ($\sqrt{s}=4.2263$ and 4.2580 GeV). This is achieved by the improvements in the optimized MC simulation model and the enhanced tracking efficiencies.

The energy-dependent cross section is fitted with different fit models for the cross section line shape, allowing for the search of resonances and the evaluation of their parameters. It was found that the structure close to 4.0 GeV is better described when using the BW function as compared to an Exp function. The $Y(4220)$ and $Y(4320)$ resonances were observed with significances larger than $10\sigma$ and their resonance parameters were estimated to be consistent with those reported in Ref. [3]. However, the presence of an additional structure around 4.5 GeV, possibly identifiable with the $\psi(4415)$, influences the evaluation of the $Y(4320)$ parameters, which are ($M,~{}\Gamma$) = ($4298\pm 12\pm 26$ MeV/$c^{2}$, $127\pm 17\pm 10$ MeV), therefore, reported with a large uncertainty.

Table 5 lists the parameters used to calculate the cross section of the $e^{+}e^{-}\rightarrow\pi^{+}\pi^{-}J/\psi$ process with the Eq. 1.