Charmonium-like resonances with J PC = 0 ++ , 2 ++ in coupled D D ¯ , D s D ¯ s scattering on the lattice
S. Prelovsek, S. Collins, D. Mohler, M. Padmanath, S. Piemonte
PPrepared for submission to JHEP
MITP/20-065
Charmonium-like resonances with J P C = 0 ++ , ++ incoupled D ¯ D , D s ¯ D s scattering on the lattice Sasa Prelovsek a,b,c
Sara Collins a Daniel Mohler d,e,f
M. Padmanath d,e
StefanoPiemonte a a Institute for Theoretical Physics, University of Regensburg, 93040 Regensburg, Germany b Faculty of Mathematics and Physics, University of Ljubljana, Slovenia c Jozef Stefan Institute, Ljubljana, Slovenia d Helmholtz-Institut Mainz, Johannes Gutenberg-Universit¨at, D-55099 Mainz, Germany e GSI Helmholtzzentrum f¨ur Schwerionenforschung, 64291 Darmstadt, Germany f Johannes Gutenberg-Universit¨at Mainz, D-55099 Mainz, Germany
E-mail: [email protected] , [email protected] Abstract:
We present the first lattice investigation of coupled-channel D ¯ D and D s ¯ D s scattering in the J P C = 0 ++ and 2 ++ channels. The scattering matrix for partial waves l = 0 , m π (cid:39) a (cid:39) .
09 fm and
L/a = 24 ,
32 are utilized. The resulting scattering matrix suggeststhe existence of three charmonium-like states with J P C = 0 ++ in the energy region rangingfrom slightly below 2 m D up to 4 .
13 GeV. We find a so far unobserved D ¯ D bound statejust below threshold and a D ¯ D resonance likely related to χ c (3860), which is believed tobe χ c (2 P ). In addition, there is an indication for a narrow 0 ++ resonance just below the D s ¯ D s threshold with a large coupling to D s ¯ D s and a very small coupling to D ¯ D . Thisresonance is possibly related to the narrow X (3915)/ χ c (3930) observed in experiment alsojust below D s ¯ D s . The partial wave l = 2 features a resonance likely related to χ c (3930).We work with several assumptions, such as the omission of J/ψω , η c η and three-particlechannels. Only statistical uncertainties are quantified, while the extrapolations to thephysical quark-masses and the continuum limit are challenges for the future. Keywords:
Charmonium, Lattice QCD, hadron scattering a r X i v : . [ h e p - l a t ] N ov ontents E cm , Riemann sheets and poles 4 D ¯ D scattering with l = 0 near threshold 128.2 D s ¯ D s scattering with l = 0 near threshold in the one-channel approximation 148.3 D ¯ D scattering with l = 2 and J P C = 2 ++ resonance 168.4 Coupled D ¯ D , D s ¯ D s scattering with l = 0 for E cm (cid:39) . − .
13 GeV 178.4.1 Analysis omitting l = 2 178.4.2 Analysis including l = 0 and l = 2 198.5 Coupled D ¯ D , D s ¯ D s scattering in a wider energy region 21 C.1 D ¯ D scattering with l = 0 near threshold 30C.2 D s ¯ D s scattering with l = 0 near threshold in the one-channel approximation 30C.3 D ¯ D scattering with l = 2 30C.4 Coupled D ¯ D , D s ¯ D s scattering with l = 0 for E cm (cid:39) . − .
13 GeV 30 D Ω( E cm ) in a wider energy region E cm (cid:39) m D − . GeV 31 – 1 –
Introduction
Since the discovery of the
J/ψ meson in 1970s a multitude of charmonium bound statesand resonances have been found with energies ranging up to almost 5 GeV. A simple c ¯ c quark model provides a reasonable description of the levels below the strong decay thresh-olds and also some of the states above, however, there are clearly too many states to fitinto this picture. Some mesons, such as the charged Z c states certainly have additionalquark content, while for other states the interpretation is not so clear. On the theory sidethe nature of these states is being explored in tetraquark, molecular, and hybrid mesonmodels, among others, while on the experimental side insight is provided by establishingtheir quantum numbers, decay modes and widths. Lattice QCD studies of the charmo-nium spectrum have a significant role to play in terms of guiding experimental searches,determining the quantum numbers of the states not well established as well as investigatingtheir internal structure.In this work we focus on the isoscalar channel I ( J P C ) = 0(0 ++ ) in the region up to4.13 GeV for which there are a number of open questions. The ground state, χ c (1 P ),found well below the D ¯ D threshold is interpreted as the P c ¯ c level of the quark modeland is the only well established state. In the energy region around 3.9 GeV, above thethreshold, one expects a corresponding excited state. So far, three hadrons have beenobserved with the possible assignment of J P C = 0 ++ : the X (3860), a broad resonancedetected by Belle [1, 2], and two narrow resonances just below the D s ¯ D s threshold — the χ c (3930) discovered in the D ¯ D channel by LHCb [3, 4] and the X (3915) observed throughit’s decay into J/ψω [5–8] (with the assignment of J P C = 0 ++ or 2 ++ ). While the latter tworesonances could be the same state, their narrowness may indicate exotic content, where X (3915) has been interpreted as c ¯ cs ¯ s meson in Ref. [9]. Predictions have also been madefor an additional, as yet unobserved, bound state just below the D ¯ D threshold [10, 11].The determination of the low lying charmonium spectrum on the lattice is relativelystraightforward, with the energy levels being directly accessed from correlation functionsmeasured on the configurations generated in the Monte-Carlo simulation. Systematics aris-ing from finite lattice spacing and simulating with unphysical light (sea) quark masses mustbe addressed by carrying out a continuum and quark-mass extrapolation. Near and abovethreshold, the analysis is considerably more challenging with information on the massesand (for resonances) also the widths being inferred from scattering amplitudes which canbe obtained from the finite volume spectra via the L¨uscher method [12–14]. Two-particleinterpolators must be included in the basis of operators for the construction of the correla-tion functions in order to reliably determine these spectra. Simulating charmonia in flightprovides additional levels with which to probe the scattering matrix, however, the identifi-cation of the continuum spin and parity quantum numbers of the levels is complicated dueto the reduced symmetry on the lattice. In addition, for the energy range of interest boththe D ¯ D and D s ¯ D s thresholds must be considered leading to a coupled-channel scatteringanalysis.So far, the coupled-channel scattering matrix has been extracted for several light-meson systems, for example, πK, ηK [15, 16], πη , K ¯ K [17] and ππ, K ¯ K, ηη [18] by the– 2 –adron Spectrum Collaboration. In the heavy sector, there has been one investigationof
Dπ, Dη, D s ¯ K scattering in isospin-1/2 [19] and a recent analysis of the Z c (3900) via D ∗ ¯ D, J/ψπ scattering [20]. The HALQCD Collaboration has also investigated the Z c (3900)using a different approach which involves solving the Schr¨odinger equation with potentialsdetermined on the lattice [21, 22]. Pioneering works such as these were limited to a singlelattice spacing and unphysical light-quark masses.The charmonium scalar channel has previously been studied by some of the authorsconsidering only D ¯ D scattering with total momentum zero [23]. Here we present a latticestudy of scattering in the coupled-channels D ¯ D and D s ¯ D s with quantum numbers I = 0 and J P C = 0 ++ , ++ . This represents the first determination of the coupled-channel scatteringmatrix from lattice QCD in the charmonium system with isospin zero. Two lattice volumesare employed for the charmonium system at rest and in flight. This analysis uses the samelattice setup as our previous article on the identification of the spin and parity of the singlehadron spectrum [24] and the investigation of single channel D ¯ D scattering for J P C = 1 −− and 3 −− [25]. While the present study represents a significant improvement on previouswork, some simplifications remain and a comparison of the results for the masses and widthswith experiment is qualitative. Within the energy range of interest, additional scatteringchannels, such as the J/ψω , η c η and those involving three particles, could in principle alsobe relevant. The effects of these channels will be investigated in the future, along withsystematics associated with finite lattice spacing and unphysical light quark masses.The remainder of the paper is organized as follows. We begin by reviewing the essentialgeneral aspects of one-channel and two-channel scattering in Section 2. The details ofthe lattice setup and methodology are given in Section 3 and the single- and two-mesoninterpolators used in the correlation functions are discussed in Section 4. Simplifyingassumptions made in this study are summarized in Section 5. The first step in extractingthe scattering amplitudes is to compute the finite-volume spectra from the correlationfunctions. Our analysis and the final spectra are presented in Section 6. An overviewof determining the scattering amplitudes from the lattice eigen-energies is provided inSection 7. Our results for the J P C = 0 ++ and 2 ++ channels are detailed in Section 8 andthe relation to states observed in experiment is discussed in Section 9. Finally, Section 10presents our conclusions. The masses and widths of strongly-decaying resonances should be inferred from the studyof scattering processes where these resonances appear. In this section, we briefly reviewrelevant concepts regarding scattering matrices, complex energy planes, pole singularities,hadron masses, and their decay widths. The first part lists definitions and notations forthe scattering amplitudes, the phase space factors, etc. . The second part discusses namingconventions for various Riemann sheets, pole singularities in the complex energy plane andtheir relation to the hadron properties. – 3 – .1 Scattering matrices for real energies
The unitary scattering amplitude S for one-channel scattering ( D ¯ D or D s ¯ D s ) of spin-lessparticles in partial wave l is generally parametrized in terms of the energy-dependent phaseshift δ ( E cm ), S = 1 + 2 i θ ( E cm − m ) ρ t = e iδ , (2.1) t − = cot δ ρ − i ρ = 2 E cm p l ˜ K − − i ρ with ˜ K − = p l +1 cot δ , where ρ ≡ p/E cm , p denotes the momentum of the scattering particles in the center-of-momentum frame and t is the scattering amplitude. The Heaviside function θ ( E cm − m )ensures that S is trivial below the threshold. The factors p − l in front of ˜ K − lead tosmooth behavior close to the threshold. In the case of t exhibiting simple Breit-Wignertype behavior, ˜ K − /E cm falls linearly as a function of E cm , t ( E cm ) = − E cm Γ( E cm ) E cm − m R + i E cm Γ( E cm ) , Γ( E cm ) = g p l +1 E cm , ˜ K − E cm = m R − E cm g . (2.2)The phase shift equals π/ E cm = m R , while the width Γ( E cm ) is parametrized in termsof the coupling g and the phase space. S , t , ˜ K and δ depend on E cm and partial wave l (the dependence on l is not written explicitly).For coupled-channel scattering of D ¯ D and D s ¯ D s in partial wave l , the scattering ma-trices S are energy-dependent 2 × S ij = δ ij + 2 i Θ ij √ ρ i ρ j t ij = (cid:40) η e iδ i i = j (cid:112) − η e i ( δ i + δ j ) i (cid:54) = j , ( t − ) ij = 2 E cm p li p lj ( ˜ K − ) ij − i ρ i δ ij , (2.3) ρ i ≡ p i /E cm = (cid:112) − (2 m i ) /E cm , i, j = 1 ( D ¯ D ) , D s ¯ D s ) . Θ ij = θ ( E cm − E th,i ) θ ( E cm − E th,j ) ensures that the elements of S differ from unity onlyabove the respective thresholds. The momenta of D and D s in the center-of-momentumframe are denoted by p and p , respectively. t is the scattering matrix and ˜ K ( E cm ) is areal symmetric matrix. We follow the definition of t by the Hadron Spectrum Collaborationand the definition of ˜ K from Ref. [26] . E cm , Riemann sheets and poles In experiment and lattice QCD simulations the scattering matrices S ( E cm ) are determinedfor real energies. The theoretical interpretation in terms of (virtual) bound states and reso-nances is conventionally performed by looking at poles of S ( E cm ), analytically continued tothe complex energy-plane. The feature that leads to interesting physics is the square root Note that unlike in Ref. [26] we do not divide our energy levels or physical quantities by the mass ofthe scattering particles. – 4 – II m I m ( E c m ) Re( E cm ) I II II III m m I m ( E c m ) Re( E cm ) Figure 1 : Sketch of the pole locations in the scattering matrix t that typically affectthe experimental rates on the physical axes (denoted by the cyan line) for one-channelscattering (left) and two coupled channels (right). The Roman numbers indicate the Rie-mann sheets where the poles are located according to Eq. (2.4). Poles immediately below athreshold, indicated by crosses, can also have observable effects on the physical axes abovethe respective threshold.branch cut related to ρ = 2 p/E cm = (cid:112) − (2 m ) /E cm starting from the threshold con-necting the physical Riemann sheet (or sheet I), conventionally chosen to have Im( ρ ) > ρ ) <
0. For a two channelsystem, there will be four Riemann sheets, such thatsheet I : Im( ρ D ) > , Im( ρ D s ) > , sheet II : Im( ρ D ) < , Im( ρ D s ) > , ( ρ i = 2 p i /E cm ) , sheet III : Im( ρ D ) < , Im( ρ D s ) < , sheet IV : Im( ρ D ) > , Im( ρ D s ) < . (2.4)Bound states, virtual bound states and resonances are related to pole singularities of t in the complex energy plane. These poles affect the physical axes, indicated by the cyan linein Fig. 1, along which the experimental measurements are made. Fig. 1 presents a schematicpicture of various pole locations in our study, that can affect scattering amplitudes/matricesalong the physical axes for one-channel and two-channel scattering. The location of thepoles are related to the masses and widths via E pcm = m − i Γ for resonances and E pcm = m for the (virtual) bound states.In the close vicinity of the pole, the scattering matrix has the energy dependence( t − ) ij ∼ c i c j ( E pcm ) − E cm for E cm (cid:39) E pcm , i = 1 ( D ¯ D ) , D s ¯ D s ) , (2.5)and the residue can typically be factorized into the couplings c i ( ), whose relative size isrelated to the branching ratios of a resonance (associated with the pole) to both channels i = 1 , We employ two ensembles generated with N f = 2 + 1 non-perturbatively O ( a ) improvedWilson dynamical fermions provided by the Coordinated Lattice Simulations (CLS) con- The couplings c i should not be confused with the coupling g (2.2) that will be used to parametrize thefull width of a resonance. – 5 – π [MeV] m K [MeV] m D [MeV] m D s [MeV] aM av [MeV]280(3) 467(2) 1927 (2) 1981(1) 3103(3) Table 1 : Hadron masses in physical units for the gauge configurations used in this project,where M av = ( m η c + 3 m J/ψ ) /
4. The hadrons containing charm quarks are from κ c =0 . m u/d,s are chosen along a trajectory that approachesthe physical point holding the average quark mass, 2 m u/d + m s , constant. The ensembles,denoted H105 and U101, have an inverse gauge coupling β = 6 /g = 3 .
4, corresponding toa lattice spacing a = 0 . N T × N L = 128 × and 96 × ,respectively [29]. Open boundary conditions in time are imposed [30] and the sources ofthe correlation functions are placed in the bulk away from the boundary. We remark thatthese correlation functions do not show any effects related to the finite time extent in thetime regions analyzed. For H105 we use replica r001 and r002 for which the issue of neg-ative strange-quark determinants described in Ref. [31] is not of practical relevance. Forour analysis we use 255 (492) configurations on two replicas for ensemble U101 (H105).The masses of the pion, kaon, D and D s mesons determined on the larger ensembleare shown in Table 1. Note that the chosen quark-mass trajectory leads to a larger thanphysical m u/d and a smaller than physical m s . This means that the splitting betweenthe D ¯ D and D s ¯ D s thresholds is smaller than in experiment, emphasizing the need for acoupled-channel analysis. We employ the charm-quark hopping parameter κ c = 0 . m c and spin-averaged 1 S -charmonium mass M av thatare slightly larger than their physical values. For estimates of the statistical uncertaintywe use the bootstrap method with (asymmetric) error bars resulting from the central 68%of the samples. Further details are collected in Appendix A. The correlation matricesare averaged over several source-time slices and momentum polarizations to increase thestatistical precision. Note that all quoted uncertainties are statistical only, and that resultsquoted in MeV have been obtained using the central value of the lattice scale withoutpropagating its statistical or systematic uncertainties into the results.For hadrons with charm quarks, non-negligible discretization effects are observed whencomputing the dispersion relation on lattices with a ≈ .
086 fm. The dispersion relationdeviates from the continuum one, however, the latter is assumed in the finite-volume anal-ysis outlined in Section 7. In order to reduce the systematic uncertainty arising fromdiscretization effects we adopt the analysis strategy described in Sect. IV.B. of Ref. [25].Here we only reiterate the most important details.First the energy shift of each interacting eigenstate with respect to a nearby non-interacting two-hadron level H ( (cid:126)p ) H ( (cid:126)p ) is computed(∆ E lat ) s = ( E lat ) s − ( E lat H ( (cid:126)p ) ) s − ( E lat H ( (cid:126)p ) ) s , (3.1)where (cid:126)p , = (cid:126)n , πL , (cid:126)p + (cid:126)p = (cid:126)P and s denotes the bootstrap sample. Here, ( E lat ) s isthe energy of the interacting two-hadron system, while ( E lat H i ( (cid:126)p i ) ) s is the energy of a single– 6 –adron (either D or D s meson in this paper) with momentum (cid:126)p i measured on the lattice.We then use ( E calc ) s = (∆ E lat ) s + (cid:16) E cont H ( (cid:126)p ) (cid:17) s + (cid:16) E cont H ( (cid:126)p ) (cid:17) s (3.2)as input to the quantization condition (see Eq. (7.1)) for each bootstrap sample s . Theenergies ( E cont H i ( (cid:126)p i ) ) s are computed from the continuum dispersion relation using the latticemomenta (cid:126)p , and the single-hadron ( D and D s ) masses at rest. The resulting energies E calc are equal to E lat in the naive continuum limit a → The main aim of this work is to investigate the coupled-channel D ¯ D - D s ¯ D s scattering ampli-tudes and cross-sections in the channel I ( J P C ) = 0(0 ++ ) in the energy range encompassingthe D ¯ D threshold up to 4 .
13 GeV. Following L¨uscher’s approach [12–14, 34], this requiresa reliable extraction of the finite-volume charmonium spectrum below 4.13 GeV on severaldifferent volumes and/or in different momentum frames. In this study, we consider thecharmonium spectrum in four different lattice irreducible representations (irreps) Λ:i) Λ P = A +1 , | (cid:126)P | = 0 , J P [ λ ] = 0 + [0] , ii) Λ = A , | (cid:126)P | = 1 , J P [ λ ] = 0 + [0] , + [0] , iii) Λ = A , | (cid:126)P | = 2 , J P [ λ ] = 0 + [0] , + [0] , ± [2] , iv) Λ = B , | (cid:126)P | = 1 , J P [ λ ] = 2 ± [2] . (4.1)The squared momenta | (cid:126)P | in the lab frame are given in units of (2 π/L ) . Chargeconjugation C = + is a good quantum number in all frames and hence is suppressedfor brevity. On the right of Eq. (4.1), we list all relevant states with quantum numbers J P [ λ ] contributing to the respective irreps. Here λ refers to the helicity of the state. Thefirst three irreps are relevant for an investigation of the J P = 0 + channel. The irrepsin the moving lab frames also receive contributions from states with J P [ λ ] = 2 + [0] and2 ± [2] within the energy range of interest . The analysis of the spectrum in the B irrepconstrains the parameters for D ¯ D scattering with l = 2. This partial wave inevitablycontributes to the finite-volume spectrum of irreps A with P >
0. We utilize a large basisof single-meson as well as two-meson interpolators in the above irreps to reliably determinethe relevant low energy spectrum.As in our previous publications [24, 25], we construct the single-meson interpolatorsfollowing the procedure in Refs. [35, 36], using up to two gauge covariant derivatives.Table 2 lists the number of single-meson operators employed in each of the finite-volumeirreps considered. The procedure discussed in Ref. [24] guides us in assigning the quantum States with other J P , such as 3 + , contribute at higher energies. We assume these higher lying statesto have negligible influence in the energy range considered. – 7 – (cid:126)P | Λ ( P ) C O N ops | (cid:126)P | Λ ( P ) C O N ops O h A ++1 ¯ cc Dic B +1 ¯ cc D (0) ¯ D (0) 2 D (2) ¯ D (1) 1 D (1) ¯ D (1) 1 Dic A +1 ¯ cc D s (0) ¯ D s (0) 2 D (2) ¯ D (0) 2 D ∗ (0) ¯ D ∗ (0) 2 D (1) ¯ D (1) 2 J/ψ (0) ω (0) 2 D (2) ¯ D (2) 1 Dic A +1 ¯ cc D (3) ¯ D (1) 1 D (1) ¯ D (0) 2 D s (2) ¯ D s (0) 1 D (2) ¯ D (1) 1 D s (1) ¯ D s (1) 1 D s (1) ¯ D s (0) 2 J/ψ (2) ω (0) 3 J/ψ (1) ω (0) 2 Table 2 : The single- and two-meson interpolators utilized in each lattice irrep Λ ( P ) C considered in this study. We use the simplified notation M ( (cid:126)p ) M ( (cid:126)p ) for the two-mesoninterpolators with the momentum (cid:126)p i of each meson ( i = 1 ,
2) given in units of 2 π/L . Thefull expressions are omitted for brevity. N ops indicates the number of operators of eachtype employed.numbers J P [ λ ] to the extracted energy levels and aids us in selecting the levels relevant forthe amplitude analysis.The D ¯ D as well as D s ¯ D s interpolators are constructed following the same procedureas in Ref. [25]. The momentum combinations implemented in this study are given inTable 2. We also include two-meson operators involving spin 1 mesons, such as J/ψω and D ∗ ¯ D ∗ (see Table 2). For non-zero momenta, the construction of such operators needsadditional care and we follow the induced representation method described in Appendix A2of Ref. [37]. In the | (cid:126)P | = 2 frame, for example, we implement three linearly independent J/ψ (2) ω (0) operators and observe three almost-degenerate eigenstates. These operatorsare not included when extracting the finite volume spectrum for the amplitude analysis,as discussed in Section 5. This study is performed using lattice gauge ensembles with two different physical volumesat a single lattice spacing and at unphysical quark masses (the resulting masses of keyhadrons are given in Table 1). As a consequence, only a qualitative comparison of theresults can be made with experiment.Unlike for light hadrons [38], scattering studies in the charmonium sector are stillat an early stage. For the physical states we are interested in, a three-particle channeland multiple two-particle channels are open and all could, in principle, be relevant. Onepossible approach is to simulate at very heavy pion (and kaon) mass, such that the numberof relevant decay modes is reduced to a few two-hadron modes, which can then be fullyexplored. This approach has the disadvantage that the quark masses are far removed– 8 –rom their physical values, making a comparison to experiment a challenge. We opt for astrategy where we simulate at a moderate pion mass of 280 MeV and take into accountthe scattering channels expected to be most relevant for the physics close to the open-charm threshold(s). Some additional channels are neglected (as discussed below), however,our assumptions about which thresholds are relevant can be relaxed successively in futurecalculations.Beyond the scattering channels investigated explicitly, our current study includes
J/ψω and D ∗ ¯ D ∗ operators in the interpolator basis. Due to the poor signal obtained for lightisoscalar mesons, the energy levels close to the non-interacting J/ψω levels are not veryprecisely determined and would not provide strong constraints on the scattering matrix.In particular, almost all energy levels dominated by
J/ψω interpolators fall within onestandard deviation of the non-interacting
J/ψω energies, and — apart from the additionalenergy levels which appear — the other finite-volume energies do not shift significantlywhen including these interpolators. Section 6 will present the finite volume spectrum up to4.13 GeV based on all the operator types in Table 2, apart from the
J/ψω operators (seeFig. 2). Note that for our lattices the
J/ψω threshold is located at approximately 3 .
95 GeV.We also neglect the η c η channel which has a threshold of around 3 .
54 GeV. We remarkthat this decay channel has not been observed for any of the experimental candidatesmentioned in the introduction (and discussed in more detail in Sect. 9). Operators withmore than two hadrons are also not implemented. The lowest three-hadron threshold is forthe decay into χ c ππ at 4 .
02 GeV. This threshold is close to the upper end of the energyregion that we consider.As in all studies of charmonium-like resonances to date, charm annihilation Wickcontractions are omitted. All the remaining contraction diagrams arising from the single-and two-meson operators in our basis (shown in Fig. 1 of Ref. [23]) are computed followingthe procedures described in our previous publications [24, 25].We stress that we determine the finite-volume spectra at a single lattice-spacing andare therefore unable to quantify the uncertainty associated with the lattice discretization.In particular, the uncertainty arising from the heavy quark discretization may be non-negligible. As discussed in the previous section, the dispersion relation deviates from thecontinuum relation in our study and spin-splittings are also likely to be affected [39, 40].In general, lattice spacing effects in heavy-light mesons and charmonium are different withthe net result that even at physical light-quark masses the open-charm thresholds can beshifted with respect to the measured charmonium states at finite lattice spacing.
This section presents the eigen-energies E n that will be used to determine the scatteringmatrices. The energies are obtained from the correlation matrices C ij ( t ) = (cid:104) O i ( t ) O † j (0) (cid:105) via the widely-used variational method. This involves solving the generalized eigenvalueproblem C ( t ) u ( n ) ( t ) = λ ( n ) ( t ) C ( t ) u ( n ) ( t ) for the eigenvalues λ ( n ) ( t ) and the eigenvectors u ( n ) ( t ) [41–43]. We use the reference time t /a = 3 or 4. The eigen-energies are extracted– 9 – / aE cm [ GeV ] , A + P =
24 323.43.53.63.73.83.94.4.1 L / aE cm [ GeV ] , A P =
24 323.43.53.63.73.83.94.4.1 L / aE cm [ GeV ] , A P =
24 323.43.53.63.73.83.94.4.1 L / aE cm [ GeV ] , B P = Figure 2 : The eigen-energies in the center-of-momentum frame ( E cm ) for the charmonium-like system with I = 0 and C = +1. Results are presented for irreducible representationsΛ P = A (+)1 , B and total momenta | (cid:126)P | = 0 , ,
2, which give information on the channelswith J P C = 0 ++ , ± + . The data points correspond to the eigen-energies obtained fromthe lattice simulation: the black circles are used to extract the coupled-channel scatteringmatrices for D ¯ D − D s ¯ D s , while the blue circles are omitted from the scattering analysis.The solid and dashed red (green) lines correspond to discrete D ¯ D ( D s ¯ D s ) eigen-energiesin the non-interacting limit: solid lines correspond to the operators that are implemented,while dashed lines correspond to the operators that are not implemented. Dotted linesrepresent thresholds. The data points indicated by the light blue circles correspond toground-state charmonia with J P C = 2 ++ and 2 − + , which appear at m (cid:39) .
56 GeV and3 .
83 GeV, respectively. Some points are shifted horizontally slightly for clarity.from 1-exponential fits to the eigenvalues λ ( n ) ( t ) = A n e − E n t with the fit range, in mostcases, starting between timeslices 10 and 12.The finite-volume spectrum of the charmonium system with isospin I = 0 and C = +1is shown in Fig. 2. We present the spectrum in the center-of-momentum (cm) frame E cm,n = [ E n − (cid:126)P ] / ( ) for irreducible representations Λ P = A (+)1 , B and total momenta | (cid:126)P | = 0 , ,
2. These irreps give information on the charmonium(like) states and D ( s ) ¯ D ( s ) scattering in the channels with J P C = 0 ++ , ± + (see Eq. 4.1). The energies indicated bythe black-circles are used to extract information on D ( s ) ¯ D ( s ) scattering. These energies arenear or above the D ¯ D threshold and are precise enough to reliably resolve the energy-shifts Here, the energy in the lattice frame E n stands for E calcn obtained according to Eq. (3.2) or Eq. (17)of Ref. [25]. – 10 –ith respect to the non-interacting energies of D ( s ) ¯ D ( s ) (indicated by the solid lines). Thelight-blue circles are the energy levels related to ground-state charmonia with J P = 2 ± . The bound states and resonances are inferred from the scattering matrices as briefly re-viewed in Section 2. The infinite-volume scattering matrix S ( E cm ) is related to the finite-volume two-hadron spectrum for real energies E cm above the threshold and somewhatbelow it through the well-known L¨uscher relation [12–14]. The eigen-energies of the cou-pled channel D ¯ D − D s ¯ D s system given in the previous section provide information on the2 × S ( E cm ) for these coupled channels via the generalization of this for-malism [34, 44, 45], considered (for other channels), for example, in Refs. [16, 17, 19]. The S matrix can be expressed in terms of a real function ˜ K ( E cm ) for one-channel scattering(2.1) and a real symmetric 2 × K ( E cm ) for two coupled channels (2.3), as detailedin the next section. ˜ K uniquely determines S , while both depend also on the partial wave l . We use the spectrum from the previous section to determine ˜ K ( E cm ) using the publiclyavailable package TwoHadronsInBox [46].The relation between discrete lattice eigen-energies E cm and ˜ K -matrix for coupled-channel scattering is referred to as the quantization condition [46]det[ ˜ K − l ; ij ( E cm ) δ ll (cid:48) − B (cid:126)P , Λ ll (cid:48) ; i ( E cm ) δ ij ] = 0 . (7.1)Both terms in the determinant are matrices in the space of partial waves l, l (cid:48) and channels i, j ( D ¯ D , D s ¯ D s or both), and the determinant is evaluated over both indices. ˜ K l ; ij δ ll (cid:48) isan unknown matrix in channel space that depends on the partial wave l ; it is diagonal in l since the good quantum numbers in continuum scattering of spin-less particles (such as D and D s ) are J , S = 0 and l = J − S = J . The B (cid:126)P , Λ ll (cid:48) ; i ( E cm ) are known box-functions [46] thatare in general non-diagonal in the partial wave index.In one-channel scattering and when only partial wave l contributes, relation (7.1) sim-plifies to ˜ K − ( E cm ) = B (cid:126)P , Λ ( E cm ), since the argument of the determinant is a 1 × K − ( E cm ) will be shown as points in figures for one-channel scattering. Fortwo coupled channels, for the case when only partial wave l contributes, the determinantequation (7.1) provides one relation between ˜ K ( E cm ), ˜ K ( E cm ) and ˜ K ( E cm ) for eachenergy level, complicating the determination of those functions. Therefore, we follow thestrategy proposed in Ref. [44], where the ˜ K ij ( E cm ) are parametrized as functions of theenergy. In this strategy, the ˜ K -matrix elements are determined by requiring that relation(7.1) is simultaneously satisfied for all relevant lattice energies E cm .We will focus on certain interesting and rather narrow energy regions, where a lineardependence on s is expected to be a good approximation˜ K − ij ( s ) √ s = a ij + b ij s , s = E cm . (7.2)Such a parametrization is equivalent to a Breit-Wigner parametrization in the resonanceregion and is also similar to the well-known effective range expansion K − ij ( s ) = c ij + d ij p – 11 –ear threshold, where p is the momentum of the scattering particles in the center-of-momentum frame. We determine the parameters a ij and b ij following the strategy discussedabove, using the determinant residual method proposed in [46], which is briefly describedin Appendix B. A posteriori, we always verify that the resulting parametrization predictsvia Eq. (7.1) the same number of eigen-energies observed in the actual simulation in therelevant energy range; this is shown in the Ω plots for some fits in Appendix C. Thisprocedure will be followed for the extraction of the coupled-channel scattering matrix aswell as for one-channel scattering.The box-function B (cid:126)P , Λ ll (cid:48) ; i ( E cm ) can have off-diagonal elements for l (cid:48) (cid:54) = l due to the lack ofrotational symmetry in a finite box. This will result in contributions from multiple partialwaves in the quantization condition Eq. (7.1) for a given lattice irrep Λ. We considerpartial waves l = 0 and l = 2 and ignore contributions from l ≥
3, which is a reasonableassumption in the energy region considered for the respective irreps. In this case the onlynon-diagonal elements B ll (cid:48) among the A ( P (cid:54) = 0) irreps that are nonzero are B P =001 ,A and B P =110 ,A . These will be taken into account in the analysis of Section 8.4.2. In this section we present our results for the scattering matrices, pole positions, masses,and widths of J P C = 0 ++ and 2 ++ charmonium(like) states in various energy regions andwith varying assumptions. The energy range from slightly below 2 m D up to 4.13 GeV isdivided into smaller intervals, where the elements of the coupled D ¯ D − D s ¯ D s scatteringmatrix are separately parametrized according to Eq. (7.2) or as a constant. The detailsof the parametrizations and the results are presented in separate subsections below, whileinformation on the energy levels considered in each case is given in Appendix C.A single description of the whole energy region requires a finite-volume analysis involv-ing many more parameters, which results in more challenging and unstable fits. Such ananalysis is beyond the scope of the current investigation. Our inferences and conclusionsare based on the finite-volume analysis of separate energy regions. To explore the appli-cability of our results collectively to a wider energy region, we also model the scatteringmatrix based on the results from separate energy regions and investigate the implicationsin finite and infinite volume. D ¯ D scattering with l = 0 near threshold The narrow energy region near the D ¯ D threshold is significantly below the D s ¯ D s thresholdand can be treated in a one-channel approach. We employ the parametrization in Eq. (7.2) p cot δ l =0 D ¯ D E cm = ( ˜ K − ) l =011 E cm = a (cid:48) + b (cid:48) E cm , (8.1)which is equivalent to the effective range expansion p cot δ = 1 /a + r p D / E cm closest to 2 m D (listed in Appendix C.1) are utilized to– 12 – .82 3.84 3.86 3.88 3.90 3.92 - E cm [ GeV ] p c o t δ / E c m ( a ) E cm [ GeV ] ρ | t ( b ) - - - - - Re ( E cm ) - D [ GeV ] I m ( E c m ) [ G e V ] ( c ) Figure 3 : D ¯ D scattering in partial wave l = 0 near threshold. The green dashed linedenotes the D ¯ D threshold in the simulation where m D (cid:39) p cot δ/E cm , related to the scattering phase shift δ , as a function ofcenter-of-mass energy E cm . The red line indicates the parametrization of Eqs (8.1) and(8.2). The orange line represents ip/E cm . A bound state is located at the energy where thered and orange curves intersect. (b) The quantity ρ | t | that is proportional to the numberof D ¯ D events in experiment N D ¯ D ∝ pσ ∝ ρ | t | ( ρ = 2 p/E cm ). (c) Position of the pole ofthe scattering matrix on sheet I: the real component corresponds to the binding energy,presented in Eq. (8.5).determine the parameters via the quantization condition (7.1). We find a (cid:48) = − . ± . b (cid:48) = (0 . ± . a , cor = (cid:34) . − . . (cid:35) , χ d.o.f. = 1 . , (8.2)where cor is the correlation matrix defined in Appendix A. The fit is shown in Fig. 3a.This scattering matrix leads to a bound state at the energy E cm = m when the scatteringmatrix t (2.1) has a pole on the real axis below threshold on sheet I t = 1 ρ (cot δ − i ) = ∞ → cot δ = i , (8.3) p cot δ/E cm = ip/E cm = −| p | /E cm , ( p = i | p | ) . (8.4)The lhs of the second equation is shown as the red line in the figure, while the rhs isindicated by the orange line. The bound state occurs at the value E cm = m , where the twocurves intersect. The slope of p cot δ at the intersection, is smaller than the slope of −| p | ,– 13 –s required for an s-wave bound state (see Section VC of [25]). The location of the polein the scattering matrix is shown in Fig. 3c. The bound state appears just below the D ¯ D threshold with the binding energy χ D ¯ Dc : m − m D = − . +3 . − . MeV . (8.5)We denote this state by χ D ¯ Dc , indicating it has J P C = 0 ++ and a strong connection to the D ¯ D threshold. This state comes in addition to the conventional χ c (1 P ), which is foundsignificantly below threshold. Experiments cannot explore D ¯ D scattering below threshold,however, a closeby bound state below threshold could be identified experimentally througha sharp increase of the rate just above threshold. Fig. 3b shows a dimensionless quantity ρ | t | related to the number of events N D ¯ D ∝ pσ ∝ ρ | t | expected in experiment. It featuresa peak above threshold, which increases much more rapidly than the phase space.Such a D ¯ D bound state was not claimed by experiments so far. A similar statewas predicted in a phenomenological models [10, 47, 48], and some indication for it wassuggested in the experimental data [11, 49] and in data from the lattice simulation of Ref.[23]. A more detailed discussion follows in the summary in Section 9. Note that it is notknown whether this bound state would also feature in a simulation with physical quarkmasses. D s ¯ D s scattering with l = 0 near threshold in the one-channel approxima-tion The D s ¯ D s channel carries the same quantum numbers as D ¯ D necessitating the consider-ation of coupled-channel scattering. In this subsection we aim to get a rough estimate of D s ¯ D s scattering in the one-channel approximation, which will also provide initial guessesfor the parameters when coupled channel scattering is considered in Section 8.4. The D s ¯ D s scattering near threshold is parametrized by p cot δ l =0 D s ¯ D s E cm = ( ˜ K − ) l =022 E cm = a + b E cm . (8.6)We employ the quantization condition (7.1) together with four lattice energies close to thisthreshold that are dominated by D s ¯ D s interpolators (listed in Appendix C.2) and obtain a = − . ± . b = (0 . ± . a , cor = (cid:34) . − . . (cid:35) , χ d.o.f. = 2 . . (8.7)The resulting fit is shown in Fig. 4a. The scattering matrix has a bound state pole at theenergy E cm = m where condition (8.3) is satisfied, see Fig. 4c. Again, the slope of p cot δ issmaller than the slope of −| p | at the position of the pole, as required for an s-wave boundstate (see Section VC of [25]).We find a shallow D s ¯ D s bound state at χ D s ¯ D s c : m − m D s = − . +3 . − . MeV (8.8)– 14 – .90 3.92 3.94 3.96 3.98 4.00 4.02 - - - - E cm [ GeV ] p c o t δ / E c m ( a ) E cm [ GeV ] ρ | t ( b ) - - - - - Re ( E cm ) - D [ GeV ] I m ( E c m ) [ G e V ] ( c ) Figure 4 : D s ¯ D s scattering in partial wave l = 0 in the one-channel approximation. Thegreen dashed line denotes the D s ¯ D s threshold in the simulation, where m D s (cid:39) p cot δ/E cm , related to the scattering phase shift δ , as a function of center-of-mass energy E cm . The red line indicates the parametrizationof Eqs (8.6) and (8.7). The orange line represents ip/E cm . A bound state is located at theenergy where the red and orange curves intersect. (b) The quantity ρ | t | that is proportionalto the number of D s ¯ D s events in experiment N D s ¯ D s ∝ pσ ∝ ρ | t | ( ρ = 2 p/E cm ). (c)Position of the pole of the scattering matrix on sheet I: the real component corresponds tothe binding energy, presented in Eq. (8.8).that we denote χ D s ¯ D s c , indicating it has J P C = 0 ++ and a strong connection to the D s ¯ D s threshold. This state is responsible for the significant increase in the D s ¯ D s rate shown inFig. 4b just above threshold. In order to search for the χ D s ¯ D s c in experiment an explorationof the D s ¯ D s invariant mass near threshold would be invaluable. In one-channel D s ¯ D s scattering, considered here, the state is decoupled from D ¯ D , while it will become a narrowresonance and acquire a small width when the coupling to D ¯ D is considered in Section 8.4.Two candidates χ c (3930) [3] and X (3915) [2] (which may correspond to the samestate) have already been observed in experiment just below the threshold 2 m exp D s (cid:39) D ¯ D and a small width. If the D s ¯ D s bound state (8.8) corre-sponds to χ c (3930) and/or X (3915), it naturally explains both features as will be discussedin Section 9. – 15 – .90 3.95 4.00 4.05 4.10 4.15 4.20 - - - - E cm [ GeV ] p c o t δ / E c m ( a ) E cm [ GeV ] ρ | t ( b ) - - Re ( E cm ) I m ( E c m ) [ G e V ] ( c ) Figure 5 : D ¯ D scattering in partial wave l = 2 in the energy region around the conventional χ c (3930) resonance. (a) The purple crosses show the quantity p cot δ/E cm and the red linerepresents a simple Breit-Wigner resonance fit. (b) The quantity ρ | t | that is proportionalto the number of D ¯ D events in experiment N D s ¯ D s ∝ pσ ∝ ρ | t | shows a resonance peak.(c) The position of the χ c (3930) resonance pole E pcm = m − i Γ of the scattering matrixon sheet II from Eq. (8.11). D ¯ D scattering with l = 2 and J P C = 2 ++ resonance This channel features charmonia with J P C = 2 ++ . It is not the main focus of our study,however, an estimate of its scattering amplitude is required to extract the l = 0 scatteringamplitude using Eq. (7.1). We consider the energy region encompassing the 2 ++ resonanceand neglect the coupling to D s ¯ D s scattering with l = 2, which we assume to be negligiblein this region. The scattering amplitude is parametrized by the Breit-Wigner form (2.2) p cot δ l =2 D ¯ D E cm = ( ˜ K − ) l =211 E cm = a + b E cm = m J − E cm g J , (8.9)and the parameters are extracted via the quantization condition (7.1). We find m J = (1 . ± . a − = (4 . ± . ,g J = (10 . ± . a = (4 . ± .
3) GeV − cor = (cid:34) . . . (cid:35) , χ d.o.f. = 1 . , (8.10)using the four lattice energies closest to the resonance region (employing irreducible repre-sentation B with total momentum | (cid:126)P | = 1 as detailed in Appendix C.3). The correspond-ing fit is shown in Fig. 5a. The mass m J corresponds to the energy where the phase-shift– 16 –eaches π/
2, which is close-to the 2 ++ resonance mass obtained from the pole positionbelow, while the coupling g J is related to its width as shown in Eq. (2.2).The position of the pole E pcm of the scattering matrix (8.9) on sheet II provides a betterway of determining the resonance mass m and width Γ. We obtain χ c (3930) : E p = (4 . +0 . − . ) − i (0 . +0 . − . ) GeV = m − i Γ . (8.11)The pole is plotted in Fig. 5c. This leads to a J P C = 2 ++ resonance with χ c (3930) : m − M av = 905 +14 − MeV , Γ = 64 +32 − MeV , g = 4 . +0 . − . GeV − , (8.12)where g parametrizes the width Γ = g p D /m . This likely corresponds to the well-established resonance χ c (3930) = χ c (2 P ) [2]; a detailed comparison with experimentis made in Section 9. The resonance mass and the coupling obtained from the pole andfrom Eq. (8.10) are consistent, which is expected for a narrow resonance. D ¯ D , D s ¯ D s scattering with l = 0 for E cm (cid:39) . − . GeV
Finally, we turn to the coupled D ¯ D − D s ¯ D s scattering. We focus on the energy region E cm (cid:39) . − .
13 GeV near the D s ¯ D s threshold and we find an indication for severalinteresting hadrons. The scattering matrix for partial wave l = 0 is parametrized as( ˜ K − ) l =0 E cm = (cid:32) a + b E cm a a a + b E cm (cid:33) ≡ (cid:32) m J − E cm g J a a a + b E cm (cid:33) , (8.13)with the off-diagonal element held constant in E cm . Of the two equivalent parametrizationsshown above, we will utilize the one on the rhs. The 5 parameters in Eq. (8.13) aredetermined using all levels of irreps A (+)1 within the energy region E cm = 3 . − .
13 GeVdisplayed in Fig. 2: there are 14 levels from three frames with P = 0 , , N L = 24 ,
32 (see the black circles in the figure). The quantizationcondition (7.1) for A (+)1 irreps depends on the scattering amplitudes for l = 0, which weaim to determine. However, it also depends on the scattering amplitudes for l = 2 when P >
0. Below we present analyses both including and excluding the contribution from the l = 2 partial wave. l = 2In the first analysis we omit the contribution of the partial-wave l = 2. This is expected tobe a fair approximation since l = 2 effects D ¯ D scattering only in the narrow 2 ++ resonanceregion that is at the upper end of the current energy range of interest. The 5 parameters ofthe scattering matrix (8.13) are determined employing 13 out of 14 eigen-energies, where– 17 – Re( E cm ) [GeV] -0.08-0.06-0.04-0.02 0.00 I m ( E c m ) [ G e V ] m D s sheet IIsheet IIIsheet IV Re( E cm ) [GeV] -0.08-0.06-0.04-0.02 0.00 I m ( E c m ) [ G e V ] m D s sheet IIsheet IIIsheet IV -0.30 0.00 0.30 0.60 0.90 1.20 Re( c ) [GeV ] -0.15-0.10-0.05 0.00 0.05 0.10 I m ( c ) [ G e V ] sheet IIsheet IIIsheet IV -0.30 0.00 0.30 0.60 0.90 1.20 Re( c ) [GeV ] -0.15-0.10-0.05 0.00 0.05 0.10 I m ( c ) [ G e V ] sheet IIsheet IIIsheet IV -2.40 -1.60 -0.80 0.00 0.80 1.60 2.40 Re( c ) [GeV ] -0.50-0.30-0.10 0.10 0.30 0.50 I m ( c ) [ G e V ] sheet IIsheet IIIsheet IV -2.40 -1.60 -0.80 0.00 0.80 1.60 2.40 Re( c ) [GeV ] -0.50-0.30-0.10 0.10 0.30 0.50 I m ( c ) [ G e V ] sheet IIsheet IIIsheet IV Figure 6 : Coupled channel D ¯ D , D s ¯ D s scattering in partial wave l = 0: the analysisomitting l = 2 (left) and including l = 2 (right). Top: Location of the poles in t ij (8.13,8.14,8.15) in the complex energy plane. The orange poles are on Riemann sheet II,the magenta poles are on sheet III and the violet poles are on sheet IV (see Eq. (2.4)).The green dashed line denotes the D s ¯ D s threshold with m D s (cid:39) c i to channels i = 1 ( D ¯ D ) and i = 2 ( D s ¯ D s ) as extracted fromthe respective pole residues according to Eq. (2.5). The symbol fillings distinguish the twoorange poles. – 18 –he level dominated by ¯ cc [ J =2] operators is omitted m J = (1 . ± . a − = (3 . ± . g J = (0 . ± . a − = (1 . ± . a = − . ± . b = (0 . ± . a a = − . ± . . − .
034 0 . .
11 0 . . .
049 0 .
041 0 . . .
12 0 . . . . ,χ d.o.f = 0 . . (8.14) The lattice energy levels and the levels predicted by this parametrization are compared inAppendix C.4, where we verify that the same number of levels is observed and predicted.The resulting scattering matrix has several poles in the energy region investigated andtheir locations ( E pcm ) in the complex energy plane are shown in the left pane of Fig. 6. Thecouplings c i to both channels are extracted from the behavior of t ij near the poles usingEq. (2.5) and are given in the same figure. l = 0 and l = 2In the following we present our main result for the coupled channel scattering in the energyregion 3 . − .
13 GeV. We fix the scattering amplitude for D ¯ D scattering in partial wave l = 2 to the values in Eq. (8.10), while D s ¯ D s scattering in partial wave l = 2 is neglected.The parameters of the coupled D ¯ D , D s ¯ D s scattering matrix (8.13) with l = 0 are thendetermined using all 14 eigen-energies in the energy region of interest (from irreps A (+)1 ) m J = (1 . ± . a − = (3 . ± . g J = (0 . ± . a − = (1 . ± . a = − . ± . b = (0 . ± . a a = − . ± . . − .
078 0 .
13 0 .
079 0 . . − . − .
14 0 . . − . − . . − . . , χ dof = 0 . . (8.15) The coupling between channels a is non-zero but small. We also performed a study wherefive parameters for l = 0 and two parameters for l = 2 are fitted simultaneously using 18levels of irreps A (+)1 and B , and we obtain consistent results.The scattering matrix has several poles in the energy region investigated and theirlocations ( E pcm ) in the complex energy plane are shown in the right pane of Fig. 6. Theelements of scattering matrix t ij (2.5) near the poles are parametrized by the couplings c i presented in the same figure. These give information on the couplings of a resonance(associated with the pole) to both channels i = 1 ,
2. The location of poles and the cor-responding c i are similar to those extracted in the analysis omitting l = 2, shown in theleft pane of Fig. 6. The experimental rates are related to the values of √ ρ i ρ j | t ij | on thephysical axes, which are presented for D ¯ D → D ¯ D , D s ¯ D s → D s ¯ D s and D ¯ D → D s ¯ D s inFig. 7. Alternatively, the unitary 2 × S (2.3) on the physical axes canbe described by the phase shifts δ D ¯ D , δ D s ¯ D s and inelasticity η , which are shown in Fig. 8. This level corresponds to the 5th excited state in the A irrep with P = 1 on the smaller volume, N L = 24, see Fig. 2. – 19 – .90 3.95 4.00 4.05 4.10012345 E cm [ GeV ] ρ D | t DD - → DD - E cm [ GeV ] ρ D s | t D s D s - → D s D s - E cm [ GeV ] ( ρ D ρ D s ) / | t DD - → D s D s - Figure 7 : Coupled channel D ¯ D , D s ¯ D s scattering in partial wave l = 0. The three plotsshow the quantity √ ρ i ρ j | t ij | ∝ p | σ | on the physical axes, that is related to the numberof events in experiment ( ρ = 2 p/E cm ) for D ¯ D → D ¯ D (left), D s ¯ D s → D s ¯ D s (middle) and D ¯ D → D s ¯ D s (right). The green lines indicate the D s ¯ D s threshold with m D s (cid:39) - - - E cm [ GeV ] δ DD - E cm [ GeV ] δ D s D s - E cm [ GeV ] η Figure 8 : Coupled channel D ¯ D , D s ¯ D s scattering in partial wave l = 0: two phase shifts δ DD , δ D s ¯ D s and inelasticity η that parametrize the 2 × S (2.3) are shownas a function of E cm . The green lines indicates the D s ¯ D s threshold with m D s (cid:39) D ¯ D and D s ¯ D s channels are not strongly coupled, which can be seenfrom the inelasticity η (cid:39) ) and from the smallness of the off-diagonal element a in Eqs. (8.15) and (8.13). Our results suggest there are two 0 ++ resonances in thisenergy region: a narrow resonance dubbed χ D s ¯ D s c just below the D s ¯ D s threshold and abroader one denoted by χ (cid:48) c . The broader 0 ++ resonance χ (cid:48) c is related to the pole indicated in red on sheet III χ (cid:48) c : E p = (3 . +0 . − . ) − i (0 . +0 . − . ) GeV = m − i Γ . (8.16)This pole affects the scattering amplitude on the physical axes above the D s ¯ D s thresholdand is responsible for a peak around 3 .
98 GeV in the D ¯ D → D ¯ D rate shown in the leftpane of Fig. 7. The presence of this pole is also reflected in the phase shift δ D ¯ D , whichincreases gradually starting from 2 m D s as is evident in the left pane of Fig. 8. The nearbypole on sheet II does not have a significant influence on the physical scattering above thesecond threshold. The pole residues indicate that this state decays predominantly to D ¯ D ,while the decay to D s ¯ D s is suppressed, as evidenced by | c | (cid:29) | c | , see in the last two rows The inelasticity is equal to one below the D s ¯ D s threshold by construction due to the Θ-function inEq. (2.3). – 20 –f Fig. 6 for the pole presented in red. The resonance parameters are χ (cid:48) c : m − M av = 880 +28 − MeV , Γ = 58 +6 − MeV , g = 1 . +0 . − . GeV , (8.17)where M av = (3 m J/ψ + m η c ), and the coupling g parametrizes the full width Γ = g p D /m .The possible relation of this state to the broad resonance χ c (3860) discovered by Belle in2017 [1, 2] is discussed in Section 9. The narrow 0 ++ resonance χ D s ¯ D s c near the D s ¯ D s threshold is related to the poleon sheet II, indicated by the top-filled orange symbols in the first row of Fig. 6. Its locationrelative to the threshold is given by χ D s ¯ D s c : E pcm − m D s = ( − . +0 . − . ) − i (0 . +2 . − . ) MeV = m − i Γ . (8.18)This resonance is related to the bound state in the analysis of D s ¯ D s -scattering in the one-channel approximation of Section 8.2. The pole on sheet II and the nearby pole on sheetIV correspond to this resonance and are mutually exclusive across the bootstrap samples.Further details on this can be found in Appendix C.4. It is clear from Figure 7 that theresonance pole leads to a sharp rise in the D s ¯ D s → D s ¯ D s and D ¯ D → D s ¯ D s rates just above2 m D s . The increased D ¯ D → D s ¯ D s rate is also responsible for a dip in the D ¯ D → D ¯ D rateat 2 m D s and all three features should be used as a signature for experimental searches ofthis state. Note that the magnitude of the D s ¯ D s → D s ¯ D s peak above 2 m D s is larger whenthe pole is closer to the threshold. χ D s ¯ D s c couples predominantly to D s ¯ D s and very weakly to D ¯ D (one can see that | c | (cid:29) | c | in Fig. 6). The mass difference of the state with respect to the threshold andits narrow total width Γ = g p D /m parametrized in terms of g are χ D s ¯ D s c : m − m D s = − . +0 . − . MeV , Γ = 0 . +2 . − . MeV , g = 0 . +0 . − . GeV . (8.19)On the experimental side, the newly discovered χ c (3930) [3] and the X (3915) [2] lie nearthe D s ¯ D s threshold and have very small or zero decay rate to D ¯ D . The indication for a D s ¯ D s state in our study explains both properties, as detailed in Section 9.The parameters of the scattering matrix obtained from the analysis including or exclud-ing l = 2 are similar, with the exception of a and b . These parametrize D s ¯ D s → D s ¯ D s and differ between the coupled-channel analysis and the one-channel approximation, how-ever, both analyses lead to a state just below the D s ¯ D s threshold on the real axis (seeFig. 4c) or slightly away from it. The conclusion that there is a near-threshold pole is ro-bust, while its exact location and the effect on physical scattering need to be investigatedin a simulation with higher statistics and a better control of systematic uncertainties. D ¯ D , D s ¯ D s scattering in a wider energy region In the previous subsections, we discussed the analysis of rather narrow energy ranges witha linear and/or constant parametrization in s for the elements ˜ K − ij ( s ) / √ s (see Eq. (8.13)).A single description of coupled D ¯ D − D s ¯ D s scattering encompassing the whole of theenergy range from 2 m D up to 4.13 GeV requires additional parameters. This is difficult– 21 – .68 1.70 1.72 1.74 1.76 - - - E cm a a K - ij / E c m Figure 9 : The modeled ˜ K − / √ s -matrix elements in lattice units for the coupled D ¯ D − D s ¯ D s scattering: red for D ¯ D → D ¯ D , black for D s ¯ D s → D s ¯ D s and blue for D s ¯ D s → D ¯ D .The green dashed lines indicate the D ¯ D and the D s ¯ D s thresholds, whereas the orangedashed lines give the bound state constraint for either channel.as, with the statistics and the number of lattice QCD ensembles available to us, the fitsbecome unstable. Below we employ an alternative approach, in which we first model theinfinite-volume scattering matrix in the wider energy range based on the parametrizationspresented in the previous subsections. We then investigate the implications of the resultingscattering matrix in the complex-energy planes and compare with the analyses above. Wealso verify that it predicts the same number of finite-volume energy levels as observed inthe actual simulation.The t -matrix elements are modeled in the energy range E cm (cid:39) m D − .
13 GeV ina piecewise fashion as shown for ˜ K − / √ s in Fig. 9 based on Eqs. (7.2) and (8.13). Werequire that these are continuous in energy, however, their derivatives with respect toenergy need not be continuous. The parameters in the high-energy region are fixed to thevalues (8.14) obtained from the coupled-channel analysis. Below we provide more detailson each t -matrix element in turn: t The energy region considered is divided into three intervals, as shown by the red linein the figure. t is parametrized as in the coupled channel analysis of Eq. (8.14) andthe single D ¯ D channel analysis of Eq. (8.2) in the high and middle energy intervals,respectively. In the lower region, t is chosen such that ˜ K − / √ s is constant in orderto prevent the occurrence of a second bound-state (this would correspond to the redand orange lines in the figure intersecting a second time). The exact form of thischoice is not important as this is beyond the region of interest. t Again the energy region is split into three intervals, see the blue line in the figure.In the upper region, starting somewhat below the D s ¯ D s threshold, t is set to theconstant value of the coupled channel analysis presented in Eq. (8.14). t is set tozero in the lower region, where the effects from D s ¯ D s channel are expected to benegligible. In between these two regions, linear behaviour ensures continuity.– 22 – Re( E cm ) [GeV] -0.06-0.05-0.04-0.03-0.02-0.01 0.00 0.01 I m ( E c m ) [ G e V ] m D s m D sheet Isheet IIsheet IIIsheet II J PC =0 ++ J PC =2 ++ J PC =0 ++ J PC =2 ++ Figure 10 : Pole singularities of the scattering amplitude/matrix in the complex energyplane, which are associated with the hadrons predicted in this work. The pole related tothe J P C = 2 ++ resonance appears on sheet II, as we have neglected the l = 2 partial wavecontribution from D s ¯ D s scattering. t This element is parametrized as in Eq. (8.14), for the entire energy region, see theblack line in the figure. Note that we ignore any crossing of the D s ¯ D s bound statecondition with the t parametrization that occurs well below the D ¯ D threshold.We observe that the number of poles with the above-designed scattering matrix is thesame as that obtained from the analysis of the separate energy regions. The pole locationson the various complex Riemann sheets and their residues are also largely unchanged.Following L¨uscher’s finite-volume analysis, we have extracted the finite-volume spectrumfrom this scattering matrix and have verified that the number of levels in the energy rangeof interest from this approach agrees with that of the actual lattice simulation. This isdemonstrated in Appendix D. Finally, we remark that our conclusions are based on theresults presented in the previous subsections, while the approach followed here provides aconsistency check of our analysis over the full energy range. Below we summarize the properties of the charmonium-like hadrons found in this simu-lation. These are denoted by the conventional names χ c (1 P ), χ c (1 P ), χ c (3930) whenthe identification with experimental states is unambiguous, while the other states foundare denoted χ (cid:48) c , χ D ¯ Dc and χ D s ¯ D s c . The subscripts c c J P C = 0 ++ and 2 ++ , respectively. The location of the poles in the complex energy planerelated to these hadrons are given in Fig. 10, while the corresponding masses are comparedto experiment in Fig. 11. Due to the unphysical quark masses employed in the simulation,the hadron masses obtained are not compared to experiment directly. Instead, the differ-ence m − E ref is utilized, where the reference energy E ref is either a nearby threshold or– 23 – ++ ++ ++ ++ ++ ++ m [ G e V ] J PC D Ds Ds D Lat : m Lat : m - E ref + E ref exp c c2 (1P) Exp c c0 (1P) c c2 (3930)X(3860) c c0 (3930)X(3915) Figure 11 : Energy spectrum of hadrons predicted in this work compared with experiment.The left and middle panes show the lattice results, where m denotes the mass obtainedfrom the lattice study. The magenta symbols correspond to hadrons extracted via thescattering analysis: diamonds represent resonances and triangles represent bound states.The blue crosses are extracted directly from the lattice energies. The reference energy inthe middle pane is E ref = 2 m D (2 m D s ) for the state closest to the D ¯ D ( D s ¯ D s ) threshold,while E ref = M av = (3 m J/ψ + m η c ) for the remaining four states. The right pane showsthe experimental spectrum [2], where χ c (3930) [3] and X (3915) [2] may be the same state.the spin-averaged charmonium mass M av = (3 m J/ψ + m η c ). The positions of the D ¯ D and D s ¯ D s thresholds are given by the masses m D (cid:39) m D s (cid:39) p l +1 evaluated for the mesonmomenta (in the cm-frame) at the resonance energy, which in turn depends on the positionof the threshold. The latter is different in the simulation and in experiment. Therefore itis customary to compare the coupling g that parametrizes the full width of a hadronΓ ≡ g p l +1 D /m with l = 0 , J P C = 0 ++ , ++ , (9.1)as g is expected to be less dependent on the quark masses than the width itself. χ c (1 P ) and χ c (1 P ) These states lie significantly below the D ¯ D threshold and their masses are extracted fromthe ground state energies m = E ( P = 0) in irreps A ++1 and E ++ , with J P C = 0 ++ and 2 ++ ,– 24 –espectively. We obtain χ c (1 P ) : m − M av = 358 ± , m = 3461 ± , (9.2) χ c (1 P ) : m − M av = 462 ± , m = 3565 ± , (9.3)which can be compared to the experimental valuesexp χ c (1 P ) : m − M av = 346 . ± . , (9.4)exp χ c (1 P ) : m − M av = 487 . ± .
14 MeV . (9.5) ++ state χ D ¯ Dc slightly below D ¯ D threshold We find a shallow D ¯ D bound state labeled χ D ¯ Dc with binding energy m − m D = − . +3 . − . MeV . (9.6)Note that it is not known whether this bound state would also feature in a simulationwith physical quark masses. Such a state has not been claimed by experiments, however,the region above the D ¯ D threshold is often populated by D ¯ D from X (3872) → D ¯ D ∗ → D ¯ Dπ (see, for example, Ref. [50]) which makes the search for a state in this region morechallenging.The existence of a shallow D ¯ D bound state dubbed X (3720) was already suggestedby an effective phenomenological model in Ref. [10] . Ref. [11] indicates that there mayalready be some evidence for such a state in the D ¯ D invariant mass distribution from Belle[49], which shows an enhancement just above threshold. The D ¯ D rate from Babar [51]also shows a hint of enhancement just above threshold (see Fig. 5 of [51]). In a molecularpicture, a 0 ++ state is expected as a partner of X (3872) via heavy-quark symmetry argu-ments [47, 48]. A similar virtual bound state with a binding energy of 20 MeV follows fromthe data of the only previous lattice simulation of D ¯ D scattering [23] . A strategy for anexperimental search for this state in B → D ¯ D K decay is proposed in Ref. [52]. ++ resonance and its relation to χ c (3930) We find a resonance with J P C = 2 ++ in l = 2 D ¯ D scattering with the following properties χ c (3930) : m − M av = 904 +14 − MeV , g = 4 . +0 . − . GeV − . (9.7)This is most likely related to the conventional χ c (3930) resonance (also called χ c (2 P ))discovered by Belle [53]exp χ c (3930) : m − M av = 854 ± , g = 2 . ± .
12 GeV − . (9.8)Here g parametrizes the width Γ = g p D /m . The masses are reasonably close, while thecoupling from lattice QCD is larger that in experiment. However, the couplings are alsonot inconsistent given the large statistical error from our simulation and the unquantifiedsystematic uncertainties discussed in Section 5. This state with m (cid:39) .
718 GeV is listed in Table 4 of Ref. [10]. The presence of this state was not mentioned in Ref. [23], as such virtual bound states were not searchedfor. – 25 – road 0 ++ resonance and its possible relation to χ c (3860) This resonance couples mostly to D ¯ D and has a very small coupling to D s ¯ D s . Its resonanceparameters are χ (cid:48) c : m − M av = 880 +28 − MeV , g = 1 . +0 . − . GeV . (9.9)These can be compared to the scalar resonance χ c (3860) discovered by Belle in 2017 [1] ,exp χ c (3860) : m − M av = 793 +48 − MeV , g = 2 . +1 . − . GeV , (9.10)based on the following arguments: The mass and coupling are reasonably consistent withexperiment, in particular, when considering the experimental errors and the systematicuncertainties in the lattice results. The mass is also close to the value obtained from theonly previous lattice study of D ¯ D scattering [23] , although the width and coupling arelarger in the present work.The experimental couplings g of χ c (3860) and χ c (3930), quoted in Eqs. (9.10) and(9.8), respectively, are similar, and their widths differ mainly because of the different phasespace. On the lattice side, the χ (cid:48) c coupling in Eq. (9.9) is smaller than the χ c (3930)coupling in Eq. (9.7), however, both have large uncertainties. Narrow 0 ++ resonance χ D s ¯ D s c and its possible relation to χ c (3930), X (3915) We find a narrow 0 ++ resonance near the D s ¯ D s threshold. It has a large coupling to D s ¯ D s and a very small coupling to D ¯ D . The latter is responsible for its small decay rate to D ¯ D and the small total width. This state corresponds to the bound state in one-channel D s ¯ D s scattering discussed in Section 8.2. We compare the resulting resonance parameters χ D s ¯ D s c : m − m D s = − . +0 . − . MeV , g = 0 . +0 . − . GeV (9.11)with two experimental states that share similar features to the state we find,exp χ c (3930) : m − m D s = − . ± . , Γ = 17 ± , g = 0 . ± .
10 GeV , exp X (3915) : m − m D s = − . ± . , Γ = 20 ± , g = 0 . ± .
10 GeV . (9.12)The χ c (3930) with J P C = 0 ++ was very recently discovered in D ¯ D decay by LHCb[3]. The X (3915) was observed by Belle [6] and BaBar [5, 7, 8] in J/ψω decay and has J P C = 0 ++ or 2 ++ , while its decay to D ¯ D was not observed [2]. They might representthe same state if X (3915) is a scalar. Both experimental states lie just below the D s ¯ D s threshold. One would naturally expect 0 ++ states with this mass to be broad, given thelarge phase space available to D ¯ D decay. Their narrow widths can only be explained iftheir decay to D ¯ D is suppressed by some mechanism. For the couplings calculated from the experimental values we vary both the mass and width by ± σ and take the maximal positive and negative deviations as the uncertainties The value to compare is m − M expav (cid:39) .
90 GeV and 0 .
93 GeV from the fits in Eqs. (6.3) and (6.7) ofRef. [23], respectively. – 26 –f the resonance found on the lattice is indeed related to X (3915)/ χ c (3930), our re-sults indicate that this state owes its existence to a large interaction in the D s ¯ D s channelnear threshold, which naturally explains why its width is small and its decay to D ¯ D issuppressed. Note that a detailed quantitative comparison of lattice and experimental re-sults in Eqs. (9.11) and (9.12) is not possible due to the unphysical masses of the quarksin the lattice study and due to the omission of decays to J/ψω and η c η , which may affectthe determination of the width. The qualitative comparison, however, suggests the exis-tence of a D s ¯ D s resonance with small coupling to D ¯ D . This could be further investigatedexperimentally by considering the D s ¯ D s invariant mass spectrum near threshold, where apeak (see Fig. 7) would be visible for a state just below threshold.The X (3915) was proposed to be a ground ¯ cc ¯ ss state within the diquark-antidiquarkapproach by Polosa and Lebed [9]. The identification ¯ cc ¯ ss was considered also in phe-nomenological studies [54, 55].
10 Conclusions
We presented a lattice study of coupled-channel D ¯ D - D s ¯ D s scattering in the J P C = 0 ++ and 2 ++ quantum channels with isospin 0. Using the generalized L¨uscher method anda piecewise parametrization of the energy region from slightly below 2 m D to 4.13 GeV,the coupled-channel scattering matrix S along the real energy axis was determined. Theresulting S was then analytically continued to search for pole singularities in the complexenergy plane that can affect the scattering amplitudes/parameters along the physical axes.Our study utilized the spectrum in three different inertial frames determined on two CLSensembles with u/d and s quarks, spatial extents ∼ ∼ ∼ χ c (1 P ), the results suggest three charmonium-like states with J P C =0 ++ below 4.13 GeV. One is a yet undiscovered D ¯ D bound state just below threshold.The second is a narrow resonance just below the D s ¯ D s threshold predominantly coupledto D s ¯ D s . This state is possibly related to the narrow resonance X (3915) /χ c (3930), whichis also below the D s ¯ D s threshold in the experiments. The third feature is a D ¯ D resonancepossibly related to the χ c (3860) observed by Belle, which is believed to be χ c (2 P ). Anoverview of the resulting pole structure of the coupled-channel D ¯ D - D s ¯ D s scattering matrixin the complex energy plane is given in Fig. 10, and the possible implications of thissingularity structure for experiments are illustrated in Figs. 7 and 8. The masses arecompared to experiment in Fig. 11 and summarized in Section 9.Turning to states with J P C = 2 ++ , the mass of the ground state χ c (1 P ) was de-termined directly from the lattice energy and is compared with the experimental value inEq. (9.3). We have assumed the 2 ++ resonance to be coupled only with the D ¯ D scatteringchannel in the l = 2 partial wave and have parametrized this with a Breit-Wigner form.The resonance parameters are extracted and compared with the experimental values of theconventional χ c (3930) in Eqs. (9.7) and (9.8). These are then fixed for the finite-volumecoupled-channel analysis discussed in Section 8.4.2. We find the estimates for positionsand residues for the poles with J P C = 0 ++ to be robust with the exclusion/inclusion of– 27 –he l = 2 partial wave contribution to the analysis. The resulting pole positions and theresidues from either study are shown in Fig. 6.In this study we worked with several simplifying assumptions (detailed in Section 5)necessary for a first investigation of this coupled-channel system. The lattice QCD ensem-bles we used have heavier-than-physical light and charm quarks and a lighter-than-physicalstrange quark. This results in a smaller-than-physical splitting between the D ¯ D and the D s ¯ D s thresholds. In future studies, it will be important to systematically improve upon thecurrent results by successively relaxing our assumptions, for example, by explicitly includ-ing η c η interpolating fields and adding this channel as well as J/ψω to the coupled-channelstudy. With regard to the pole structure observed in this work, it would be particularlyinteresting to investigate how our observations, such as the shallow D ¯ D bound state andour X (3860) candidate evolve when simultaneously approaching the limit of physical quarkmasses and the continuum limit. Acknowledgments
We thank G. Bali, V. Baru, T. Gershon, F.-K. Guo, D. Johnson, C. B. Lang, J. Nieves, E.Oset and A. Sch¨afer for useful discussions. We thank the authors of Ref. [46] for makingthe TwoHadronsInBox package public and Ben H¨orz and C. B. Lang for contributions tothe computing codes we used. We use the multigrid solver of Refs. [56–59] for the inversionof the Dirac operator. Our code implementing distillation is written within the frameworkof the Chroma software package [60]. The simulations were performed on the RegensburgiDataCool and Athene2 clusters, and the SFB/TRR 55 QPACE 2 [61] and QPACE 3 ma-chines. We thank our colleagues in CLS for the joint effort in the generation of the gaugefield ensembles which form a basis for the computation. The Regensburg group was sup-ported by the Deutsche Forschungsgemeinschaft (collaborative research centre SFB/TRR-55), the European Union’s Horizon 2020 Research and Innovation programme under theMarie Sklodowska-Curie grant agreement no. 813942 (ITN EuroPLEx), and the STRONG-2020 project under grant agreement no. 824093. M. P. acknowledges support from the EUunder grant no. MSCA-IF-EF-ST-744659 (XQCDBaryons). S. Prelovsek was supported bySlovenian Research Agency ARRS (research core funding No. P1-0035 and No. J1-8137).We are grateful to the Mainz Institute for Theoretical Physics (MITP) for its hospitalityand its partial support during the course of this work.
A Error treatment
All quoted estimates are obtained from the ensemble average of hadron correlation mea-surements and the errors are based on N b = 999 bootstrap samples. The 1 σ standard errorformulae for a Gaussian-distributed quantity Q provides the range that capture central68% of the samples. We present resonance masses, widths/couplings, pole positions, phaseshifts and | t | with the asymmetric errors ¯ Q + σ + − σ − , where the interval [ ¯ Q − σ − , ¯ Q + σ + ] con-tains the central 68% of the bootstrap samples. The remaining quantities such as energies,parameters of the scattering matrices are provided with symmetric errors σ = √ cov ii . Here– 28 –ov is the modified correlation matrix defined ascor ij = cov ij √ cov ii cov jj , cov ij = MN b (cid:88) b (cid:48) ( Q bi − ¯ Q i )( Q bj − ¯ Q j ) , M = . . (A.1)The sum indicated by a prime runs over the bootstrap samples b in which Q bi and Q bj arenot among the 16% of the values excluded on either end. cov ii in Eq. (A.1) is equal tothe standard covariance N b (cid:80) N b b =1 ( Q bi − ¯ Q i ) for Gaussian-distributed quantities. cov ij alsocoincides with the standard covariance for completely correlated or uncorrelated Gaussian-distributed quantities i and j . We have verified that the standard covariance and covrender almost identical errors on the energies and most of the parameters of the scatteringmatrix. The advantage of the modified correlation matrix is the exclusion of outliers fornon-Gaussian distributions with long tails. Such distributions might occur for the bootstrapsamples of scattering-matrix parametrizations due to the highly non-linear nature of thebox-functions B ( E cm ) in Eq. (7.1). B Fitting the parameters of the scattering matrix
The parameters of the matrix ˜ K are determined from the energies presented in Fig. 2via the quantization condition (7.1) following the determinant residual method proposedin Ref. [46]. In this method, one determines the parameters such that the zeros of theΩ( E cm ) function (which are identical to the zeros of the determinant in Eq. (7.1))Ω( E cm ) = det( A )det(( µ + AA † ) / ) , A ( E cm ) = ˜ K − ( E cm ) − B ( E cm ) , (B.1)are as close as possible to E cm extracted from the lattice. Our results are based on fits with µ = 1 and the parameters change negligibly for values of µ in the interval [0 . , E cm ) as a function of E cm for the resulting parameters of the scattering matrix aregiven in Figs. 12 and 13. The values of E cm where Ω( E cm ) crosses zero are indeed nearthe positions of the observed lattice energies. We have verified that the number of zero-crossings equals the number of observed eigenstates in the energy-regions relevant for allour fits. C More details on the analysis of the scattering channels
Below we specify the energy levels used for each analysis, referring to levels E cm,n plottedin Fig. 2 counting from the lowest state with n = 1. Their values, their errors and covari-ance matrices can be obtained from the authors on request. The masses of the scatteringparticles are m D a = 0 . m D s a = 0 . L/a = 32 ensemble and m D a = 0 . m D s a = 0 . L/a = 24 ensemble. For simplicity, weuse the phase space factor from the full ensemble average during the pole search for allbootstrap samples. The central values are not affected by this procedure.– 29 – .1 D ¯ D scattering with l = 0 near threshold This analysis employs four energy levels closest to the 2 m D threshold in the A (+)1 irrepsshown in Fig. 2: these are levels n = 2(3) from | (cid:126)P | = 0(1) on both volumes. An analysisalso including levels n = 3 , P = 2 on the N L = 32 ensemble(indicated by the gray points in Fig. 3 below E cm = 3 . p cot δ/E cm just fails to cross the orange line in Fig. 3a. An additional 6% ofthe bootstrap samples render a virtual bound state - this corresponds to the bootstrapsfor which p cot δ/E cm crosses ip/E cm = | p | /E cm rather than −| p | /E cm in Fig. 3a. Both ofthese scenarios happen in extreme cases and end up within the 32% of bootstrap samplesthat are excluded when computing the errors via (A.1). C.2 D s ¯ D s scattering with l = 0 near threshold in the one-channel approxima-tion This analysis employs only those eigenstates whose overlaps are dominated by D s ¯ D s op-erators and that do not have significant overlap with D ¯ D operators. These are the fourlevels in Fig. 2 near the D s ¯ D s threshold in the A (+)1 irreps: levels n = 3 , | (cid:126)P | = 0 , N L = 24 ensemble and levels n = 4 , | (cid:126)P | = 1 , N L = 32. Here 97% of thebootstrap samples result in a bound-state pole, while 2.3% result in a virtual bound-statepole and 0.7% do not render any poles on the real axis – the latter two cases end up amongthe extremal 32% of bootstrap samples. C.3 D ¯ D scattering with l = 2 D ¯ D scattering in partial wave l = 2 is the not the main focus of our study. It was consideredin order to investigate and constrain its contribution to the A irreps we have studied. Theextraction of the phase shift in Eq. (8.9) employs four lattice levels in the B irrep with | (cid:126)P | = 1: these are levels n = 3 , n = 1 , J P C = 2 ++ and 2 − + , respectively). C.4 Coupled D ¯ D , D s ¯ D s scattering with l = 0 for E cm (cid:39) . − . GeV
Fig. 12 shows an example of Ω( E cm ) (B.1) for the parameters of the coupled-channelscattering matrix given in Eq. (8.14). The values of E cm at which Ω crosses zero are indeednear the observed eigen-energies (indicated by the black circles). The number of crossingsagrees with the number of observed levels in the relevant energy ranges.In relation to the poles near 2 m D s in Fig. 6, 70% of the samples have a pole on sheetII, while a pole appears at a similar location on sheet IV for the rest. Hence the results forthe pole location related to χ D s ¯ D s c given in Eq. (8.18) are based on the samples for whichthis pole appears on sheet II, whereas the errors on the rates and phase shifts presented inFigs. 7 and 8 are computed from the entire bootstrap samples.– 30 – Ω ( E cm ) E cm [ GeV ] P =
0, L / a = - Ω ( E cm ) E cm [ GeV ] P =
1, L / a = - Ω ( E cm ) E cm [ GeV ] P =
2, L / a = - Ω ( E cm ) E cm [ GeV ] P =
0, L / a = - Ω ( E cm ) E cm [ GeV ] P =
1, L / a = - Ω ( E cm ) E cm [ GeV ] P =
2, L / a = Figure 12 : The function Ω( E cm ) (defined in Eq. (B.1)) for the resulting scattering matrixof the coupled channels D ¯ D − D s ¯ D s (8.14) is given by the orange line. The observedeigen-energies are given by the circles: the black levels are employed to fit the parameters(8.14), while the blue circles are not. D Ω( E cm ) in a wider energy region E cm (cid:39) m D − . GeV
In Fig. 13, we present the Ω( E cm ) function (defined in Eq. (B.1)) for the coupled D ¯ D − D s ¯ D s scattering matrix as modeled in Section 8.5. The points indicate the energy levelsobserved in the actual simulation. The zeros of the Ω( E cm ) are the predictions for thefinite-volume spectrum derived from the scattering matrix. It is clear from the figure thatthe predicted spectrum agrees qualitatively with the observed spectrum. References [1] Belle, K. Chilikin et al. , Phys. Rev. D , 112003 (2017), [arXiv:1704.01872]. – 31 – Ω ( E cm ) E cm [ GeV ] P =
0, L / a = - Ω ( E cm ) E cm [ GeV ] P =
1, L / a = - Ω ( E cm ) E cm [ GeV ] P =
2, L / a = - Ω ( E cm ) E cm [ GeV ] P =
0, L / a = - Ω ( E cm ) E cm [ GeV ] P =
1, L / a = - Ω ( E cm ) E cm [ GeV ] P =
2, L / a = Figure 13 : The Ω( E cm ) function (defined in Eq. (B.1)) for the scattering matrix of thecoupled channels D ¯ D − D s ¯ D s in the wider energy region ( E cm (cid:39) m D − . [2] Particle Data Group, P. Zyla et al. , PTEP , 083C01 (2020).[3] LHCb, R. Aaij et al. , 2009.00026.[4] LHCb, R. Aaij et al. , 2009.00025.[5] BaBar, B. Aubert et al. , Phys. Rev. Lett. , 082001 (2008), [arXiv:0711.2047].[6] Belle, S. Uehara et al. , Phys. Rev. Lett. , 092001 (2010), [arXiv:0912.4451].[7] BaBar, P. del Amo Sanchez et al. , Phys. Rev. D , 011101 (2010), [arXiv:1005.5190].[8] BaBar, J. Lees et al. , Phys. Rev. D , 072002 (2012), [arXiv:1207.2651].[9] R. F. Lebed and A. D. Polosa, Phys. Rev. D , 094024 (2016), [arXiv:1602.08421]. – 32 –
10] D. Gamermann, E. Oset, D. Strottman and M. Vicente Vacas, Phys. Rev. D , 074016(2007), [arXiv:hep-ph/0612179].[11] D. Gamermann and E. Oset, Eur. Phys. J. A , 189 (2008), [arXiv:0712.1758].[12] M. L¨uscher, Commun. Math. Phys. , 153 (1986).[13] M. L¨uscher, Nucl. Phys. B354 , 531 (1991).[14] M. L¨uscher, Nucl. Phys.
B364 , 237 (1991).[15] Hadron Spectrum, J. J. Dudek, R. G. Edwards, C. E. Thomas and D. J. Wilson, Phys. Rev.Lett. , 182001 (2014), [arXiv:1406.4158].[16] D. J. Wilson, J. J. Dudek, R. G. Edwards and C. E. Thomas, Phys. Rev. D , 054008(2015), [arXiv:1411.2004].[17] Hadron Spectrum, J. J. Dudek, R. G. Edwards and D. J. Wilson, Phys. Rev. D , 094506(2016), [arXiv:1602.05122].[18] R. A. Briceno, J. J. Dudek, R. G. Edwards and D. J. Wilson, Phys. Rev. D , 054513(2018), [arXiv:1708.06667].[19] G. Moir, M. Peardon, S. M. Ryan, C. E. Thomas and D. J. Wilson, JHEP , 011 (2016),[arXiv:1607.07093].[20] CLQCD, T. Chen et al. , Chin. Phys. C , 103103 (2019), [arXiv:1907.03371].[21] HAL QCD, Y. Ikeda et al. , Phys. Rev. Lett. , 242001 (2016), [arXiv:1602.03465].[22] HAL QCD, Y. Ikeda, J. Phys. G , 024002 (2018), [arXiv:1706.07300].[23] C. Lang, L. Leskovec, D. Mohler and S. Prelovsek, JHEP , 089 (2015), [arXiv:1503.05363].[24] M. Padmanath et al. , Phys. Rev. D , 014513 (2019), [arXiv:1811.04116].[25] S. Piemonte, S. Collins, D. Mohler, M. Padmanath and S. Prelovsek, Phys. Rev. D ,074505 (2019), [arXiv:1905.03506].[26] R. Brett et al. , Nucl. Phys. B , 29 (2018), [arXiv:1802.03100].[27] M. Bruno et al. , JHEP , 043 (2015), [arXiv:1411.3982].[28] RQCD, G. S. Bali, E. E. Scholz, J. Simeth and W. S¨oldner, Phys. Rev. D94 , 074501 (2016),[arXiv:1606.09039].[29] M. Bruno, T. Korzec and S. Schaefer, Phys. Rev.
D95 , 074504 (2017), [arXiv:1608.08900].[30] M. L¨uscher and S. Schaefer, Comput. Phys. Commun. , 519 (2013), [arXiv:1206.2809].[31] D. Mohler and S. Schaefer, Phys. Rev. D , 074506 (2020), [arXiv:2003.13359].[32] Hadron Spectrum, M. Peardon et al. , Phys. Rev.
D80 , 054506 (2009), [arXiv:0905.2160].[33] C. Morningstar et al. , Phys. Rev.
D83 , 114505 (2011), [arXiv:1104.3870].[34] R. A. Brice˜no, Phys. Rev.
D89 , 074507 (2014), [arXiv:1401.3312].[35] J. J. Dudek, R. G. Edwards, M. J. Peardon, D. G. Richards and C. E. Thomas, Phys. Rev.D , 034508 (2010), [arXiv:1004.4930].[36] C. E. Thomas, R. G. Edwards and J. J. Dudek, Phys. Rev. D , 014507 (2012),[arXiv:1107.1930]. – 33 –
37] J. J. Dudek, R. G. Edwards and C. E. Thomas, Phys. Rev. D , 034031 (2012),[arXiv:1203.6041].[38] R. A. Briceno, J. J. Dudek and R. D. Young, Rev. Mod. Phys. , 025001 (2018),[arXiv:1706.06223].[39] A. X. El-Khadra, A. S. Kronfeld and P. B. Mackenzie, Phys. Rev. D , 3933 (1997),[arXiv:hep-lat/9604004].[40] M. B. Oktay and A. S. Kronfeld, Phys. Rev. D , 014504 (2008), [arXiv:0803.0523].[41] C. Michael, Nucl. Phys. B , 58 (1985).[42] M. L¨uscher and U. Wolff, Nucl. Phys. B , 222 (1990).[43] B. Blossier, M. Della Morte, G. von Hippel, T. Mendes and R. Sommer, JHEP , 094(2009), [arXiv:0902.1265].[44] M. D¨oring, U.-G. Meissner, E. Oset and A. Rusetsky, Eur. Phys. J. A , 139 (2011),[arXiv:1107.3988].[45] M. T. Hansen and S. R. Sharpe, Phys. Rev. D , 016007 (2012), [arXiv:1204.0826].[46] C. Morningstar et al. , Nucl. Phys. B , 477 (2017), [arXiv:1707.05817].[47] C. Hidalgo-Duque, J. Nieves, A. Ozpineci and V. Zamiralov, Phys. Lett. B , 432 (2013),[arXiv:1305.4487].[48] V. Baru et al. , Phys. Lett. B , 20 (2016), [arXiv:1605.09649].[49] Belle, P. Pakhlov et al. , Phys. Rev. Lett. , 202001 (2008), [arXiv:0708.3812].[50] LHCb, R. Aaij et al. , JHEP , 035 (2019), [arXiv:1903.12240].[51] BaBar, B. Aubert et al. , Phys. Rev. D , 092003 (2010), [arXiv:1002.0281].[52] L. Dai, J.-J. Xie and E. Oset, Eur. Phys. J. C , 121 (2016), [arXiv:1512.04048].[53] Belle, S. Uehara et al. , Phys. Rev. Lett. , 082003 (2006), [arXiv:hep-ex/0512035].[54] J. F. Giron and R. F. Lebed, Phys. Rev. D , 014036 (2020), [arXiv:2005.07100].[55] W. Chen, H.-X. Chen, X. Liu, T. Steele and S.-L. Zhu, Phys. Rev. D , 114017 (2017),[arXiv:1706.09731].[56] S. Heybrock et al. , Lattice QCD with Domain Decomposition on Intel Xeon PhiCo-Processors, in The International Conference for High Performance Computing,Networking, Storage, and Analysis: SC14: HPC matters (SC) New Orleans, LA, USA,November 16-21, 2014 , 2014, [1412.2629].[57] S. Heybrock, M. Rottmann, P. Georg and T. Wettig, PoS
LATTICE2015 , 036 (2016),[arXiv:1512.04506].[58] D. Richtmann, S. Heybrock and T. Wettig, PoS
LATTICE2015 , 035 (2016),[arXiv:1601.03184].[59] P. Georg, D. Richtmann and T. Wettig, PoS
LATTICE2016 , 361 (2017),[arXiv:1701.08521].[60] SciDAC, LHPC, UKQCD, R. G. Edwards and B. Joo, Nucl. Phys. Proc. Suppl. , 832(2005), [arXiv:hep-lat/0409003], [,832(2004)].[61] P. Arts et al. , PoS
LATTICE2014 , 021 (2015), [arXiv:1502.04025]., 021 (2015), [arXiv:1502.04025].