Cross Section Measurement of $e^+e^- \rightarrow \pi^+\pi^-\psi(3686)$ from $\sqrt{s}=4.0076$~GeV to 4.6984~GeV

Using data samples with a total integrated luminosity of $20.1~\rm fb^{-1}$ collected by the BESIII detector operating at the BEPCII collider, the cross section of the process $e^+e^- \rightarrow \pi^+\pi^-\psi(3686)$ is measured at center-of-mass energies between 4.0076 and 4.6984 GeV. The measured cross section is consistent with previous results, but with much improved precision. A fit to the measured energy-dependent cross section, which includes three Breit-Wigner functions and a non-resonant contribution, confirms the existence of the charmonium-like states $Y(4220)$, $Y(4390)$, and $Y(4660)$. This is the first observation of the $Y(4660)$ at the BESIII experiment.


I. INTRODUCTION
Over the past several decades, a series of charmoniumlike states with J PC ¼ 1 −− , referred to as the Y states, have been discovered and confirmed by numerous experiments. However, the internal structure of the Y states remains controversial. Possible interpretations include hybrid mesons, meson molecules, hadrocharmonium, and tetraquarks, but none of them has so far proved conclusive [1]. Therefore, comprehensive studies of production and decay patterns, as well as the measurement of resonance parameters from different experiments, are desirable to provide information that will help probe the nature of the Y states.
Dedicated studies of Y states were initially triggered by the discovery of the Yð4260Þ in e þ e − → γ ISR π þ π − J=ψ using the initial-state-radiation (ISR) approach by the BABAR experiment [2], which was subsequently confirmed by CLEO [3,4] and Belle [5,6]. The latest results from the BESIII experiment on the channel e þ e − → π þ π − J=ψ reveal with unprecedented precision that the structure previously identified as the Yð4260Þ consists of two structures referred to as Yð4220Þ and Yð4320Þ [7]. This observation is also confirmed in the process e þ e − → π 0 π 0 J=ψ [8]. Furthermore, similar structures to the Yð4220Þ also appear in the processes e þ e − → ωχ c0 [9], π þ π − h c [10], π þ D 0 D Ã− [11], and ηJ=ψ [12]. However, the parameters of the Yð4220Þ observed in different processes are different. Further studies both theoretically and experimentally are still needed to understand the internal structure of the Yð4220Þ.
Analogously, studies of the process e þ e − → γ ISR π þ π − ψð3686Þ have been performed at the Belle and BABAR experiments, where two Y states, namely the Yð4360Þ and Yð4660Þ [13][14][15] were observed and confirmed. BESIII measured the cross sections of the process e þ e − → π þ π − ψð3686Þ with much improved precision at center-of-mass (c.m.) energies from 4.0076 to 4.5995 GeV [16]. The existence of the Yð4360Þ was confirmed, but with a mass close to 4.39 GeV=c 2 , so we refer to it here as the Yð4390Þ. Furthermore, the Yð4220Þ was reported for the first time in the final state π þ π − ψð3686Þ, but the Yð4660Þ was not studied due to the lack of data in the corresponding energy region. It is worth noting that a Yð4390Þ with close mass is also reported in the processes e þ e − → π þ π − h c [10], and ηJ=ψ [12] by BESIII. The differences in parameters of the Y states of different experiments require more studies to understand, and the e þ e − → π þ π − ψð3686Þ process is an appropriate channel to provide more information on the Y states.
In this paper, an analysis of e þ e − → π þ π − ψð3686Þ is presented using an analysis approach similar to that of Ref. [16], with old data sets used in Ref. [16] and new data sets collected by BESIII since 2017. Additional data samples between 4.23 and 4.36 GeV and new data samples between 4.60 and 4.70 GeV are used to study the properties of the Y states, especially the Yð4660Þ. The c.m. energy and the corresponding luminosity [8,17,18] of each of the samples used in this analysis can be found in Table I. Following Ref. [16], the ψð3686Þ is reconstructed via its charged decay mode ψð3686Þ → π þ π − J=ψ (Mode I) and the neutral decay mode ψð3686Þ → neutrals þ J=ψ (Mode II), where "neutrals" refers to π 0 ; π 0 π 0 or η from ψð3686Þ and γγ from the cascade decay ψð3686Þ → γχ cJ → γγJ=ψ. Only the neutral decays of π 0 and η are considered in TABLE I. The Born cross section σ B at individual c.m. energies. The subscript "cha" or "neu" denotes the charged mode (Mode I) or the neutral mode (Mode II), respectively. The uncertainties on the signal yields and the Born cross sections for individual decay modes are statistical only. The first uncertainties for the combined Born cross sections are statistical and the second are systematic. Mode II. The J=ψ is reconstructed via the decay modes J=ψ → l þ l − ðl ¼ e=μÞ.

II. THE BESIII EXPERIMENT AND THE DATA SETS
The BESIII detector [19] records symmetric e þ e − collisions provided by the BEPCII storage ring [20], which operates with a peak luminosity of 1 × 10 33 cm −2 s −1 in the c.m. energy range from 2.0 to 4.946 GeV. BESIII has collected large data samples in this energy region [21]. The cylindrical core of the BESIII 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 GeV=c is 0.5%, and the dE=dx 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 multigap resistive plate chamber technology, providing a time resolution of 60 ps [22].
Simulated data samples produced with a GEANT4based [23] Monte Carlo (MC) package, which includes the geometric description of the BESIII detector and the detector response, are used to determine detection efficiencies and to estimate background contributions. The simulation models beam energy spread and ISR effects in the e þ e − annihilations with the generator KKMC [24]. The signal process e þ e − → π þ π − ψð3686Þ is simulated by incorporating the line shape of the cross section from the previous study [16]. The inclusive MC sample includes the production of open charm processes, the ISR production of vector charmonium(-like) states, and the continuum processes incorporated in KKMC [24]. The known decay modes are modeled with EVTGEN [25] using branching fractions taken from the Particle Data Group (PDG) [26], and the remaining unknown charmonium decays are modeled with LUNDCHARM [27]. The QED final-state radiation (FSR) correction from charged final-state particles is incorporated using the PHOTOS package [28].

III. EVENT SELECTION
The topology of a signal event includes either six charged tracks (Mode I) or four charged tracks and at least two photons (Mode II). The good charged tracks detected in the MDC are required to be within a polar angle (θ) range of j cos θj < 0.93, where θ is defined with respect to the z axis. For charged tracks, the distance of closest approach to the interaction point must be less than 10 cm along the z axis, jV z j < 10 cm, and less than 1 cm in the transverse plane, jV xy j < 1 cm. Photon candidates are identified using showers in the EMC. The deposited energy of each shower must be more than 25 MeV in the barrel region (j cos θj < 0.80) and more than 50 MeV in the end cap region (0.86 < j cos θj < 0.92). To exclude showers that originate from charged tracks, the angle between the position of each shower in the EMC and the closest extrapolated charged track must be greater than 10 degrees. To suppress electronic noise and showers unrelated to the event, the difference between the EMC time and the event start time is required to be within (0, 700) ns.
To improve the efficiency, candidate events from Mode I are required to have five charged tracks with AE1 net charge, allowing for one missed charged pion, or six charged tracks with zero net charge. Candidate events from Mode II are required to have four charged tracks with zero net charge and at least two good photon candidates. For the signal candidates of both Modes, the pions and leptons are kinematically well separated, and thus the charged tracks with momenta above 1.0 GeV=c are assigned to be leptons, and those with momenta below 0.65 ð0.80Þ GeV=c are assigned to be pions for data samples with c.m. energies below (above) 4.465 GeV. Electron and muon separation is performed using the deposited energy (E) in the EMC, i.e., both electrons must satisfy E=p > 0.7, while both muons E < 0.45 GeV, where p is the momentum of charged tracks. Signal candidates are required to have a lepton pair of the same flavor and opposite charge. A vertex fit is applied to charged tracks from Mode I and Mode II, respectively. The results measured with data samples from J=ψ → e þ e − and J=ψ → μ þ μ − agree well with each other.
In Mode I, a four-constraint (4C) kinematic fit imposing energy-momentum conservation under the hypothesis e þ e − → π þ π − π þ π − l þ l − is implemented for the candidate events with six charged tracks, while a one-constraint (1C) kinematic fit constraining the missing mass to the pion mass is applied for those events with five charged tracks. The corresponding χ 2 is required to satisfy χ 2 4C < 60 and χ 2 1C < 15, respectively. Afterwards, a clear J=ψ signal is observed as shown in Fig. 1(a). Mðl þ l − Þ is calculated with the four-momentum of the lepton pair after the 4C (1C) kinematic fit. Candidate events are further required to satisfy 3.05 < Mðl þ l − Þ < 3.15 GeV=c 2 , shown as the region between the close red arrows in Fig. 1(a). To improve the mass resolution of ψð3686Þ, Mðl þ l − Þ is constrained to the PDG [26] value of the J=ψ mass, resulting in a kinematic fit with either five constraints (5C) or two constraints (2C). The events with Mðl þ l − Þ in the signal region are subjected to these 5C or 2C fits and the resulting kinematic variables are used in the subsequent analysis. Finally, the ψð3686Þ signal is extracted from the invariant mass of the l þ l − π þ π − system, Mðl þ l − π þ π − Þ, where there are four π þ π − combinations in an event. The one with Mðl þ l − π þ π − Þ closest to the ψð3686Þ nominal mass [26] is selected.
In Mode II, the number of photons in the final state is not fixed due to the fact that several ψð3686Þ decay modes are being included, and no kinematic fit is performed. The π þ π − opening angle is required to satisfy cos θ π þ π − < 0.9 in order to remove radiative Bhabha and dimuon background events in which the radiative photon converts into an e þ e − pair and is misidentified as a π þ π − pair. To remove the background due to e þ e − → π þ π − J=ψ, the recoiling mass of the π þ π − l þ l − system M rec ðπ þ π − l þ l − Þ, which corresponds to the invariant mass of the system of all neutral particles in the final state, is used. Based on exclusive MC   simulation, the veto is set at 3 times the width, and events with jM rec ðπ þ π − l þ l − Þj > 63 MeV=c 2 are accepted. To eliminate the background due to e þ e − → neutrals þ ψð3686Þ with the subsequent decay ψð3686Þ → π þ π − J=ψ, the ψð3686Þ mass with a J=ψ resolution correction where Mðψð3686ÞÞ is the ψð3686Þ mass from the PDG [26]. A requirement of jMðγγπ þ π − Þ − MðηÞj > 50 MeV=c 2 is used to eliminate background events from e þ e − → ηJ=ψ with a subsequent decay η → π þ π − π 0 , where γγ are the two photons with the largest energy and MðηÞ is the nominal η mass from the PDG [26]. The J=ψ signal is extracted from the Mðl þ l − Þ distribution by requiring 3.05 < Mðl þ l − Þ < 3.15 GeV=c 2 as indicated in Fig. 1(b). From the study of data sets with ffiffi ffi s p ¼ 4.1784, 4.2263, 4.2580, 4.4156 GeV, the small peak around 3.0 GeV=c 2 in Fig. 1(b) is determined to be due to statistical fluctuation. Finally, the ψð3686Þ signal is extracted from the invariant mass recoiling against the π þ π − system, M rec ðπ þ π − Þ. After applying all of the event selection criteria mentioned above, the distributions of Mðl þ l − π þ π − Þ for Mode I and M rec ðπ þ π − Þ for Mode II are shown in Fig. 2 for the data sample with a c.m. energy at 4.3583 GeV. The ψð3686Þ signals are clearly visible with small background contributions. Potential backgrounds without a J=ψ are studied using events in the J=ψ sideband region, where the J=ψ sideband regions are defined as 2.97 < Mðl þ l − Þ < 3.02 GeV=c 2 and 3.17 < Mðl þ l − Þ < 3.22 GeV=c 2 , as shown by the green arrows in Fig. 1. In Mode I, the 5C or 2C kinematic fit is not performed for the sideband events. The corresponding distributions of J=ψ sideband events are shown in Fig. 2, where no peak around the ψð3686Þ signal is observed. Potential background events including a J=ψ in the final state are studied using inclusive as well as dedicated exclusive MC samples, and no peak around the ψð3686Þ signal is found.

IV. EXTRACTION OF THE BORN CROSS SECTION
The Born cross section of the process is determined by where N obs is the observed signal yield, L int is the integrated luminosity, (1 þ δ) is the ISR correction factor obtained from the KKMC generator, 1 j1−Πj 2 is the vacuum polarization factor from QED calculations [29], The measured cross section for each c.m. energy is obtained by performing simultaneous unbinned maximum likelihood fits to the corresponding Mðl þ l − π þ π − Þ and M rec ðπ þ π − Þ spectra, where the two decay modes share the same cross section values. The signal shape is described with a MC-simulated shape convolved with a Gaussian function. The Gaussian function represents the mass resolution difference between data and MC simulation, and its parameters are allowed to float during the fit procedure. The background shape is a linear function. The fit curves for the data sample with c.m. energy 4.3583 GeV are shown in Fig. 2. The signal yields for Mode I and Mode II are also given by individual fits. The obtained cross section results for the two ψð3686Þ decay modes are consistent with each other. The cross section results and signal yields are summarized in Table I for the different c.m. energies. A comparison shown in Fig. 3 validates the consistency of cross section results among BESIII [16], Belle [14] and BABAR [15]. The BESIII results have the best precision.
Several sources of systematic uncertainties are considered in the measurement of the Born cross section and they can be classified into two categories: the shared uncertainties between the two ψð3686Þ decay modes and the individual ones. The uncertainties for the combined results of both ψð3686Þ decay modes are obtained by the following approaches: the shared uncertainties among the two ψð3686Þ decay modes, which are listed with the symbol " †" in Table II, are directly propagated to the combined   The shared uncertainties include those related to the luminosity, the vacuum polarization, the radiative corrections, the lepton tracking efficiency and the residual sources. The uncertainty of the integrated luminosity is 1% as obtained by analysing the large-angle Bhabha scattering events [30]. The uncertainty of the vacuum polarization factor is taken as 0.5% from the QED calculation in Ref. [31]. In Eq. (1), (1 þ δ) and ϵ depend on the input line shape of the cross section due to the effect of ISR. Thus, the measurement of the cross section is performed iteratively until the last variation of ð1 þ δÞϵ becomes negligible, where the measured cross section from this iteration is used as the input of the next iteration. Consequently, the uncertainty related to the ISR correction is dominated by the uncertainty of the input line shape of the Born cross section. To estimate the corresponding uncertainty, the ð1 þ δÞϵ is evaluated 100 times varying the input cross section line shape parameters within the uncertainties obtained from the nominal result of this analysis. The standard deviation of the obtained ð1 þ δÞϵ distribution is considered as the systematic uncertainty. The parametrization of the input cross section will be discussed in the next section. The precision of the ISR calculation in the KKMC generator, 0.5%, is taken as an additional uncertainty related to the ISR correction. The uncertainty in the tracking efficiency of leptons is 1.0% per track [32]. Therefore, two leptons in the final state result in an uncertainty of 2.0%. The residual uncertainty from other sources, e.g., lepton separation, trigger efficiency and FSR is small and is conservatively estimated to be 1.0%.
The individual uncertainties include contributions from the branching fractions, the tracking efficiency of pions and photons, the J=ψ mass window, the kinematic fit, the requirement to suppress backgrounds, the fit procedure, and the MC model. The uncertainties of the branching fractions of the ψð3686Þ and other intermediate states are taken from the PDG [26]. The uncertainty in the tracking efficiency for pions is 1.0% per track [32]. A weighted value of 3.5% is assigned to Mode I since the yields ratio between the three-and four-pion cases is about 1, and 2.0% for Mode II with two pions in the final state. The uncertainty in the photon detection efficiency is 1.0% per photon [33]. For Mode II with at least two photons, an alternative measurement of the cross sections with a 2.0% efficiency change is performed to determine the relative differences, which are assigned as the uncertainties from the photon detection efficiency. The uncertainty associated with the J=ψ mass requirement is estimated by smearing the Mðl þ l − Þ distribution of MC samples according to the resolution discrepancy between data and MC, and the resulting uncertainties in detection efficiencies are about 0.1% for Mode I and 0.8% for Mode II. In Mode I, the uncertainty for the kinematic fit is estimated by tuning the helix parameters of charged tracks according to data as in Ref. [34], and the case without the helix parameters correction is taken as the nominal case. In Mode II, the uncertainties associated with the requirements of cos θ π þ π − , M rec ðπ þ π − l þ l − Þ, M corr ðψð3686ÞÞ and Mðγγπ þ π − Þ are evaluated by comparison to the alternative conditions cos θ π þ π − < 0.8,  [16], Belle [14] and BABAR [15]. The open circles represent the results of this work. The full triangles are from previous BESIII work. The open squares are from BABAR work. The full squares are from Belle work. TABLE II. Summary of the systematic uncertainties ð%Þ for different data sets. The sources marked with "Ã" are shared systematic uncertainties for different data sets, and the sources marked with " †" are shared systematic uncertainties between the two ψð3686Þ decay modes. The "Á Á Á" means that the value is very small and can be neglected.   jM rec ðπ þ π − l þ l − Þj > 105 MeV=c 2 , jM corr ðψð3686ÞÞ − Mðψð3686ÞÞj > 10 MeV=c 2 and jMðγγπ þ π − Þ − MðηÞj > 60 MeV=c 2 , respectively. The uncertainties related to the fit procedure are investigated by changing the fit range, replacing the linear function by a quadratic function for the background description and by varying the width of the convolved Gaussian function for the signal shape. The selection efficiency is obtained with the signal MC sample generated according to the PWA results. To estimate the corresponding uncertainty of the MC model, 100 sets of signal MC samples are generated to obtain the detection efficiency distribution, and the resulting standard deviation is taken as the contribution to the systematic uncertainty. In each set, the MC sample is generated by varying all the PWA parameters randomly according to a multivariate Gaussian function, where the mean and width are the nominal value and error of the parameters with correlation considered.
The uncertainties for the combined results of the data samples with large statistics are summarized in Table II. For those data samples with low statistics, the uncertainties are set as the values of the closest data sets in Table II. Assuming all sources of systematic uncertainties to be independent, the total uncertainties are obtained by adding the individual values in quadrature and are found to be in the range of 4.5% to 5.1%.

V. FIT TO THE CROSS SECTION
To study possible Y states in the process e þ e − → π þ π − ψð3686Þ, a binned χ 2 fit is performed to the dressed cross sections according to the method in Ref. [35] by incorporating both statistical and systematical uncertainties and considering both the correlated and uncorrelated terms according to the formula where ΔX is the difference between the measured dressed cross section and the expected value calculated by the probability density function (PDF) for each signal c.m. energy. M is the covariance matrix of elements where the index iðjÞ represents the iðjÞth data set, Δ stat i is the asymmetric statistic uncertainty for the ith data set and the value is chosen following the strategy in Ref. [10], Δ sys i is the total systematic uncertainty, and σ i ðσ j Þ is the measured cross section. ϵ s is the relative correlated systematic uncertainty, which is the quadratic sum of terms marked with "*" in Table II. Because some correlated uncertainties are different for different data sets, to be conservative, the largest value is used for the line shape fit.
The PDF used in the fit is parametrized as a coherent sum of the possible resonant components together with a phase space component for the continuous contribution, i.e., where BW is the Breit-Wigner function with the three-body phase space (PHSP) factor representing the possible resonant structures Yð4220Þ, Yð4390Þ and Yð4660Þ. ϕ k and ϕ cont are the corresponding phases relative to that of the Yð4390Þ. The continuous part is parametrized as In Eq. (5), M k , Γ tot k , Γ ee k and B k are the parameters representing mass, full width, partial width coupling to e þ e − , and the branching fraction to the π þ π − ψð3686Þ final state. Φð ffiffi ffi s p Þ is the standard three-body PHSP factor [26].
In Eq. (6), a and n are free parameters. A nominal fit including Yð4220Þ, Yð4390Þ and Yð4660Þ as well as the continuous components is shown in Fig. 4, and the results are summarized in Table III TABLE III. Results of the fit to the e þ e − → π þ π − ψð3686Þ cross section for the case when the continuous part is described by ψ cont in Eq. (6). The uncertainties involve statistical and systematic ones propagated from the cross section measurement in Table I Table IV. Following Ref. [37], the continuous component is also described with a Fano-like interference function, where g 1 and g 2 are free parameters and the u is the same as that in Eq. (7). The resulting fit quality is χ 2 =d:o:f: ¼ 40.0=20, so this fit method is not accepted.
To examine the significance of the three Y states, the fits are repeated using only two of the three Y states. The significances for the Yð4220Þ, Yð4390Þ, and Yð4660Þ are 4.0σ, 16.1σ, and 8.1σ, respectively. Also, a fit is performed by fixing the resonant parameters of the Yð4220Þ to be those of the ψð4160Þ from the PDG [26] and the resulting fit quality is χ 2 =d:o:f: ¼ 55.5=22. The fit quality indicates that the contribution from the Yð4220Þ cannot be replaced by the ψð4160Þ. It is worth noting that the width of the Yð4220Þ in the nominal results is much smaller than TABLE IV. Results of the fit to the e þ e − → π þ π − ψð3686Þ cross section when the continuous part is described by σ NY in Eq. (7). The uncertainties involve statistical and systematic ones propagated from the cross section measurement in Table I [26], e.g., ψð4160Þ and ψð4415Þ, or floating mass and width are performed. However, none of the extra resonances are observed with a significance larger than 3σ. The systematic uncertainties for the resonant parameters in the cross section fit are discussed below. The uncertainty of the c.m. energy is 0.8 MeV, which is obtained from a measurement of dimuon events [38]. To estimate its effects on the cross section fit, the uncertainty of the c.m. energy is included in the construction of the χ 2 , and the fit is repeated. The resulting changes with respect to the nominal values are taken as the systematic uncertainties. The c.m. energy spread of BEPCII is 1.6 MeV, which is determined by the Beam Energy Measurement System [39]. Thus, the fit is repeated by convolving a Gaussian function with a width of 1.6 MeV to the PDF of the resonant structures, and the resulting differences are taken as the systematic uncertainties. The uncertainties associated with the PDF modeling are considered to be the differences of the results between the alternative fit with ffiffiffiffiffiffiffi ffi σ NY p in Eq. (7) for continuous components and the nominal fit. Table V summarizes the systematic uncertainties for the mass and total width of the resonances Yð4220Þ, Yð4390Þ and Yð4660Þ, where the total uncertainties are the quadratic sums of the individual values.

VI. SUMMARY
Using data with an integrated luminosity of 20.1 fb −1 collected by the BESIII detector operating at the BEPCII collider, the cross section of the process e þ e − → π þ π − ψð3686Þ is measured at c.m. energies between 4.0076 and 4.6984 GeV. The fit results confirm the existence of the Yð4220Þ, Yð4390Þ and Yð4660Þ as well as the contribution of a continuous component. The nominal fit yields a χ 2 =d:o:f: ¼ 31.0=20 and the masses and widths of the Y states are determined to be M ¼ 4234. 4  respectively, where the first uncertainties involve statistical and systematic ones propagated from the cross section measurement and the second are systematic from the line shape fit procedure. This is a supersession of the previous BESIII result [16] with improved precision.
A comparison of the masses and widths of the Yð4220Þ and Yð4390Þ among different processes investigated at the BESIII experiment is shown in Fig. 5(a), where the parameters of the Yð4390Þ are consistent among different decay modes within uncertainties. However, the width of the Yð4220Þ in this analysis is smaller than results from others [7][8][9][10][11][12] and the previous measurement [16]. The differences of the parameters require more study from not only experiment but also theory. The masses and widths of the Yð4360Þ, Yð4390Þ and Yð4660Þ observed in the process e þ e − → π þ π − ψð3686Þ are compared among BESIII, Belle [14] and BABAR [15], as shown in Fig. 5(b). The mass of the Yð4390Þ observed by BESIII is larger than that of the Yð4360Þ observed by Belle and BABAR, while the masses of the Yð4660Þ are in good agreement among the three experiments. Owing to the large scan data sets collected in BESIII, finer structures are observed in the Yð4360Þ energy region. From Fig. 3, the Yð4390Þ observed by BESIII can be considered as the same resonant structure as the Yð4360Þ observed by Belle and BABAR in the same process. The improvement of size and energy sampling on the BESIII data set as well as the interference of finer structures likely account for the mass difference. These results may be further improved in the future by additional energy points and a better understanding of possible intermediate decay modes.  [7][8][9][10][11][12]. (b) Masses versus widths of resonances observed in e þ e − → π þ π − ψð3686Þ from this work and other experiments [14,15].