Charged charmoniumlike structures in the $e^+ e^- \to \psi(3686) \pi^+ \pi^-$ process based on the ISPE mechanism

In 2017, the BESIII Collaboration announced the observation of a charged charmonium-like structure in the $\psi(3686)\pi^\pm$ invariant mass spectrum of the $e^+ e^- \to \psi(3686) \pi^+ \pi^-$ process at different energy points, which enables us to perform a precise study of this process based on the initial single pion emission (ISPE) mechanism. In this work, we perform a combined fit to the experimental data of the cross section of $e^+ e^- \to \psi(3686) \pi^+ \pi^-$, and the corresponding $\psi(3686)\pi^\pm$ and dipion invariant mass spectra. Our result shows that the observed charged charmonium-like structure in $e^+ e^- \to \psi(3686) \pi^+ \pi^-$ can be well reproduced based on the ISPE mechanism, and that the corresponding dipion invariant mass spectrum and cross section can be depicted with the same parameters. In fact, it provides strong evidence that the ISPE mechanism can be an underlying mechanism resulting in such novel a phenomenon.


I. INTRODUCTION
Since 2003, study on the charmoniumlike XYZ states has become a hot spot of hadron physics and particle physics. Since it has the close relation to non-perturbative behavior of quantum chromodynamics (QCD), the relevant investigation of the charmoniumlike XYZ states is helpful to deepen our understanding of strong interaction. In the past sixteen years, there have been extensive discussions about this issue (see reviews [1][2][3][4] for more details).
To solve the puzzling phenomena existing in the hiddenbottom dipion decays of Υ(10860), the initial single pion emission (ISPE) mechanism was proposed in Refs. [5][6][7], where the charged Z b (10610) and Z b (10650) structures can be qualitatively reproduced. In 2012, considering the similarity between the bottomonium and charmonium families, the ISPE mechanism was applied to study the hidden-charm dipion decays of higher charmonia and charmoniumlike state Y(4260), where the predicted charged charmoniumlike structures near DD * or D * D * threshold may exist in the corresponding J/ψπ + , ψ(3686)π + , and h c π + invariant mass spectra [8]. These predictions given by the ISPE mechanism have provided valuable information to search for such charged charmoniumlike structures.
In 2013, the BESIII and Belle collaborations indeed discovered a charged charmonium-like Z c (3900) in the process e + e − → J/ψπ + π − around √ s = 4.26 GeV [9,10], which was quickly and further confirmed by CLEO-c in the same process but at √ s = 4.17 GeV [11]. Later, another charged charmo-niumlike structure named Z c (4020) was observed in the h c π + invariant mass spectrum of the e + e − → π + π − h c process [12]. Additionally, at the center-of-mass energy of 4.26 GeV, the BESIII Collaboration reported two charged charmoniumlike structures, Z c (3885) in the (D * D ) ± invariant mass spectrum of e + e − → π ∓ (D * D ) ± process [13] and Z c (4025) in the (D * D * ) ± invariant mass spectrum of the e + e − → π ∓ (D * D * ) ± process [14]. The story of finding charged charmoniumlike structures is continued. In 2017, the BESIII experiment made a progress on observing charged charmoniumlike structure in the ψ(3686)π + invariant mass spectrum by analyzing the data of cross section of the e + e − → ψ(3686)π + π − process at different energy points √ s = 4.226, 4.258, 4.358, 4.387, 4.416, and 4.600 GeV [15]. In fact, this observation deepens our confidence again. As shown in Fig. 1, the predicted charged charmoniumlike structures in the J/ψπ + , ψ(3686)π + , and h c π + invariant mass spectrum have been reported in experiments until now.
The BESIII's experimental measurement of the e + e − → ψ(3686)π + π − process forces us to carry out a further study of this process based on the ISPE mechanism. Before the present work, in Ref. [16], the similar idea was once adopted to study the e + e − → J/ψπ + π − process combined with the BE-SIII data of the J/ψπ ± invariant mass spectrum and the corresponding dipion invariant mass spectrum of e + e − → J/ψπ + π − at one energy point √ s = 4.26 GeV, which shows that the Z c (3900) structure can be reproduced well by the ISPE mechanism. Different from Ref. [16], the present work will perform a combined fit to the released data of the ψ(3686)π ± and π + π − invariant mass spectrum from the e + e − → ψ(3686)π + π − process at different energy points. Comparing with the former work [16], the present work is obviously a tough task and is full of challenge. In this work, we will reproduce the structures observed in the e + e − → ψ(3686)π + π − process step by step based on the ISPE mechanism. It also reveals the relation between the observed charged charmoniumlike structures and the ISPE mechanism, which will be helpful to deepen our understanding on such a novel phenomenon.
This work is organized as follows. After introduction, we perform an analysis of the mechanisms working in the e + e − → ψ(3686)π + π − , and then fit the present model with the experimental data. In Sec. III, we present our fit results to the cross sections for e + e − → ψ(3686)π + π − and the mass spectra of ψ(3686)π ± and dipion invariant mass spectra. Sec. IV is devoted to conclusions and discussion.
II. MECHANISMS OCCURRING IN e + e − → ψ(3686)π + π − In Ref. [16], we once reproduced the Z c (3900) structure in the e + e − → J/ψπ + π − process based on the ISPE mechanism. In the present work, we will apply the same idea to discuss e + e − → ψ(3686)π + π − . There exist serval different decay mechanisms occurring in the process e + e − → ψ(3686)π + π − shown in Fig. 2. Fig. 2 (1) represents the nonresonance contribution, where the final states ψ(3686)π + π − directly couple to the photon produced by the electron-positron annihilation. Fig. 2 (2)-(4) correspond to the processes with an intermediate charmonium resonance, where the virtual photon first couples to a vector charmonium ψ R , which decays into ψ(3686)π + π − . There are three different kinds of decay mechanisms in the higher charmonium dipion decays: i) As shown in Fig. 2 (2), the first is the one in which the dipion is produced by the gluon hadronization with gluons emitted by a charm quark. ii) In the second decay mechanism, the dipion is produced through a scalar meson decay, while the scalar meson and ψ(3686) are connected to the initial ψ R by a charmed meson loop. This is presented in Fig. 2 (3). iii) The last decay mechanism is the ISPE mechanism, which is presented in Fig. 2 (4). We will list resonaces R inserted into Fig. 2 later in Table I. In the following, we will present the amplitudes corresponding to different mechanisms. For Fig. 2 (1), its amplitude is expressed as where g NoR is a coupling constant, e = √ 4πα denotes an electric charge with α = 1/137 being the fine structure constant, and F NoR (s) is a form factor. In the present work, we adopt the following form factor [16] F NoR (s) = e −a NoR ( where a NoR is a phenomenological parameter, which should be fixed when fitting the data. Σm f = m ψ + 2m π ± is the sum of the masses of the involved particles appearing in the final states. For Figs. 2 (2)-(4), a general amplitude, with a help of the vector meson dominance (VMD) ansatz [17,18], can be written as where f R is a decay constant of an intermediate R and p 0 denotes the momentum of the intermediate R, which satisfies p 2 0 = (k 1 + k 2 ) 2 = s. A µν indicates the decay amplitude of ψ R → ψ(3686)π + π − . Here, we also introduce the form factor F R to balance the otherwise over-increased decay width with an increased phase space, which isgiven by, where a R is a parameter obtained when fitting the data.
For Fig. 2 (2), the dipion is produced directly by the gluon emission. Hereafter, such a process is named as a direct production process and the decay amplitude is parameterized as [19], where the terms in the first and second square brackets indicate that the orbital momenta of the dipion are S -wave and D-wave, respectively. ∆M = m R − m ψ is the mass splitting of the initial ψ R and final ψ(3686). F Dir and κ Dir are parameters necessary for fitting and α is the angle between ψ R and π − in the π + π − rest frame. As shown in Fig. 2 (3), the dipion is produced through a scalar meson, and the scalar meson as well as ψ(3686) are connected to ψ R via a charmed meson loop. Such a kind of a meson loop has been proven to be important in the dipion transitions between heavy quarkonia [20,21]. Due to the kinetic limit of the involved process, we take the σ meson into consideration. Moreover, in the present work, we mainly aim at the ISPE mechanism. Thus, we will not directly estimate the triangle diagrams but parameterize their contributions in the form [16], where f σ and g σ are the S and D wave coupling constants between ψ R and ψ(3686)S, respectively, and ϕ σ is the phase angle between the S and D wave couplings. We should notice that the mass and width of the σ meson are comparable, and thus, we adopt a momentum dependent width, which is, where P(m ππ ) = m 2 ππ /4 − m 2 π is the pion momentum in the dipion rest frame, while | P(m σ )| is the pion momentum with the on-shell σ meson.
As for the ISPE mechanism, we directly estimate the triangle diagrams. To fit with the experimental data, we reduce the amplitudes of the ISPE mechanism in the form, where the coefficients C i j , {(i, j) = (0, 1, 2, 3)} are obtained by the loop integrals of the amplitudes in Appendix A.
With above preparations, we get the total amplitude, which could be the coherent sum of the amplitudes corresponding to different mechanisms, and reads as with φ (i) R being the phase angles in front of different amplitudes for each resonance. For convenience, we use φ Dir , φ σ and φ ISPE to refer to these phase angles in the amplitudes of direct production process, σ production process and ISPE processes, respectively. The differential cross section is [22] where M Tot 2 indicates an average over the spin of initial electron and positron and a sum over the spin of the final state ψ(3686). dΦ 3 is the phase space factor, which is expressed as where m 23 corresponds to m ψ π + , m 12 is m π + π − , and the definitions of angles θ and η are shown in Fig. 3. . The definitions of angles θ and η that appear in Eq. (11). In this figure, we also give the relative positions of k 1 , k 2 , and p i (i = 1, 2, 3). Since p 0 = k 1 + k 2 = ( √ s, 0, 0, 0), p 0 lies on the origin of the xyz frame.

III. NUMERICAL RESULTS AND DISCUSSION
In the differential cross section in Eq. (10), there are two free parameters, g NoR and a NoR , in the nonresonance amplitude. As for the intermediate charmonium contributions, there are 12 additional free parameters for each resonance in Table I, which are a R : a parameter in the resonance form factor, φ Dir : a phase angle of the direct production amplitude, F Dir , κ Dir : parameters in the direct production amplitude, φ σ : a phase angle of the σ meson production amplitude, f σ , g σ : S and D wave coupling constants in the σ meson production amplitude, ϕ σ : a relative phase angle between S and D wave terms in the σ meson production amplitude, φ ISPE : a phase angle of the ISPE amplitude, g RD ( * )D( * ) π : coupling constants of the four-particle interactions in the ISPE mechanism.
In the following, we need to fit our results with the experimental data of e + e − → ψ(3686)π + π − , which include the measured total cross section, the ψ(3686)π ± , and π + π − invariant mass spectra released in Ref. [15]. For the e + e − → ψ(3686)π + π − process, there exist contributions from intermediate vector charmonia as shown in Fig. 2 (2)-(4). Thus, it is a key point how to select intermediate vector charmonia. In the present work, we adopt two scenarios, which will be addressed in Sec. III A and Sec. III B.
In the present section, we first try to fit the experimental data by considering ψ(4160), ψ(4220), and ψ(4415) as intermediate state contributions to e + e − → π + π − ψ(3686). Their resonance parameters are collected in Table I.
After integrating over m π ± ψ(3686) or m π + π − in Eq. (10) at each fixed √ s, one gets the differential cross sections of e + e − → π + π − ψ(3686) for each invariant mass squared distribution. The cross section depending on √ s can be obtained by integrating over both m ψ(3686)π ± and m π + π − in Eq. (10). Here, we perform a combined fit of the invariant mass squared dependent cross sections to the experimental data of the differential cross sections at a fixed √ s. In the three-resonance scenario, there are 38 parameters, which should be determined by the present fit. Our fitted values for all the parameters are presented in Table II. With these parameters, χ 2 /d.o.f is estimated to be 3.704. The fitted differential cross sections for e + e − → π + π − ψ(3686) at a fixed √ s in three-resonance scenario are presented in Fig. 4 (black dashed curves), where √ s is fixed at 4.226, 4.258, 4.358, 4.387, 4.416, and 4.600 GeV, respectively. Most of the differential cross sections can be reproduced. In particular, both π ± ψ(3686) and π + π − invariant mass spectra are well reproduced for √ s = 4.226 GeV. As for √ s = 4.258 GeV, the π ± ψ(3686) invariant mass spectrum is a bit different from the experimental data. In the vicinity of m 2 π ± ψ(3686) = 15 GeV 2 , the experimental data show a jump, but such a feature is not reproduced in three-charmonium scenario. On the other hand, the fitted curve of m π + π − distribution is slightly larger than the experimental data in the low and high m π + π − energy range. Similar to the case for the fit to the data at 4.358 GeV, the fitted curve of the m 2 π + ψ(3686) distribution is a bit lower on the whole than the experimental data, while for m 2 π + π − distribution, the fitted curve is a little bit larger than the experimental data in the lower m π + π − energy range. For •    Our combined fit for the ψ(2S )π ± and π + π − invariant mass spectra lying on 6 different center-of-mass energies which are reported by BESIII in Ref. [15]. Here, the black dashed and red solid curves are the fitted results under three-charmonium and four-charmonium scenarios, respectively.

Events
the cases of √ s = 4.387 GeV and √ s = 4.416 GeV, the fitted curves for both m 2 π ± ψ(3686) and m 2 π + π − well reproduce the experimental data. As for the data at 4.600 GeV, the fitted curve obtained under three-charmonium scenario is lower than the experimental data for the m 2 π ± ψ(3686) invariant mass distribution, while for m 2 π + π − distribution, since the errors of experimental data are very large, thus, the fitted curve is roughly consistent with the data. The fitted cross section for e + e − → π + π − ψ(3686) is presented in Fig. 5 (black dashed curve). The structure in the vicinity of 4.2 GeV and the broad enhancement near 4.4 GeV are well reproduced. However, one can find that the fitted curve is above the experimental data in the vicinity of 4.3 GeV. Such a phenomenon indicates that there may exist other intermediate charmonium contributions to the e + e − → π + π − ψ(3686) process, which inspires us to introduce a fourcharmonium scenario to study the e + e − → π + π − ψ(3686) process, which will be discussed in the next subsection. In Fig. 6, we present the individual contributions from the nonresonance background and three charmonia. each of which is the coherent sum of direct contributions, the scalar productions, and ISPE processes for each resonance. One can find the contribution from ψ(4160) could reach up to nearly 100 pb, while the contributions from ψ(4220) and ψ(4415) are about 70 pb. From Fig. 6, one can find that resonances are still in a Breit-Wigner form. Thus, we simulate each lineshape in Fig. 6 with a Breit-Wigner resonance, which is in the form 12πΓ e + e − B(R → π + π − ψ(3686)) With the dilepton decay widths listed in Table I, one can roughly estimate the center values of the branching ratios for ψ R → π + π − ψ(3686), where ψ R = (ψ(4160), ψ(4220), ψ(4415)). As shown in Table III, one can find the center values of branching ratios are of the order of ∼ 10 −3 for ψ(4220) and ψ(4415) and ∼ 10 −2 for ψ(4160).  B. Four-charmonium scenario to e + e − → ψ(3686)π + π − As discussed in the three-charmonium scenario in Fig. 5, the experimental data in the vicinity of 4.3 GeV cannot be well reproduced. In Ref. [28], the mass spectra and decay properties of higher charmonia have been investigated, and a charmonium ψ(4380) as a partner of ψ(4220) was predicted, which has the resonance parameters m = 4384 MeV and Γ = 84 MeV. This progress stimulates us to introduce the fourcharmonium scenario to study e + e − → ψ(3686)π + π − , namely, four charmonia ψ(4160), ψ(4415), ψ(4220), and ψ(4380) are considered when fitting the data. All the resonance parameters are listed in Table I. Indeed, we find that the details of the cross section of e + e − → ψ(3686)π + π − around √ s = 4.3 GeV can be depicted.
The fitted parameters are presented in Table IV. With the center values of these parameters, one can get the differential cross sections for each invariant mass squared distribution and the cross section depending on √ s, which are presented in Fig. 4 (red curves) and Fig. 5 (red curve), respectively. After including ψ(4380), one finds that almost all the differential cross sections can be quantitatively fitted better than the three resonance-scenario. Especially, the jump near at 15 GeV 2 in m 2 π ± ψ(3686) invariant mass spectrum at √ s = 4.258 GeV can be well reproduced and m 2 π + π − invariant mass spectrum at this energy point is well fitted. The fits to differential cross sections at √ s = 4.358 GeV are also much better than the one in the three-resonance scenario. As for the √ s dependent cross section, the four-resonance scenario can well reproduce experimental data, especially the data near at 4.3 GeV. Quantitatively, the χ 2 /d.o.f is reduced to 2.77 in the four-charmonium scenario. Similar to the case of the threecharmonium scenario, we present the individual contributions of the non-resonance background and resonances in Fig. 7. One can find that the resonance contributions are smaller than the corresponding one in the three-charmonium scenario. The fitted branching ratios of ψ R → π + π − ψ(3686) are presented in Table III, where the central values of these branching ratios for ψ(4160), ψ(4220), and ψ(4415) are of the order of ∼ 10 −3 , while the one for ψ(4380) is of the order of ∼ 10 −2 .
However, before definitely identifying these charged charmomium-like structures as exotic tetraquark states, we need to exhaust all the possibilities of explaining them in the conventional theoretical framework. Here, we must mention our previous prediction of charged charmonium-like structures in our theoretical work [8], where these predicted charged charmonium-like structures by the ISPE mechanism may exist in the J/ψπ ± , h c π ± , and ψ(3686)π ± invariant mass spectra of higher charmonia, and in the charmonium-like state Y(4260) decays into J/ψπ + π − , h c π + π − and ψ(3686)π + π − [8]. After two years of our paper, the BESIII and Belle collaborations discovered a charged charmonium-like Z c (3900) in the process e + e − → J/ψπ + π − around √ s = 4.26 GeV [9,10], and CLEO-c confirmed it in the same process but at √ s = 4.17 GeV [11]. In 2013, BESIII observed another charged charmoniumlike structure named as Z c (4020) in the h c π + invariant mass spectrum of the e + e − → π + π − h c process [12]. These experimental observations make us believe that the ISPE mechanism indeed plays an important role in producing these novel phenomena. What is more important is that these measured experimental data make us possible to study them with higher precision. In Ref. [16], we fitted the experimental result of the charged Z c (3900) observed in the J/ψπ ± invariant mass spectrum [9,10] and indicated that the charged Z c (3900) structure can be well established based on the ISPE mechanism. This study further enforces theorist's confidence to explain these experimental results in the conventional theoretical framework. Similar studies were proposed in Refs. [56][57][58][59][60]. Later, the Lattice QCD calculation [61,62] also supports such a scenario.
In 2017, the BESIII Collaboration again observed charged charmonium-like structures in the ψ(3686)π ± invariant mass spectrum of the e + e − → ψ(3686)π + π − process at different energy points √ s = 4.226, 4.258, 4.358, 4.387, 4.416, 4.600 GeV [15]. This new observation inspires our interest in further testing the ISPE mechanism combined with the present precise experimental data. A crucial task of this work is to examine whether the cross section of the e + e − → ψ(3686)π + π − process, and the corresponding ψ(3686)π ± and π + π − invari-ant mass spectra at different √ s values can be reproduced by the ISPE mechanism. As illustrated by our calculation, we find that the reported charged charmonium-like structures in the ψ(3686)π ± invariant mass spectrum of the e + e − → ψ(3686)π + π − process at different energy points can be depicted in a unified framework, which sheds light on the properties of the observed charged charmonium-like structures. We need to emphasize that our study also further supports the existence of a new higher charmonium ψ(4380) which was predicted in Ref. [28], since the fitting result under the four-charmonium scenario is better than that under the three-charmonium scenario. We strongly suggest the BE-SIII and BelleII collaborations to search for this predicted missing charmonium ψ(4380), especially by analyzing the e + e − → ψ(3686)π + π − process in near future.
In the following years, studies on charmonium-like XYZ states will still be an interesting research topic at the typical facilities like BESIII, LHCb, and BelleII. Since the observed charmonium-like XYZ states can be a good candidate of an exotic tetraquark matter, it is a main task for both experimentalists and theorists how to identify them as a genuine exotic multiquark state. To achieve this aim, we need to check whether the charmonium-like XYZ states can be assigned into a conventional charmonium family or can be settled down in a conventional theoretical framework. The present work is the effort along this line. Besides theoretical investigation, the Lattice QCD study about charmonium-like XYZ will provide valuable information to shed light on the underlying properties of the XYZ states. As emphasized in a recent review article [3], a joint effort from phenomenological method, experimental analysis, and Lattice QCD calculation should be paid more attentions on, which must promote our knowledge of how to form charmonium-like XYZ states. sider the diagrams caused by isospin symmetry and initial π + emission. For isospin symmetry, the diagrams can be obtained by the replacement {D ( * )+ , D ( * )0 ,D ( * )0 } → {D ( * )0 , D ( * )− , D ( * )+ }, which are the same contributions as for a π − emission. For contributions of an initial π + emission, the typical diagrams of the ISPE mechanism can be obtained by applying charge conjugation on each Fig. 8 (i) (i=a∼f). Thus, the total decay amplitude can be expressed as where the factor 2 reflects the isospin symmetry mentioned above, and A (i) µζ[p 2 p 3 ] means that the contributions of an initial π + emission can be obtained by applying p 2 p 3 transformation on A (i) µζ . To write the decay amplitudes A (i) µζ in Eq. (A1), the effective Lagrangian approach is adopted. For the four particles interaction, the effective Lagrangian is [5,63] L ψD ( * ) D ( * ) π = −ig RDDπ ε µναβ ψ µ ∂ ν D∂ α π∂ βD + g RD * Dπ ψ µ (DπD * µ + D * µ πD) −ig RD * D * π ε µναβ ψ µ D * ν ∂ α πD * β − ih RD * D * π ε µναβ ∂ µ ψ ν D * α πD * β .
We want to note here that, to reduce the number of free parameters in our fit, we assume h RD * D * π = g RD * D * π in this work.
With above effective Lagrangians, the expressions for decay amplitudes A (i) µζ can be written as ×[gD * D * ψ g ζξ q κ − g ζκ k 2ξ − g κξ (q ζ − k 2ζ ) ] ×[gD * D * ψ g ζξ q κ − g ζκ k 2ξ − g κξ (q ζ − k 2ζ ) ] , In A (i) µζ , a monopole form factor F (q 2 ) = q 2 −Λ 2 m 2 E −Λ 2 is introduced to reflect the structure effect at each vertex and the offshell effect of the exchanged D ( * ) meson. In this form factor, m E is the mass of the exchanged D ( * ) meson, Λ is parameterized as Λ = m E + α Λ Λ QCD with Λ QCD = 220 MeV. Since Ref. [8] has shown that the ISPE mechanism is weakly dependent on α Λ , thus in this work we also set α Λ = 1 to present our results.