Vector and scalar charmonium resonances with lattice QCD

We perform an exploratory lattice QCD simulation of $D \bar D$ scattering, aimed at determining the masses as well as the decay widths of charmonium resonances above open charm threshold. Neglecting coupling to other channels, the resulting phase shift for $D \bar D$ scattering in p-wave yields the well-known vector resonance $\psi(3770)$. For $m_\pi = 156$ MeV, the extracted resonance mass and the decay width agree with experiment within large statistical uncertainty. The scalar charmonium resonances present a puzzle, since only the ground state $\chi_{c0}(1P)$ is well understood, while there is no commonly accepted candidate for its first excitation. We simulate $D \bar D$ scattering in s-wave in order to shed light on this puzzle. The resulting phase shift supports the existence of a yet-unobserved narrow resonance with a mass slightly below 4 GeV. A scenario with this narrow resonance and a pole at $\chi_{c0}(1P)$ agrees with the energy-dependence of our phase shift. Further lattice QCD simulations and experimental efforts are needed to resolve the puzzle of the excited scalar charmonia.


Introduction
Charmonia cc well below open-charm threshold DD are among the best understood hadrons.Their spectra and selected transition matrix elements are successfully described by lattice QCD simulations and QCD motivated models.Recent lattice calculations have performed the necessary extrapolations and considered spectra [1,2] as well as certain radiative transitions [3,4].For states well-below open charm threshold, the main remaining uncertainty is the neglect of charm-annihilation Wick contractions in lattice simulations.
The most interesting charmonium and charmonium-like states lie near or above open charm thresholds.During the past decade a plethora of states that can likely not be interpreted as conventional cc have been discovered in experiment (for a review see for example [5,6]).These states have been treated theoretically making simplifying assumptions and reliable quantitative results for those hadrons are not available.In particular, all of the lattice simulations so far have ignored the strong decay of the charmonium resonances to a pair charmed mesons cc → D( * ) D ( * ) , which typically represents the main decay mode.Except in a few simulations [7][8][9][10], the effect of the threshold on the near-threshold states has been neglected.The most extensive spectrum of charmonia has been obtained in simulations with N f = 2 + 1 dynamical flavors at m π ≃ 400 MeV [11], but the determination neglects the unstable nature of the states and relies on extracting the energy levels with only quark-antiquark interpolating fields, which may lead to unphysical results close to multi-hadron thresholds [12][13][14].
Here we present a lattice QCD simulation of the vector (J P C = 1 −− ) and scalar (0 ++ ) charmonium resonances above DD threshold, taking into account their strong decay to DD for the first time.The lowest vector resonance above open charm threshold, the ψ(3770) is well established in experiment [15], and we extract its width by simulating DD scattering in p-wave.In contrast to that, the experimental and theoretical status of scalar charmonia is puzzling: the only well-established state is the ground state χ c0 (1P ), while there is no commonly accepted candidate for its first excitation χ c0 (2P ).We present a study of DD scattering in s-wave, aiming to address this open problem.We also consider possible effects of the DD threshold on the vector ψ(2S) and scalar χ c0 charmonia, which lie below threshold.
2 Open questions for charmonia of interest

Vector charmonia
The ψ(3770) with M = 3773.15± 0.33 MeV and Γ = 27.2 ± 1.0 MeV is located only ≃ 45 MeV above DD threshold [15,16].We focus here on its dominant decay mode ψ(3770) → DD in p-wave with branching fraction 0.93 +8 −9 [15].It is a well-established experimental resonance and is generally accepted to be predominantly the conventional 2s+1 nL J = 3 1D 1 cc state [17][18][19].There is an ongoing discrepancy between results from BES-II [20] and Cleo [21] regarding the non-DD part of the branching fraction which may be connected to neglecting interference effects in the BES-II analysis.Significant non-DD decays into light hadrons can occur if there is non-negligible mixing with the ψ(2S) [22].For our analysis we neglect disconnected contributions that would cause decay into light hadrons and treat the decay into DD as elastic, neglecting the decays into J/ψ ππ and J/ψ η that have tiny branching fractions [15].Our aim is to perform the first determination of the ψ(3770) resonance mass and ψ(3770) → DD decay width using a lattice simulation for DD scattering in p-wave.
We also investigate whether the DD threshold has any effect on ψ(2S), which is the first radial excitation of J/ψ and is situated ≃ 42 MeV below threshold.Such a possibility was discussed in relation to the Fermilab-MILC preliminary results [23] where a simple analysis of the spin-averaged 2S state appeared high with respect to experiment, although large systematic uncertainties related to excited state contaminations were observed.A more recent HPQCD study [2] finds no significant discrepancy.The mixing of the vector charmonia with a pair of two charmed mesons was first simulated in [7], where only D 1 D in s-wave was considered and the width of ψ(3770) was not extracted.

Scalar charmonia
The only well established scalar charmonium state is the ground state χ c0 (1P ), interpreted as the 3 1P 0 cc and located well below the open charm threshold.A further known resonance, the X(3915) with Γ = 20 ± 5 MeV is seen only in J/ψ ω and γγ decay channels [15].BaBar has determined its J P quantum numbers to be 0 + [24] which would only allow J P C = 0 ++ .This spin-parity determination by BaBar assumes that a J P = 2 + resonance would be produced in the helicity 2 state, which might not be justified for an exotic meson [5].As a consequence, the PDG recently assigned X(3915) to be χ c0 (2P ) [15], but a number of convincing reasons given by Guo & Meissner [25] and Olsen [26] raise serious doubts about this assignment: • The dominant decay mode of scalar charmonium above open charm threshold is expected to be a "fall-apart" mode into DD that would lead to a relatively broad resonance.In particular the width into DD is expected to be much larger than for the well-established χ c2 (2P ) [15], which decays to DD in d-wave.Yet m D D invariant mass spectra of several experiments show no evidence for X(3915) → D D. This also indicates that the DD width extracted from the present lattice simulation cannot be compared to X(3915).
• The spin-splitting m χ c2 (2P ) − m χ c0 (2P ) within this assignment seems too small compared to • The partial width for the OZI suppressed X(3915) → ωJ/ψ seems too large [25], which is translated to two contradicting limits for this decay in [26].
The intriguing χ c0 (2P ) was related to the broad structures in DD invariant mass in the same references [25,26].The process γγ → DD from BaBar [27] and Belle [28] leads Guo&Meissner to1  Obviously the spectrum of scalar charmonia beyond the ground state presents an open question.Our aim is to shed light on this issue by simulating DD scattering in s-wave on the lattice, and look for possible resonances in the extracted scattering matrix.Preliminary results based on the same simulation and only one ensemble have been presented in Ref. [30].

Lattice setup and charm-quark treatment
The simulation is performed on two lattice ensembles with the parameters listed in Table 1.Both ensembles have rather low m π L but this is not a serious issue for charmonia and Ensemble (1) Ensemble ( 2)  1) are from [31,32].Those of ensemble (2) are provided by the PACS-CS collaboration [33].N L and N T denote the number of lattice points in spatial and time directions, N f the number of dynamical flavors and a the lattice spacing.
DD scattering in this simulation, where pions do not enter explicitly.Further details about the ensembles and our implementation of charm quarks may be found in [12,14,31,32] for ensemble (1) and in [33,34] for ensemble (2).
To minimize heavy-quark discretization effects at finite lattice spacing the Fermilab method [35,36] is used for the charm quarks.The corresponding dispersion relation [37] for a meson M containing charm quarks is where p = 2π L q and q ∈ N 3 .On both ensembles the charm quark hopping parameter κ c is tuned [14,34]  which is the relevant reference mass for our spectra of charmonium.The M 1,2,4 for the spinaveraged charmonium were determined based on the lattice data from the lattice dispersion relation (3.1), setting W 4 to zero.Then κ c was fixed by tuning the kinetic mass M 2 to mexp .The corresponding values for the spin-averaged M 1 are given in Table 2.
To investigate the DD scattering we need the dispersion relation E D (p) for D mesons, which is also given by Eq. (3.1) with parameters M 1,2,4 in Table 2.The common feature of spectra in the scalar and vector charmonium channel are two-particle states DD that have a discrete spectrum on the finite lattice.In the absence of interactions, D(q) D(−q) have energies according to (3.1) which will be shifted due to the interaction.Within the Fermilab approach, the rest masses have large discretization effects but mass differences are expected to be close to physical [38,38]   experiment.In order to compare the splitting E lat − mlat with E exp − mexp , we will sometimes plot and compare it with E exp .An important quantity is the position of the DD threshold with respect to our reference mass.The splitting 2m D − m for ensemble (2) is very close to the experimental value 2m exp D − mexp ≃ 0.666 GeV, while it is a bit larger for ensemble (1) due to the heavier pion mass and larger discretization effects (see Table 2).
Our charm quark treatment has been verified on ensemble (1) for low-lying charmonia, D meson resonances [14] and D s mesons [34,39], where reasonable agreement with experiment was found.The spectrum for D s mesons and some other hadrons containing charm quarks were also determined on ensemble (2) [34,39] with even better agreement due to the lower pion mass and smaller discretization effects.

Analysis details
Interpolating fields O are used to create and annihilate the physical system with J P C = 1 −− or 0 ++ , isospin I = 0 and total momentum zero.All quark fields in the interpolators are smeared according to the distillation method q ≡ Nv k=1 v (k) v (k) † q point [40,41].We use N v = 192 eigenvectors of the lattice laplacian v (k) for ensemble (2) and N v = 96 or 64 for ensemble (1).The distillation method is convenient for calculating a variety of Wick contractions.The full distillation method [40] is employed on ensemble (1) with a smaller volume and details of the implementation are given in [12,14].The stochastic version [41] is used on ensemble (2) with larger volume and details of our implementation are provided in [34].
where i denotes polarization, while Q ijk and R ijk = R jik are listed in [42].The cc interpolators O cc 1−14 for vector channel T −− 1 are listed in Table X of [14].The momentum is projected for each D meson separately, so that the O DD couple to p-wave.For ensemble (1) N v = 64 is used for O DD 2 , and N v = 96 for the remaining interpolators.
The irreducible representation T −− 1 contains J P C = 1 −− states of interest, and also ψ 3 states with J P C = 3 −− coupling due to the broken rotational symmetry on the lattice.In the continuum limit, O cc 1−14 contain only 1 −− , while O cc 15,16 contain 1 −− and 3 −− [42], which will help us to identify the spin 3 admixture related to ψ 3 .

Scalar channel
DD in s-wave and J/ψ ω are the dominant two-meson states in the energy region of interest E ≤ 4 GeV.Seven cc, four DD, and two J/ψ ω interpolating fields are used in the relevant irreducible representation A ++ 1 : O cc 1−7 are listed in Table X of [14].The momenta are projected for each meson separately in O DD and O J/ψ ω .For ensemble (1) , and N v = 96 for the remaining interpolators.
The irreducible representation A ++ 1 contains J P C = 0 ++ states of interest, and in general also states with J ≥ 4, which appear at energies beyond our interest.
The interpolator O DD 4 is not used for ensemble (1) since D(2)D(−2) appears above 4 GeV.The O J/ψ ω are not used on ensemble (2) since the results from ensemble (1) indicate that J/ψ ω is almost decoupled from the rest of the system.We omit contractions where the charm quark annihilates.A red solid line represents a c quark, while the black dashed line represents a u or d quark.

Towards the spectrum
The correlation matrix contains the information on energies E n and the overlaps Z n j ≡ Ω|O j |n .We evaluate all Wick contractions for O ≃ cc, (qc)(cq), (cc)(qq) (4.1,4.3)shown in Fig. 1.We omit Wick contractions where charm quark annihilates as in almost all previous lattice simulation of charmonia; these induce mixing with I = 0 decay channels containing only light quarks u, d, s, they are Okubo-Zweig-Iizuka suppressed and present a challenge for current lattice simulations.It is noteworthy that these decays might be important to clarify the experiment puzzle with regard to non-DD hadronic decays [18,22].
The energies and overlaps are extracted from the correlation matrix using the generalized eigenvalue method [43][44][45][46] where λ (n) (t) ∝ e −Ent at large t.Correlated two or one-exponential fits to λ (n) (t) are used and t 0 = 2, 3.The errors-bars correspond to statistical errors obtained using singleelimination jack-knife.3, while (b,d) utilize just O cc from the same sets.
5 Results for the vector channel

Discrete spectrum
The energy levels in the vector channel are shown in Fig. 2a and 2c together with the experimental masses.The full set of operators gave noisier signals than suitable subsets, and the chosen subsets are listed in Table 3.The circles denote the energy levels that are related to J P C = 1 −− states J/ψ, ψ(2S), ψ(3770), D(1)D(−1) (from bottom to top), while D(0) D(0) does not appear for p-wave.The diamond indicates a level related to the J P C = 3 −− state ψ 3 , that is present in representation T −− 1 due to the broken rotational symmetry on the lattice.
The highest state (n = 5) has largest overlap with O DD and disappears when these interpolators are excluded from the basis, as shown in Fig. 2b and 2d.Each energy level in addition to D(1) D(−1) indicates the presence of a bound state or a resonance.Good resemblance with the experimental spectrum is indeed confirmed in Fig. 2. The J/ψ is significantly below threshold and no effect from threshold is expected.The ψ(2S) is situated   Figure 3.The overlaps Z n j = Ω|O j |n for the vector channel show the matrix elements of interpolators O j between the vacuum Ω| and the physical eigenstate |n on the lattice.We present the overlap ratios Z n j /max m Z m j on ensemble (1) (top) and on ensemble (2) (bottom).The denominator is the maximal |Z m j | at given operator number j.These ratios are independent on the normalization of the interpolators O j .Levels n = 1, .., 5 are ordered from lowest to highest E n in Figs.2a and 2c for both ensembles, respectively.The order of interpolators j on the abscissa (listed in caption of Fig. 3) is the same as in the list (4.1).
≃ 42 MeV below threshold in experiment, and the corresponding finite volume energy on the lattice does not depend (within uncertainties) on whether DD interpolators are used or not (see Fig. 2).The appearance of levels n = 3 and 4 is related to the ψ(3770) resonance and to the spin 3 admixture and the corresponding ψ 3 resonance.Level n = 4 is related to ψ 3 due to smaller overlaps O  in the continuum limit only to 1 −− (which is responsible for small O cc 1−14 |ψ 3 at finite a), while O cc 15,16 couple to 1 −− and 3 −− .Further support is given by the near-degeneracy with the energies from the irreducible representation A −− 2 where a 3 −− state comes as the ground state (see Table 4).For the ψ(3770) the avoided-level crossing scenario suggests E 3 in the energy region m ± Γ, which is reasonably satisfied by comparing to experiment.In order to really determine the resonance mass and width for ψ(3770) one needs to consider the phase shifts for DD scattering in p-wave.

DD scattering in p-wave
We assume that DD scattering in p-wave near the resonance ψ(3770) is elastic, which is a good approximation since Br[ψ(3770) → D D] = 93 ± 9%, while the remaining part goes mainly to J/ψππ.In the elastic case, the scattering phase shift δ is given by Lüscher's relation [47,48] which applies for the total momentum zero employed in our case.The momentum p of D mesons is extracted from the measured energy levels E lat n = 2E D (p) using the dispersion relation (3.1).The resulting momenta and phase shifts for all eigenstates except for the spin 3 admixture and for the finite volume state related to J/ψ are collected in Table 3.The large absolute value of p 3 cot δ corresponds to feeble scattering, while small p 3 cot δ is related to significant scattering.
We fit our data in two ways: (i) A resonance ψ(3770): The scattering matrix in the vicinity of a resonance has a Breit-Wigner form The width is parametrized in terms of the phase space for p-wave decay and the ψ(3770) → D D coupling g.It is expected that g depends only mildly on m u/d , and that the leading dependence of Γ on m u/d is captured by phase space.Equations (5.2, 5.3) lead to Ensemble ( 1)  where p R is the D meson momentum at the resonance peak.The values of g and p R follow from the linear fit (5.4) through the energy levels n = 3, 5 in the vicinity of the resonance, where the Breit-Wigner form applies (level 4 is omitted since it is attributed to ψ 3 as discussed above).The fit is shown in Fig. 4, while the resulting resonance parameters are given in Table 5.The resonance mass m R corresponding to the p R on the lattice is given by inserting the Fermilab dispersion relation (ii) A resonance ψ(3770) and a bound state ψ(2S): In addition to the Breit-Wigner form (5.4), which is linear in p 2 , we make use also of the square form in p 2 which in general has a longer range of applicability.It aims to capture also the DD scattering in the vicinity of ψ(2S): there the (imaginary) phase shift in Table 3  The results from both fit (i) and fit (ii) on ensemble (2) are compatible with the experimental data 3 within large statistical uncertainties (see Table 5).Note that the higher-lying ψ(4040) resonance does not influence the results for this ensemble, since it lies significantly higher than the relevant energy levels.
Confronting the results for ψ(3770) from fit (i) on ensemble (1) with experiment gives a smaller resonance momentum p R than in experiment, which we attribute to the unphysical threshold on ensemble (1) at m π ≃ 266 MeV and the finite lattice spacing.The resonance mass m R calculated as in Eq. 3.4 compares favorably.The coupling constant from fit (i) is to small compared to experiment which which is likely related to to the closeness of the ψ(4040) resonance neglected in the analysis.The assumption that the resonance ψ(4040) does not affect the energy level related to D(1) D(−1) is probably not justified on ensemble (1), where energy level lies higher (and closer to ψ(4040)) than on ensemble (2).Roughly estimating the effect by comparing the one-resonance and two-resonance scenarios, estimating g and p R for ψ(3770) and ψ(4040) from available experimental data [49], the coupling we observe is consistent with this interpretation 4 .Given the possibly large influence from the ψ(4040) we can not conclude that fit (ii) is better than fit (i) on this ensemble.
The resulting 1 −− spectrum is summarized and compared to experiment in Fig. 5.
3 Since we work in the isospin-symmetric limit we measure the sum of the neutral and charged decay modes; therefore we compare to the experimental value [15].Notice also that averaging the results from recent experiment resonance mass determinations for the ψ(3770) leads to a value of m exp R = 3778.1(1.2),much larger than the fit by the PDG (which relies on an experiment neglecting interference with non-resonant background) and consistent with the most recent results in [16]. 4The maximal effect of ψ( 4040 6 Results for the scalar channel

Discrete spectrum
The energy levels in the scalar channel are shown in Figs. 6.The only experimentally well established state is χ c0 (1P ).The triangles represent the intriguing experimental candidates for χ c0 (2P ), none of which is commonly accepted (see Section 2).The spectrum from a lattice simulation consists both of energy levels that have large overlap with qq operators as well as energy levels with dominant overlap to DD operators.The latter appear near their non-interacting energies E n.i.DD of Eq. (3.3), which are denoted by dashed lines in Figs.6a,d.On ensemble (1) levels n = 2, 4 appear near the noninteracting D(0) D(0) and D(1) D(−1) (cf.Fig. 6a).Levels n = 2, 3, 4 on ensemble (2) have dominant overlap to DD scattering operators and are close to non-interacting D(0) D(0), D(1) D(−1) and D(2) D(−2) energies (cf.Fig. 6d).
Each energy level in addition to the number of expected D(q) D(−q) scattering levels is related to the presence of a bound state or a resonance.There are two such states, that cannot be attributed to D(q) D(−q) for both ensembles.The ground state is related to χ c0 (1P ) and is close to its experimental mass.The second of these two levels appears above threshold and corresponds to n = 3 for ensemble (1) and n = 5 for ensemble (2), as shown in Figs.6a and 6d.The avoided level crossing scenario suggests that an additional level appears somewhere in the range E ≃ m ± Γ, which suggests the existence of a resonance suggested in [25,26].The dashed lines shown energies of non-interacting D(q) D(−q) with q = 0, 1, , while dot-dashed line represents m J/ψ + m ω .Interpolators used in (a,d) are given in Table 6, ( in addition, while (c,e) are based only on O cc 1, roughly at m ≃ 3.9 − 4.0 GeV (naive estimate from E n ) .( This is close to the first excitation obtained using just O cc interpolators in Figs.6c and 6e.Such a basis gives a rough estimate of resonance masses but is not well suited to capture two-particle states or resonances and bound states close to threshold [12][13][14]. The spectrum including J/ψ ω interpolating fields is shown in Fig. 6b for ensemble (1).An energy level related to J/ψ(0) ω(0) appears at roughly m J/ψ + m ω while the energies of all the other levels remain unaffected with respect to Fig. 6a.We have verified also that the overlaps for the remaining levels are not affected if O J/ψ ω are in the basis or not.This indicates that the J/ψ ω channel is decoupled from DD channel to a good approximation.

DD scattering in s-wave
We now study DD scattering assuming that J/ψ ω channel is decoupled as argued in the previous paragraph.The energy shifts of the extracted E lat n with respect to E n.i.DD give the size of the s-wave scattering phase shift δ according to (5.1).On ensemble (1) we observe statistically significant energy shifts with respect to the dashed lines in Fig. 6a.The energies yield D-meson momenta p via E lat n = 2E D (p), and the corresponding phase shifts   δ(p) via Eq.(5.1).These are provided for all levels in Table 6 and plotted in Figs.7 and  8.
The uncertainties on the energies E n=2,3,4 are rather large for ensemble (2) and they are within errors compatible with non-interacting energies E n.i.DD (3.3).This implies that we are not able to reliably determine the energy shifts, and the resulting errors on the scattering matrix will be large, as illustrated in Fig. 7.If E lat n ≃ E n.i.DD within errors, this implies δ ≃ 0 modulo π and cot δ ≃ ±∞ within errors.The extracted p cot δ from n = 2, 3, 4 have large errors, which allow almost all p cot δ expect for small |p cot δ|.For n = 2, 3, 4 we plot central values f (p 2 ) with the ranges [f (p 2 − σ p 2 ), f (p 2 + σ p 2 )] where f = p cot δ or p cot δ/ √ s.The error for all other levels (on both ensembles and both channels) is the usual jack-knife.
The resulting p cot δ/ √ s for ensemble (1) has a puzzling behavior and we are going to confront it with various hypothesis, collected in Fig. 8.The errors of p cot δ/ √ s on ensemble (2) are large and do not allow reliable fits.We will still compare data from ensemble (2) with fits that are based on ensemble (1), and plot them as function of p 2 on the same figure.Note that the slightly different positions of thresholds 2m D on the two ensembles do not play a crucial role in the scalar channel, since the anticipated resonances and bound states are not expected very close to the threshold (in contrast to ψ(3770) in the vector channel).
(i) A narrow resonance: In the vicinity of a Breit-Wigner resonance (5.2) one expects and the zero gives the position of the resonance.The upper three points in Fig. 8a however do not fall onto one line, so our results cannot be reconciled with a single Breit-Wigner resonance in the region between 2m D and 4 GeV.The highest two points support the existence of a narrow resonance between them and a linear fit (6.2) over two levels shown in Fig. 8(i,a) renders (5.5) m R = 3.966(20) GeV, g = 1.26 (18) GeV, ( p lat R = 0.614 (33) GeV, Γ lat = 62 (17) MeV .
The scattering data from ensemble (1) therefore suggest the existence of a yet unobserved scalar state called χ ′ c0 .Note that its mass is within the range of the naive estimate (6.1).It has a width Γ lat in our simulation, while the corresponding width g 2 p/m 2 R in experiment would be modified due to a different phase space via p = [m 2 R /4 − (m exp D ) 2 ] 1/2 , leading to Γ exp predict = 67 (18) MeV .(6.4) We have assumed that g and m R do not depend on the pion mass here.It is unlikely that this state corresponds to X(3915) since the DD decay channel was not observed for this state.A narrow resonance is roughly consistent also with the result from ensemble (2) within huge errors (see Fig. 8(i,b)), however there must be some additional interaction between D and D near the threshold according to ensemble (1).
(ii) A narrow resonance and a bound state χ c0 (1P): Our next hypothesis assumes that χ c0 (1P ) represents a pole in DD scattering on the first Riemann sheet, leading to p cot δ ≃ i|p B |i = −|p B | at the position of the bound state.The negative value of p cot δ below threshold might be a possible reason why p cot δ at threshold is smaller than expected based on narrow resonance (6.2,6.3).In this case the value of p cot δ at threshold is influenced by the resonance and a bound state.To investigate this situation, we attempted several fits over all four levels on ensemble (1).A form that is motivated by our data is presented in Fig. 8(ii,a), where A = 0.13 (15), B = 0.66(18)/GeV The bound state is attributed to χ c0 (1P ) and its mass is very close to the one obtained from the ground state energy.The hypothesis also supports a narrow resonance at p R = 0.668 (35) GeV where function (6.5) crosses zero, and m R = 4.002 (24) GeV, g = 0.85(65) GeV, ( Γ lat = 30 (45) MeV , Γ exp predict = 32 (48) MeV .
This is roughly consistent with the χ ′ c0 in (6.3).This hypothesis based on ensemble (1) is consistent also with the result from ensemble (2) within huge errors in Fig. 8(ii,b).An interesting feature of this hypothesis is the large p cot δ or equivalently small cross-section at p 2 ≃ D, which corresponds to √ s ≃ 4.0 GeV.This feature seems to be present also in the experimental data from Belle [26] where a dip seems to appear at similar invariant mass.
(iii) A broad resonance: The broad resonances (2.1,2.2) proposed by Meissner&Guo [25] or Olsen [26] are confronted with our lattice data in Fig. 8(iii).This shows a Breit-Wigner shape (6.2) with p R and g extracted from the experimental data (2.1,2.2).Although they are roughly compatible with our scattering results near threshold, they cannot be reconciled with it in the region above threshold where our data indicates either a much narrower resonance or a more complicated situation.
(iv) Two resonances: Since neither one narrow or one broad resonance describe our scattering data near and above threshold, we next try an hypothesis with two elastic resonances p cot δ(s) With this parametrisation there are two resonance poles in the scattering amplitude, separated by a zero.Figure 8(iv,a) shows an example with g A = 2.1 GeV, p R A = 0.23 GeV, g B = 1.0 GeV and p R B = 0.57 GeV that is consistent with the upper three scattering points for ensemble (1) 5 .This hypothesis however predicts another energy level near p 2 ≃ 0.1 GeV 2 where the model (6.8) crosses with the Lüscher curve.Another energy level is expected in the two-resonance scenario also according to naive reasoning that each resonance or bound state leads to a level in addition to DD.Such an additional energy level at p 2 ≃ 0.1 GeV 2 is not observed in ensemble (1) indicating that this hypothesis is not supported by our data.An analogous conclusion is reached when confronting this hypothesis with the data from ensemble (2): the hypothesis predicts five energy levels in the region p 2 = [−0.1,0.5] GeV and we observe four levels only.

Conclusions and outlook
We performed a lattice QCD simulation of DD scattering in s-wave and p-wave to study vector and scalar charmonium resonances on two rather different ensembles.This is a first simulation aimed at determining the strong decay width of charmonium resonances above open charm threshold.Ensemble (1) has N f = 2 and m π = 266 MeV, while ensemble (2) has N f = 2 + 1 and m π = 156 MeV.Several cc and D D interpolating fields were used in both channels, where the (stochastic) distillation method was used to evaluate the Wick contractions.
In the vector channel, the well known ψ(3770) resonance is present just above DD threshold.The phase shift for D D scattering in p-wave was determined using the Lüscher formalism.We performed a Breit-Wigner fit in vicinity of the ψ(3770) to obtain its resonance mass at 3.784 (7) (10) GeV and 3.786(56) (10) GeV for ensembles (1) and ( 2), respectively.Our determination of its decay width might be affected by the Ψ(4040) on ensemble (1).Ensemble (2) does not suffer from this issue, and the determination of the resonance parameters is more reliable, but its statistical accuracy is poor.The resulting spectrum in the vector channel, including also J/ψ and ψ(2S), is compared to experiment in Figure 5.This work presents a first step towards a determination of the ψ(3770) resonance parameters from lattice QCD.Future lattice studies in multiple volumes (and or using multiple momentum frames) with a larger ensemble size should provide precision results for this resonance.
In the scalar channel, only the ground state χ c0 (1P ) is understood and there is no commonly accepted candidate for its first excitation χ c0 (2P ).Guo & Meissner [25] as well as Olsen [26] argued that the higher lying X(3915) can probably not be identified with the χ c0 (2P ).They suggest that a broad structure observed in D D invariant mass represents χ c0 (2P ).This posed a particular motivation to extract the phase shift for D D scattering in s-wave in the present work.The resulting scattering data on the ensemble with m π = 156 MeV is unfortunately noisy.The simulation at m π = 266 MeV renders the scattering phase shift only at a few values of the D D invariant mass, which also does not allow a clear answer to the puzzles in this channel.We obtain the χ c0 (1P ) and our data provides an indication for a yet-unobserved narrow resonance slightly below 4 GeV with Γ[χ ′ c0 → D D] below 100 MeV.A scenario with this narrow resonance and a pole in the D D scattering matrix at χ c0 (1P ) agrees with the energy-dependence of our phase shift.We checked the consistency of other three scenarios with our data: just one narrow resonance, just one broad resonance proposed in Guo & Meissner [25] and Olsen [26], or one narrow and one broad resonance.None of these three scenarios agree with our data in the whole energy region we probed.For the scalar channel this leaves us with a situation where puzzles remain, both from theory and experiment.To clarify the situation, further experimental and lattice QCD efforts are required to map out the s-wave D D scattering in more detail.

4. 1
Vector channel DD in p-wave is the dominant two-meson contribution for E ≤ 4 GeV, while D 1 D appears higher.Sixteen cc and two DD interpolating fields are used in the relevant irreducible representation T −− 1 :

Figure 1 .
Figure 1.Wick contractions computed for the correlation matrix (4.4) with interpolators (4.1,4.3).We omit contractions where the charm quark annihilates.A red solid line represents a c quark, while the black dashed line represents a u or d quark.

Figure 2 .
Figure 2. The energies E (see Eq. (3.4)) in the vector channel on both ensembles, together with the experimental masses.The circles represent J P C = 1 −− states, while the diamond represents a 3 −− admixture present in the irreducible representation T −− 1 and related to the ψ 3 .The dashed lines show the non-interacting energy of D(1) D(−1) (3.3), and the dotted line represents the threshold 2m D .The D(0) D(0) state does not appear for p-wave.Interpolators used in (a,c) are given in Table3, while (b,d) utilize just O cc from the same sets.

Figure 4 . p 3
Figure 4. p 3 cot δ/ √ s versus p 2 for DD scattering in p-wave in the region of the ψ(2S) bound state and the ψ(3770) resonance.The p denotes the momentum of D meson.We show the Breit-Wigner fit (i) and the extended fit (ii), which aims to capture also the behavior around ψ(2S).

Figure 5 .
Figure 5.The comparison of the final 1 −− spectrum to the experiment.The magenta diamond denotes ψ(3770) resonance mass from the Breit-Wigner fit (i) or extended fit (ii), given in Eqs.(5.4) and (5.6), respectively.The magenta triangle denotes ψ(2S) obtained as a pole in D D channel.The blue triangles denote masses of J/ψ and ψ(2S) extracted as energy levels in the finite box.The statistical and scale setting errors have been summed in quadrature.

Figure 7 .
Figure 7.We show p cot δ and p cot δ/ √ s versus p 2 for DD scattering in s-wave, where p denotes the momentum of the D meson.The circles with (sizable) errors denote the lattice data, while the solid lines show p cot δ = 2Z 00 /(L √ π) according to Lüscher's relation (5.1).When the momentum is compatible with the non-interacting momentum p = 2πq/L (q ∈ N 3 ), one has δ = 0 and | cot δ| = ∞, which is responsible for the huge errors on p cot δ on ensemble (2).

Figure 8 .
Figure 8.The p cot δ/ √ s versus p 2 for DD scattering in s-wave, where p denotes the momentum of the D meson.The lattice data (blue and red circles) is confronted with p cot δ/ √ s based on various hypothesis (dashed lines) described in Section 6 of the main text.The thin dot-dashed lines in the plots at the bottom denote p cot δ = 2Z 00 /(L √ π) (5.1).

Table 1 .
The gauge configurations of ensemble (

Table 2 .
The parameters in the dispersion relation(3.1)forDmesons and spin-averaged charmonium14 (m ηc + 3m J/ψ ).The last line is in GeV, others in lattice units.

Table 4 .
The energy of ψ 3 with J P C = 3 −− from the ground state in the A −−

Table 5 .
[15]meters of the resonance ψ(3770) and bound state ψ(2S) from fits (i)(5.4) and (ii)(5.6).The ψ(3770) → D D width Γ = g 2 p 3 /(6πs) is parametrized in terms of the coupling g and compared the value of the coupling derived from experiment[15].The p R denotes Dmeson momenta at the peak of the resonance and |p B | the binding momentum.The first errors are statistical and the second errors (where present) are from the scale setting uncertainty.The experimental data and errors are based on PDG values.Errors on experimental p R/B are suppressed as they are very small.
nearly satisfies the condition for the bound state cot δ ≃ i on the physical Riemann sheet p B = i|p B |, leading to p 3 cot δ ≃ |p B | 3 .The fit (5.6) through levels n = 2, 3, 5 in Fig. 4 assumes that the ψ(2S) state still affects the DD scattering.It renders (A, B, C) ≃ at this zero, lead to the parameters of ψ(3770) resonance in Table5.This model also leads to a bound state ψ(2S) at p B = i|p B | where the scattering amplitude T (5.2) has a pole and cot δ(p B ) = i.The bound state mass m B in Table5is indeed close to experimentally measured ψ(2S).
2, C = 0.028(63) GeV 2 and D = 0.513(77) GeV 2 are obtained from the fit.This hypothesis supports a bound state at p B = i|p B | which corresponds to a pole in T (5.2) or equivalently cot δ(p B ) = i,