Amplitude analysis and the nature of the Zc(3900)

The microscopic nature of the XYZ states remains an unsettled topic. We show how a thorough amplitude analysis of the data can help constraining models of these states. Specifically, we consider the case of the Zc(3900) peak and discuss possible scenarios of a QCD state, virtual state, or a kinematical enhancement. We conclude that current data are not precise enough to distinguish between these hypotheses, however, the method we propose, when applied to the forthcoming high-statistics measurements should shed light on the nature of these exotic enhancements.

the exchange of a cross-channel D 1 (2420), which is a good candidate to create an anomalous cusp. The prediction of the model was compared with the ππ and J/ψ π spectra of the J/ψ π + π − decay mode. The authors conclude that the cusp alone is not sufficient to describe the Z c (3900) peak and argue in favor of a resonance, although no quantitive measures are given. Numerous other works on cusps and/or poles typically assume a particular scenario for producing peaks and compare model predictions to a subset of available data [16][17][18][33][34][35][36].
Given that the available published data are not corrected for acceptance or efficiencies, and there is no polarization information, it is difficult to make a case for a systematic fit of all the datasets. Nevertheless, we will attempt such an analysis. On the theoretical side, we use several parametrizations of the amplitudes which focus on the role of various singularities, without entering into the details of which model would be able to describe their microscopic origin.

Amplitude model
Consider the three-body decay A → BCD. Under special kinematic conditions [37], a cusp in the mass distribution of BC can be generated, if there is another available direct channel and if a resonance occurs in one of the two crossed channels near the physical region [38,39]. In the absence of a coupled channel, the crossed channel resonances lead to an enhancement in the Dalitz plot, which cancels out upon mass projection [40]. Such cusps are part of the production amplitude, aka left hand side branch points of partial waves. In addition to this, partial waves have direct channel (right hand) singularities, like threshold branch points, or virtual or resonance poles. The definition of the channels relevant to this analysis is given in Figure 1. The peak at √ s ∼ 3900 MeV may thus originate from a true s-channel resonance pole (the Z c ), a virtual state, the left hand branch point, or a combination of both. The best candidate to produce a triangle cusp is the D 1 (2420) resonance in the t-channel process Y D → πD * . We also consider other possible exchanges, like the D 0 (2400) in Y D * → πD, and the f 0 (980) and σ in Y J/ψ → ππ, but the induced s-channel singularities are further away from the √ s ∼ 3900 MeV region, and give little contribution to the peak. If instead the peak is due to a pole singularity, the amplitude analysis can provide insights into the (phenomenological) microscopic nature of the Z c (3900). Consider the schematic plot in Figure 2. The poles related to compact QCD states are expected to become narrower and narrower (i.e. approach the real axis of the complex s-plane) if the coupling to the open channels is made weaker and weaker (for example, in the large N c limit [41][42][43][44][45]). Thus, they are expected to be on the sheet closest to the physical axis, the II sheet if below theDD * threshold (blue dot in the figure, reached from the physical axis with path a), or the III sheet if above (red dot in the figure, reached from the physical axis with path b). A bound state generated by inter-hadron forces would also migrate to the real axis upon switching off of the coupling to the lighter channel, and it is likely to lay on the II sheet as well. On the other hand, poles on the IV sheet (green dot) are too far from the physical axis, and would likely stay on the unphysical sheet [46]. The latter case can thus be interpreted as a virtual state, i.e. meson-meson configuration for which the attractive interaction is not strong enough to bind the constituents, but nevertheless provides an enhancement in the scattering amplitude, with a typical cusp-like shape.
The information about the angular distributions is scarce. The only published plots confirm that the Y → πZ c (→ DD * ) decay is dominated by the S -wave [25,26]. In absence of this, there is no point of considering spin. Thus we treat all particles as spinless interacting in the S -wave, at the same time we use physical masses and widths. In doing so, the constraint on the Y → πZ c angular distribution is automatically fulfilled. One may ask if this approximation is valid for the D 1 meson as well, which is known to decay into D * π in D-wave [47]. However, the D-wave barrier factor (function of t) does not affect the s-channel projection we are interested in, and can be effectively absorbed in the coupling. Moreover, the helicity distribution (s dependence) turns out to be rather flat if all the external spins are included, and if the Y →DD 1 decay is dominated by the S -wave (the weak sensitivity to this distribution was already commmented in [48]). In this respect, the spinless S -wave approximation gives a more realistic description with respect to a spinless D-wave treatment, of the D 1 . This choice reduces the number of free parameters, and only turns out in a poorer description of the reflected peak in the J/ψ π channel at ∼ 3.45 GeV.
We denote by f i (s, t, u) the scalar amplitudes for the two reactions shown in Figure 1, with i = 1 referring to Yπ →DD * and i = 2 to Yπ → J/ψ π. These are given by sums over a finite number of isobar amplitudes in the various channels [49], with z x being the cosine of the scattering angle in the center of mass frame of the x = s, t, u channel. We consider all the exchanges to happen in S -wave, the higher waves being kinematically suppressed, a (x) l,i = 0, for l > 0. The s-channel partial wave amplitudes are given by By construction, the isobars a 0,i contain right hand singularities only, whereas the projections of the crossed channels isobars induce left hand singularities in the b 0,i amplitudes. Unitarity determines the discontinuity ∆ f l, with t i j the 2 × 2 S -wave scattering matrix, and ρ j the phase space in the j channel, i.e. ρ j (s) = λ 1/2 s, m 2 j,1 , m 2 j,2 /s. The solution to Eq. (3a) is given by the well known Omnès representation [39], with s j the threshold of channel j. We ignored possible contributions from left hand singularities in the scattering matrix. We subtract the integral once to improve its convergence, and to take into account any other short-range exchange.
The original projections b l,i (s) are thus modified by an additional term describing the final state interactions. If this happens only for a finite number of partial waves (only S -waves in this model, Eq. (4)), the partial wave series can be summed back and it simply reconstructs the original isobars in the crossed channel, where the subtraction constants c j are explicitly shown. We do not expect the fits to the Dalitz plot projections to be sensitive to details of the lineshapes in the crossed channel, so we parametrize the isobar amplitudes as simple with x = t, u. With these, we define the dispersed projections, cf. Eq. (5) by Figure 2: Schematic representation of the scattering amplitude t i j as a function of complex s. The zig-zag lines represent the unitarity (right-hand) cuts. The physical axis connects to the I sheet right on top of the unitarity cut. Below the (heavier)DD * threshold, the closest unphysical sheet is II (see path a). A pole on the II sheet below threshold (blue dot), if close enough to the real axis, will produce a Breit-Wigner-like lineshape. Similarly, above theDD * threshold, the closest unphysical sheet is the III, (see path b) and similarly a pole on the III sheet above threshold, if close enough to the real axis, will result in a Breit-Wigner-like lineshape. On the other hand, poles on the III sheet below threshold (red open circle), or poles in the IV sheet (green disk) are further from the physical region, but can still give rise to a cusp-like peak on the physical axis, if they are close to theDD * threshold.
with r = D 0 , D 1 , f 0 , σ referring to the various cross channel exchanges that we take into account. The thresholds s r and the phase space ρ r are related to the channel the exchanged resonance appears in, namely channel 1 for D 0 , D 1 and channel 2 for f 0 , σ. The final expressions for the amplitudes are given by and and the expression for the Dalitz projections being and similarly for the projections in the √ t or √ u variables. Until now, we have not given any detail on the nature of the t i j scattering matrix. We use a K matrix parametrization, This parametrization contains spurious left hand cuts, which we remove by approximating [50] Alternatively, we have considered a Chew-Mandelstam phase space, but this choice has very little impact on the fits, and we will not discuss it any further. We consider four different scenarios: 1. III: In this case we consider the parametrization which is as close as possible to the one used in the original experimental analyses, even if it violates unitarity. To wit, for the K matrix we use the Flatté parametrization, i.e. K i j = g i g j /(M 2 − s). This choice produces poles in the closest sheet to the real axis, i.e. the III sheet above theDD * threshold, or the II sheet below threshold. The former case might be interpreted as a genuine QCD state, the latter could be a QCD state of a hadron molecule. We artificially remove the influence of triangle singularities by imposing H = 0, thus breaking unitarity. 2. III+tr.: The K matrix is as in case "III", but we reinstate the correct value for H, which gives rise to a triangle singularity close to the physical region. That is, the S -waves in the s-channel near the physical region, can have both the resonance pole and the logarithmic branch point from the triangle singularity.
3. IV+tr.: In this case we choose for K a symmetric constant matrix. This choice can produce poles in the IV Riemann sheet that can be interpreted as virtual states with respect to the heavierDD * channel. This would be more likely due to hadron-hadron interactions. 4. tr.: The K matrix is as in case "IV+tr.' except that we force the possible pole in t to be far from the J/ψ π threshold. We do this by imposing a penalty on the χ 2 which linearly decrease with the distance of the pole from the point s 0 = 15 GeV 2 , which corresponds to the position of the peak, and vanishes outside the disk of radius 10 GeV 2 . With this model we can assess whether the triangle singularity alone is able to generate the observed structure.
Similar analyses have been perfomed in [36,51], considering scenarios comparable with III+tr. and IV+tr., although the parametrizations employed are different.

Description of the dataset
All the relevant mass distributions that we discuss are not corrected by acceptance or efficiency. This prevents us from giving the absolute normalization of our amplitudes, or to quote physical values for the couplings. Most of the experimental analyses include some reducible incoherent background from sidebands or MonteCarlo (MC) simulations, which should be subtracted before comparing with our amplitude prediction. However, in the analyses we consider, this background seems to be rather small (∼ 15% of events) and flat [23,26,29]. Thus, since we do not give absolute normalizations of the amplitudes and this background does not affect the shape of the distributions we simply neglect it. The only exception is the neutralDD * π 0 channel [28], where the mis-reconstructed events are a large fraction of the Dalitz plot and have a nontrivial shape. A curve parametrizing this background, obtained from MC simulations, is shown in [28], but with no associated uncertainties. In our analysis we use the same shape to subtract from the signal and assume Poissonian uncertainties.
As discussed in the introduction, the Z c (3900) has been observed in e + e − → J/ψ π + π − by BESIII at fixed beam energy of E CM = 4260 MeV [23]. We include in the fit the three projections of the Dalitz plot m(J/ψ π + ), m(J/ψ π − ), and m(π + π − ) quoted in the paper, ignoring their correlations. Belle published a similar analysis, but the final state is produced in association with an undetected ISR photon [24]. The systematics which affect this observation mode are rather different from those by BESIII, and since this dataset is smaller we do not use it. Similarly, we do not consider the low statistics analysis of the CLEO-c data [27]. BESIII also reported the observation in the neutral channel, e + e − → J/ψ π 0 π 0 [29]. The paper shows only the m(J/ψ π 0 ) projection, at the energies E CM = 4230, 4260, and 4360 MeV. The distributions at 4230 and 4360 MeV are shown only for m(J/ψ π 0 ) > 3650 MeV. We include in the fit the two datasets at 4230 and 4260 MeV. To match the charged data, the 4260 MeV dataset has been rescaled by isospin, binning, and efficiency (with the values quoted in [23,29]).
For the open charm channel, we consider the double-tag analysis by BESIII [26] of e + e − →D 0 D * + π − and e + e − → D * 0 D + π − . The paper quotes the mass projection m(DD * ) only, at the energies E CM = 4230 and 4260 MeV. We include in the fit all four datasets. The previous BESIII single-tag analysis [25] is somehow statistically independent from the latter, but the data are affected by larger incoherent backgrounds from mis-reconstructed D ( * ) mesons, and we do not include them in the fit. Instead, we consider the four m(DD * ) distributions of the neutral channel e + e − →D 0 D * 0 π 0 and e + e − → D − D * + π 0 , at energies E CM = 4230 and 4260 MeV. We subtract the incoherent background from these data, then we rescale to match the number of events in the charged channel, to take into account the unquoted different efficiencies.
Our full dataset has 566 experimental points. We work in the isospin symmetric limit, so we use m π = (2m π + + m π 0 )/3, and m D ( * ) = (m D ( * )0 + m D ( * )+ ) /2. Four points happen to be below the iso-symmetricDD * threshold, and are removed from the fit. The values of masses and widths of the final state mesons, and of the intermediate D 1 (2420) and D 0 (2400) are taken from the Particle Data Group (PDG) [52]. Since we are not directly interested in the m(ππ) distribution, we parametrize the ππ resonances with "effective" f 0 (980) and σ, whose masses and widths are M f 0 = 920 MeV, Γ f 0 = 223 MeV, and M σ = 112 MeV, Γ σ = 906 MeV, respectively, as determined from a preliminary fit to the ππ distribution only. The data points at m(ππ) > 1 GeV are not well described by this choice, and since they do not affect the J/ψ π andDD * distribution we are interested in, we remove these points from the fit. This parametrization breaks unitarity in the t-channel, and does not include either the KK channel (important to have a good description of the f 0 ), or the Adler zero (which generates the σ pole). Nevertheless is good enough to describe the the projection in  [23]; red points are the J/ψ π 0 data [29], rescaled as described in the text. As expected, the fit does not reproduce the peaking structure at 3.45 GeV, which is the reflection of the peak at the right, and would require to take spins properly into account. the s-channel. We remark that, because of these approximations, the Breit-Wigner parameters we quote for the σ and f 0 can by no means be compared with the right ones extracted with more refined techniques. Higher statistics will require a thorough parametrization of the ππ scattering [53,54]. We consider the datasets at the different center-of-mass energies as independent samples, that is couplings at different E CM are independent fit parameters. Thus our model has 15 fit parameters: for each one of the two centerof-mass energies, the amplitudes in Eqs. (8) and (9) have one coupling for each one of the 4 exchanged resonances, and the two short-range coefficients (subtraction constants). Both center-of-mass energies share the same K matrix, which is parametrized with three constants K 11 , K 12 , and K 22 in the scenarios IV+tr. and tr., and by two couplings and a mass (g 1 , g 2 , M) in the scenarios III and III+tr.. Considering this, and the number of points removed from the dataset, the four fits have 532 degrees of freedom.

Fit results
We perform a minimum χ 2 fit using minuit [55]. In Figures 3, 4, 5, and 6 we show the results of the fits for the four scenarios. The starting values of the fit parameters have been set by looking for the best χ 2 of O(10 4 ) preliminary fits with randomly chosen initial parameters. The mean value and uncertainty of the fitted curve have been computed using the bootstrap technique [56][57][58], which allows us to take into account correlations among fit parameters and to properly propagate the uncertainties not only to all the observables but also to any quantity that can be extracted from the amplitude (e.g. the pole positions). Specifically, for each one of the four models, we generate 2000 datasets by randomly sampling the experimental points according to Gaussian distributions. For each pseudo-dataset, we perform an independent minimum χ 2 fit, using as initial conditions the ones of the original fit. We can thus select the best 68% fits (1σ confidence level). The χ 2 s of the best fits are reported in Table 1.
All the models have a χ 2 /DOF ∼ 1.3, and give a rather good description of the dataset, as can be seen from Figures 3, 4, 5, and 6. The peak at √ s 3.4 GeV is due to the reflection of the structure at right, and cannot be reproduced properly if spins are neglected. To show separately the contribution of the triangle singularity and of the pole in the scattering matrix, in Figure 7 and 8 we show the magnitude and the phase of the t 12 amplitude, of the unitarized term c 1 + H(s, D 1 ) + H(s, D 0 ), and the product of the two. We also plot separately the contributions of the D 0 and D 1 exchanges. We show that the latter only is in the kinematical regime to generate a sharp peak close to the physical region.
The best χ 2 is obtained with model III+tr., but the difference with the other models does not seem significant. To properly compare the quality of one model fit vs another we use the ∆χ 2 estimator on a number of MC generated datasets [59,60]. More specifically, to give the significance, say of model A with respect to model B, we generate 2000 pseudodata samples according to either one of the two models, with the fit parameters obtained from the best fit discussed above. Each data point is generated according to a Poissonian distribution, whose mean value is given by the value of the theoretical model at the center of the bin. Each generated dataset is fitted again with both models, and we fill a histogram with the ∆χ 2 = χ 2 (A) − χ 2 (B) estimator 2 . We can thus compare the distribution of ∆χ 2 of the datasets generated according to A (which is expected to peak at negative values of ∆χ 2 ), with the distribution of ∆χ 2 of the datasets generated according to B (which is expected to peak at positive values of ∆χ 2 ). These distributions are used to calculate the fraction of samples in which ∆χ 2 has a value larger (for A) or smaller (for B) than the one obtained 2 which is equivalent to a likelihood-ratio test, if one assumes Gaussian errors.         from data, which can be translated into Gaussian significance. The ∆χ 2 histograms are shown in Figure 9, and the significances are listed in Table 2. We can appreciate the peculiar behaviour of the ∆χ 2 distributions for the tr. model, which peaks at ∆χ 2 0 and exhibits a long tail towards negative values. This is due to the penalty introduced in the χ 2 to push the pole far into the complex plane, thus affecting the pure statistical meaning of the χ 2 . The significances relative to the tr. scenario have thus to be considered as mere indications (that is why in Table 2 we report them under quotation marks). Anyway, we note that all the significances are never greater than 3σ. These are going to be even more diluted if we were to consider the systematic uncertainties. We conclude that present statistics prevents us from drawing any strong statements, but the robustness of the tools we have discussed here will allow us to distinguish the different phenomenological models, when new data will be available.

Pole searches
The existence of a Z c state is equivalent to the appearance of a pole in the unphysical sheets of the scattering amplitude. As discussed in Section 1, the Riemann sheet where the pole appears can give hints on its microscopic origin. For each one of the three scenarios that allow for the presence of a pole, we can calculate the pole position, and estimate its statistical uncertainty according to the bootstrap analysis we discussed in previous section. In Figure 10 we show the pole position according to the 68% fraction of best χ 2 obtained in the bootstrap analysis. This can be translated into the 1σ region where the pole is expected to occur. The results are summarized in Table 3, and the main observations are as follows: 1. III: The pole appears above theDD * threshold, on the III sheet (the closest to the physical region), and the width is Γ 50 MeV. This is marginally compatible with the value quoted in the PDG, M = 3886.6 ± 2.4 MeV, Scenario III+tr. IV+tr. tr. III 1.5σ (1.5σ) 1.5σ (2.7σ) "2.4σ" ("1.4σ") III+tr. − 1.5σ (3.1σ) "2.6σ" ("1.3σ") IV+tr.
− − "2.1σ" ("0.9σ") Table 2: Significance of each model versus another. The number in the cell AB indicates the probability for the ∆χ 2 generated according to A (B) to be greater (smaller) than the ∆χ 2 obtained from the fit to the real data. The significances relative to the tr. option are affected by the penalty in the χ 2 and should be considered as mere indications.   Γ = 28.1 ± 2.6 MeV [52]. The reasons for this slight discrepancy are twofold: i) in the fits performed in the experimental analysis the sum of the signal (Breit-Wigner) and background (phase-space shaped) is performed incoherently, which tends to provide narrower values for the width; ii) in particular for the J/ψ π 0 π 0 data, we cannot disentangle the Breit-Wigner width from the experimental resolution, effectively giving a slightly larger width to the resonance. 2. III+tr.: The presence of the logarithmic branching point close to the physical region allows for the pole to be slightly deeper in the complex plane, with a width Γ 90 MeV. The mass is still safely above threshold. 3. IV+tr.: In this case the peak is generated by the combination of the logarithmic branching point with the virtual state pole on the IV sheet. Given the presence of the triangle singularity, the position of the pole is not well constrained. The width, defined in analogy to the Breit-Wigner case as 2 Im √ s P , is broader than in the other scenarios, Γ 250 MeV, but the mass is unchanged, albeit with errors of ∼ 100 MeV. 4. tr.: By construction, this scenario does not allow for poles close to the physical region.  Table 3: Mass and width of the Z c (3900) according to the scenarios which allow for the presence of a pole. The error quoted is the 1σ statistical uncertainty obtained with the bootstrap analysis of Section 3.

Conclusions
The literature on XYZ states abounds with discussions about their microscopic nature. In this letter we show how a thorough amplitude analysis can help in constraining the various different phenomenological models. We tested four different scenarios, corresponding to pure QCD states, virtual states, or purely kinematical enhancements. The best fit is obtained for a compact QCD state, but the rejection of the other scenarios is not significant. We conclude that given the present data, specifically mass projections, it is not possible to distinguish between the different hypotheses. Future high-statistics measurements and the study of the full Dalitz plot, thus including angular correlations, will improve the discrimination power of our analysis, in particular by constraining the contribution of the D 1 exchange. This new information, together with a combined analysis of other reactions, e.g. Y(4260) → h c ππ or photoproduction off protons, will allow us to shed more light on the nature of the exotic charmonium sector.
All material will be gathered onto an interactive website which will available online [61,62].
shop. This material is based upon work supported in part by the U.