The BESIII collaboration recently presented results from a high-precision measurement of the exclusive process e þ e -→ D D [1]. Their report included a fit to the data with a coherent sum of eight Breit-Wigner amplitudes, one of which is a possible new resonance labeled as Rð3900Þ. Although the authors point out that the parameters of this purported resonance depend upon the chosen fit model, they nevertheless extract the values m ¼ 3872 AE 14 AE 3 MeV and Γ ¼ 180 AE 14 AE 7 MeV. This ostensible state has already been interpreted as a possible D DÃ P-wave molecular resonance [2].
The measured e þ e -→ D D cross section is compatible with earlier results from BABAR [3,4] and Belle [5]. Although these experiments see a clear peaking structure near ffiffi ffi s p ¼ 3900 MeV, sometimes denoted Gð3900Þ, the collaborations do not claim the presence of a new resonance, chiefly because this behavior was predicted more than 40 years ago in a classic model calculation by the Cornell group [6]. The model computes charmonium masses with nonrelativistic kinematics and an assumed "Coulomb-plus-linear" cc potential. The authors then extend the model by coupling pure quarkonium states to continuum channels (in practice, D ðÃÞ DðÃÞ and D ðÃÞ s DðÃÞ s ), thereby permitting the computation of e þ e -→ hadrons in the charm region [7]. Their result in the D D channel shows a clear enhancement near 3900 GeV, which the authors attribute to the opening of the D Dà þ D à D channel, as well as to the presence of a node in the ψð3SÞ wave function.
Here we examine this mechanism in a general framework to test its viability for explaining the recent BESIII data. We find that a K-matrix fit can explain the enhancement, and that its primary cause is a combination of interference and the opening of the D Ã D channel. We also find that momentum dependence of the ψð4040Þ resonance coupling to the D Ã D channel over hadronic scales is sufficient to drive the reaction rate down, thereby creating the enhancement that was fit as a Breit-Wigner resonance by the BESIII collaboration. Our result implies that the appearance of a radial node in the ψð4040Þ wave function is not required to explain the steep dip observed in the cross section just above 4 GeV. Nevertheless, we find that fits with a modeled node do improve fit quality and conclude that, although not required to explain the data, there is weak evidence for the presence of a node in the decays of ψð4040Þ. Extracted pole positions are similar to those reported by the Particle Data Group (PDG) [8], but the width we obtain for ψð4040Þ is Γ ¼ 130 AE 30 AE 125 MeV (the second uncertainty reflecting variation over models), somewhat larger than the 80 AE 10 MeV reported by the PDG. We have also determined the e þ e -coupling of ψð4040Þ to be 180 AE 100 AE 170 eV, to be compared with the value 860 AE 170 eV reported by the PDG.
Previous work related to the questions addressed here includes a perturbative treatment of ψð2SÞ -ψð3770Þ mixing in an effective Lagrangian approach [9], wherein the authors treat Gð3900Þ as a resonance but find that it can also be explained by the D Ã D threshold. Du et al. [10] examine charm production in a similar mass range in a Lippmann-Schwinger approach, using contact interactions that are constrained by heavy-quark symmetry, and find a pole near 3900 MeV. Uglov et al. [11] perform a K-matrix analysis of Belle data in the energy range 3.4-4.7 GeV. This model describes the channels D D, D Ã D, D Ã D Ã , and D Dπ using 40 parameters, but no bare pole near 3900 MeV is found to be required.
Cao and Lenske [12] examine the line shape of ψð3770Þ in a coupled-channel T-matrix approach, finding good agreement with the data, although with a small pole mass of 3716 (30) MeV. They also conclude that the Gð3900Þ enhancement is due to a distortion of the ψð3770Þ tail by the D DÃ threshold.
The K-matrix formalism has recently been applied to the ψð3770Þ line shape [13]. This work concludes that BES(I,II,III) data are consistent with ψð3770Þ being due to a single pole, and that its non-D D branching fraction is less than 6%, [the first version of this paper reported a larger branching fraction because the authors used an uncorrected σðe þ e -→ D DÞ cross section].
Nakamura et al. have performed a partially unitary coupled-channel analysis of charm production over the energy range 3.75-4.7 GeV, finding evidence for two Z c poles [14].
Finally, von Detten et al. [15] examine charmonium production near 4.2 GeV with a combination of perturbative and Dyson-equation methods, obtaining evidence for a resonance with a pole at 4227 -25i MeV, which they interpret as a D 1 D molecule (although with a dataset that is not as complete as the one used here).
In Sec. II we describe the K-matrix fit model (introduced in Ref. [16]) used here, and we present the results of the fits to e þ e -→ D D data in Sec. III. Section IV summarizes our findings and concludes. Corrections to the measured cross sections for σðe þ e -→ D DÞ are presented in the Appendix.
It is generally recognized that fitting complex data with an amplitude model that is a sum of Breit-Wigner resonances is inadequate in the case of overlapping resonances. Indeed, an amplitude that is a naive sum of Breit-Wigner profiles does not respect unitarity. This type of model is also inadequate when nonresonant behavior is important, such as when nearby thresholds can cause dramatic changes in event rates. In view of charmed thresholds at DD Ã (3871.7 MeV for D0 D Ã0 and 3880.0 MeV for D -D Ãþ ) and D þ s D - s (3936.7 MeV), it is prudent to construct an amplitude model that incorporates coupled channels and enforces unitarity. We choose to employ the K-matrix model as implemented in Ref. [16] for this purpose, presenting only the most crucial details here.
The K-matrix model couples resonances R to two-body continuum channels μ via couplings g R∶μ , and permits scattering between continuum channels via a contact-like term f μ∶ν , as follows:
Guided by constituent quark models and phenomenology, the resonance couplings are written as
Here, β is a typical universal hadronic scale that is either fit to the data or fixed in different K-matrix models. The relative momentum in the channel μ (with component masses m 1;2 ) is
and l μ is the dominant orbital angular momentum in the channel μ.
In view of the claims of the Cornell group, we have also fit the data with a coupling model that permits a radial node in the simplest possible way:
where the node location k 0 is an additional fit parameter. In a similar fashion, the nonresonant couplings f are parametrized by
This form was chosen in Ref. [16] for simplicity, consistency with the coupling model, and because the inclusion of high-energy damping was found to be advantageous for fit robustness in the case of B ðÃÞ ðsÞ BðÃÞ ðsÞ and ϒ production. The model was not applied to the e þ e -channel, where a hadronic form factor is unnecessary.
In general, one expects the couplings g R∶μ to be isopin independent, since the resonances themselves are isoscalar. However, the continuum couplings f μ;ν for different charge states can vary substantially since they receive contributions from both the I ¼ 0 and I ¼ 1 components of the photon, which can induce substantial isospin asymmetry through mixing with the threshold strong-interaction processes. We will examine the effect of varying isospin constraints in the following discussion.
These model choices are inserted into the modelindependent K-matrix formalism. The scattering amplitude is written as
where C is the Chew-Mandelstam function that introduces the imaginary part of the scattering amplitude required for unitarity. The specific form chosen for the Chew-Mandlestam function in this case [17][18][19] is designed to improve the analytic properties of the amplitude in the vicinity of thresholds. The imaginary part of C is given by
where ρ is a diagonal phase-space matrix in channel space, with elements given by
Finally, Aitchison's P-vector formalism [20] is used to implement production via the e þ e -channel:
where P ν ¼ K ν;ee , while K and Ĉ are defined in the restricted channel space that excludes the initial channel.
Our fit models incorporate between 24 and 30 parameters, as reported in Table I.
Our choice for datasets to be fit, the specific implementation of the K-matrix model, and the fit procedure employed are described in the following subsections.
We choose to model the region ffiffi ffi s p ¼ 3.7-4.2 GeV with two bare charmonium resonances, corresponding to ψð3770Þ and ψð4040Þ [nominally the ψð1DÞ and ψð3SÞ
states of the Cornell model, respectively]. No additional or novel charmonia are assumed to contribute. We truncate the data at 4200 MeV to avoid complications from known exotic states like ψð4230Þ and from the effects of pairs of orbitally excited D mesons such as D 1 ð2420Þ. We have found it necessary to separate the D 0 D0 and D þ D -data, presumably because the thresholds at 3730 and 3739 MeV, respectively, are sufficiently different to drive the differences visible in the cross section data between the two channels, as discussed further below. We also employ (and report in the Appendix) precise unpublished BESIII data for D D cross sections in the vicinity of ψð3770Þ [21]. It is useful to remove initial-state radiation and vacuum polarization (VP) effects from the data when fitting with the K-matrix model because the model is designed to extract hadronic information from the specific final state. To this end, the observed cross
Model 1 plus node [Eq. ( 4)] 3 0 section σ obs of e þ e -→ D D has been related to the dressed cross section σ dre by
where the initial-state radiation correction factor 1 þ δ is obtained via an iterative procedure, following Ref. [22]. The Born cross section σ B , which is the one used in our fits, is then calculated from the expression
where the VP factor 1=j1 -Πj 2 is taken from Ref. [23]. Figure 1 displays the observed, dressed, and Born cross sections for the neutral and charged channels, respectively; uncertainties are statistical only. A common systematic uncertainty of 3.22% [21]. The cross sections are compiled in Tables III and IV in the Appendix.
In view of its claimed importance, the D DÃ þ D Ã D channel is also included in our model fit. BESIII data [24] have been supplemented with data from Belle [25] and CLEO-c [26,27], because BESIII measurements in this channel are sparsely spaced in energy.
Finally, if D D ↔ D DÃ þ D Ã D rescattering does in fact turn out to be important, it is prudent to also consider the nearby D Ã DÃ channel in the coupled-channel fit model. The data employed here are readily available from BESIII [24], Belle [25], and CLEO-c [26,27].
Previous work on bottom production has shown that incorporating the inclusive reaction e þ e -→ b b into the fit model is beneficial both because of the high quality of the data and the additional constraints it imposes upon the fit [16]. We have therefore incorporated an inclusive e þ e -→ cc channel in this work. The cross section was obtained by subtracting light-quark contributions to R in the relevant mass region using Eqs. (9.7)-(9.9) of the PDG review on quantum chromodynamics [8], with α s ðsÞ calculated using RunDec [28], and multiplying by σðe þ e -→ μ þ μ -Þ. While useful, a complication introduced by this method is that the sum of exclusive channels typically does not saturate the inclusive cross section. This deficit is accounted for by introducing a "dummy" channel-meant to represent all the neglected channelsthat absorbs the difference in the rates, as discussed further in the next section.
We note that the D þ s D - s channel opens at 3937 MeV, which provides a another potential source for a significant threshold effect. However, we have chosen not to include this channel in the fit because of the paucity of available data. A very recent BESIII measurement of the D þ s D - s cross section [29] shows a clear ψð4040Þ peak and will be considered in future analyses.
With these choices, the task is to fit a two-resonance, five-channel model to
(the latter two isospin-averaged), and the total cross section data, along with the dummy channel. In general, this task is a difficult one, and we have found it useful to implement the fit by progressively increasing the number of included channels. Fits thus start with data exclusively from the D D channels, which forms the starting point for a three-channel fit to D D and D DÃ þ D Ã D data, and proceeds to the full five-channel model.
The need for a dummy channel to account for missing cc production has been discussed. We have chosen to model this channel in two ways that appear most relevant: as a three-body channel J=ψππ treated as a two-body channel J=ψðππÞ in a relative S wave, or a D þ s D - s channel in a relative P wave.
The stability of the fit was tested by generating 200 pseudo data sets, fitting to these, and then using the resulting minimum of each fit as a new starting point for fits to the measured data. In general, all model parameters are unconstrained; however, we have chosen to impose the resonance isospin constraint g R∶D þ D -¼ g R∶D 0 D 0 , which is expected to hold to good accuracy. A similar constraint might be expected to hold for the cross-channel nonresonant couplings f D þ D -∶ μ and f D 0 D 0 ∶μ , as well as for the elastic couplings
In view of these possibilities, we have tested model sensitivity by fitting five models, varying the constraints and dummy channels as shown in Table I.
Finally, the hadronic scale is set to β ¼ 1 GeV, which was found to be optimal in terms of fit quality and stability in the case of bottom production [16].
All fit models reproduce the data fairly well in all channels, with the pull distributions indicating that the fairly large values of χ 2 =ndf are to a certain degree driven by mild inconsistencies between the experimental datasets (see, e.g., Fig. 4, left). In the two D D channels, where data solely from the BESIII experiment is used, the fit quality is excellent. Also note that we neglect mass differences between neutral and charged D ðÃÞ in D à D þ D Dà and D à Dà , which can affect the fit quality for near-threshold data points.
Results for all the models are very similar, with the exception that the fits for the D Ã D channel (Fig. 3, left) that include a dummy D þ s D - s channel (models 3 and 4) exhibit a stronger rise near 3.95 GeV that better accommodates the high data points in this region. As it gives the smallest overall χ 2 for those models not containing a node, the fit results for model 1 are shown in Figs. 234. The red area in Figs. 34represents 68% confidence levels, while green indicates 90% confidence levels. We also note that, while the data in the range up to 4.2 GeV is described very well, the predictive power of the K-matrix model is very limited. Extending the model to higher c.m. energies will certainly require additional channels and bare resonances.
Open data points are the Born cross section values based on observed cross sections, as reported in Ref. [21]; closed data points are from Ref. [1]. Right: e þ e -→ D Ã DÃ . The red region indicates the 68% confidence level, while green is the 90% confidence level. Black data points are from BESIII [24], red data is from CLEO-c [26,27], blue data are from Belle [25].
We have extracted resonance parameters by analytically continuing the fit model and numerically searching for poles on the closest physical sheets [defined by negative imaginary part of the relative momentum k μ ðzÞ for channels below the real part of the pole mass, and positive for channels above]. Partial widths are obtained by determining the residue at the pole in the relevant channel with the aid of Cauchy's theorem. Our results are presented in Table II and Fig. 4 (right). The latter shows 200 pole locations for each of the five models, as determined by repeating the fits over 200 sets of pseudodata.
We have remarked that the Cornell group has suggested that the enhancement at 3.9 GeV is due to a combination of the D Ã D threshold and a node in the decay amplitudes of the ψð4040Þ [7]. We have tested this by setting D Ã D couplings to zero (without refitting) and find surprisingly little effect on the fit, thus disfavoring the Cornell explanation. Alternatively, Nakamura et al. [14] observe that the enhancement occurs due to a combination of interferences between the ψð3770Þ and the D Ã D threshold and the ψð4040Þ and nonresonant effects. In view of this, we have probed the mechanism further by setting the ψð4040Þ continuum couplings to zero, and find that the enhancement remains. Thus, it appears that the K-matrix method is general enough to create enhancements via multiple mechanisms. In fact, we find that interference between resonances and the background or between different continuum channels can also lead to a similar effect. We therefore conclude that establishing the specific cause of the enhancement is likely to depend on the details of the fit model, although we can confidently conclude that a Gð3900Þ resonance is not required.
We have also examined the claim that a node in the decay amplitudes of ψð4040Þ can affect open-charm production in the mass range considered here. The good fit quality of models 1-4 indicates that such a node is not required.
We nevertheless consider a fifth model that replicates model 1 with modified resonance couplings, as shown in Eq. ( 4). We find that the resulting fit quality is somewhat improved, with χ 2 =ndf decreasing from 2.20 to 1.88. Thus, we find weak evidence in favor of a node in the decay amplitudes of ψð4040Þ. 1 We find it significant that our bestfit node parameter is k 0 ¼ 620 MeV, remarkably close to the value of 750 MeV anticipated by the Cornell group [7].
We have performed a coupled-channel K-matrix fit to the D 0 D0 , D þ D -, D Ã Dþ c.c., D Ã DÃ , and inclusive hiddencharm cross sections from e þ e -for 3.7-4.2 GeV, and find a good fit in all channels for a variety of models. The fit only assumes two bare poles, corresponding to ψð3770Þ and ψð4040Þ. The excellent quality of all the model fits is a strong indication that an additional bare pole near 3900 MeV is not needed to explain the data, in keeping with an old prediction from the Cornell group [7]. This observation highlights the need for coupled-channel analyses when interpreting and fitting data in regions where thresholds are expected.
Although a node in the decay amplitude for ψð4040Þ → D D is expected on general grounds, we have found that it is not required to explain the enhancement at 3900 MeV likely because the D Ã D threshold enhancement in combination with hadronic couplings that decrease over typical hadronic scales is sufficient to create the observed line shape. Nevertheless, a model with ψð4040Þ couplings including a node does produce a fit with a somewhat better χ 2 =ndf, providing weak evidence for this mechanism and its effects on hadronic observables. Left: e þ e -→ hadrons as fit using model 1. Data is from Refs. [30] (BES, black closed), [31] (BESII, black open), [32] (SPEAR, orange), [33] (SPEAR, green), and [34] (SPEAR, light blue). Right: pole locations for model 1 (red), model 2 (blue), model 3 (green), model 4 (yellow), as well as model 1 including a node (model 5, pink). Each point corresponds to a pole location in a fit to one of 200 sets of pseudodata.
A node in the amplitude for ψð4040Þ → D D implies similar nodes in ψð4040Þ → D ðÃÞ DðÃÞ under very general assumptions [35].
Resonance parameters are reported in Table II for the different models. The rows labeled "fit summary" are provided for convenience, as a summary of this work. It would be inappropriate to perform a statistical model averaging of our results, so we simply report our preferred result (model 1) with an additional "systematic uncertainty." We find that the hadronic and electronic widths of ψð3770Þ are in most cases very close to those reported by the PDG. It is interesting, however, that our preferred mass is 5 MeV above the PDG value-a significant difference, in view of the small uncertainties quoted. Indeed, we have found that forcing the nonresonant couplings to be equal (models 2 and 4) raises the ψð3770Þ pole location substantially.
Similarly, the preferred ψð4040Þ mass and hadronic widths lie close to PDG values; however, the extracted electronic partial width (for models 1, 3, and 4) is substantially smaller. This discrepancy is not due to a statistical fluctuation; rather, the PDG value is obtained from analyses that assume all of the structure in the total hadronic cross section to arise from resonances (e.g., from a model consisting of a sum of Breit-Wigner amplitudes). Such a choice, which ignores nonresonant contributions and rescattering, necessarily overestimates the electronic width. In our model, substantial structure is generated by nonperturbative effects and by nonresonant scattering, naturally leading to a smaller residue and partial width. Thus, we recommend that Γ ee ½ψð4040Þ be updated to a smaller value, as obtained here. We remark that similar problems affect bottomonium vectors, as reported in Ref. [16].
We note that models 2 and 5 yield electronic widths that greatly deviate from the other models. This discrepancy appears to be related to bootstrap fits in which the two bare pole positions that ultimately produce ψð3770Þ and ψð4040Þ unexpectedly move close to each other, or in which the scattering couplings become anomalously large (perhaps indicating that the models are sufficiently complex as to overfit the data).
In summary, a unitary coupled-channel analysis of charm production near 4.0 GeV finds strong evidence for the wellestablished ψð3770Þ and ψð4040Þ resonances, with no need for additional poles. The extracted resonance parameters are in broad agreement with PDG values; however, we find a slightly higher ψð3770Þ pole mass at 3778.7(7)(50) MeV, and a lower electronic width for the ψð4040Þ of 180(100) (170) eV. Our results are in agreement with the predictions of the Cornell group, and affirm the importance of performing coupled-channel analyses when fitting data with overlapping resonances or nearby thresholds.
Details and numerical results for the new BESIII measurement of the observed, dressed, and Born cross sections for e þ e -→ D D are presented here.
We remark that statistical uncertainties in the dressed and Born cross sections are not independent. To account for this issue, 10,000 pseudo cross sections were generated at each energy point. An iterative method (Ref. [27]) was then employed to extract the dressed and Born cross sections for each of the pseudodata sets. The resulting average cross section and the root of the diagonal elements of the covariance matrix are reported in the tables as the central value and the uncertainty for each entry. In practice, we found it unnecessary to include the full covariance matrix in the fits, since off-diagonal correlations are negligible.