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

# Charmonium resonances with Jpc=1−− and 3−− from ¯DD scattering on the lattice

Stefano Piemonte Institute for Theoretical Physics, University of Regensburg, 93040 Regensburg, Germany    Sara Collins Institute for Theoretical Physics, University of Regensburg, 93040 Regensburg, Germany    Daniel Mohler Helmholtz-Institut Mainz, 55099 Mainz, Germany Johannes Gutenberg-Universität Mainz, 55099 Mainz, Germany    M. Padmanath Institute for Theoretical Physics, University of Regensburg, 93040 Regensburg, Germany    Sasa Prelovsek Institute for Theoretical Physics, University of Regensburg, 93040 Regensburg, Germany Faculty of Mathematics and Physics, University of Ljubljana, 1000 Ljubljana, Slovenia Jozef Stefan Institute, 1000 Ljubljana, Slovenia
July 2, 2019
###### Abstract

We present a lattice QCD study of charmonium resonances and bound states with and near the open-charm 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 non-perturbatively 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 charm-quark 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 charmonium-like states.

preprint: MITP/19-029

## 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 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 non-perturbative calculations, which can be performed using lattice QCD. Ab-initio 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 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 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 single-hadron approach, taking into account the strong decays of charmonium resonances above the open-charm threshold.111The 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 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 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 non-zero total momentum) can provide additional information on the relevant two-meson 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 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 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 well-established 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 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 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 finite-sized box with periodic boundary conditions.

We are interested in scattering, i.e. single-channel scattering of spin-zero particles with equal masses, in inertial frames with zero and non-zero total momenta . In this case the scattering matrix can be expressed in terms of a single phase shift in the partial-wave , 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 spin-less 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

 det[~K−1l(Ecm) δl′l−B→P;Λl′l(Ecm)]=0. (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

 (ρ tl)−1=2i (Sl−1)−1=p−2l−1~K−1l−i, (2)

with , or equivalently

 ~K−1l(Ecm)=p2l+1cotδl(p). (3)

The kinematic variable is the momentum of each meson in its center-of-momentum frame

 p2=s4−m2D, (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)).222However, 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 box-matrix , 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 off-diagonal entries in the partial wave indices and , 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 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 well-known Lüscher relation Lüscher (1986)

 pcotδ1(p)−2Z00(p2)√πL=0, (5)

for zero total momentum.

In the case when the effects of higher partial waves are non-negligible, 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 near-threshold 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 Wilson-clover action  Sheikholeslami and Wohlert (1985) is non-perturbatively 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 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 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 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 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 ,

 Mav=14(mηc+3mJ/ψ). (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 spin-average mass is equal to 3103(3) MeV for and 2820(3) MeV for . The energy gap between and 2 is

 Mav−2mD = −750(3)~{}MeV    κc=0.12315 Mav−2mD = −703(3)~{}MeV    κc=0.12522 Mav−2mD = −665~{}MeV    ~{}~{}~{}~{}experimental

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 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. 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 time-source/sink. Therefore, we use stochastic distillation Morningstar et al. (2011) for the light quark propagators on the H105 ensemble only. The perambulator, defined as

 τ(t′,μ′,i′;t,μ,i)=⟨vμ′i′(t′)|(D−1|vμi(t)⟩), (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

 kn(t,i,μ)=1√2(±1±i), (8)

using full-time dilution. The stochastic source given to the multigrid solver is constructed for every time-slice as

 |wn(t)⟩=∑j,νkn(t,j,ν)|vνj(t)⟩. (9)

The perambulator can be approximated as

 τ(t′,μ′,i′;t,μ,i)=1Ns∑nk†n(t,i,μ)⟨vμ′i′(t′)|(D−1|wn(t)⟩), (10)

which is equal to the exact expression of Eq. (7) in the limit , since

 1NsNs∑n=1k†n(t,i,μ)kn(t,j,ν)≃δijδμν. (11)

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 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 two-point correlation functions

 Cij(t)=⟨Oi(t)O†j(0)⟩=∑nZniZn∗je−tE{lat}n, (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) ΛPC=T−−1 in the rest frame |→P|2=0, ii) ΛC=A−1 in the moving frame |→P|2=1 and iii) ΛC=A−1 in the moving frame |→P|2=2,

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 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. 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 single-meson 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 two-meson interpolators are built using the projection

 O¯DD(→p1+→p2) = ∑R∈G(→p1+→p2)TΛr,r(R)∗{¯D(R→p1)D(R→p2) −¯D(R→p2)D(R→p1)}, (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 two-meson 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 two-meson interpolators in the basis brings additional Wick contractions. We compute all relevant Wick contraction diagrams for the single- and two-meson 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)

 Cij(t)vin(t) = λn(t)Cij(t0)vin(t), λn(t) ∝ e−tE{lat}n(1+O(e−ΔEnt)) (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 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. 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 non-negligible deviations of the dispersion relation of the -meson from its continuum counterpart

 EcontD(→p)=√m2D+|→p|2. (15)

As is evident from Table 2, there are non-negligible 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 non-interacting state

 (ΔE)s=(Elat)s−(ElatD(→p1))s−(ElatD(→p2))s, (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,333The Lüscher -functions, that enter the box matrix in the quantization condition, have poles at energies with and , while in the Lüscher-type 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

 (Ecalc)s=(ΔElat)s+(EcontD(→p1))s+(EcontD(→p2))s (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 phase-shift when in Eq. (16) is zero, even if non-negligible heavy-quark discretization effects are present in our simulation.

## V Results for charmonia with Jpc=1−− and 3−−

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 near-threshold 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 charm-annihilation.

Our second aim is to study the lowest charmonium resonance. A candidate for this state referred to as the 444The 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 “double-pole”

 p3cot(δ1)√s=(G21m21−s+G22m22−s)−1, (18)

 p3cot(δ1)√s=A+Bs+Cs2, (19)

where . Other parameterizations were investigated including, for example, linear and triple-pole 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 spin-identified 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 low-lying partial wave excitation on the rest of the discrete spectrum, we parametrize it with a simple Breit-Wigner form given by

 p7cot(δ3)√s=m23−sg23. (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

 p2l+1cot(δ)√s=⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩([0.63(33)]2[1.4966(30)]2−s+[3.69(37)]2[1.5457(32)]2−s)−1l=1[1.568(11)]2−s[0.07(3)]2l=3 (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

 p2l+1cot(δ)√s=⎧⎪⎨⎪⎩−0.60(+20−24)+0.55(+20−17)s−0.127(+38−44)s2  l=1[1.5639(+83−14)]2−s[0.058(+94−55)]2  l=3 (23)

with a , while for

 p2l+1cot(δ)√s=⎧⎪⎨⎪⎩−0.08(+11−13)+0.11(+10−8)s−0.028(+15−17)s2  l=1[1.6869(+33−32)]2−s[0.042(+53−22)]2  l=3 (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,

 Ω(μ,A(s))=det(A)det((μ2+AA†)12), (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

 tl(s)=1ρ cot(δl)−i ρ, (26) ρ=2p√s=√1−4m2Ds.

The two Riemann sheets I and II, from the square root branch cut are defined to have Im and Im, respectively555The square-root 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 Chew-Mandelstam phase space), either of which results in pole structure deep below the energy region being studied.. The product

 ρcotδl=2p2l [p2l+1 cotδl√s] (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

 p2l+1cot(δl) = −(p2)l√−p2  (bound state) (28) p2l+1cot(δl) = (p2)l√−p2  (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.

### v.3 Consistency check for a bound state in partial wave l

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 pre-factor in front of the pole to be

 S(p)=−(−1)l i β2p−iκB, (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

 p2l+1cot(δl)−[−(p2)l√−p2] (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

 p2l+1cot(δ)−[−(p2)l√−p2] = −2(−1)l(−p2)l(κB√−p2+p2)(−1)lβ2−κB+√−p2. (32)

Taking the derivative of the previous expression with respect to and evaluating it at , we get

 ddp2(p2l+1cot(δ)−[−(p2)l√−p2])∣∣∣p2=−κ2B=−κ2lBβ2<0. (33)

The above condition implies that the slope of in terms of has to be smaller666Note 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 sanity-check 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 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 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 ψ(3770)

The physical parameters for bound states and resonances are based on our fits of 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.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 double-pole 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

 m=mlat−Mlatav+Mexp% av, (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.

The parameters for the resonance are extracted from the linearized (Breit-Wigner) behavior of Eq. (21) in the resonance region

 p3cotδ1√s|s≃m2=6πg2(m2−s) ,Γ=g2 p36π s. (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.

The mass of the bound state is obtained from of Eq. (21) by finding a pole for real energies below threshold according to Eq. (28): it corresponds to the crossing of the red and blue curves in Fig. 9(a). The resulting mass in Table 3 is about MeV smaller than the experimental mass.

• ( 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 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  MeV in Table 3 is determined according to Eq. (28) and corresponds to the near-threshold 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 bound-state condition of Eq. (28), has a large error due to the huge uncertainty of the coupling in the double-pole parametrization in Eq. (22). This mass is consistent within errors with the value MeV obtained directly from the first-excited energy-level 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 charm-quark masses considered. The mass of the is about 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 -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.