The $Z_c(3900)$ peak does not come from the"triangle singularity"

We compare contributions from the triangle diagram and the $D\bar D^*$ bubble chain to the processes of $e^{+}e^{-}\rightarrow J/\psi\pi^{+}\pi^{-}$ and $e^{+}e^{-}\rightarrow (D\bar D^*)^\mp\pi^{\pm}$. By fitting the $J/\psi\pi$ maximum spectrum and the $D\bar D^*$ spectrum, we find that the triangle diagram cannot explain the new data from BESIII Collaboration at center of mass at 4.23GeV and 4.26GeV, simultaneously. On the contrary, the molecular assignment of $Z_c(3900)$ gives a much better description.


Introduction
The charged charmonium-like state Z c (3900) was observed in J/ψπ ± mass spectrum by BES III Collaboration in e + e − → J/ψππ process [1], and confirmed by Belle [2] and CLEO [3] Collaborations in the same processes. Afterwards, it was also observed in the (DD * ) ± invariant mass spectrum in the process of e + e − → DD * π ∓ , and the quantum number of Z c (3900) was determined to be I(J P ) = 1(1 + ) by angular distribution analysis on πZ c (3900) system [4]. The experimental discovery has stimulated a lot of discussions because of the unique nature of Z c (3900), as it could be the first unambiguous candidate of the long wanted tetra-quark state.
In a recent paper [5], we have made a detailed comparison between the DD * molecule picture and the "elementary" picture, and concluded that Z c (3900) is of DD * molecular nature, using the pole counting method [6].
However, it is also found in the literature another possible mechanism, called anomalous triangle singularity (ATS), to explain the singularity structure at Z c (3900). ATS refers to a branch cut in a three point loop function other than the normal threshold. The study of ATS can be traced back to about 60 years ago. In Ref. [7] S. Mandelstam worked out the ATS branch point and discussed its effects on the deuteron electromagnetism form factor, and in Ref. [8] L. D. Landau applied Landau equations to triangle diagrams to analyse ATS. Extensive studies on the triangle singularity using dispersion techniques can also be found in Ref. [9]. Especially in the paper by Lucha, Melikhov and Simula of Ref. [9], a detailed dispersive analysis is given on different variables.
ATS has attracted renewed interests recently because it may contribute to peaks in some certain invariant-mass spectrums. In other words, some so-called "exotic hadron states" could be just the ATS peak rather than real particles; or even if real exotic hadron states exist, there may be some non-negligible contributions from ATS. For example, it is suggested in Refs. [10][11][12] that the singularity structure of the triangle diagram (see Fig. 1), which contains both the normal threshold effect and anomalous threshold effect, may lead to the peak at 3900 MeV. In Refs. [13,14] it is emphasized that the anomalous triangle singularity may have significant impact in understanding the nature of the near-threshold state.
This paper devotes to the study of triangle diagram contribution to the Z c (3900) peak. In Sec. 2 we give a pedagogical analysis on general three point loop functions using Feynman parameter representation that can be found in most textbooks, and discuss the properties of the ATS. In Sec. 3 we calculate the triangle diagram corresponding to e + e − → J/ψπ + π − and e + e − → (DD * ) ∓ π ± processes and fit the experimental data to test whether the Z c (3900) peak comes from triangle diagram. In Sec. 4 some conclusions are drawn. Basically, it is found that the new experimental results from Ref. [15] play a crucial role in clarifying the issue on triangle diagram contribution: the experimental data indicate that the peak at 4.23 GeV is higher than that at 4.26 GeV, whereas the triangle diagram predicts an opposite behavior. Our analysis reveals that Z c (3900) peak cannot be explained from the triangle diagram contribution from Fig. 1. Hence, combining with our previous analysis in Ref. [5], the molecular nature of Z c (3900) is firmly established.

Theoretical framework
To start let us first look upon a general triangle diagram as shown in Fig. 1. The amplitude for such a triangle diagram contains the following scalar three-point function: where k 1 = q − p 2 , k 2 = p 1 − q, k 3 = q, s i = p 2 i (i = 1, 2, 3), and ǫ → 0 + is a small positive parameter. This scalar three-point function could be evaluated using the standard method as follows. First one can rewrite it in terms of Feynman parameters x, y and z After straightforward calculation, the function can be presented as where One can perform a momentum translation q = l − xp 2 + yp 1 , then To evaluate that integration, Wick rotation is done 1 : Then the integration in momentum space can be worked out anyway regardless of the sign of ∆(x, y) 2 Further, the integration of y is then calculated analytically (for simplicity we omit the −iǫ in the following discussion) where Now it's not difficult to analyse the singularity structures of T (s 1 , s 2 , s 3 ) based on Eq. 8 and 9, and for simplicity we only study the singularities of s 3 variable. There are two "arctan" terms in Eq. 8 ; the singularity of the first term is given by the ±i branch points of arctan function, i.e.
Its solution with respect to the integration variable x is If λ(s 3 , m 2 1 , m 2 2 ) → 0 − , pinch singularity happens. Especially when s 3 → (m 1 + m 2 ) 2 − 0 + , two singularities from Eq. 12 x = m 2 m 1 + m 2 ± i0 + will pinch the integration interval [0, 1] while a(s 3 , x) stays finite, hence it's reasonably concluded that s 3 = (m 1 + m 2 ) 2 is a singularity of T (s 1 , s 2 , s 3 ) 3 . In fact, that result is nothing but the well known normal threshold of the three point function. One can do analogous analyses on the second arctan term of Eq. 8, and it's found that the solution of x is independent of s 3 , so the second term has no contribution to the singularities of T (s 1 , s 2 , s 3 ) with respect to variable s 3 .
On the other hand, if the singularities occur in the denominator a(s 3 , x) = 0, then the solution for x variable is 2 When ∆(x, y) < 0, the −iǫ term guarantees the validity of this integration, see for example, M. E. Pesking and D. V. Schroeder, An Introduction to Quantum Field Theory, Westview Press, 1995, Page 808. 3 In fact there is another solution corresponding to λ(s 3 , m 2 1 , m 2 2 ) → 0 − , that is, s 3 = (m 1 − m 2 ) 2 , which is called pseudo-threshold and only appears on the un-physical sheet. So it's less relevant to our discussion.
When s approaches the following two points, a pinch singularity occurs in T (s 1 , s 2 , s 3 ). However, that singularity may not always appear on the physical Riemann sheet (i.e. the first sheet). Actually there exist two cases leading the singularity to be on the un-physical (the second) sheet: firstly, when the pinch points in Eq. 13 lie off the integral interval [0, 1]; secondly, when the numerator of the integrand in Eq. 8 approaches 0 simultaneously as a(s 3 , x) → 0, giving a well defined value of the integrand. We will meet examples of these two cases later.
To proceed, we in the following discuss two situations, one is that all particles are stable, the other is that some particles are unstable like the case we meet in the real situation (X(4260) → J/Ψππ, DD * π).
For the first situation, . When we substitute the s + in Eq. 14 for s 3 in Eq. 13, it will be found that x is not in the integral interval, thus s + locates on the second sheet. But the situation is more complicated when we focus on s − , which is usually called "anomalous threshold" or "anomalous triangle singularity". Since the pinch points in Eq. 13 stay in [0, 1] when s → s − , one has to test whether the numerator of the integrand in Eq. 8 approaches 0. We have the following observations: , the s − locates on the second sheet and is below the normal threshold; , the s − just rides on the normal threshold; , the s − locates on the first sheet and is below the normal threshold.
To obtain above results, one needs to analyse the behaviour of the numerator of the integral in Eq. 8 at s − by looking upon the arguments of the two arctan functions. For example, when using arctan(−∞) = −π/2, one finds the two arctan terms cancel each other in Eq. 8 when s 3 → s − , so in this case s − is not on the physical sheet. But when the first argument becomes +∞, leading the numerator to be nonzero thus s − appears on the first sheet. In above discussions the behavior of s − is actually the classical example discovered long time ago by Mandelstam in Ref. [7], being used to explain the long tail of deuteron wave function.
In reality, however, the stability condition in obtaining the above results may not hold. For example, one may consider the following kinematics (which corresponds to the kinematics of Fig. 4): After some analyses similar to the situation of stable particles, we find that when s 1 satisfies s − locates on the physical sheet, and is above the normal threshold. Otherwise, s − would be on the un-physical sheet.
In addition, to understand the dependence of s − on s 1 , one gives s 1 a small positive imaginary part, s 1 → s 1 + i0 + . Then s − can be expressed as s − (s 1 + i0 + ) = s − (s 1 ) + i ∂s − ∂s1 0 + derived from Eq. 14. As s 1 increases, the near-threshold trajectory of s − both in classical stable case (see Fig. 2(a)) and under the kinematics of Eq. 15 (see Fig. 2(b)) can thus be drawn. The aim of this paper is to investigate whether Z c (3900) peak is mainly from triangle singularity or not, so the processes X(4260) → π + π − J/ψ, DD * π are considered, with respect to the triangle diagrams shown in Fig. 4 as suggested by Refs. [10] - [14]. We set s 1 to be the square of X(4260) 4− momentum, and the pertinent masses to the masses of those particles, then according to Eq. 16, it's found that the ATS lies on the second sheet when √ s 1 lies between 4230MeV and 4260MeV, hence the Z c (3900) peak cannot be a direct manifestation of the anomalous threshold. We plot the modulus-square of the amplitude in Eq. 8 with different center of mass energies √ s 1 , as shown in Fig. 3.
It is however found that the location and effect of the anomalous threshold are very sensitive to the energy of X(4260). When the anomalous threshold is on the second sheet as shown in Fig.  3(a), the closer it is to normal threshold the more influence it has on the amplitude. Since the anomalous threshold (on the second sheet) can be rather close to the normal threshold, one still needs to check whether the (anomalous and normal) threshold effects can cause the experimentally observed Z c (3900) peak. The discussions above are only brief qualitative analyses aiming at studying the dependence of ATS peak on s 1 variable, and what we have calculated above is not the whole amplitude of that process to fit experimental data, since in the full amplitude there exists X(4260) production process, X(4260) propagator, spin structures, as well as derivative couplings. The detailed formulae can be found in Appendices. A and B, with which the numerical discussions are made in the next section.

Numerical Analyses and Discussions
3.1 Fit to the data of Ref. [1,4] As pointed out in Ref. [13], the ATS contribution to the decay of X(4260) to J/ψππ (see Fig. 3(b)) may have great dependence on the center-of-mass energy. However Ref. [13] just made a rough numerical discussion.
On the other side, the result in Ref. [5] supports molecular state interpretation rather than "elementary" state explanation. Since the aim of this paper is to test whether the triangle diagram can provide the Z c (3900) peak, we make three independent fits to compare the fit quality: • Fit I: X(4260) decays to final states only through triangle diagram as depicted in Fig. 4.
• Fit II: the final states are only produced by DD * rescattering, and the Feynman diagram is shown in Fig. 5.
• Fit III: the mixed situation by combining the triangle diagram and DD * bubble chains, as shown in Fig. 6.
Recall that both the bubble and the triangle diagrams are ultraviolet divergent. We use MS − 1 scheme of dimensional renormalization method to deal with the divergences, which leaves another somewhat arbitrary parameter: the renormalization scale µ (see Appendix. B). Considering the physical process we are studying, it is reasonable to expect that the µ parameter should be around or not much differ from the mass of X(4260). Indeed, we find that when we set the renormalization scale at a reasonable value µ = 5GeV, the bubble chains give a satisfactory fit result as shown in the following text, and it's verified the fit quality is not sensitive to the variation of µ parameter. On the contrary, the triangle diagram does not give a good description to the data when µ is set at 5GeV. The fit results shown in the following text correspond to an unreasonably small µ (in getting the Fig. 7 and Fig. 8). 4 In the numerical fit two sets of data including J/ψπ maximum invariant mass spectrum [1] and DD * invariant mass distribution [4] are taken into account.  The Lagrangian we use is taken from Ref. [5] with three additional pieces and where X µν = ∂ µ X ν − ∂ ν X µ . The corresponding Feynman rules can be found in Appendix. A. With these preparations, it is possible to make a combined fit on the J/ψπ maximum spectrum and DD * mass distribution [1,4]. Except for various coupling parameters, two parameters for the DD * incoherent background and two normalization constants are further introduced. In total, there are 8 and 10 free parameters for Fit I and Fit II, respectively.
Since the value of the center-of-mass energy of X(4260)/γ * severely influences the ATS contribution as discussed previously, we in the fit also carefully analyze the effect of energy resolution both in X(4260) channel and in Z c (3900) channel. However it is found that the effect of the energy resolution does not obviously improve the fit quality, since the energy resolution parameters are much smaller than the particle widths. 5 The fit results are shown in Fig. 7. The χ 2 /d.o.f of Fit I and Fit II are 2.5 and 0.96, respectively. In Fit II we find a bound state pole of DD * , √ s = 3.8747 ± 0.0148iGeV, which is consistent with our previous result in Ref. [5]. Further, Fit III gives very similar χ 2 to Fit II with pole located at √ s = 3.8749 ± 0.0145iGeV, which suggests that the triangle diagram plays only a minor role as comparing with the bubble chain contribution. In Fit II and Fit III, the denominators of amplitudes of the processes e + e − → J/ψππ(DD * π) take the form where the λ 1 represents the DD * D D * contact coupling constant, and c 0 simulates the contribution from other lighter channels to the Z c (3900) width. The function D(l) is the DD * meson loop integral. The parameters λ 1 and c 0 which decide the pole positions are listed in Table 1.  Although the bubble chains (the molecule picture) fit data better than the triangle diagram, the latter still cannot be firmly excluded here. Generally speaking, one believes that the cusp effects are much weaker singularities than poles, and the triangle diagram does not generate poles except branch point singularities. However, the undetermined overall normalization factor in the present fit makes up the defect of branch cuts (i.e., being weak in general), and hence prevents us from excluding the triangle diagram mechanism.

Fit to the new data
Fortunately, the new data from BESIII Collaboration [15] indicates that there are more events in Z c (3900) peak at √ s 1 = 4.23GeV (with a integrated luminosity of L = 1092 pb −1 ) than at √ s 1 = 4.26GeV (with L = 827 pb −1 ), after background subtraction (see Fig. 8). On the contrary, the magnitude of the triangle diagram in Fig. 4 at the Z c peak is smaller when √ s 1 = 4.23GeV comparing with the magnitude when √ s 1 = 4.26GeV. 7 It is noticed that the s 1 dependence of 5 G. Y. Tang, private communications. 6 Other coupling constants and the normalization constants are multiplied to each other and are not quite interesting physically, so we do not list them here. 7 Since the latter is closer to theDD 1 threshold. The results are shown in Fig. 8. The χ 2 dof = 5.3 for pure triangle diagram (Fit I') and the χ 2 dof = 1.6 for pure bubble resummation (Fit II'). The pole of Fit II' locates at √ s = 3.8804 ± 0.0150iGeV. Hence the triangle diagram gives a much worse fit comparing with the bubble chain diagram, and hence can be ruled out. Apparently, the new data [15] are crucial in supporting of the DD * molecule explanation of Z c (3900). We also fit the new data of Fig. 8 in the mixed situation (Fit III'), however it does not obviously improve the total χ 2 comparing with Fit II'. I.e., the χ 2 d.o.f. is almost the same as Fit II', with the pole location √ s = 3.8822 ± 0.0119iGeV, see Table 2. Hence we may draw the conclusion that the bubble chain mechanism plays a dominant role in reproducing the experimentally observed peak structure.
Here we should mention that in the above fit using triangle diagrams, the renormalization scale µ runs to a ridiculously small number, ∼ 10 −7 GeV. If we fix µ at 5 GeV, the fit quality of triangle diagrams gets even worse.

Conclusion
To summarize, we have investigated whether the triangle singularity mechanism proposed in the literature can be responsible for the experimentally observed Z c (3900) peak. It is found that, though the triangle diagram could barely explain the line shape, up to an arbitrary normalization constant, it fails to explain the dependence of the process e + e − → J/ψπ + π − on the center of mass energy, not to mention the weird value of the renormalization scale it requires. Therefore, we conclude that Z c (3900) peak is dominantly contributed by the pole of the DD * molecular state.
The authors are grateful to Chang-zheng Yuan for bringing Ref. [15] to our attention, and Guang-Yi Tang for valuable discussions. We also would like to thank Qiang Zhao and Feng-Kun Guo for helpful discussions. This work is supported in part by National Nature Science Foundations of China (NSFC) under Contract Nos. 10925522, 11021092.
Note added: when this paper is being completed, we became aware of a recent paper [19], where the authors also attacked the same problem using four different amplitude parameterizations. They reach a conclusion that at this stage they can not have a preference on one of these parameterizations, which contradicts our conclusion. We point out here that they were only able to use the old data (at c.m. energy 4.23GeV) in the neutral channel of J/ψπ 0 [20], which is much worse in statistics comparing with the data of [15]. Hence we urge the authors of Ref. [19] redo their analysis with the new data incorporated and compare with our result.

Appendix A The pertinent Feynman rules
For charged final state DD * or J/ψπ, Feynman vertices in Fig. 4 are given as follows.
iV DD * µ DD * ν = 2iλ 1 g µν , where µ, ν are Lorentz indices of spin−1 particles, f π is the pion decay constant, the flat spacetime metric is g µν = diag{1, −1, −1, −1}, and p A denotes the four-momentum of particle A in the vertex. As for propagators of D * , D 1 and J/ψ, we adopt the standard formula of Proca fields with four-momentum k mass M P , i.e.
When one calculates the modulus-square of the amplitudes, the polarization summing of the final state D * and J/ψ is needed, and the formula of Proca field physical polarization summing is as follows: with ε(k) being the final state Proca field polarization vector and k being the four-momentum of the Proca field. The initial e + e − → γ * (µ) QED process is written as v(p + , s)(−ieγ µ )u(p − , r) , with the initial positron 4− momentum being p + and polarization being s, and the electron's p − and r. Moreover, the γ * (µ) and X(4260) (ν) two-point coupling is iV γ * µ X ν = 2i(g µν p 2 1 − p µ 1 p ν 1 )c γ , with p 1 = p + + p − . That form ensures QED Ward identity.
B Amplitudes of the total process In this section we denote the 4− momentum of X(4260) as p 1 , the 4− momentum of final state π in the 3− vertex as p 2 , and 4− momentums of final state π and J/ψ in the 4− vertex as p π and p J/ψ , respectively, and similar for the final state D and D * (see Fig. 9). The total amplitudes with triangle diagrams of DD * π or J/ψππ final states can be written as iM =v(p + , s)γ α u(p − , r) 2iec 1 (p 2 1 g αµ − p 1α p 1µ ) while the triangle parts: iT µσ J/ψπ = − with