Charmonium resonances with and from scattering on the lattice
Abstract
We present a lattice QCD study of charmonium resonances and bound states with and near the opencharm threshold, taking into account their strong transitions to . 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 scattering amplitudes for partial waves and are extracted on the lattice by means of the Lüscher 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 nonperturbatively improved Wilson quarks are employed with MeV, a single lattice spacing of fm and two lattice spatial extents of and . Two values of the charm quark mass are considered to examine the influence of the position of the threshold on the hadron masses. For the lighter charm quark mass we find the vector resonance with mass MeV and coupling (related to the width), both consistent with their experimental values. The vector appears as a bound state with MeV. The charmonium resonance with is found at MeV, consistent with the recently discovered by LHCb. At our heavier charmquark mass the as well as the are bound states and the 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 charmoniumlike 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 opencharm decay threshold are found to fit within a quarkantiquark picture, the properties of several resonances above the threshold, generally referred to as states, call for exotic interpretations. Various suggestions have been made based on phenomenological approaches, effective field theories, potential models, etc. These include compact diquarkantidiquark 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. Abinitio determinations of ground state charmonia, such as in Refs. Burch et al. (2010); Donald et al. (2012); Bečirević and Sanfilippo (2013); Briceño et al. (2012); Yang et al. (2015); Galloway et al. (2014); Mathur et al. (2018); DeTar et al. (2019), illustrate the power of these nonperturbative 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 charmoniumlike 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 of type and assumes that this qualitatively captures the single meson spectrum in a finite volume Dudek et al. (2008); Bali et al. (2011); Liu et al. (2012); Mohler et al. (2013); Cheung et al. (2016). So far, only one lattice study has gone beyond the singlehadron approach, taking into account the strong decays of charmonium resonances above the opencharm threshold.^{1}^{1}1The strong decay of the resonance with isospin 1 has been studied in the HALQCD approach Ikeda et al. (2016); Ikeda (2018). In that study an exploratory determination of the scattering amplitude in the and partial waves was performed Lang et al. (2015). The lowlying 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 Dudek et al. (2008); Bali et al. (2011); Liu et al. (2012); Mohler et al. (2013); Prelovsek and Leskovec (2013); Padmanath et al. (2015); Lang et al. (2015); Cheung et al. (2016). However, lattice calculations in the moving frames (i.e. with nonzero total momentum) can provide additional information on the relevant twomeson scattering amplitudes (for a review, see Ref. Briceño et al. (2018)). The challenge in spectroscopy using lattice techniques is that, due to the reduced symmetry on the lattice, several 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 (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 Padmanath et al. (2019). Such a spin identified singlehadron 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 of and mesons in the partial wave, which features charmonium resonances and bound states with , to determine the respective masses as well as the couplings with the scattering channel. We also consider the scattering in partial wave yielding information on the lowest charmonium resonance with , for which LHCb discovered a candidate named Aaij et al. (2019).
The scattering channel in partial wave is the dominant hadronic decay mode of the (with branching ratio ). It is a wellestablished vector resonance generally accepted to be predominantly a conventional state. We therefore assume elastic scattering, neglecting the influence of all other allowed decay modes of the and of higher thresholds. The determination of the relevant scattering amplitudes and resonance parameters, such as the mass and decay width of the , will serve as a demonstration of our realization of Lüscher’s finite volume treatment. In contrast to the only previous study Lang et al. (2015), the scattering amplitude in 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 . In this way we explore the dependence of the spectra on the position of the threshold.
The recent LHCb discovery of the charmonium with likely quantum numbers Aaij et al. (2019) in part motivates our investigation of this narrow state with MeV. This resonance is also considered since it inevitably affects the spectrum of vector excitations in a finite hypercubic 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 twoparticle states above the strong decay threshold. Lüscher’s finite volume formalism and its extensions Lüscher (1986, 1991a, 1991b); Briceño (2014) (see the review Ref. Briceño et al. (2018) to find more references) show that the energy dependence of infinite volume scattering amplitudes determines the discrete energy spectrum in a finitesized box with periodic boundary conditions.
We are interested in scattering, i.e. singlechannel scattering of spinzero particles with equal masses, in inertial frames with zero and nonzero total momenta . In this case the scattering matrix can be expressed in terms of a single phase shift in the partialwave , where is the momentum of each hadron in the center of momentum frame. The spin of the system is zero, since the scattering particles are spinless and the total angular momentum equals the orbital angular momentum ().
The quantization condition connecting the infinite volume scattering amplitudes with the finite volume discrete energy of the eigenstates can be expressed in our case as
(1) 
Here the determinant is computed over all partial waves and is the eigenenergy in the center of momentum frame. The real function parametrizes the infinite volume scattering amplitude and the transition amplitude in the partial wave as
(2) 
with , or equivalently
(3) 
The kinematic variable is the momentum of each meson in its centerofmomentum frame
(4) 
with . We follow the definitions of in Ref. Brett et al. (2018) and the definition of used by the Hadron Spectrum Collaboration (see e.g. Ref. Dudek et al. (2014)).^{2}^{2}2However, unlike in Ref. Brett et al. (2018), we do not divide our energy levels or physical quantities by the mass of the scattering particles. The factor in Eq. (2) ensures a smooth kinematic behaviour of at the threshold.
The boxmatrix , as introduced in Ref. Morningstar et al. (2017), depends on , the lattice spatial size and the lattice irreducible representation (irrep) under investigation. Its entries are customarily written in terms of the known functions defined in Ref. Lüscher (1986). Note that also has offdiagonal entries in the partial wave indices and , due to the reduced rotational symmetries in a hypercubic 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 to mesons scattering in wave. In the simplest approximation, in which the influence of higher partial waves are neglected, the quantization condition in Eq. (1) allows one to compute directly from a given measured energy level in a finite box. In this case, the determinant condition reduces to the wellknown Lüscher relation Lüscher (1986)
(5) 
for zero total momentum.
In the case when the effects of higher partial waves are nonnegligible, the properties of the scattering amplitudes are accessible only after finding and fitting a suitable parametrization of all the relevant , 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 Rummukainen and Gottlieb (1995); Dudek et al. (2013), as in our case. Hence the contamination of higher partial waves in the excitations due to broken rotational symmetry can occur only from odd partial waves . The contributions from the 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 and (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 the matrix to our energy levels, following the “determinant residual” method Morningstar et al. (2017).
Iii Lattice setup
The results presented in this work have been determined from two ensembles labeled U101 and H105 generated by the Coordinated Lattice Simulation (CLS) consortium Bruno et al. (2015); Bali et al. (2016). The inverse gauge coupling is and the lattice spacing is fm Bruno et al. (2017). The lattice volume is for H105 and for U101. The Wilsonclover action Sheikholeslami and Wohlert (1985) is nonperturbatively improved to remove the lattice artefacts, with the clover term set equal to . The strange and light hopping parameters are and , respectively, leading to pion and the kaon masses of MeV and MeV. Note that the light quark mass () is heavier than its physical value while the strange quark mass () 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 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.
iii.1 Choosing the charm quark mass
The bare charm quark mass, set by the hopping parameter for Wilsontype actions, is a parameter to be determined with care. As we aim to study the strong decay of nearthreshold charmonium resonances, the particles relevant for tuning the charm quark mass are the lowest states, the and the , and the decay products of their excited states, such as the and 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 and the , rather than the mass of an individual charmonium state. We aim to explore how charmonia near and above the opencharm 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 given by 0.12522 and 0.12315, corresponding to the meson mass being below and above its physical value, see Table 3. The main results for the position of the resonances in physical units will be presented relative to the spin average mass of the and the ,
(6) 
This is motivated by the fact that some of the discretization effects (which can be significant at a lattice spacing of fm) are canceled when computing energy differences. The spinaverage mass is equal to 3103(3) MeV for and 2820(3) MeV for . The energy gap between and 2 is
Note that the relative position of the threshold with respect to is closer to the experimental value for than for .
iii.2 Setup of the distillation algorithm
We employ the distillation method to compute our correlation functions Peardon et al. (2009), 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 quarklines 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. Peardon et al. (2009)). The cost of full distillation for the determination of light quark propagators on the H105 ensemble would have been prohibitively expensive, as inversions of the Dirac operators on a large volume would have been required for each flavor and each timesource/sink. Therefore, we use stochastic distillation Morningstar et al. (2011) for the light quark propagators on the H105 ensemble only. The perambulator, defined as
(7) 
requires to compute the inverse of the Dirac matrix on a source which is zero everywhere except for the timeslice and spin index , where it is equal to the Laplacian eigenvector . The number of inversions required can be reduced by considering random sources. We introduce complex stochastic noise vectors
(8) 
using fulltime dilution. The stochastic source given to the multigrid solver is constructed for every timeslice as
(9) 
The perambulator can be approximated as
(10) 
which is equal to the exact expression of Eq. (7) in the limit , since
(11) 
Thanks to the use of fulltime 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 fulldilution for the purely charm contributions and employ full distillation to compute the correlation functions between type operators, that need to be estimated precisely for the study of charmonium excited states we are interested in.
Iv Energy spectrum
iv.1 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 is extracted from Euclidean time twopoint correlation functions
(12) 
constructed from a basis of interpolators {} with the desired quantum numbers. Here refers to the operator state overlap.
For the purpose of investigating meson scattering in wave, we determine the charmonium spectrum in three lattice irreducible representations, each in a different inertial frame:
i)  
ii)  
iii) 
where the values of the square of the total momentum are given in units of . Note that charge conjugation is a good quantum number in all frames, whereas parity 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 {} that includes singlemeson as well as twomeson operators in all the irreps is required.
For the singlemeson interpolators, we follow the procedure proposed in Refs. Dudek et al. (2010); Thomas et al. (2012) 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 singlemeson interpolators and, following the spin identification procedure discussed in Ref. Padmanath et al. (2019), the continuum quantum numbers 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 naive description that enter into our phase shift analysis.
The twomeson interpolators are built using the projection
(13) 
where and are the momenta of the scattering particles such that , is a group element of the point symmetry group of the inertial frame with momentum , is the representation for the group element in the irrep and is the row of the irrep . As is the predominant decay mode of the resonance, we include interpolators for all momentum combinations within the relevant low energy region. In Table 1, we show the twomeson interpolators considered in our calculation in each of the lattice irreps being analyzed. In this study, we neglect the decay into final states including three or more hadrons, such as , and , that have extremely small branching fractions.
The inclusion of twomeson 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. Lang et al. (2015)). 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 Michael (1985); Lüscher and Wolff (1990); Blossier et al. (2009) by solving a generalized eigenvalue problem (GEVP)
(14) 
for the correlation matrices in Eq. (12). We set in our analysis. The eigenvalues are fitted to a single or a double exponential function and the quality of the fits is compared to estimate the best fitting intervals. 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 ) in all the three frames studied for and , 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 the levels indicates the quantum numbers. Blue levels are related to the bound state and the green levels are identified to have quantum numbers . The remaining levels in red arise from charmonia with and from 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 level in the spectrum. Taking the naive energy level for 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 charmonium state observed by Belle Bhardwaj et al. (2013) and BESIII Ablikim et al. (2015). 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, is an unnatural quantum number and is an irrelevant partial wave. This level is absent when we omit the relevant 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.
iv.2 Approach to dispersion relations and charmquark discretization errors
lat  0.7732(9)  0.8132(10)  0.8516(11)  1.2318(20)  

cont    0.8163(9)  0.8573(8)    
lat  0.7711(6)  0.7942(5)  0.8166(6)  1.2348(13)  
cont    0.7957(9)  0.8296(8)    
lat  0.8453(13)  0.8814(10)  0.9164(65)  1.3559(20)  
cont    0.8849(13)  0.9228(12)    
lat  0.8433(6)  0.8642(5)  0.8844(6)  1.3587(14)  
cont    0.8659(8)  0.8878(8)   
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. Morningstar et al. (2017)), 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 and , since, in dimensionless units, even on many currently available lattices. At finite lattice spacing , there are nonnegligible deviations of the dispersion relation of the meson from its continuum counterpart
(15) 
As is evident from Table 2, there are nonnegligible deviations from particularly for the smaller physical volume (for which the lattice momentum is larger), and hence care is required when extracting the scattering matrix in the channel. In this section we describe our approach to this issue.
First, we determine the energy shift of each interacting eigenstate with total momentum with respect to the nearest noninteracting state
(16) 
where and . All three energies on the r.h.s. are extracted from their corresponding eigenvalues of the GEVP, separately. Here, is the energy of the interacting system, presented in Section IV.1, while is the energy of a single meson with momentum measured on the lattice. All are extracted on a given bootstrap or jackknife sample . If the scattering matrix is equal to the identity, the 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,^{3}^{3}3The Lüscher functions, that enter the box matrix in the quantization condition, have poles at energies with and , while in the Lüschertype relations is a sum of terms with functions in the denominator. Therefore at those energies . we ensure this important constraint is satisfied by using the energies
(17) 
as input to Eq. (1). The energies and are computed from Eq. (15) using the lattice momenta and the meson mass at rest (also determined from the bootstrap or jackknife sample ). Note that, the energies are equal to in the continuum limit by construction (Eqs. (16) and (17)). Furthermore, the above also clearly provides a zero phaseshift when in Eq. (16) is zero, even if nonnegligible heavyquark discretization effects are present in our simulation.
V Results for charmonia with and
The ground state in the vector channel is the , while the first excited state, the , lies below the open charm decay threshold. Beyond these two states, there is a relatively narrow resonance, the with MeV, and a broader resonance, the , both which can decay to in wave.
The nearthreshold charmonium resonance is the main target in this investigation. In particular, we aim to extract the pole position of the and its coupling to , which is related to the 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 more complex cases where several decay channels are involved (and for which the experimental situation is unclear). We assume that the main contribution to the phase shift of the is coming from the elastic scattering of mesons. Experimentally, the decays into with a branching ratio of Tanabashi et al. (2018). In this work, we neglect the effects from other hadronic decay modes, such as , and , which collectively have a branching ratio below 0.2%. We also neglect decays into light hadrons through charmannihilation.
Our second aim is to study the lowest charmonium resonance. A candidate for this state referred to as the ^{4}^{4}4The LHCb collaboration has not identified the quantum numbers of this state but assumes as the mass and width of the fit model expectations. was recently discovered by LHCb in decay and has width MeV Aaij et al. (2019).
v.1 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. The ground state level in each frame, corresponding to , is excluded from our fits. We include the energy level corresponding to the 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 “doublepole”
(18) 
and the “quadratic” form
(19) 
where . Other parameterizations were investigated including, for example, linear and triplepole forms. However, the former did not describe the data well while for the latter we found additional data at large momenta would be required to enable a third pole to be resolved.
Although many partial wave amplitudes in principle contribute to the finite volume spectrum in any given irrep starting at the 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 spinidentified finite volume spectrum for charmonium from our previous investigation Padmanath et al. (2019) indicates a low lying resonance, that contributes to the partial wave. In order to investigate the effects of such a lowlying partial wave excitation on the rest of the discrete spectrum, we parametrize it with a simple BreitWigner form given by
(20) 
Our fits include all the fourteen excited energy levels we have been able to determine on the H105 and the U101 ensembles. The “double pole” fit is presented for the channel in Fig. 3(a) for and in Fig. 3(b) for . For this parametrization (and Eq. (20)), we have been able to fit all our data with
(21) 
for with a , while we have
(22) 
for with a .
The quadratic fit is presented for the channel in Fig. 4(a) for and in Fig. 4(b) for . Compared to the double pole fit form, the minimum of the correlated 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 is equal to
(23) 
with a , while for
(24) 
with a .
The quality of our fits for the combined channel case is presented in Figs. 5, 6, 7 and 8, in terms of the function,
(25) 
as introduced in Ref. Morningstar et al. (2017) to minimize the in the “determinant residual method”. Here the matrix is the argument of the determinant in Eq. (1) and is a regularization parameter. Each crossing of the function with the vertical axis represents a predicted energy level from our parametrization, which can be compared with our observed energy spectrum. It is interesting to observe that our fitted function has a double zero in case that there are two approximately degenerate energy levels, for instance when . For both , our fits are more constrained by the data from the H105 ensemble. We aimed to capture the phase shift dependence only in the energy region between the lowest and highest lattice levels. Note that there is an additional crossing with zero in certain irreps for , see Fig. 5. but below the energy region of interest.
The parameters for the resonance are quite unstable across the bootstrap samples, resulting in large uncertainties, especially for the coupling . To check the stability of the fits in the vector channel, we have also performed an alternative analysis: the levels related to are identified (following the methods utilized in Ref. Padmanath et al. (2019)) and excluded from the fits. The resulting couplings and masses for the resonance are found to be compatible with the values quoted above and we can therefore conclude that the influence of the resonance on the is negligible in our simulation.
v.2 Bound states, virtual bound states and resonances
Experiments and lattice QCD determine the scattering amplitudes for real energies. The physical interpretation in terms of (virtual) bound states and resonances is however conventionally performed by looking at poles of the in Eq. (2), continued to the complex plane
(26)  
The two Riemann sheets I and II, from the square root branch cut are defined to have Im and Im, respectively^{5}^{5}5The squareroot branch cut for is chosen in the conventional sense, such that it runs from the threshold to along the positive real axis. We neglect the effects from the left hand cut from s (singularity of at ) and from the freedom in choosing the real part of in Eq. (26) (such as using a ChewMandelstam phase space), either of which results in pole structure deep below the energy region being studied.. The product
(27) 
is an analytic function of only, where our lattice results for were expressed in terms of in Section V.1.
The scattering amplitude in Eq. (26) has a pole at a given if . A bound state appears as a pole of the matrix at , i.e. on the real axis below the scattering threshold on the first Riemann sheet. A virtual state appears as a pole at , i.e. on the real axis on the second Riemann sheet. In both cases denotes a real positive number. It is convenient to express these pole conditions in terms of the quantity , which is considered in the Section V.1 for . Such poles appear if the conditions
(28)  
(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.
v.3 Consistency check for a bound state in partial wave
The matrix for a bound state has a pole at , as discussed in the previous section. The analyticity of the matrix near the real energy axes implies the sign of the prefactor in front of the pole to be
(30) 
where is a real and positive number. This functional form is derived from the general form of the solutions for continued to complex in Ref. Sitenko (1971) (see Eq. (7.60)) and applied for wave in, for example, in Ref. Iritani et al. (2017). Equation (30) must be satisfied for all bound states, but it does not apply to virtual bound states.
Let us express the condition in Eq. (30) in terms of in Eq. (28) which was fitted from our data. For this purpose we consider the dependence of the quantity
(31) 
on near the bound state pole, where it equals zero in Eq. (28). Inserting in Eq. (31) with from Eq. (30), and expressing with as , we obtain
(32) 
Taking the derivative of the previous expression with respect to and evaluating it at , we get
(33) 
The above condition implies that the slope of in terms of has to be smaller^{6}^{6}6Note that this applies for the values of these quantities, not for absolute values of these quantities. than the slope of at the position of the bound state in partial wave , 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 wave was referred to as a sanitycheck in Ref. Iritani et al. (2017).
Let us now turn to the fitted phase shifts in Section V.1 and investigate which results satisfy the condition of Eq. (30) or equivalently the consistency check of Eq. (33). The slope of (the red curves in Fig. 9) has to be smaller than the slope of (the blue curve in Fig. 9) where these two curves cross. This is satisfied for the doublepole 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 has to be larger than the absolute value of slope for .
v.4 Bound states in the vector channel, and mass and coupling of the
The physical parameters for bound states and resonances are based on our fits of scattering in pwave. 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 doublepole fits in Section V.1 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 (Eqs. (26) and (27)) in the complex energy plane is shown in Figs. 10 and 11 for and , respectively. The location of the poles are related to (virtual) bound states and resonances according to the general criteria from Section V.2. The final masses of various states will be quoted according to
(34) 
where the mass splitting with respect to of Eq. (6) is determined from the lattice. Figure 12 and Table 3 summarize our results for various states, which are detailed below:

( MeV): For our lighter charm quark mass there is a complex conjugated pair of poles in the second sheet corresponding to the resonance, see Fig. 10(b). Figure 10(a) shows a single pole of the matrix on the real axis in the first Riemann sheet, corresponding to the bound state . In addition, there is a virtual bound state on the real axis on the second sheet.
Figure 10: The amplitude modulus of Eq. (2) for scattering in the vector channel plotted in the complex energy plane for . The bound state () is the pole on the real axis of the first Riemann sheet, while the resonance () appears on the second Riemann sheet as poles off the real axis above the scattering threshold. In the same position, a shoulder appears on the first sheet above threshold. The parameters for the resonance are extracted from the linearized (BreitWigner) behavior of Eq. (21) in the resonance region
(35) The resonance width can not be directly compared to experiment since the phase space (proportional to ) for our unphysical meson mass is different to that in experiment. Therefore we extract the dimensionless coupling that parametrizes the width. This coupling agrees within errors with the experimental value obtained from . The resonance mass equals of Eq. (18). The fitted value is given in Eq. (21) and (utilizing Eq. (34)) this yields a final mass MeV, which is consistent with the experimental value of MeV.

( MeV): 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 nearthreshold 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 MeV in Table 3 is determined according to Eq. (28) and corresponds to the nearthreshold crossing of the red and blue lines in Fig. 3(b). The cannot strongly decay into , but we still extract the coupling , which parametrizes the slope of near the real axes according to Eq. (35).
The corresponds to the lighter of the two bound states. Its mass MeV, determined from the boundstate condition of Eq. (28), has a large error due to the huge uncertainty of the coupling in the doublepole parametrization in Eq. (22). This mass is consistent within errors with the value MeV obtained directly from the firstexcited energylevel at and . 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 mass at this .
Our results for the masses of the and states are summarized and compared to experiment in Fig. 12 and Table 3. The mass and the coupling for the agree with the experimental values within errors, which applies to both charmquark masses considered. The mass of the is about below the experimental value, which is not unreasonable given the lack of continuum and quarkmass 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 subleading 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 meson mass, we see that our U101 data have a tendency to prefer a bound state (unlike the H105 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 H105 data alone. We find the mass and the coupling from such an analysis to be MeV and , which are in agreement with the results from our combined fits.
lat (present work)  lat (present work)  exp  lat Lang et al. (2015)  
[MeV]  MeV  
[MeV]  
[MeV]  
[MeV]  280  280  MeV  266  
resonance  bound st.  resonance Tanabashi et al. (2018)  resonance  
g  13.2(1.2)  
[MeV]  711(7)  707(7)  704.25(35)  715(7)  
[MeV]  9(7)  43(8)  38.52(35)  
[MeV]  3780(7)  3776(7)  3773.13(35)^{7}^{7}7We consider PDG fit estimate as the experimental value for the mass of the resonance throughout this article.  3784(7)  
bound st.  bound st.  bound st. Tanabashi et al. (2018)  bound st.  
[MeV]  597(10)  596(9)  617.347(25)  605(6)  
[MeV]  105(11)  154(10)  48.383(25)  
[MeV]  3666(10)  3665(9)  3686.097(25)  3674(6)  
resonance  resonance  resonance Aaij et al. (2019)  
[MeV]  762()  754()  773.9(2)  
[MeV]  59()  4()  108.2(2)  
[MeV]  3831  3822  3842.7(2) 