Simultaneous description of the $e^+e^- \to J/ \psi \, \pi \pi\, (K \bar{K})$ processes

In this work, we provide a simultaneous and accurate description of the $\pi^+\pi^-$ and $\pi^{\pm} J/\psi$ invariant mass distributions of the recent BESIII data on $e^+ e^- \to J/\psi \; \pi^+ \pi^-$ together with the $e^+ e^- \to J/\psi \; K^+K^-$ cross sections at $e^+e^-$ center-of-mass energies $q=4.23$ GeV and $q=4.26$ GeV. The rescattering effects between pions in the S and D-waves are taken into account through the Muskhelishvili-Omn\`es formalism. Since the physical region of the $\pi\pi$ invariant mass extends above 1 GeV, the important $K\bar{K}$ intermediate state in the S-wave is implemented through coupled-channel unitarity. For the left-hand cuts, we account for the well established charged exotic state $Z_c(3900)$ in $t$- and $u$-channels, while the other contributions are absorbed in the subtraction constants. For the $e^+ e^- \to J/\psi\, K \bar{K}$ we provide the prediction of the two-kaon invariant mass distribution. The constructed amplitudes serve as an essential framework to interpret the present and forthcoming measurements by the BESIII and Belle II Collaborations.


I. INTRODUCTION
The charged exotic charmonium-like state Z c (3900) was discovered simultaneously by the BESIII and Belle Collaborations in 2013 both in direct production [1] and using initial-state radiation [2], in the process e + e − → J/ψπ + π − and soon confirmed using the CLEO-c data [3]. In 2015, the neutral partner was observed by the BESIII Collaboration in the same reaction with neutral pions e + e − → J/ψπ 0 π 0 [4]. Recently, the D0 Collaboration, using proton-antiproton collisions, has found a signal of Z c (3900) in non-prompt semi-inclusive weak decays of b-flavored hadrons [5,6]. Furthermore, in recent years, BESIII has observed Z c (3900) in the e + e − → (DD * ) ∓ π ± process using a single-tag analysis [7], a double-tag analysis [8], and also by analyzing the neutral channel e + e − → (DD * ) 0 π 0 [9]. The most precise data so far has been reported in [10], where an updated BESIII analysis of e + e − → J/ψπ + π − allowed to determine the spin-parity J P = 1 + assignment of the Z c (3900).
From the theory side, the nature of Z c (3900) is still a puzzle [11][12][13][14]. Most likely it corresponds to a pole in the unphysical Riemann sheet, which could be a hadrocharmonium [15,16], molecular state [13,17] or a virtual state [18,19]. The peak at the Z c (3900) position has also been interpreted through a kinematic effect [20][21][22]. The most popular scenarios correspond to the triangle singularity associated with D * 0 (2300)D * D [21] or D 1 (2420)D * D [22] loops. In both cases, the left-hand cut branch point stays relatively far away from the physical region, either due to the large width of D * 0 (2300), or due to the off-shellness of D 1 (2420) for q = 4.23 GeV and q = 4.26 GeV, and only theD * D threshold cusp gets enhanced. However, as it was pointed out in [23], the recent BESIII data [10] indicate that the Z c (3900) peak is more enhanced for q = 4.23 GeV compared to q = 4.26 GeV, in contrast to what one expects from the threshold cusp enhancement mechanism due to the triangle singularity. Additionally, the contribution from the rescattering process has to be accounted for, which typically smooths out kinematic singularities. To shed further light on this puzzle, it will be very helpful to observe the Z c (3900) in other decay modes [24]. Besides, it is important to clarify if there exists a possible strange partner of Z c (3900), the so-called Z cs , which can show up in the KJ/ψ distribution of the e + e − → J/ψK + K − process. So far, Belle [25] and BESIII [26] have not seen a clear structure in the KJ/ψ mass distribution, and future high statistics measurements are necessary.
The purpose of the present work is to demonstrate a dispersive amplitude analysis, which can be applied in the experimental works to describe the whole Dalitz plot with minimum assumptions about the nature of the charged Z c state. Our work is a continuation of the previous work [27], where for the first time, a dispersive amplitude analysis was applied to describe e + e − → ψ(2S)π + π − Dalitz plot projections [28,29]. In our current analysis, the recent BESIII [10] data on e + e − → J/ψπ + π − play the central role. We present a simultaneous description of the π + π − and π ± J/ψ invariant mass distributions by providing rigorous dispersive treatment of the ππ final state interactions. We account for Z c (3900) as an explicit degree of freedom in the t-and u-channels and unitarize the ππ final state interaction on the base of the Muskhelishvili-Omnès formalism. Other possible left-hand cut contributions are absorbed in the subtraction constants which we determine from a combined fit to the e + e − → J/ψπ + π − Dalitz plot data and the total cross-section data for e + e − → J/ψK + K − . Due to the relatively large physical region of the ππ invariant mass, we also extend our previous analysis of [27] to the coupled-channel in the ππ S-wave and include the Dwave. Allowing for a minimum number of parameters, which enter in the form of subtraction constants, and assuming the absence of Z cs at q = 4.23 GeV and q = 4.26 GeV 1 , we provide a prediction for the invariant mass dis-tribution for the process e + e − → J/ψ K + K − .
In our analysis we do not aim at a description of the full e + e − → J/ψπ + π − cross-section and instead implement the q 2 dependence model independently, by applying our formalism for each q-value independently. The study of the two possible resonance structures seen in the e + e − → J/ψπ + π − total cross-section [31] is beyond the scope of this paper. Rather, we want to use the available Dalitz plot projection data to make a simultaneous description of both π + π − and π ± J/ψ invariant mass distributions and obtain a prediction of the K + K − and K ± J/ψ invariant mass distributions. This is different from the analysis performed in Ref. [32], which focused only on the π + π − invariant mass distribution to get insights into the structure of the Y (4260) state from the light-quark perspective. Though the analysis of the ππ final state interaction is similar in spirit to ours, there are several technical differences, which we will point out below.

II. KINEMATICS
The double differential cross section for the e − (p 1 ) e + (p 2 ) → γ * (p γ * ) → J/ψ(p ψ ) π + (p π + ) π − (p π − ) process can be written as where we have neglected the electron mass compared to the e + e − center of mass (CM) energy q = p 2 γ * . In Eq.(1) the helicity amplitudes H λ1λ2 are defined in the usual way, where λ 1 (λ 2 ) denote the γ * (J/ψ) helicities, respectively. For the process γ * → J/ψ π + π − the following Mandelstam variables are chosen, the strange partner of Zc(3900) would have a heavier mass (in particularly, Ref. [30] predicts a mass of 3.97 ± 0.08 GeV) and therefore Zcs cannot be seen as peak in the KK invariant mass distribution for q = 4.23 − 4.26 GeV.
which satisfy In the following we use the kinematics in the CM frame of the di-pion system, and define z ≡ cos θ s as the cosine of the angle between the p π + and the p ψ momenta. Thus, in this frame the following relations hold where and λ being the Källen function. Consequently, z can be written in terms of t and u as

III. DISPERSIVE FORMALISM
In this section, we briefly describe the dispersive formalism that we adopt to account for the rescattering between two pions (kaons), which generates the most important singularities at low energies in the s-channel. The partial wave (p.w.) expansion reads where I is the isospin, Λ = λ 1 − λ 2 and d (J) Λ,0 is the Wigner rotation function. For better readability, below we will consistently suppress the isospin indices, and retrieve them at the beginning of Sec.IV. On account of causality, the p.w. amplitudes should have contributions from the left-and right-hand cuts, where the branch cut due to the two-pion interaction starts at s = 4 m 2 π . We note that the amplitudes h λ1λ2 (s) are subject to kinematical constraints, which in principle have to be removed before application of dispersion relations. The hadron tensor H µν of γ * → J/ψ ππ can be decomposed into a suitable set of Lorentz structures given in [27] (see also [33][34][35][36][37]), with F i the corresponding invariant amplitudes. One can then show that for the S-wave the p.w. helicity amplitudes are correlated at the kinematic points while for the D-wave the kinematic correlations between different p.w. helicity amplitudes are more complicated and can be found in [37]. As it will be shown in the next section, for the considered kinematics most of these constraints have a negligible impact on the results, since the sum in Eq.(1) in the physical region can be written in terms of H ++ only, i.e.
Under this approximation it is enough to take into account only the so-called centrifugal barrier factor for , which comes from the properties of the Legendre polynomials entering p.w. expansion in Eq. (9). We note, however, while Eq. (14) is exact for s = 4m 2 π , a zero at s = (q − m ψ ) 2 is only approximate and typically a few MeV away. This is related to the approximation made in Eq. (13), which we will discuss further on.
The discontinuity across the branch cut in the schannel is given by which can be straightforwardly extended to the case of two cuts (coupled-channel case) in the S-wave The two-body phase space ρ(s) is given by where σ αα (s) = λ 1/2 (s, m 2 α , m 2 α )/s, with α = π or K. The {ππ, KK} coupled-channel scattering amplitude t(s) is normalized as . We note, that in the p.w. expansion of the γ * (q) → J/ψ KK process we include an extra factor 1/ √ 2 in contrast to γ * (q) → J/ψ ππ in order to match our normalization for the hadronic p.w. amplitudes, which ensure the same unitarity relations for the identical and non-identical particles. For the Swave the standard Muskhelishvili-Omnès representation for the left-hand cut subtracted p.w. amplitude is given by (modulo subtractions) where the coupled-channel Omnès function (with 1 = ππ and 2 = KK) satisfies the following unitarity relation Disc Since the tail of the f 2 (1270) resonance could overlap with the physical region, we include D-wave singlechannel ππ-rescattering. As discussed previously, we factor out the known threshold factor and write a dispersion relation for h where under the dispersive integral we slightly adjusted a zero of γ(s ) at s = (q − m ψ ) 2 to match exactly a zero of h (2),L ++ (s ), which is few MeV away. One can notice, that the overall threshold factor γ(s) is also needed to compensate the singularities of z = cos θ s (see Eq. (8)) of the full amplitude H R ++ (s, t) at the borders of the Dalitz plot (i.e. at s = 4m 2 π and s = (q − m ψ ) 2 ). This is different from [32] where in the dispersive representation no threshold factors were taken into account in the Dwave.
In our formalism, we are accounting for the ππ rescattering effects only in S-and D-waves, and beyond that (for J > 2) the p.w. amplitudes in Eq.(10) are approximated by the first term, h (J),L λ1λ2 (s). In other words, we keep the cross channel p.w. expansion to all orders. That is crucial to get the description of the full Dalitz plot, where there are peaks structures in both ππ and πJ/ψ systems. The final result for the total helicity amplitude can be written as 2 where the sum goes only over even J values due to Bose symmetry of two pions and C-parity conservation. 2 We note the difference between Eq. (22) and the reconstruction theorem written in [27]. The latter is correct only for the scalar particles and needs to be modified for the particles with spin [38].
Since we only considered rescattering effects in the s-channel, and h (0),t (22), this has no effect on the results in [27].

A. Left-hand cuts
The cuts associated with the crossed channel exchange terms, i.e. h (J),L λ1λ2 (s), are approximated by the charged Z c exchanges, motivated by the experimental data [10], where the Z c (3900) axial-vector state and its kinematic reflection show up as clear peaks in the πJ/ψ projection for both e + e − -CM energies q = 4.23 GeV and q = 4.26 GeV. According to the mechanism γ * (q 2 ) → π ∓ +(Z ± c → J/ψ + π ± ), the helicity amplitude can be expressed in a general form as follows where S νµ (Q z ) is the axial meson propagator. We use the following vertices [39], where Q z = (p γ * − p π ), C Zcψπ the coupling among Z c , J/ψ and π and F γ * πZc (q 2 ) is the corresponding transition form factor. For the present analysis, the latter should in principle encode for two resonances, as observed in the data [31]. In our formalism we will perform two independent analyses at q = 4.23 GeV and q = 4.26 GeV, without any assumption for F γ * πZc (q 2 ) to avoid possible model dependence.
By inspecting Eq.(23) for our particular kinematics, we observe that the helicity amplitudes, H Zc ++ and H Zc 00 give the main contribution compared to other helicity amplitudes. Furthermore, H Zc ++ and H Zc 00 turn out to be numerically very close to each other |H Zc ++ | ≈ |H Zc 00 |. Therefore, one can write the sum in Eq.(1) in terms of only H Zc ++ (see Eq. (13)) and this approximation has less than 1% error in the physical region. A similar observation was also made in Refs. [40,41] based on the heavy-quark nonrelativisitic expansion.
The expression of the helicity amplitude H Zc ++ in terms of the invariant amplitudes F Zc i (s, t) is given by where Due to the polynomial ambiguity of the p.w. amplitudes, we will consider only the pole contribution. Based on the fixed-s Mandelstam representation one can show that the pole contribution corresponds to fixing t = m 2 Z and u = m 2 Z in the numerators of Eq. (26). This procedure is in line with the definition of the on-shell transition form factor F γ * πZc (q 2 ) and does not change the amplitude in the physical region.

B. Triangle Singularities
For the three-body decays, it is frequent that the lefthand cut overlaps with the right-hand cut and requires special treatment in the dispersive formalism. In our previous analysis of γ * (q 2 ) → ψ(2S)ππ [27] such an overlap required a distortion of the integration path which was performed by including an additional anomalous piece [42][43][44]. For the processes considered in the present paper, the overlap of the left and right-hand cuts does not introduce anomalous thresholds, but still require the proper analytical continuation for the energy variable q 2 → q 2 + i [45,46] due to the presence of the so-called triangle singularity [11] associated with Z c ππ loop. Indeed, for q = 4.23 GeV and q = 4.26 GeV the exchange Z c (3900) state in the triangle loop can be on-shell, satisfying the Coleman-Norton conditions q 2 > (m Z + m π ) 2 and m 2 Z > (m ψ +m π ) 2 [47]. This implies that the branch point s − associated with the left-hand cut is located just above the two pion threshold but infinitesimally below the real axis [22]. We note, that the q 2 → q 2 + i continuation guarantees that the branching point never crosses the unitarity cut and the dispersive representations of Eqs. (18) and (21) are correct. Due to the finite resonance width, however, the effect of the triangle singularity smears out, since the singular point is shifted further away from the physical region.
There are several ways of accounting for the width of the Z c (3900) state. The proper implementation requires modeling the propagator using a spectral representation [46], i.e., it should have sound analyticity properties, such as pole on the unphysical Riemann sheet and the righthand cuts starting at πJ/ψ and DD * thresholds. This analysis is beyond the scope of our paper due to the lack of experimental information. Since the width of the Z c (3900) meson is relatively small (Γ Z = 28.2 MeV) [48], we follow here a pragmatic approach by implementing the finite width in the denominators of Eq. (26). In this case, it is possible to cross-check our dispersive implementation on an example of a toy model of scalar fields with a constant (equal to unity) interaction between pions. As one can see in Fig.1, the result of the dispersive calculation and the calculation via Feynman parameters in perturbation theory give the same results. For illustrative purpose, we also show in Fig.1 the result based on using the spectral representation of the Z c propagator [46], but accounting for just one channel πJ/ψ as it was done in [32]. As expected the difference is negligible. Due to the narrowness of the Z c (3900) state one can also observe in Fig.1 that the peak is still relatively sharp. However, the inclusion of the unitarization through the Muskhelishvili-Omnès representation smears it out in the Dalitz plot.

C. Omnès functions
For the S-wave isospin I = 0 amplitude, we use the coupled-channel Omnès function from a dispersive summation scheme [49,50] which implements constraints from analyticity and unitarity. The method is based on the N/D ansatz [51], where the set of coupled-channel integral equations for the N -function are solved numerically with the input from the left-hand cuts which we present in a model-independent form as an expansion in a suitably constructed conformal mapping variable. These coefficients in principle can be matched to χPT at low energy [52]. Here we use a data-driven approach, and determine these coefficients directly from fitting to Roy analyses for ππ → ππ [53], ππ → KK [54] and existing experimental data for these channels. After solving the linear integral equation for N (s), the D-function (inverse of the Omnès function) is computed. The obtained coupledchannel Omnès function has already been successfully applied for the photon-fusion reactions γ ( * ) γ ( * ) → ππ in [37,[55][56][57]. The Omnès function for the D-wave (I = 0) is constructed directly from the ππ phase shift [53] and given by since the inelasticity around f 2 (1270) peak is suppressed [48].

IV. RESULTS AND DISCUSSION
For the S-wave contribution we write a twicesubtracted dispersive representation. Due to the coupledchannel there are in total four subtraction constants. For the D-wave we allow for one subtraction. Even though the dispersive integrals are formally convergent with less subtractions 3 they acquire significant corrections from the integration over large s. Therefore, we implement over-subtracted dispersion relations in order to reduce the sensitivity to the high energy region and the effects of additional unknown left-hand cuts, such as possible Dmeson loops or contact interaction [32,41]. To check on the physical importance of the latter, we will also compare in the following for the S-wave contribution the fitted subtraction constants with the sum rule result which one would obtain from a once-subtracted dispersive formalism. Since the dispersion relations in Eqs. (18) and (21) are written for I = 0 we need to encode the transformation coefficients between isospin and the physical amplitudes 3 For the S-wave both h  Table I). The BESIII data is taken from [10], which was normalized to the total cross section given in [31].
Therefore, for e + e − → J/ψ π + π − one obtains + Ω Below we perform a simultaneous fit to the π + π − and π ± J/ψ invariant mass distributions [10] together with the total cross-section data for σ(J/ψK + K − ) [26]. To ensure that the e + e − → J/ψK + K − total cross-section constraint is well accounted for and contributes realistically to the total χ 2 , we re-scale its error by the amount of experimental data points above the KK threshold in the ππ distributions. In our fits we therefore minimize where The Due to an overlap of left-and right-hand cuts, the subtraction constants (a, b, c, d, e) can in principle be complex, which together with the product F γ * πZ C Zψπ leaves us with eleven parameters for each e + e − center-of-mass energy to describe the data. We definitely do not want to over-fit the data and describe some variations in the data that could just be statistical noise. Therefore, we decided to start with the most economical fit in which we fit four parameters as described in the following, and    (31) which were adjusted to reproduce the empirical ππ and πJ/ψ invariant mass distributions together with the cross-section σ(J/ψK + K − ) at e + e − center-of-mass energies q = 4.23 GeV and q = 4.26 GeV. Tildes on top of a subtraction constants indicate that they are given relative to the couplings constants entering h (J),Zc 0,++ , for instanceã ≡ a/(Fγ * πZ C Zψπ ). For easier comparison of the fits with real subtraction constants (φi = 0) and the fits with complex subtraction constants (φi = 0), we restricted φi in the region (−π/2, π/2), i.e. allowing to have ± signs in front of the absolute value. Errors on fit parameters are shown in brackets.
will then compare it with our best fit which has seven parameters. The summary of the fit results is given in Table I.
We start with the case when all the subtraction constants are real in the S-wave while for the D-wave we use an unsubtracted dispersive representation. It turns out that the fitted value of the c parameter is consistent with zero and therefore it is justified to ignore it for this initial fit. This leaves us with four real parameters. Even though this parameterization is not perfect, it provides a good description of the data as shown in Fig.2. In particular, in the ππ mass distribution the dip structure around the KK threshold comes out naturally in our formalism due to the f 0 (980) resonance. In addition, the PDG [48] averaged mass and the width of Z c (3900), m Z = 3.8872(23) GeV and Γ Z = 28.2(2.6) MeV, seem to be well in agreement with the data for the πψ mass distribution. Furthermore, it is worth mentioning an interesting observation: if we fit only the πψ invariant mass distribution for q = 4.23 GeV or q = 4.26 GeV, the post-diction for the ππ distribution reproduces very well the major features of the data 4 . This implies that our framework has the correct ingredients in the simultaneous description of the data. As seen from the parameter values of Fit 1 in Table I, we also find that they not vary much between q = 4.23 GeV and q = 4.26 GeV. This is in accordance with our expectation since the consid- 4 The opposite is not true, because by fitting only the ππ distribution it is hard to constrain well the parameters of the Zc state and the post-diction of the πψ distribution is then only qualitative.
ered e + e − center-of-mass energies are different only by 30 MeV. Therefore, the parameters of Fit 1 determine the starting values of our improved fit. A significant improvement over Fit 1 can be obtained by adding a phase to the parameter a and to a lesser extent also to the parameter b, since the subtraction constants a and b are mainly responsible for the description of the data below the KK threshold. The region above KK threshold is a bit more complicated since the c and d parameters play a significant role there. Due to the absence of the KK mass distribution data, we keep the subtraction constants c and d real. From the analysis of different fits we found that for q = 4.23 GeV a small non-zero value of the phase φ b allows to improve the fit more, while for q = 4.26 GeV the parameter c plays a more prominent role. In addition, we allow for one subtraction in the D-wave, which may differ from the unsubtracted sum rule value. As a result, we decided to limit ourselves to the following "best" fit scenario with seven parameters: for q = 4.23 GeV we consider the product F γ * πZ C Zψπ , a, φ a , b, φ b , d, and e as fit parameters, while for q = 4.26 GeV we consider the product F γ * πZ C Zψπ , a, φ a , b, c, d, and e as fit parameters.
The resulting parameters and χ 2 are collected under the Fit 2 in Table I and shown in Fig.3. We see that our results are in very good agreement with the data. As a conservative error estimate we show in Fig.3 the spread between the Fit 1 (our most economical fit) and Fit 2 (our best fit) results. We found that the remaining parameters, beyond the seven parameters considered, have a rather small effect on the ππ and πψ distributions and can be determined only when very precise data will be available. The black curves are the total fit results. The individual contributions from the ππ re-scattering in the S, D waves and Zc(3900) are indicated by red, green and blue curves, respectively. The shaded bands indicate the spread between Fit 1 (thin curves) and Fit 2 (thick curves) results. The BESIII data is taken from [10], which was normalized to the total cross section given in [31].
It is instructive to compare the fitted values of theb,d andẽ parameters from Table I  which forb andd are approximately 20 times smaller in magnitude (and even more forẽ) than the fitted values. This implies that besides the direct production of the Z c (first term in Eq.(30)), which is responsible for the peak regions in the πψ distribution, the two pions are predominantly produced directly in the transition from the Y state to the J/ψ state through a contact term and subsequently rescatter. Our analysis thus shows that the rescattering of the two pions happens predominantly without going through the Z c (3900) state.
Since we obtained a simultaneous and accurate description of the BESIII data for the π + π − and π ± J/ψ invariant mass distributions, we find it justified to predict the K + K − mass distribution. As one can see in Fig.4, the obtained shape has a rapid rise just above the threshold, which is quite different from the pure phase space, i.e. when K ++ (s, t) is replaced by a constant. This behavior is due to f 0 (980) resonance and we expect to see it in future experimental measurements. For completeness, we also provide the prediction to K ± J/ψ mass distribution, which is just a pure phase space in our approximation.

V. SUMMARY
In this work, we provided a quantitative and simultaneous description of the π + π − and π ± J/ψ invariant mass distributions of the recent BESIII data on e + e − → J/ψ π + π − together with the total cross-sections σ(J/ψK + K − ) at e + e − center-of-mass energies q = 4.23 GeV and q = 4.26 GeV. A crucial element of our analysis is the well established charged exotic state Z c (3900), which we account for explicitly in t-and u-channels. The final state interaction of the two pions in the S-and Dwaves is treated using the dispersion theory. For the S-wave, we consider coupled-channel unitarity since the kinematical region goes beyond the inelastic channel KK and the effect from the f 0 (980) resonance impacts significantly the observables. On the other hand, for the Dwave a single-channel Omnès approach is adopted since the lowest resonance in that channel, the f 2 (1270) tensor resonance, decays predominantly into two pions. The final amplitudes depend on a set of subtraction constants, which have been fitted to the BESIII data. A simultaneous description of π + π − and π ± J/ψ mass distributions together with the cross-sections σ(J/ψK + K − ) is achieved through a four-parameter fit. We showed that the latter can be further improved by adding phases to the subtraction constant and allowing for one subtraction in the D-wave contribution. We found that the resulting seven parameter fit yields a very good description of the π + π − and π ± J/ψ mass distributions together with the cross-sections σ(J/ψK + K − ). Our dispersive formalism shows that besides the direct production of the Z c , responsible for the peak regions in the πψ distributions, the two pions are predominantly produced through a contact term in the transition from the Y state to the J/ψ state and subsequently rescatter. For the e + e − → J/ψ KK we provided the first theoretical prediction for the twokaon invariant mass distribution, which is significantly different from the pure phase space.
The constructed amplitudes provide powerful tools to analyze future data by the BESIII and Belle II Collaborations. It can also be readily applied to study e + e − annihilation into Υ(nS)π + π − , where charged bottomonia like Z ± b states have been observed.