Charmonium resonances with $J^{PC}=1^{--}$ and $3^{--}$ from $\bar DD$ scattering on the lattice

We present a lattice QCD study of charmonium resonances and bound states with $J^{PC}=1^{--}$ and $3^{--}$ near the open-charm threshold, taking into account their strong transitions to $\bar DD$. Vector charmonia are the most abundant in the experimentally established charmonium spectrum, while recently LHCb reported also the first discovery of a charmonium with likely spin three. The $\bar DD$ scattering amplitudes for partial waves $l=1$ and $l=3$ are extracted on the lattice by means of the L\"uscher formalism, using multiple volumes and inertial frames. Parameterizations of the scattering amplitudes provide masses and widths of the resonances, as well as the masses of bound states. CLS ensembles with 2+1 dynamical flavors of non-perturbatively $O(a)$ improved Wilson quarks are employed with $m_\pi\simeq 280$ MeV, a single lattice spacing of $a\simeq0.086$ fm and two lattice spatial extents of $L=24$ and $32$. Two values of the charm quark mass are considered to examine the influence of the position of the $\bar{D}D$ threshold on the hadron masses. For the lighter charm quark mass we find the vector resonance $\psi(3770)$ with mass $m=3780(7)$ MeV and coupling $g=16.0(^{+2.1}_{-0.2})$ (related to the width), both consistent with their experimental values. The vector $\psi(2S)$ appears as a bound state with $m=3666(10)$ MeV. The charmonium resonance with $J^{PC}=3^{--}$ is found at $m=3831(^{+10}_{-16})$ MeV, consistent with the $X(3842)$ recently discovered by LHCb. At our heavier charm-quark mass the $\psi(2S)$ as well as the $\psi(3770)$ are bound states and the $X(3842)$ remains a resonance. We stress that all quoted uncertainties are only statistical, while lattice spacing effects and the approach to the physical point still need to be explored. This study of conventional charmonia sets the stage for more challenging future studies of unconventional charmonium-like states.


I. INTRODUCTION
The charmonium spectrum determined from experiment challenges our theoretical understanding of the internal dynamics of mesons containing charm quarks. While charmonia below the open-charm decay threshold are found to fit within a quark-antiquark picture, the properties of several resonances above the threshold, generally referred to as XY Z states, call for exotic interpretations. Various suggestions have been made based on phenomenological approaches, effective field theories, potential models, etc. These include compact diquark-antidiquark states, mesonic molecules, hybrid mesons with gluonic excitations, or resonances generated through coupled channel scattering. The theoretical approaches are motivated by phenomenology and approximate certain regimes of the strong interaction. However, a complete description of the observed charmonium spectrum from Quantum ChromoDynamics (QCD) also requires complementary first principles nonperturbative calculations, which can be performed using lattice QCD. Ab-initio determinations of ground state charmonia, such as in Refs. [1][2][3][4][5][6][7][8], illustrate the power of these non-perturbative approaches in precisely determining masses and splittings. Lattice calculations of other conventional charmonium bound states and resonances aim to satisfactorily describe well understood excitations and demonstrate the efficacy of the methods employed. This gives confidence in investigations of unconventional charmonium-like states and of those for which the experimental situation is less clear.
To date, several lattice QCD calculations of the excited charmonium spectrum with isospin zero have been performed within the single hadron approach, where one employs only mesonic interpolators ofcc type and assumes that this qualitatively captures the single meson spectrum in a finite volume [9][10][11][12][13]. So far, only one lattice study has gone beyond the single-hadron approach, taking into account the strong decays of charmonium res-onances above the open-charm threshold. 1 In that study an exploratory determination of theDD scattering amplitude in the l = 1 and l = 0 partial waves was performed [16]. The low-lying vector and scalar charmonium spectra were calculated in a single lattice volume in the center of momentum frame (CMF). The authors extracted the resonance and bound state parameters in the vector channel, while the conclusions in the scalar channel indicated different possible interpretations and left several unresolved puzzles for both theory and experiment.
Note that all previous lattice studies of charmonium were limited to the CMF [9][10][11][12][13][16][17][18]. However, lattice calculations in the moving frames (i.e. with non-zero total momentum) can provide additional information on the relevant two-meson scattering amplitudes (for a review, see Ref. [19]). The challenge in spectroscopy using lattice techniques is that, due to the reduced symmetry on the lattice, several J can contribute (in the CMF) to any given lattice irrep. In the moving frames, the situation is much worse as the finite volume spectrum in any lattice irrep gets a contribution from states with different J P (in their rest frame), as they are not good quantum numbers already in the infinite volume continuum and also on the lattice. As the first step, we extracted the excited charmonium spectrum within the single hadron approach on the lattice in multiple inertial frames and identified the continuum quantum numbers [20]. Such a spin identified single-hadron spectrum tells us where to expect effects from resonances in different partial waves, which in turn provides valuable information for the parameterization of scattering amplitudes. Thus this study allows one to fully exploit the data in moving frames in the present and related future studies.
The current article presents a lattice QCD investigation of charmonium resonances taking into account their most important strong decay, which goes beyond the single hadron approach. We focus on the scattering of a pair ofD and D mesons in the l = 1 partial wave, which features charmonium resonances and bound states with J P C = 1 −− , to determine the respective masses as well as the couplings with theDD scattering channel. We also consider theDD scattering in partial wave l = 3 yielding information on the lowest charmonium resonance with J P C = 3 −− , for which LHCb discovered a candidate named X(3842) [21].
TheDD scattering channel in partial wave l = 1 is the dominant hadronic decay mode of the ψ(3770) (with branching ratio Br = 93 ± 9%). It is a well-established vector resonance generally accepted to be predominantly a conventionalcc state. We therefore assume elasticDD scattering, neglecting the influence of all other allowed decay modes of the ψ(3770) and of higher thresholds. The determination of the relevant scattering amplitudes and resonance parameters, such as the mass and decay width of the ψ(3770), will serve as a demonstration of our realization of Lüscher's finite volume treatment. In contrast to the only previous study [16], theDD scattering amplitude in p-wave is determined utilizing the finite volume spectrum from two lattice ensembles with different physical volumes, each in three different inertial frames. The study is performed for two charm quark masses, one below and one above the physical m c . In this way we explore the dependence of the spectra on the position of theDD threshold.
The recent LHCb discovery of the charmonium X(3842) with likely quantum numbers J P C = 3 −− [21] in part motivates our investigation of this narrow state with Γ 3 MeV. This resonance is also considered since it inevitably affects the spectrum of vector excitations in a finite hyper-cubic volume. While we are able to determine the mass of this charmonium resonance, the finite volume spectrum in our lattice setup is not dense enough to reliably estimate the very narrow width of this state.
The layout of this paper is as follows. In Section II, the relevant quantization condition relating the infinite volume scattering amplitudes to the finite volume energy spectrum is discussed. The details of the lattice setup and related technical information are briefly outlined in Section III. In Section IV, our procedure for the determination of the excited charmonium spectrum and our approach for dealing with discretization errors in hadron observables involving charm quarks is discussed. We summarize our findings in Section V and conclude in Section VI.

II. SCATTERING IN A FINITE BOX: LÜSCHER FORMALISM
In a Euclidean box of finite size, the QCD spectrum is discrete. Energy levels corresponding to bound states are affected by exponentially suppressed finite (spatial) volume effects, while a power law dependence on the size of the box is expected for energy levels related to two-particle states above the strong decay threshold. Lüscher's finite volume formalism and its extensions [22][23][24][25] (see the review Ref. [19] to find more references) show that the energy dependence of infinite volume scattering amplitudes determines the discrete energy spectrum in a finite-sized box with periodic boundary conditions.
We are interested inDD scattering, i.e. single-channel scattering of spin-zero particles with equal masses, in inertial frames with zero and non-zero total momenta P . In this case the scattering matrix S l (p) = e 2iδ l (p) can be expressed in terms of a single phase shift δ l (p) in the partial-wave l, where p is the momentum of each hadron in the center of momentum frame. The spin of the system is zero, since the scattering particles are spin-less and the total angular momentum equals the orbital angular momentum (J = l).
The quantization condition connecting the infinite vol-ume scattering amplitudes with the finite volume discrete energy E of the eigenstates can be expressed in our case as Here the determinant is computed over all partial waves l and E cm is the eigenenergy in the center of momentum frame. The real functionK l parametrizes the infinite volume scattering amplitude S l and the transition amplitude t l in the partial wave l as with ρ = 2p/E cm , or equivalentlỹ The kinematic variable p is the momentum of each D meson in its center-of-momentum frame with s = E 2 cm . We follow the definitions ofK in Ref. [26] and the definition of t l used by the Hadron Spectrum Collaboration (see e.g. Ref. [27]). 2 The factor p −2l−1 in Eq. (2) ensures a smooth kinematic behaviour ofK at the threshold.
The box-matrix B P ;Λ l l (E cm ), as introduced in Ref. [28], depends on E cm , the lattice spatial size L and the lattice irreducible representation (irrep) Λ under investigation. Its entries are customarily written in terms of the known Z-functions defined in Ref. [22]. Note that B also has offdiagonal entries in the partial wave indices l and l, due to the reduced rotational symmetries in a hyper-cubic box. The latter leads to the finite volume spectrum for any lattice irrep containing information from scattering in multiple partial waves.
The main focus of the present work are the resonances and bound states of charmonium in the vector channel that are related toDD mesons scattering in p-wave. In the simplest approximation, in which the influence of higher partial waves l > 1 are neglected, the quantization condition in Eq. (1) allows one to compute δ 1 directly from a given measured energy level in a finite box. In this case, the determinant condition reduces to the well-known Lüscher relation [22] for zero total momentum. In the case when the effects of higher partial waves are non-negligible, the properties of the scattering amplitudes S l are accessible only after finding and fitting a suitable parametrization of all the relevantK l , so that the determinant condition in Eq. (1) is solved for each of the energy levels in the finite volume spectrum. The finite volume mixing of even and odd partial waves is not allowed if the scattering particles have equal masses [29,30], as in our case. Hence the contamination of higher partial waves in the l = 1 excitations due to broken rotational symmetry can occur only from odd partial waves l = 3, 5, . . . . The contributions from the l > 5 partial waves are expected to be strongly suppressed in the nearthreshold region explored here and it can be neglected. We will therefore investigate partial waves l = 1 and l = 3 (as well as their mutual influence) through the finite volume quantization condition of Eq. (1). We use the package TwoHadronsInBox to compute the determinant in Eq. (1) when fitting the parametrization of thẽ K-matrix to our energy levels, following the "determinant residual" method [28].

III. LATTICE SETUP
The results presented in this work have been determined from two N f = 2 + 1 ensembles labeled U101 and H105 generated by the Coordinated Lattice Simulation (CLS) consortium [31,32]. The inverse gauge coupling is β = 6/g 2 = 3.4 and the lattice spacing is a = 0.08636(98)(40) fm [33]. The lattice volume is T × L 3 = 96 × 32 3 for H105 and 128 × 24 3 for U101. The Wilson-clover action [34] is non-perturbatively improved to remove the O(a) lattice artefacts, with the clover term set equal to c sw = 1.986246. The strange and light hopping parameters are κ l = 0.13697 and κ s = 0.13634079, respectively, leading to pion and the kaon masses of m π = 280(3) MeV and m K = 467(2) MeV. Note that the light quark mass (m l ) is heavier than its physical value while the strange quark mass (m s ) is lighter. This is because the U101 and H105 ensembles lie on a trajectory where the physical point is approached keeping the average quark mass 2m l + m s fixed. The gauge and the fermion fields fulfill open boundary conditions in the time direction and periodic boundary conditions in the spatial directions. The total number of configurations analyzed is 255 for the U101 ensemble and 492 for H105, with one configuration taken every 20 and 16 MDUs, respectively, from two different replica in each case.

A. Choosing the charm quark mass
The bare charm quark mass, set by the hopping parameter κ c = 1/(2m c + 8) for Wilson-type actions, is a parameter to be determined with care. As we aim to study the strong decay of near-threshold charmonium resonances, the particles relevant for tuning the charm quark mass are the lowestcc-states, the J/ψ and the η c , and the decay products of their excited states, such as the D and D s mesons. Away from the physical point, the masses of all these particles cannot be tuned simultaneously to their experimental values. However, the most critical quantity for our current investigation is the position of the strong decay threshold relative to the J/ψ and the η c , rather than the mass of an individual charmonium state. We aim to explore how charmonia near and above the open-charm threshold depend on its position. We also want to investigate the systematics coming from the tuning of the charm quark mass. To this end we employ two different values of κ c given by 0.12522 and 0.12315, corresponding to the D meson mass being below and above its physical value, see Table III. The main results for the position of the resonances in physical units will be presented relative to the spin average mass of the J/ψ and the η c , This is motivated by the fact that some of the discretization effects (which can be significant at a lattice spacing of a ≈ 0.09 fm) are canceled when computing energy differences. The spin-average mass is equal to 3103(3) MeV for κ c = 0.12315 and 2820(3) MeV for κ c = 0.12522. The energy gap between M av and 2m D is Note that the relative position of theDD threshold with respect to M av is closer to the experimental value for κ c = 0.12522 than for κ c = 0.12315.

B. Setup of the distillation algorithm
We employ the distillation method to compute our correlation functions [35], i.e. the propagation of the quarks is computed in the distillation space of the eigenvectors of the 3D Laplacian. As a consequence, we are able to separate the construction of the interpolators from the computation of the light and charm propagators without resorting to storing full propagators. Only at a later time, the quark-lines of the Feynman diagrams are contracted and summed to obtain the entries of the correlation matrix. We employ 90 and 150 eigenvectors of the lattice Laplacian for the U101 and the H105 gauge field ensembles, respectively. We use full distillation for U101. Full distillation is also used on H105 for the charm propagators and the φ-matrices (defined in Eq. (12) of Ref. [35]). The cost of full distillation for the determination of light quark propagators on the H105 ensemble would have been prohibitively expensive, as 150×4 inversions of the Dirac operators on a large volume would have been required for each flavor and each time-source/sink. Therefore, we use stochastic distillation [36] for the light quark propagators on the H105 ensemble only. The perambulator, defined as Energy spectrum E lat n in various lattice irreducible representations for κc = 0.12522 for the U101 (L = 24) and H105 (L = 32) ensembles. The blue energy levels correspond to the naive ψ(2S) energy levels, while the green points have been identified as having J P C = 3 −− . The red points correspond to the remaining energy levels, which arise from charmonia with J P C = 1 −− and fromDD. Black solid lines show the energies of the non-interactingD(p 2 1 )D(p 2 2 ) lattice energy levels, where p 2 1,2 = 0, 1 or 2, respectively in units of (2π/L) 2 . The black dashed line in the rest frame indicates theDD threshold, while the grey dashed lines refer to the DsDs and DD * thresholds.
requires to compute the inverse of the Dirac matrix D on a source |v µ i (t) which is zero everywhere except for the timeslice t and spin index ν, where it is equal to the Laplacian eigenvector j. The number of inversions required can be reduced by considering random sources. We introduce N s complex Z 2 stochastic noise vectors using full-time dilution. The stochastic source given to the multigrid solver is constructed for every time-slice t as The perambulator can be approximated as which is equal to the exact expression of Eq. (7) in the limit N s → ∞, since Thanks to the use of full-time dilution in our stochastic distillation scheme, there is no bias for contributions to the correlation matrix involving light annihilation diagrams. However, we need to consider two independent sets of random vectors for diagrams involving light quarks propagating at different timeslices. We therefore use two sets of 20 noise vectors, and we average over the two only for light annihilation diagrams. Following this approach, we are able to reduce the number of the expensive inversions of the Dirac operator for the light quarks by a factor ten on our largest volume. At the same time, we retain the advantages of full-dilution for the purely charm contributions and employ full distillation to compute the correlation functions betweencc-type operators, that need to be estimated precisely for the study of charmonium excited states we are interested in.

IV. ENERGY SPECTRUM
A. Determination of the finite volume spectrum Many excited energy levels need to be accurately determined on several different volumes and/or in different momentum frames in order to extract the scattering matrix according to Lüscher's method. The excited energy spectrum E lat n is extracted from Euclidean time two-point correlation functions constructed from a basis of interpolators {O i } with the desired quantum numbers. Here Z n i = 0|O i |n refers to the operator state overlap.
Our main purpose is to investigateDD-meson scattering in partial wave l = 1. We determine the charmonium spectrum in three lattice irreducible representations that contain l = 1 (they all also contain partial wave l = 3): 1 in the moving frame | P | 2 = 1 and iii) Λ C = A − 1 in the moving frame | P | 2 = 2 , where each irreducible representation corresponds to a different inertial frame, with total momentum P . The values of | P | 2 are given in units of (2π/L) 2 . Note that charge conjugation C is a good quantum number in all frames, whereas parity P is a good quantum number only in the rest frame. In order to reliably extract the relevant low energy spectrum, a large basis of interpolators {O i } that includes single-meson as well as two-meson operators in all the irreps is required.
For the single-meson interpolators, we follow the procedure proposed in Refs. [37,38] and construct all interpolators with up to two gauge covariant derivatives. As a first step in the analysis, the spectrum is extracted using a basis of these single-meson interpolators and, following the spin identification procedure discussed in Ref. [20], the continuum quantum numbers J P of the levels in the energy region of interest are identified. In this way, we can understand the quantum numbers of the energy levels that are related to a naivecc description that enter into our phase shift analysis.
The two-meson interpolators are built using the projection where p 1 and p 2 are the momenta of the scattering particles such that P = p 1 + p 2 , R is a group element of the point symmetry group G of the inertial frame with momentum P , T Λ r,r (R) is the representation for the group element R in the irrep Λ and r is the row of the irrep Λ. AsDD is the predominant decay mode of the ψ(3770) resonance, we includeDD interpolators for all momentum combinations within the relevant low energy region.
In Table I, we show the two-meson interpolators considered in our calculation in each of the lattice irreps being analyzed. In this study, we neglect the p-wave decay J/ψη and in addition the final states including three or more hadrons, such as J/ψπ + π − and J/ψπ 0 π 0 .
The inclusion of two-meson interpolators in the basis brings additional Wick contractions. We compute all relevant Wick contraction diagrams for the single-and twomeson operators (see Fig. 1 of Ref. [16]). We exclude diagrams with disconnected charm quark contractions, that lead to the decay of charmonium into light hadrons. The correlation functions constructed from the single and two meson operators are averaged over the spin and momentum polarization to increase the statistical precision. We also bin our data in blocks of size equal to one, two and four to correctly address the autocorrelation of our configurations. For further details, see Appendix A.
The extraction of the energy levels proceeds via a variational approach [39][40][41] by solving a generalized eigenvalue problem (GEVP) for the correlation matrices in Eq. (12). We set t 0 = 2 in our analysis. The eigenvalues λ n (t) are fitted to a single or a double exponential function and the quality of the fits is compared to estimate the best fitting intervals. We remark that increasing t 0 leads to consistent results, however, the errors also increase. Corrections to these fit forms arise when using open boundary conditions in the temporal direction and the source or sink timeslice of the correlation function is close to one of the boundaries. However, all our measurements are made in the bulk of the lattice, 28 time slices away from either boundary, and no such corrections are observed. In Figs. 1 and 2, we present the finite volume excited energy spectrum (above the J/ψ) in all the three frames studied for κ c = 0.12522 and κ c = 0.12315, respectively. The statistical uncertainty is smaller on the H105 ensemble, in part due to the larger number of configurations analyzed compared to U101. The color coding of For brevity we use the simplified notation D(p 2 1 )D(p 2 2 ), which refers to two-meson interpolators, with good charge conjugation C, built from aD and D meson operator with momentum p1 and p2, respectively in units of 2π/L. The full expressions for these interpolators are much longer and are omitted for brevity.
the levels indicates the quantum numbers. Blue levels are related to the ψ(2S) bound state and the green levels are identified to have quantum numbers J P C = 3 −− . 3 . The remaining levels in red arise from charmonia with J P C = 1 −− and fromDD scattering levels shifted due to interaction. We extract up to 4 or 5 excited states for each lattice irrep considered.
Within the energy range explored, we also observe a 2 −− level in the A − 1 (| P | 2 = 2) spectrum. The helicity two components of the 2 −− state appear in this irrep as can be seen from Table II of [38] and Table I of [20]. Taking the naive energy level for κ c = 0.12522 and using the continuum dispersion relation, we determine the mass of this state in its rest frame to be 3825(8) MeV, which is consistent with the lowest 2 −− charmonium state observed by Belle [42] and BESIII [43]. The coupling of this excitation with the rest of the spectrum is expected to be negligible as, in a study of vector channel scattering of two equal mass pseudoscalar mesons, J P = 2 − is an unnatural quantum number and l = 2 is an irrelevant partial wave. This level is absent when we omit the relevant 2 −− interpolators from the operator basis, while the rest of the spectrum remains unaffected. Hence we exclude this level from the amplitude analysis and also do not show it in the figures.
The thresholds for the channels D s D s andDD * , which are omitted in our simulations, are shown in Fig. 1 and 2. These thresholds are lower than some of the employed energy levels on the N L = 24 ensemble, however we expect that the effect of the omitted channels is very small because they appear in p-wave and they open at relatively high energy. Note that there are no close-by two-particle levels of the omitted channels because they appear in p-wave. The influence of the omitted channels on the results from the N L = 32 ensemble is expected to be negligible 4 . aED(p 2 = 0) aED(p 2 = 1) aED(p 2 = 2) aMav κc = 0.12522, NL = 24 lat 0.7732 (9) 0.8132 (10) (15). Above p 2 are given in units of (2π/L) 2 . The spin-averaged charmonium mass of Eq. (6) is also given.

B. Approach to dispersion relations and charm-quark discretization errors
The lattice energy levels presented in Figs. 1 and 2 are inputs to the quantization condition given in Eq. (1) (realized using the TwoHadronsInBox package of Ref. [28]), which is based on continuum dispersion relations. However, lattice energy levels are subject to discretization effects and can deviate from continuum expectations. Such lattice artifacts are larger in magnitude for charm than for the light quarks u, d and s, since, in dimensionless units, 1 ∼ am c am u/d,s even on many currently available lattices. At finite lattice spacing a, there are nonnegligible deviations of the dispersion relation E lat D (p) of the D-meson from its continuum counterpart As is evident from Table II, there are non-negligible deviations from E cont D ( p) particularly for the smaller physical volume (for which the lattice momentum is larger), and hence care is required when extracting the scattering matrix in theDD channel. In this section we describe our approach to this issue.
First, we determine the energy shift of each interacting eigenstate with total momentum P with respect to the nearest non-interacting stateD( where p 1,2 = n 1,2 2π L and p 1 + p 2 = P . All three energies on the r.h.s. are extracted from their corresponding eigenvalues of the GEVP, separately. Here, (E lat ) s is the energy of the interactingDD system, presented in Section IV A, while (E lat D( p1) ) s is the energy of a single D meson with momentum p 1 measured on the lattice. All are extracted on a given bootstrap or jackknife sample s. If the scattering matrix is equal to the identity, the monia from the N L = 32 ensemble alone are compatible with the ones based on combined fits of all volumes. phase shift is equal to zero and hence the energy shifts in Eq. (16) should be equal to zero. Given that the quantization condition is based on the continuum dispersion relation, 5 we ensure this important constraint is satisfied by using the energies as input to Eq. (1). The energies (E cont D( p1) ) s and (E cont D( p2) ) s are computed from Eq. (15) using the lattice momenta p 1,2 and the D-meson mass at rest (also determined from the bootstrap or jackknife sample s). Note that, the energies E calc are equal to E lat in the continuum limit a → 0 by construction (Eqs. (16) and (17)). Furthermore, the above E calc also clearly provides a zero phase-shift when ∆E in Eq. (16) is zero, even if non-negligible heavy-quark discretization effects are present in our simulation.
The ground state in the vector channel is the J/ψ, while the first excited state, the ψ(2S), lies below the open charm decay threshold. Beyond these two states, there is a relatively narrow resonance, the ψ(3770) with Γ exp ≈ 27 MeV, and a broader resonance, the ψ(4040), both which can decay toDD in p-wave.
The near-threshold charmonium resonance ψ(3770) is the main target in this investigation. In particular, we aim to extract the pole position of the ψ(3770) and its coupling toDD, which is related to the ψ(3770) decay width. The study of the vector channel provides a benchmark for our simulations enabling us to assess systematic uncertainties which will also arise when considering 5 The Lüscher Z-functions, that enter the box matrix B in the quantization condition, have poles at energies E = E cont D( p 1 ) + E cont D( p 2 ) with p 1,2 = n 1,2 2π L and p 1 + p 2 = P , while tan δ in the Lüscher-type relations is a sum of terms with Z-functions in the denominator. Therefore δ = 0 at those energies E. more complex cases where several decay channels are involved (and for which the experimental situation is unclear). We assume that the ψ(3770) resonance is only coupled with theDD scattering channel and hence its parameters can be completely determined from the elastic DD scattering amplitude. Experimentally, the ψ(3770) decays intoDD with a branching ratio of 93 −9 +8 % [44]. In this work, we neglect the effects from other hadronic decay modes, such as J/ψη, J/ψπ 0 π 0 and J/ψπ + π − , which collectively have a branching ratio below 0.2%. We also neglect decays into light hadrons through charm-

annihilation.
Our second aim is to study the lowest 3 −− charmonium resonance. A candidate for this state referred to as the X(3842) 6 was recently discovered by LHCb inDD decay and has width Γ exp ≈ 2.8 MeV [21].

A. Fits of the phase shift
We perform a combined fit of the energy levels calculated on the U101 and H105 ensembles using bootstrapping for the determination of the uncertainties (see Appendix A for details). The ground state level in each frame, corresponding to J/ψ, is excluded from our fits. We include the energy level corresponding to the ψ(2S) to use the maximal information from our spectrum and to correctly describe the spectrum in the vicinity of the threshold. The fitting forms we have explored for the vector channel are the "double-pole" inspired by two states in the region of interest, the ψ(2S) and the ψ(3770) 7 ; and the "quadratic" form where p 2 = s/4 − m 2 D . Other parameterizations were investigated. A linear form did not describe the data well; a triple-pole form would require additional data at large momenta to enable a third pole to be resolved; a double pole form with an additional constant is described in Appendix C.
Although many partial wave amplitudes in principle contribute to the finite volume spectrum in any given irrep starting at theDD threshold, the effect of higher partial waves should be small, unless there is a resonance contributing to a higher partial wave in the energy region of interest. The spin-identified finite volume spectrum for charmonium from our previous investigation [20] indicates a low lying J P C = 3 −− resonance, that contributes to the l = 3 partial wave. In order to investigate the effects of such a low-lying l = 3 partial wave excitation on the rest of the discrete spectrum, we parametrize it with a simple Breit-Wigner form given by Our fits include all the fourteen excited energy levels we have been able to determine on the H105 and the U101 ensembles. For both κ c , the fits are more constrained by the data from the H105 ensemble. The "double pole" fit is presented for the 1 −− channel in Fig. 3(a) for κ c = 0.12522 and in Fig. 3 for κ c = 0.12315 with a χ 2 /d.o.f. = 0.3. The quadratic fit is presented for the 1 −− channel in Fig. 4(a) for κ c = 0.12522 and in Fig. 4(b) for κ c = 0.12315. Compared to the double pole fit form, the minimum of the correlated χ 2 is more unstable across different bootstrapping resamples, resulting in an asymmetric distribution of the fitted parameters with long tails. We therefore quote asymmetric uncertainties for this parametrization. The fitted quadratic functional form for κ c = 0.12522 is equal to with a χ 2 /d.o.f. = 1.32, while for κ c = 0.12315 7 Here cot(δ 1 ) = 0 at s = m 2 1 and s = m 2 2 . Between these two energies, cot(δ 1 ) has a pole and T has a zero at s =  The quality of our fits for the combined channel case is presented in Figs. 5, 6, 7 and 8, in terms of the Ω function, as introduced in Ref. [28] to minimize the χ 2 in the "determinant residual method". Here the matrix A(s) is the argument of the determinant in Eq. (1) and µ = 8 is a regularization parameter, see Appendix A for further details on the choice of µ. Each crossing of the Ω function with the vertical axis Ω = 0 represents a predicted energy level from our parametrization, which can be compared with our observed energy spectrum. The number of energy levels is indeed equal to the number of solutions of Ω = 0 in the energy region between E m ψ(2S) (the blue energy level) and the highest level. In particular, there are two approximately degenerate energy levels when the fitted Ω-function has a double zero (this occurs for 3 −− when P 2 = 2). We remark that our scattering amplitude captures the energy dependence only in the energy region where we have data and is not applicable outside of this region. There is an additional crossing of Ω = 0 in certain irreps for κ c = 0.12315 in Fig. 5, however, it occurs for E < m ψ(2S) . The amplitude is not constrained and is, therefore, not applicable there. 8 The parameters for the 3 −− resonance are quite unstable across the bootstrap samples, resulting in large uncertainties, especially for the coupling g 3 . To check the stability of the fits in the vector channel, we have also performed an alternative analysis: the levels related to J P C = 3 −− are identified (following the methods uti- . 7. Ω-function in the three lattice irreducible representations in various frames for the vector channel resulting from the fit of all energy levels using the "quadratic" fit ansatz for the phase shift for κc = 0.12315. lized in Ref. [20]) and excluded from the fits. The resulting couplings and masses for the 1 −− resonance are found to be compatible with the values quoted above and we can therefore conclude that the influence of the 3 −− resonance on the ψ(3770) is negligible in our simulation.

B. Bound states, virtual bound states and resonances
Experiments and lattice QCD determine the scattering amplitudes t l for real energies. The physical interpretation in terms of (virtual) bound states and resonances is however conventionally performed by looking at poles of the t l (s) in Eq. (2), continued to the complex s-plane The two Riemann sheets I and II, from the square root branch cut are defined to have Im(ρ) > 0 and Im(ρ) < 0, . 8. Ω-function in the three lattice irreducible representations in various frames for the vector channel resulting from the fit of all energy levels using the "quadratic" fit ansatz for the phase shift for κc = 0.12522. respectively 9 . The product is an analytic function of s only, where our lattice results for p 2l+1 cot δ l / √ s were expressed in terms of s in Section V A. The scattering amplitude t l in Eq. (26) has a pole at a given p if cot δ(p) = i. A bound state appears as a pole of the t l matrix at p = iκ B , i.e. on the real axis below the scattering threshold on the first Riemann sheet. A virtual state appears as a pole at p = −iκ B , i.e. on the real axis on the second Riemann sheet. In both cases κ B denotes a real positive number. It is convenient to express these 9 The square-root branch cut for ρ is chosen in the conventional sense, such that it runs from the threshold s = 4m 2 D to +∞ along the positive real axis. We neglect the effects from the left hand cut and from the freedom in choosing the real part of iρ in Eq. (26) (such as using a Chew-Mandelstam phase space), either of which results in pole structure deep below the energy region being studied.  (28) and virtual bound state (29) conditions. The intersection of the red and blue curves corresponds to a bound state, while the intersection of the red and orange curves corresponds to a virtual bound state. pole conditions in terms of the quantity p 2l+1 cot δ, which is considered in the Section V A for l = 1, 3. Such poles appear if the conditions p 2l+1 cot(δ l ) = −(p 2 ) l −p 2 (bound state) (28) p 2l+1 cot(δ l ) = (p 2 ) l −p 2 (virtual bound state) (29) are fulfilled. A resonance corresponds instead to a pair of poles away from the real axis above threshold on the second sheet, influencing the physical region in the first sheet.
C. Consistency check for a bound state in partial wave l The S(p) = e 2iδ(p) matrix for a bound state has a pole at p = iκ B , as discussed in the previous section. The analyticity of the S(p) matrix near the real energy axes implies the sign of the pre-factor in front of the pole to be where β 2 is a real and positive number. This functional form is derived from the general form of the solutions for S(k) continued to complex p in Ref. [45] (see Eq. (7.60)) and applied for s-wave in, for example, in Ref. [46]. Equation (30) must be satisfied for all bound states, but it does not apply to virtual bound states. We have explicitly verified that the relation (30) with positive β 2 is satisfied for s-wave and p-wave bound states in 3D nonrelativistic quantum mechanics for various forms of the potential (square well, Yukawa, Wood-Saxon, and several other forms of the type V (r) = 1 r n exp (−Ar m )). Let us express the condition in Eq. (30) in terms of p 2l+1 cot(δ l ) in Eq. (28) which was fitted from our data. For this purpose we consider the dependence of the quantity on p 2 near the bound state pole, where it equals zero in Eq. (28). Inserting in Eq. (31) cot δ = i(S + 1)/(S − 1) with S from Eq. (30), and expressing p with p 2 as p = i|p| = i −p 2 , we obtain Taking the derivative of the previous expression with respect to p 2 and evaluating it at p 2 = −κ 2 B , we get The above condition implies that the slope of p 2l+1 cot δ in terms of p 2 has to be smaller 10 than the slope of −(p 2 ) l −p 2 at the position of the bound state in partial wave l, while it does not apply for virtual bound states. Note that the position of the bound state is where these two curves cross. The condition for s-wave was referred to as a sanity-check in Ref. [46]. Let us now turn to the fitted phase shifts in Section V A and investigate which results satisfy the condition of Eq. (30) or equivalently the consistency check of Eq. (33). The slope of p 2l+1 cot(δ) (the red curves in Fig. 9) has to be smaller than the slope of −(p 2 ) l −p 2 (the blue curve in Fig. 9) where these two curves cross. This is satisfied for the double-pole fits, while it is not satisfied for the quadratic fit as can be seen from Fig. 9. Note that all these slopes are negative, and in this case the absolute value of slope for p 2l+1 cot(δ) has to be larger than the absolute value of slope for −(p 2 ) l −p 2 . The physical parameters for bound states and resonances are based on our fits of DD scattering in p-wave. The first step is to ensure that our fitted parametrizations fulfill all the physical requirements, and in particular the consistency check given by the condition in Eq. (33). The double-pole fits in Section V A satisfy this condition, while the quadratic fits do not, as explained in the previous subsection. Therefore, in the following, in terms of the physics conclusions we focus on the doublepole fits.
The absolute value of the resulting amplitude t l=1 (s) (Eqs. (26) and (27)) in the complex energy plane is shown in Figs. 10 and 11 for κ c = 0.12522 and 0.12315, respectively. The locations of these poles on both Riemann sheets are collected in Fig. 12. These poles are related to (virtual) bound states and resonances according to the general criteria from Section V B. The final masses of various states will be quoted according to where the mass splitting [m−M av ] lat with respect to M av of Eq. (6) is determined from the lattice. Figure 13 and Table III summarize  real axis in the first Riemann sheet, corresponding to the bound state ψ(2S). In addition, there is a virtual bound state on the real axis on the second sheet.
The parameters for the ψ(3770) resonance are extracted from the linearized (Breit-Wigner) behavior of Eq. (21) in the resonance region The resonance width Γ can not be directly compared to experiment since the phase space (proportional to p 3 ) for our unphysical D meson mass is different to that in experiment. Therefore we extract the dimensionless coupling g = 16.0( +2.1 −0.2 ) that parametrizes the ψ(3770) →DD width. This coupling agrees within errors with the experimental value g exp = 18.7(9) obtained from Γ exp . The resonance mass equals m 2 of Eq. (18). The fitted it corresponds to the crossing of the red and blue curves in Fig. 9(a). The resulting mass in Table   11 The phase space p 2 /s in (35)   III is about 20 MeV smaller than the experimental mass.
• κ c = 0.12315 (m D = 1927(2) MeV): ψ(3770) is a shallow bound state for our heavier charm quark mass, in contrast to experiment where it is a strongly decaying resonance. It corresponds to a near-threshold pole in the first sheet of Fig. 11(a), arising from the fact that the crossing of our fits with the real axis occurs on the left of the axis origin (below threshold) in Fig. 3(b). The mass of the bound state m = 3776(7) MeV in Table  III is determined according to Eq. (28) and corresponds to the near-threshold crossing of the red and blue lines in Fig. 3(b). The ψ(3770) cannot strongly decay intoDD, but we still extract the coupling g = 18.9( +0.8 −0.7 ), which parametrizes the slope of p 3 cot δ 1 / √ s near the real axes according to Eq. (35).
The ψ(2S) corresponds to the lighter of the two bound states. Its mass m 3687 MeV, determined from the bound-state condition of Eq. (28), has a large error due to the huge uncertainty of the coupling G 2 in the double-pole parametrization in Eq. (22). This mass is consistent within errors with the value m = 3665(9) MeV obtained directly from the first-excited energy-level at P = 0 and L = 32. This agreement is expected as the influence of the threshold on a state which lies about 150 MeV below is very small. We quote the latter value as the final result for the ψ(2S) mass at this κ c .
Our results for the masses of the ψ(3770) and ψ(2S) states are summarized and compared to experiment in Fig. 13 and Table III. The mass m and the coupling g for the ψ(3770) agree with the experimental values within errors, which applies to both charm-quark masses considered. The mass of the ψ(2S) is about 2σ below the experimental value, which is not unreasonable given the lack of continuum and quark-mass extrapolations to the physical values.
The above results (that are based on a combined fit to two different ensembles) are subject to a technical subtlety. Both ensembles have been generated with the same bare lattice parameters, but with different lattice sizes. Hence sub-leading exponentially suppressed finite volume effects, which are generally neglected in the Lüscher finite volume formalism, can sometimes lead to different physical situations for different box sizes. Such a delicate situation can happen when a pole singularity in the complex energy plane is expected to be very close to the threshold: a small volume may exhibit different physics from that of a larger volume. For example, in our case for the lighter D-meson mass, we see that our N L = 24 data have a tendency to prefer a bound state (unlike the N L = 32 data), although the data are compatible with a resonance within the large errors. In order to confirm that our combined fit does not lead to conclusions that are different from the physical situation in the larger volume (where the neglected exponential corrections are small) we performed an analysis of the N L = 32 data alone. We find the mass and the coupling from such an analysis to be m 2 = 3782(7) MeV and g = 15( +2 −1 ), which are in agreement with the results from our combined fits.
The LHCb collaboration recently reported the discovery of the X(3842) charmonium state in theDD invariant mass [21] with a mass m exp 3842 MeV, a very narrow width Γ exp 2.8 MeV and likely quantum numbers J P C = 3 −− . The presence of such a 3 −− resonance was explored by fitting the parametrization of theK matrix in Eq. (1) to our energy levels including the Breit-Wigner form in Eq. (20) for l = 3.
The resulting mass for κ c = 0.12522 (see Eq. (21)) is The mass obtained from our lattice simulation is consistent with the mass of the discovered X(3842) state. This adds confidence to the interpretation of the X(3842) as a 3 −− charmonium resonance. The results for both values of the charm quark mass are compared to experiment in Table III. Our simulation is not sensitive to the width of   Table III. this very narrow resonance since there are noDD eigenstates in the very narrow energy region m ± Γ. 12 In particular, the finite volume lattice energy by itself would lead to compatible results for the mass of the X(3842). For this reason we can not compare the width or the X(3842) →DD coupling with experiment.

F. Comparison to previous results
ElasticDD scattering and its relation to the ψ(2S) and ψ(3770) resonances have previously been investigated in a lattice study in Ref. [16] for m π 156 MeV and 266 MeV. In this subsection we compare our current study to these older results, focusing on those from the m π 266 MeV ensemble, since the pion mass is similar and the errors are smaller. The most significant improvement over the calculation in Ref. [16] lies in the use of two lattice volumes at the same lattice spacing and light-quark masses, combined with calculations in moving frames. The latter is facilitated by a new approach to dealing with discretization effects described in Section IV B along with the spin identification techniques pioneered in Ref. [49] and employed on our ensemble(s) in J P C lat (present work) lat (present work) exp lat [16] [16]) which are the appropriate masses for the Fermilab approach [47,48] pursued in Ref. [16]. These differ from the rest masses by discretization effects.
Ref. [20]. The use of two volumes and several moving frames results in a significant increase in the number of energy levels obtained, which renders more information on the dependence of the l = 1, 3 scattering amplitudes on E cm .
One important consequence is that the J P C = 3 −− resonance, which contributes to the same irreducible representations on a hyper-cubic lattice, can now be included in our fits. In Ref. [16] it had to be assumed that after omission of an energy level related to the presence of this resonance, its effect on the remaining finite volume energies was negligible. For our current ensembles this assumption can instead be tested, and our results agree with the mass of the recently observed X(3842) [21].
The l = 1 scattering amplitude in Ref. [16] was constrained only at three values of E cm and the results for ψ(2S) and ψ(3770) were summarized in Table V of Ref. [16] for two types of fits. Fit (ii) in Eq. (5.6) was similar to the "quadratic fit" used in the present work and does not satisfy the consistency check of Eq. (30) for the ψ(2S) bound state (this constraint on the physical parameterizations was not checked in Ref. [16]). Fit (i) was the Breit-Wigner fit of the energy region above threshold and rendered the ψ(3770) resonance parameters summarized in the right column of Table III; the mass is consistent with our present result while the coupling is somewhat smaller. The ψ(2S) mass m = 3674(6) MeV obtained from the naive energy level (second level in Table III of Ref. [16]) agrees within errors with our present result. The two-pole fit with four parameters could not be preformed using just three energy levels in Ref. [16].
Another advantage of the present study is the use of the CLS gauge ensemble library, which in principle enables us to extend the current study to further ensembles with different lattice spacings and light quark masses, while maintaining our current setup.

VI. CONCLUSIONS
We determined the elasticDD scattering amplitude in the energy region of the ψ(2S) bound state and the ψ(3770) resonance from lattice QCD simulations using the finite volume Lüscher method. The results we obtain are from simulations with 2+1 flavors of dynamical light quarks at unphysical masses and at a single lattice spacing. We investigate different parametrizations of the p-wave scattering amplitude and impose a consis-tency check on the extracted bound state poles to ensure the physical bound state constraint in Eq. (33). Using a "double pole" parametrization defined in Eq. (18), the bound state mass of the ψ(2S) and the resonance parameters of the ψ(3770) have been extracted from a simultaneous fit to the energy spectrum from three different moving frames and two lattice volumes.
The lattice study of a near-threshold charmonium state away from the physical point requires a careful choice of the charm quark mass for the strong decay toDD to be allowed. To this end, on each gauge field ensemble, we employ two charm quark mass values leading to D-meson masses 100 MeV below and 60 MeV above the physical value. The masses of the ψ(2S) and ψ(3770), measured as splittings from the spin-averaged ground-state mass, are in good agreement with experiment for these values of the heavy quark mass. The result for the coupling g of the ψ(3770) resonance with theDD scattering channel is compatible with the experimental value, even for the charm quark mass where the ψ(3770) is a bound state. The results from both heavy quark masses employed in this work, the results from previous existing work and the corresponding experimental values are summarized in Table III.
We also investigate theDD scattering in partial wave l = 3 that contains 3 −− resonances. One of the aims is to investigate how this partial wave contributes to the finite volume spectra. Including the lattice energy levels related to the lowest 3 −− state in the amplitude analysis and parametrizing the correspondingDD scattering amplitude with a Breit-Wigner form (Eq. 20), we are able to extract the mass of this state. Nice agreement is found with the mass of the recently discovered X(3842), which is interpreted as a 3 −− charmonium resonance. However, we are unable to reliably determine its coupling with thē DD scattering channel as the finite volume spectrum of DD scattering channels are not sufficiently dense to be sensitive to its decay width. We also find our results for the vector channel remain unchanged with the exclusion of the l = 3 partial wave from the amplitude analysis and hence conclude the influence of the low lying 3 −− resonance on the lattice estimates for the ψ(3770) resonance parameters is negligible in our calculation.
The main challenge for future determinations of scattering amplitudes in the charmonium spectrum will be the presence of multiple coupled channels for the study of conventional and exotic resonances away from theDDthreshold. A first preliminary analysis of the scalar channel has been presented in Ref. [50], and we plan to complete the study in a forthcoming publication. In addition it would be desirable to obtain a denser set of points by using more and/or larger volumes, to reduce the discretization effects by calculations on CLS ensembles with a smaller lattice spacing, and to investigate the approach to the physical point by utilizing ensembles with lighter pion masses. In particular, discretization errors are an important source of systematic uncertainty of our current work. Using the CLS trajectory with physical strange-quark mass [32] would be particularly attractive, as the splitting between the D and D s mesons on these ensembles is closer to the physical splitting and therefore results in a larger energy window for elasticDD-scattering. ondary quantities extracted or fitted from the raw data. Each bootstrapping resample consists of 255 and 492 randomly chosen configurations from the U101 and H105, respectively. The first step in the analysis is the determination of the covariance matrix for the correlated χ 2 minimization of our K-matrix parametrization. This is estimated simply as the standard correlation between two variables as given by bootstrapping. The minimum of the correlated χ 2 , defined from the Ω-function of Eq. (25) as given in Eqs. (38) and (44) of Ref. [28], is recomputed for each bootstrapping resample. We extract from here the distribution of secondary derived quantities such as the parameters of the K-matrix or the positions of the poles in the complex plane. The Ω-function depends on an additional regularization parameter µ, and we have verified that varying µ in the range between 2 and 16 does not change the central value nor the error of the fitted parameters significantly. The crossing with the vertical axis, determining the expected spectrum, does not change significantly if µ is varied in the above interval, see Fig. 16. Note however that the value of µ determines the slope and the curvature of the Ω-function at the crossing. Given that our fits combine two different ensembles with a different number of configurations, single elimination jackknife is not an optimal choice. However, it is worth to compare bootstrapping and jackknife if we restrict the fits to a single ensemble, say U101. Looking to Fig. 15, we can see that jackknife (for which the spread of the resamples is smaller by a factor √ N conf compared to bootstrapping) determines the error as a linear approximation of the Lüscher Z-functions around the central value of the phase shift. It might happen however that a given energy level lies very close to a turning point of the Z-function, as in the case of the four rightmost points of Fig. 15. In this situation, the strong non-linear nature of the phase-shift analysis could result in a systematic underestimation of the error when using jackknife. For the same reason, the bootstrap distribution of the phase shift is not symmetric and has long tails. We therefore conclude that bootstrapping provides a more reliable determination of the uncertainty compared to jackknife for phase-shift studies.
The error on the fitted parameters, as well as on the couplings and masses of the charmonium states, is determined as the standard deviation of the bootstrap samples, except for the case when the distribution is asymmetrical or has long tails. In this situation, we quote asymmetric uncertainties as determined by the interval that fits 68 % of the bootstrap samples, i.e. by an interval that excludes 16 % of the bootstrap samples from the left and right tail.

Appendix B: Correlation matrices of the fitted parameters
In this appendix we include the correlation matrix C between the parameters κ i of the "double pole" fits described in Sec. V A which fulfill the consistency checks of Sec. V C. The correlation matrix is defined as where the expectation values appearing in the numerator and the denominator are computed directly from the bootstrapping distributions andκ i = κ i . For the "double pole" fit form, the six parameters κ i are given in the order m 2 , m 1 , m 3 , G 2 , G 1 and g 3 . The correlation matrices are The "double pole" fit form is able to reproduce the charmonium spectrum including the ψ(2S), as discussed in Sec. V A. The presence of an additional crossing of the Ω-function below the ψ(2S) is an interesting feature that suggests the possibility of describing the full charmonium spectrum including even the J/ψ. For example, by using a slightly modified "double pole" fit form with an additional linear parameter ρ The corresponding Ω-function determined from the fits are presented in Fig. 17. While the couplings and the masses G 1 , G 2 , m 1 and m 2 are consistent within error with Eq. (21) and (22) after the inclusion of the J/ψ, the parameter ρ fluctuates strongly across the bootstrapping resamples. Given that the spectrum below the ψ(2S) is well separated from the two resonances we are investigating, we can safely exclude the J/ψ from our main analysis. We emphasize again that our focus is the scattering amplitude in the energy region above the ψ(2S) (the consistency check of Sec. V C would not be satisfied for the J/ψ with the parameterization C1). lattice irreducible representation in rest frames for the vector channel resulting from the fit of all energy levels including the J/ψ using the modified "double pole" of Eq. (C1) fit ansatz.