Charmonium-like resonances with $J^{PC}=0^{++},2^{++}$ in coupled $D\bar D$, $D_s\bar D_s$ scattering on the lattice

We present the first lattice investigation of coupled-channel $D\bar D$ and $D_s\bar D_s$ scattering in the $J^{PC}=0^{++}$ and $2^{++}$ channels. The scattering matrix for partial waves $l=0,2$ and isospin zero is determined using multiple volumes and inertial frames via L\"uscher's formalism. Lattice QCD ensembles from the CLS consortium with $m_{\pi}\simeq280$ MeV, $a \simeq 0.09 $ fm and $L/a=24,~32$ are utilized. The resulting scattering matrix suggests the existence of three charmonium-like states with $J^{PC}=0^{++}$ in the energy region ranging from slightly below $2m_D$ up to 4.13 GeV. We find a so far unobserved $D\bar D$ bound state just below threshold and a $D\bar D$ resonance likely related to $\chi_{c0}(3860)$, which is believed to be $\chi_{c0}(2P)$. In addition, there is an indication for a narrow $0^{++}$ resonance just below the $D_s\bar D_s$ threshold with a large coupling to $D_s\bar D_s$ and a very small coupling to $D\bar D$. This resonance is possibly related to the narrow $X(3915)$/$\chi_{c0}(3930)$ observed in experiment also just below $D_s\bar D_s$. The partial wave $l=2$ features a resonance likely related to $\chi_{c2}(3930)$. We work with several assumptions, such as the omission of $J/\psi\omega$, $\eta_c\eta$ and three-particle channels. Only statistical uncertainties are quantified, while the extrapolations to the physical quark-masses and the continuum limit are challenges for the future.


Introduction
Since the discovery of the J/ψ meson in 1970s a multitude of charmonium bound states and resonances have been found with energies ranging up to almost 5 GeV. A simple cc quark model provides a reasonable description of the levels below the strong decay thresholds and also some of the states above, however, there are clearly too many states to fit into this picture. Some mesons, such as the charged Z c states certainly have additional quark content, while for other states the interpretation is not so clear. On the theory side the nature of these states is being explored in tetraquark, molecular, and hybrid meson models, among others, while on the experimental side insight is provided by establishing their quantum numbers, decay modes and widths. Lattice QCD studies of the charmonium 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 investigating their internal structure.
In this work we focus on the isoscalar channel I(J P C ) = 0(0 ++ ) in the region up to 4.13 GeV for which there are a number of open questions. The ground state, χ c0 (1P ), found well below the DD threshold is interpreted as the 3 1P 0 cc level of the quark model and is the only well established state. In the energy region around 3.9 GeV, above the threshold, one expects a corresponding excited state. So far, three hadrons have been observed with the possible assignment of J P C = 0 ++ : the X(3860), a broad resonance detected by Belle [1,2], and two narrow resonances just below the D sDs threshold -the χ c0 (3930) discovered in the DD channel by LHCb [3,4] and the X(3915) observed through it's decay into J/ψω [5][6][7][8] (with the assignment of J P C = 0 ++ or 2 ++ ). While the latter two resonances could be the same state, their narrowness may indicate exotic content, where X(3915) has been interpreted as ccss meson in Ref. [9]. Predictions have also been made for an additional, as yet unobserved, bound state just below the DD threshold [10,11].
The determination of the low lying charmonium spectrum on the lattice is relatively straightforward, with the energy levels being directly accessed from correlation functions measured on the configurations generated in the Monte-Carlo simulation. Systematics arising from finite lattice spacing and simulating with unphysical light (sea) quark masses must be addressed by carrying out a continuum and quark-mass extrapolation. Near and above threshold, the analysis is considerably more challenging with information on the masses and (for resonances) also the widths being inferred from scattering amplitudes which can be obtained from the finite volume spectra via the Lüscher method [12][13][14]. Two-particle interpolators must be included in the basis of operators for the construction of the correlation functions in order to reliably determine these spectra. Simulating charmonia in flight provides additional levels with which to probe the scattering matrix, however, the identification of the continuum spin and parity quantum numbers of the levels is complicated due to the reduced symmetry on the lattice. In addition, for the energy range of interest both the DD and D sDs thresholds must be considered leading to a coupled-channel scattering analysis.
So far, the coupled-channel scattering matrix has been extracted for several lightmeson systems, for example, πK, ηK [15,16], πη, KK [17] and ππ, KK, ηη [18] by the Hadron Spectrum Collaboration. In the heavy sector, there has been one investigation of Dπ, Dη, D sK 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ödinger equation with potentials determined on the lattice [21,22]. Pioneering works such as these were limited to a single lattice spacing and unphysical light-quark masses.
The charmonium scalar channel has previously been studied by some of the authors considering only DD scattering with total momentum zero [23]. Here we present a lattice study of scattering in the coupled-channels DD and D sDs with quantum numbers I = 0 and J P C = 0 ++ , 2 ++ . This represents the first determination of the coupled-channel scattering matrix from lattice QCD in the charmonium system with isospin zero. Two lattice volumes are employed for the charmonium system at rest and in flight. This analysis uses the same lattice setup as our previous article on the identification of the spin and parity of the single hadron spectrum [24] and the investigation of single channel DD scattering for J P C = 1 −− and 3 −− [25]. While the present study represents a significant improvement on previous work, some simplifications remain and a comparison of the results for the masses and widths with experiment is qualitative. Within the energy range of interest, additional scattering channels, such as the J/ψω, η c η and those involving three particles, could in principle also be relevant. The effects of these channels will be investigated in the future, along with systematics associated with finite lattice spacing and unphysical light quark masses.
The remainder of the paper is organized as follows. We begin by reviewing the essential general aspects of one-channel and two-channel scattering in Section 2. The details of the lattice setup and methodology are given in Section 3 and the single-and two-meson interpolators used in the correlation functions are discussed in Section 4. Simplifying assumptions made in this study are summarized in Section 5. The first step in extracting the scattering amplitudes is to compute the finite-volume spectra from the correlation functions. Our analysis and the final spectra are presented in Section 6. An overview of determining the scattering amplitudes from the lattice eigen-energies is provided in Section 7. Our results for the J P C = 0 ++ and 2 ++ channels are detailed in Section 8 and the relation to states observed in experiment is discussed in Section 9. Finally, Section 10 presents our conclusions. More details are given in several Appendices.

Generalities on scattering matrices, poles, hadron masses and widths
The masses and widths of strongly-decaying resonances should be inferred from the study of scattering processes where these resonances appear. In this section, we briefly review relevant concepts regarding scattering matrices, complex energy planes, pole singularities, hadron masses, and their decay widths. The first part lists definitions and notations for the scattering amplitudes, the phase space factors, etc.. The second part discusses naming conventions for various Riemann sheets, pole singularities in the complex energy plane and their relation to the hadron properties.

Scattering matrices for real energies
The unitary scattering amplitude S for one-channel scattering (DD or D sDs ) of spin-less particles in partial wave l is generally parametrized in terms of the energy-dependent phase shift δ(E cm ), where ρ ≡ 2p/E cm , p denotes the momentum of the scattering particles in the center-ofmomentum frame and t is the scattering amplitude. The factors p −2l in front ofK −1 lead to smooth behavior close to the threshold. In the case of t exhibiting simple Breit-Wigner type behavior,K −1 /E cm falls linearly as a function of E 2 cm , (2. 2) The phase shift equals π/2 at E cm = m R , while the width Γ(E cm ) is parametrized in terms of 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 DD and D sDs in partial wave l, the scattering matrices S are energy-dependent 2 × 2 unitary matrices, 3) The momenta of D and D s in the center-of-momentum frame are denoted by p 1 and p 2 , respectively. t is the scattering matrix andK(E cm ) is a real symmetric matrix. We follow the definition of t by the Hadron Spectrum Collaboration (e.g. [15]) and the definition of K from Ref. [26] 1 .

Continuation to complex E cm , Riemann sheets and poles
In experiment and lattice QCD simulations the scattering matrices S(E cm ) are determined for real energies. The theoretical interpretation in terms of (virtual) bound states and resonances is conventionally made via the poles in the t-matrix, analytically continued to the complex s-plane. The feature that leads to interesting physics is the square root branch cut related to ρ = 2p/E cm = 1 − (2m) 2 /E 2 cm starting from the threshold connecting the physical Riemann sheet (or sheet I), conventionally chosen to have Im(ρ) > 0, to the I II 2m 2 Im(E cm ) Re (E cm )   I  II   II  III   2m 1  2m 2   2 Im(E   cm ) Re(E cm ) Figure 1: Sketch of the pole locations in the scattering matrix t that typically affect the experimental rates on the physical axes (denoted by the cyan line) for one-channel scattering (left) and two coupled channels (right). The Roman numbers indicate the Riemann sheets where the poles are located according to Eq. (2.4). Poles immediately below a threshold, indicated by crosses, can also have observable effects on the physical axes above the respective threshold.
unphysical Riemann sheet (or sheet II), which has Im(ρ) < 0. For a two channel system, there will be four Riemann sheets, such that Bound states, virtual bound states and resonances are related to pole singularities of t in the complex s-plane. These poles affect the physical axes, indicated by the cyan line in Fig. 1, along which the experimental measurements are made. Fig. 1 presents a schematic picture of various pole locations in our study, that can affect scattering amplitudes/matrices along the physical axes for one-channel and two-channel scattering. The location of the poles are related to the masses and widths via E p cm = m − i 2 Γ for resonances and E p cm = m for the (virtual) bound states.
In the close vicinity of the pole, the scattering matrix has the energy dependence 5) and the residue (c i c j ) can typically be factorized into the couplings c i ( 2 ), whose relative size is related to the branching ratios of a resonance (associated with the pole) to both channels i = 1, 2.   [29]. Open boundary conditions in time are imposed [30] and the sources of the correlation functions are placed in the bulk away from the boundary. We remark that these correlation functions do not show any effects related to the finite time extent in the time regions analyzed. For H105 we use replica r001 and r002 for which the issue of negative strange-quark determinants described in Ref. [31] is not of practical relevance. For our 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 ensemble are shown in Table 1. Note that the chosen quark-mass trajectory leads to a larger than physical m u/d and a smaller than physical m s . This means that the splitting between the DD and D sDs thresholds is smaller than in experiment, emphasizing the need for a coupled-channel analysis. We employ the charm-quark hopping parameter κ c = 0.12315 corresponding to a charm-quark mass m c and spin-averaged 1S-charmonium mass M av that are slightly larger than their physical values. For estimates of the statistical uncertainty we 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 matrices are averaged over several source-time slices and momentum polarizations to increase the statistical precision. Note that all quoted uncertainties are statistical only, and that results quoted in MeV have been obtained using the central value of the lattice scale without propagating its statistical or systematic uncertainties into the results.
For hadrons with charm quarks, non-negligible discretization effects are observed when computing the dispersion relation on lattices with a ≈ 0.086 fm. A comparison of the finite momentum lattice energies and the continuum dispersion relation for the D meson on the two ensembles utilised in this work is given in Table II of Ref. [25]. The deviations found are small but statistically significant. A similar picture is observed for the D s meson. Note that, these deviations may spoil the finite-volume analysis outlined in Section 7, which assumes the continuum dispersion relation. In particular, it is important to ensure that if the energy shifts observed with respect to nearby non-interacting two hadron levels are zero then the resulting phase shift arising from the finite-volume analysis is also zero. In order to achieve this and mitigate the affect of the discretisation effects we adopt the analysis strategy described in Sect. IV.B. of Ref. [25]. Below we reiterate the most important details of the method.
First the energy shift of each interacting eigenstate with respect to a nearby noninteracting two-hadron level where p 1,2 = n 1,2 2π L , p 1 + p 2 = P and s denotes the bootstrap sample. Here, (E lat ) s is the energy of the interacting two-hadron system, while (E lat H i ( p i ) ) s is the energy of a single hadron (either D or D s meson in this paper) with momentum p i measured on the lattice. We then use (E calc ) s = (∆E lat ) s + E cont as input to the quantization condition (see Eq. (7.1)) for each bootstrap sample s. The energies (E cont H i ( p i ) ) s are computed from the continuum dispersion relation using the lattice momenta p 1,2 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 → 0 by construction. The non-interacting levels are chosen via an analysis of the overlap factors by identifying those levels that are dominated by the corresponding two hadron interpolators. 3 In the case where more than one suitable nearby level was identified, we found the results obtained for E calc were consistent. A comparison of E calc with E lat is presented in Appendix B.
For further details of the lattice methodology, in particular of the setup for computing the quark propagators with the (stochastic) distillation method [32,33] we refer the reader to our previous papers [24,25].

Interpolators
The main aim of this work is to investigate the coupled-channel DD-D sDs scattering amplitudes and cross-sections in the channel I(J P C ) = 0(0 ++ ) in the energy range encompassing the DD threshold up to 4.13 GeV. Following Lüscher's approach [12][13][14]34], this requires a reliable extraction of the finite-volume charmonium spectrum below 4.13 GeV on several different volumes and/or in different momentum frames. In this study, we consider the charmonium spectrum in four different lattice irreducible representations (irreps) Λ: The squared momenta | P | 2 in the lab frame are given in units of (2π/L) 2 . Charge conjugation C = + is a good quantum number in all frames and hence is suppressed for 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. The first three irreps are relevant for an investigation of the J P = 0 + channel. The irreps The single-and two-meson interpolators utilized in each lattice irrep Λ (P )C considered in this study. We use the simplified notation M 1 ( p 2 1 )M 2 ( p 2 2 ) for the two-meson interpolators with the momentum p i of each meson (i = 1, 2) given in units of 2π/L. The full expressions are omitted for brevity. N ops indicates the number of operators of each type employed.
in the moving lab frames also receive contributions from states with J P [λ] = 2 + [0] and 2 ± [2] within the energy range of interest 4 . The analysis of the spectrum in the B 1 irrep constrains the parameters for DD scattering with l = 2. This partial wave inevitably contributes to the finite-volume spectrum of irreps A 1 with P > 0. We utilize a large basis of single-meson as well as two-meson interpolators in the above irreps to reliably determine the relevant low energy spectrum.
As in our previous publications [24,25], we construct the single-meson interpolators following 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-volume irreps considered. The procedure discussed in Ref. [24] guides us in assigning the quantum numbers J P [λ] to the extracted energy levels and aids us in selecting the levels relevant for the amplitude analysis.
The DD as well as D sDs interpolators are constructed following the same procedure as in Ref. [25]. The momentum combinations implemented in this study are given in Table 2. The two operators for D (s) (0)D (s) (0) differ in terms of the gamma matrices employed: γ 5 or γ t γ 5 for each meson. Similarly, for D * (0)D * (0) and J/ψ(0)ω(0), two operators are constructed by employing γ i or γ t γ i for the spin structure. Only one eigenstate related to J/ψ(0)ω(0) or D * (0)D * (0) is expected in the non-interacting limit. 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 needs additional care and we follow the induced representation method described in Appendix A2 of Ref. [37]. In the | P | 2 = 2 frame, for example, we implement three linearly independent J/ψ(2)ω(0) operators and observe three almost-degenerate eigenstates. These operators are not included when extracting the finite volume spectrum for the amplitude analysis, as discussed in Section 5.

Assumptions and simplifications in the present study
This study is performed using lattice gauge ensembles with two different physical volumes at a single lattice spacing and at unphysical quark masses (the resulting masses of key hadrons are given in Table 1). As a consequence, only a qualitative comparison of the results can be made with experiment.
Unlike for light hadrons [38], scattering studies in the charmonium sector are still at an early stage. For the physical states we are interested in, a three-particle channel and multiple two-particle channels are open and all could, in principle, be relevant. One possible approach is to simulate at very heavy pion (and kaon) mass, such that the number of relevant decay modes is reduced to a few two-hadron modes, which can then be fully explored. This approach has the disadvantage that the quark masses are far removed from their physical values, making a comparison to experiment a challenge. We opt for a strategy where we simulate at a moderate pion mass of 280 MeV and take into account the scattering channels expected to be most relevant for the physics close to the opencharm threshold(s). Some additional channels are neglected (as discussed below), however, our assumptions about which thresholds are relevant can be relaxed successively in future calculations.
Neglecting certain scattering channels in our study is relevant in two different ways, which could be seen as two different assumptions. Some channels are already neglected in constructing the correlator matrices. This implicitly assumes that the neglected multihadron correlators would simply yield additional energy levels rather than significantly modifying the extracted spectrum. Additionally, we assume that the resulting energy levels can be analyzed with the (coupled-channel) formalism for the channels we deem to be dominant, which might fail if there is significant coupling to neglected channels. Beyond the scattering channels investigated explicitly, our current study includes J/ψω and (some) D * D * operators in the interpolator basis. Due to the poor signal obtained for light isoscalar mesons, the energy levels close to the non-interacting J/ψω levels are not very precisely determined and would not provide strong constraints on the scattering matrix. In particular, almost all energy levels dominated by J/ψω interpolators fall within one standard deviation of the non-interacting J/ψω energies, and -apart from the additional energy levels which appear -the other finite-volume energies do not shift significantly when including these interpolators. Section 6 will present the finite volume spectrum up to 4.13 GeV based on all the operator types in Table 2, apart from the J/ψω operators (see Fig. 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 remark that this decay channel has not been observed for any of the experimental candidates mentioned in the introduction (and discussed in more detail in Sect. 9). Operators with more than two hadrons are also not implemented. The lowest three-hadron threshold is for the decay into χ c0 ππ at 4.02 GeV. This threshold is within the energy region we consider, close to the upper end.
The analysis of DD scattering with l = 2 assumes that the coupling to the channel D sDs with l = 2 is negligible in the analyzed energy region and hence is omitted. We also neglect the coupling to DD * with l = 1, which contributes to irrep B 1 . The DD * threshold opens at 4.0 GeV, while the lowest non-interacting level D * (2)D(1) would appear at E cm 4.2 GeV and 4.1 GeV on the N L = 24 and 32 ensembles, respectively, which is at the upper limit of the analyzed region (see Fig. 2). We also assume negligible effect of the D * D * channel with threshold at 4.1 GeV.
As in all studies of charmonium-like resonances to date, charm annihilation Wick contractions are omitted. All the remaining contraction diagrams arising from the singleand two-meson operators in our basis (shown in Fig. 1 of Ref. [23]) are computed following the procedures described in our previous publications [24,25].
We stress that we determine the finite-volume spectra at a single lattice-spacing and are therefore unable to quantify the uncertainty associated with the lattice discretization. In particular, the uncertainty arising from the heavy quark discretization may be nonnegligible. As discussed in the previous section, the dispersion relation deviates from the continuum 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 with the net result that even at physical light-quark masses the open-charm thresholds can be shifted with respect to the measured charmonium states at finite lattice spacing.

Determination of the finite-volume spectrum
This section presents the eigen-energies E n that will be used to determine the scattering matrices. The energies are obtained from the correlation matrices C ij (t) = O i (t)O † j (0) via the widely-used variational method. This involves solving the generalized eigenvalue problem C(t)u (n) (t) = λ (n) (t)C(t 0 )u (n) (t) for the eigenvalues λ (n) (t) and the eigenvectors u (n) (t) [41][42][43]. We use the reference time t 0 /a = 3 or 4. The eigen-energies are extracted from 1-exponential fits to the eigenvalues λ (n) (t) = A n e −Ent with the fit range, in most cases, starting between timeslices 10 and 12.
The finite-volume spectrum of the charmonium system with isospin I = 0 and C = +1 is shown in Fig. 2. We present the spectrum in the center-of-momentum (cm) frame 1 , B 1 and total momenta | P | 2 = 0, 1, 2. These irreps give information on the charmonium(like) states and D (s)D(s) scattering in the channels with J P C = 0 ++ , 2 ±+ (see Eq. 4.1). The energies indicated by the black-circles are used to extract information on D (s)D(s) scattering. These energies are near or above the DD threshold and are precise enough to reliably resolve the energy-shifts with respect to the non-interacting energies of D (s)D(s) (indicated by the solid lines). The light-blue circles are the energy levels related to ground-state charmonia with J P = 2 ± .   Figure 2: The eigen-energies in the center-of-momentum frame (E cm ) for the charmoniumlike system with I = 0 and C = +1. Results are presented for irreducible representations 1 , B 1 and total momenta | P | 2 = 0, 1, 2, which give information on the channels with J P C = 0 ++ , 2 ±+ . The data points correspond to the eigen-energies obtained from the lattice simulation: the black circles are used to extract the coupled-channel scattering matrices for DD − D sDs , while the blue circles are omitted from the scattering analysis. The solid and dashed red (green) lines correspond to discrete DD (D sDs ) eigen-energies in the non-interacting limit: solid lines correspond to the operators that are implemented, while dashed lines correspond to the lowest-lying energies from operators that are not implemented. Dotted lines represent thresholds. The data points indicated by the light blue circles correspond to ground-state charmonia with J P C = 2 ++ and 2 −+ , which appear at m 3.56 GeV and 3.83 GeV, respectively. Some points are shifted horizontally slightly for clarity.

Determining scattering matrices from lattice finite-volume energies
The bound states and resonances are inferred from the scattering matrices as briefly reviewed in Section 2. The infinite-volume scattering matrix S(E cm ) is related to the finitevolume two-hadron spectrum for real energies E cm above the threshold and somewhat below it through the well-known Lüscher relation [12][13][14]. The eigen-energies of the coupled channel DD − D sDs system given in the previous section provide information on the in 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 determineK(E cm ) using the publicly available package TwoHadronsInBox [46].
The relation between discrete lattice eigen-energies E cm andK-matrix for coupledchannel scattering is referred to as the quantization condition [46] Both terms in the determinant are matrices in the space of partial waves l, l and channels i, j (DD, D sDs or both), and the determinant is evaluated over both indices.K l;ij δ ll is an 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 P ,Λ ll ;i (E cm ) are known box-functions [46] that are in general non-diagonal in the partial wave index.
In one-channel scattering and when only partial wave l contributes, relation (7.1) simplifies toK −1 (E cm ) = B P ,Λ (E cm ), since the argument of the determinant is a 1 × 1 matrix. The values of K −1 (E cm ) will be shown as points in figures for one-channel scattering. For two coupled channels, for the case when only partial wave l contributes, the determinant equation (7.1) provides one relation betweenK 11 (E cm ),K 22 (E cm ) andK 12 (E cm ) for each energy level, complicating the determination of those functions. Therefore, we follow the strategy proposed in Ref. [44], where theK ij (E cm ) are parametrized as functions of the energy. In this strategy, theK-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 linear dependence on s is expected to be a good approximatioñ Such a parametrization is equivalent to a Breit-Wigner parametrization in the resonance region and is also similar to the well-known effective range expansion K −1 ij (s) = c ij + d ij p 2 near threshold, where p is the momentum of the scattering particles in the center-ofmomentum frame. We determine the parameters a ij and b ij following the strategy discussed above, using the determinant residual method proposed in [46], which is briefly described in Appendix C. A posteriori, we always verify that the resulting parametrization predicts via Eq. (7.1) the same number of eigen-energies observed in the actual simulation in the relevant energy range; this is shown in the Ω plots for some fits in Appendix D. This procedure will be followed for the extraction of the coupled-channel scattering matrix as well as for one-channel scattering.
This study is based on the parametrization in Eqn. (7.2). Alternatively, one could parametrizeK ij (s) itself with common pole terms in both channels, such as those tabulated in Table IV of Ref. [47]. We have performed fits with different parametrizations ofK ij (s) (single pole, double pole, triple pole, poles with polynomial terms, etc.). We find that fits (for coupled DD − D sDs scattering) with a single-pole in the higher energy region are not consistent with our data. Including two or more poles/resonances leads to fits with six or more parameters. We observed that the data used in this work is insufficient to accommodate such a large number of parameters and hence such an analysis is beyond scope of this work. An investigation of the model-independence of the findings presented here requires extending the lattice calculation to include a larger set of ensembles with high statistics.
The box-function B P ,Λ ll ;i (E cm ) can have off-diagonal elements for l = l due to the lack of rotational symmetry in a finite box. This will result in contributions from multiple partial waves in the quantization condition Eq. (7.1) for a given lattice irrep Λ. We consider partial waves l = 0 and l = 2 and ignore contributions from l ≥ 3, which is a reasonable assumption in the energy region considered for the respective irreps. In this case the only non-diagonal elements B ll among the A 1 (P 2 = 0) irreps that are nonzero are B P =001,A 1 02 and B P =110,A 1 02 . These will be taken into account in the analysis of Section 8.4.2.

Results for various channels and energy regions
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 and with varying assumptions. The energy range from slightly below 2m D up to 4.13 GeV is divided into smaller intervals, where the elements of the coupled DD − D sDs scattering matrix are separately parametrized according to Eq. (7.2) or as a constant. The details of the parametrizations and the results are presented in separate subsections below, while information on the energy levels considered in each case is given in Appendix D.
A single description of the whole energy region requires a finite-volume analysis involving many more parameters, which results in more challenging and unstable fits. Such an analysis is beyond the scope of the current investigation. Our inferences and conclusions are based on the finite-volume analysis of separate energy regions. Similar parametrizations to those employed for the separate energy regions, are employed collectively to a wider energy range in Appendix E as an additional consistency check.
8.1 DD scattering with l = 0 near threshold The narrow energy region near the DD threshold is significantly below the D sDs threshold and can be treated in a one-channel approach. We employ the parametrization in Eq. (7.2) which is equivalent to the effective range expansion p cot δ 0 = 1/a 0 + r 0 p 2 /2 near threshold. Four lattice energy levels with E cm closest to 2m D (listed in Appendix D.1) are utilized to determine the parameters via the quantization condition (7.1). We find  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 scattering matrix t (2.1) has a pole on the real axis below threshold on sheet I The lhs of the second equation is shown as the red line in the figure, while the rhs is indicated by the orange line. The bound state occurs at the value E cm = m, where the two curves intersect. The slope of p cot δ at the intersection, is smaller than the slope of −|p|, as required for an s-wave bound state (see Section VC of [25]). The location of the pole in the scattering matrix is shown in Fig. 3c. The bound state appears just below the DD threshold with the binding energy We denote this state by χ DD c0 , indicating it has J P C = 0 ++ and a strong connection to the DD threshold. This state comes in addition to the conventional χ c0 (1P ), which is found significantly below threshold. Experiments cannot explore DD scattering below threshold, however, a closeby bound state below threshold could be identified experimentally through a sharp increase of the rate just above threshold. Fig. 3b shows a dimensionless quantity ρ|t| 2 related to the number of events N DD ∝ pσ ∝ ρ|t| 2 expected in experiment. It features a peak above threshold, which increases much more rapidly than the phase space.
Such a DD bound state was not claimed by experiments so far. A similar state was predicted in phenomenological models [10,48,49], and some indication for it was suggested in the experimental data [11,50] and in data from the lattice simulation of Ref. [23]. A more detailed discussion follows in the summary in Section 9.
Details of the fit (8.2) and some variations thereof are provided in Appendix D.1. In these fits, the ensemble average of the data gives rise to a bound state, while a very small proportion of the bootstrap samples instead produce a virtual bound state. This indicates that our lattice results, at the employed quark masses, favour the existence of a bound state. However, with the present statistical accuracy, one cannot completely rule out the existence of a virtual bound state. The robust conclusion is that we observe a significant DD interaction near threshold, leading to one state just below threshold. Such a state leads to an increase of the DD rate above threshold irrespective of whether it is a bound or a virtual bound state. Note that it is not known whether this state would also feature in a simulation with physical quark masses.
8.2 D sDs scattering with l = 0 near threshold in the one-channel approximation The D sDs channel carries the same quantum numbers as DD necessitating the consideration of coupled-channel scattering. In this subsection we aim to get a rough estimate of D sDs scattering in the one-channel approximation, which will also provide initial guesses for the parameters when coupled channel scattering is considered in Section 8.4. The D sDs scattering near threshold is parametrized by We employ the quantization condition (7.1) together with four lattice energies close to this threshold that are dominated by D sDs interpolators (listed in Appendix D.2) and obtain The resulting fit is shown in Fig. 4a. The scattering matrix has a bound state pole at the energy E cm = m where condition (8.3) is satisfied, see Fig. 4c. Again, the slope of p cot δ is smaller than the slope of −|p| at the position of the pole, as required for an s-wave bound state (see Section VC of [25]  that we denote χ DsDs c0 , indicating it has J P C = 0 ++ and a strong connection to the D sDs threshold. This state is responsible for the significant increase in the D sDs rate shown in Fig. 4b just above threshold. In order to search for the χ DsDs c0 in experiment an exploration of the D sDs invariant mass near threshold would be invaluable. In one-channel D sDs scattering, considered here, the state is decoupled from DD, while it will become a narrow resonance and acquire a small width when the coupling to DD is considered in Section 8.4. Two candidates χ c0 (3930) [3] and X(3915) [2] (which may correspond to the same state) have already been observed in experiment just below the threshold 2m exp Ds 3936 MeV; they have a small coupling to DD and a small width. If the D sDs bound state (8.8) corresponds to χ c0 (3930) and/or X(3915), it naturally explains both features as will be discussed in Section 9.  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 scattering amplitude using Eq. (7.1). We consider the energy region encompassing the 2 ++ resonance and neglect the coupling to D sDs scattering with l = 2, which we assume to be negligible in this region. The scattering amplitude is parametrized by the Breit-Wigner form (2.2)  Fig. 5a. The mass m J2 corresponds to the energy where the phase-shift reaches π/2, which is close-to the 2 ++ resonance mass obtained from the pole position below, while the coupling g J2 is related to its width as shown in Eq. (2.2).
The position of the pole E p cm of the scattering matrix (8.9) on sheet II provides a better way of determining the resonance mass m and width Γ. We obtain The pole is plotted in Fig. 5c. This leads to the lowest J P C = 2 ++ resonance above DD threshold with where g parametrizes the width Γ = g 2 p 5 /m 2 . This likely corresponds to the wellestablished resonance χ c2 (3930) = χ c2 (2P ) [2]; a detailed comparison with experiment is made in Section 9. The resonance mass and the coupling obtained from the pole and from Eq. (8.10) are consistent, which is expected for a narrow resonance. The next higher 2 ++ charmonium is estimated 6 to be near E cm 4.2 GeV, which is above our region of interest. We assume it to be narrow and to have a negligible effect on the analysis of the lower-lying 2 ++ resonance. Finally, we turn to the coupled DD − D sDs scattering. We focus on the energy region E cm 3.93 − 4.13 GeV near the D sDs threshold and we find an indication for several interesting hadrons. The scattering matrix for partial wave l = 0 is parametrized as with the off-diagonal element held constant in E cm . Of the two equivalent parametrizations shown above, we will utilize the one on the rhs. The 5 parameters in Eq. (8.13) are determined using all levels of irreps A (+) 1 within the energy region E cm = 3.93 − 4.13 GeV displayed in Fig. 2: there are 14 levels from three frames with P 2 = 0, 1, 2 and from two spatial volumes N L = 24, 32 (see the black circles in the figure). The quantization condition (7.1) for A (+) 1 irreps depends on the scattering amplitudes for l = 0, which we aim 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.

Analysis omitting l = 2
In the first analysis we omit the contribution of the partial-wave l = 2. This is expected to be a fair approximation since l = 2 effects DD scattering only in the narrow 2 ++ resonance region that is at the upper end of the current energy range of interest.   The coupling between channels a 12 is non-zero but small. We also performed a study where five parameters for l = 0 and two parameters for l = 2 are fitted simultaneously using 18 levels of irreps A   3) on the physical axes can be described by the phase shifts δ DD , δ DsDs and inelasticity η, which are shown in Fig. 8.
We find that the DD and D sDs channels are not strongly coupled, which can be seen from the inelasticity η 1 in Fig. 8 ( 8 ) and from the smallness of the off-diagonal element a 12 in Eqs. (8.15) and (8.13). Our results suggest there are two 0 ++ resonances in this energy region: a narrow resonance dubbed χ DsDs c0 just below the D sDs threshold and a broader one denoted by χ c0 .
The broader 0 ++ resonance χ c0 is related to the pole indicated in red on sheet III This pole affects the scattering amplitude on the physical axes above the D sDs threshold and is responsible for a peak around 3.98 GeV in the DD → DD rate shown in the left pane of Fig. 7. The presence of this pole is also reflected in the phase shift δ DD 0 , which increases gradually starting from 2m Ds as is evident in the left pane of Fig. 8. The nearby pole on sheet II does not have a significant influence on the physical scattering above the second threshold. The pole residues indicate that this state decays predominantly to DD, while the decay to D sDs is suppressed, as evidenced by |c 1 | |c 2 |, see in the last two rows of Fig. 6 for the pole presented in red. The resonance parameters are (8.17) where M av = 1 4 (3m J/ψ +m ηc ), and the coupling g parametrizes the full width Γ = g 2 p D /m 2 . The possible relation of this state to the broad resonance χ c0 (3860) discovered by Belle in 2017 [1,2] is discussed in Section 9.
The narrow 0 ++ resonance χ DsDs c0 near the D sDs threshold is related to the pole on sheet II, indicated by the top-filled orange symbols in the first row of Fig. 6. Its location relative to the threshold is given by This resonance is related to the bound state in the analysis of D sDs -scattering in the onechannel approximation of Section 8.2. The pole on sheet II and the nearby pole on sheet IV correspond to this resonance and are mutually exclusive across the bootstrap samples. Further details on this can be found in Appendix D.4. It is clear from Figure 7 that the resonance pole leads to a sharp rise in the D sDs → D sDs and DD → D sDs rates just above 2m Ds . The increased DD → D sDs rate is also responsible for a dip in the DD → DD rate at 2m Ds and all three features should be used as a signature for experimental searches of this state. Note that the magnitude of the D sDs → D sDs peak above 2m Ds is larger when the pole is closer to the threshold. χ DsDs c0 couples predominantly to D sDs and very weakly to DD (one can see that |c 2 | 2 |c 1 | 2 in Fig. 6). The mass difference of the state with respect to the threshold and its narrow total width Γ = g 2 p D /m 2 parametrized in terms of g are On the experimental side, the newly discovered χ c0 (3930) [3] and the X(3915) [2] lie near the D sDs threshold and have very small or zero decay rate to DD. The indication for a D sDs state in our study explains both properties, as detailed in Section 9. The parameters of the scattering matrix obtained from the analysis including or excluding l = 2 are similar, with the exception of a 22 and b 22 . These parametrize D sDs → D sDs and differ between the coupled-channel analysis and the one-channel approximation, however, both analyses lead to a state just below the D sDs threshold on the real axis (see Fig. 4c) or slightly away from it. The conclusion that there is a near-threshold pole is robust, while its exact location and the effect on physical scattering need to be investigated in a simulation with higher statistics and a better control of systematic uncertainties.
2m Ds 2m D sheet I sheet II sheet III sheet II J PC = 0 + + J PC = 2 + + J PC = 0 + + J PC = 2 + + Figure 9: Pole singularities of the scattering amplitude/matrix in the complex energy plane, which are associated with the hadrons predicted in this work. The pole related to the J P C = 2 ++ resonance appears on sheet II, as we have neglected the l = 2 partial wave contribution from D sDs scattering.
the identification with experimental states is unambiguous, while the other states found are denoted χ c0 , χ DD c0 and χ DsDs c0 . The subscripts c0 and c2 indicate the assignment of J P C = 0 ++ and 2 ++ , respectively. The location of the poles in the complex energy plane related to these hadrons are given in Fig. 9, while the corresponding masses are compared to experiment in Fig. 10 The resonance decay widths depend on the phase space p 2l+1 evaluated for the meson momenta (in the cm-frame) at the resonance energy, which in turn depends on the position of the threshold. The latter is different in the simulation and in experiment. Therefore it is customary to compare the coupling g that parametrizes the full width of a hadron Γ ≡ g 2 p 2l+1 D /m 2 with l = 0, 2 for J P C = 0 ++ , 2 ++ , (9.1) as g is expected to be less dependent on the quark masses than the width itself. Note that it is not known whether this bound state would also feature in a simulation with physical quark masses. Such a state has not been claimed by experiments.
The existence of a shallow DD bound state dubbed X(3720) was already suggested by an effective phenomenological model in Ref. [10] 9 featuring also exchanges of vector mesons. Ref. [11] indicates that there may already be some evidence for such a state in the DD invariant mass distribution from Belle [50], which shows an enhancement just above threshold. The DD rate from Babar [51] also shows a hint of enhancement just above threshold (see Fig. 5 of [51]). In a molecular picture, a 0 ++ state is expected as a partner of X(3872) via heavy-quark symmetry arguments [48,49]. A similar virtual bound state with a binding energy of 20 MeV follows from the data of the only previous lattice simulation of DD scattering [23] 10 .
2 ++ resonance and its relation to χ c2 (3930) We find a resonance with J P C = 2 ++ in l = 2 DD scattering with the following properties This is most likely related to the conventional χ c2 (3930) resonance (also called χ c2 (2P )) discovered by Belle [57] exp χ c2 (3930) : m − M av = 854 ± 1 MeV , g = 2.65 ± 0.12 GeV −1 . (9.8) Here g parametrizes the width Γ = g 2 p 5 D /m 2 . The masses are reasonably close, while the coupling from lattice QCD is larger that in experiment. However, the couplings are also not inconsistent given the large statistical error from our simulation and the unquantified systematic uncertainties discussed in Section 5.
Broad 0 ++ resonance and its possible relation to χ c0 (3860) This resonance couples mostly to DD and has a very small coupling to D sDs . Its resonance parameters are  Table 4 of Ref. [10]. 10 The presence of this state was not mentioned in Ref. [23], as such virtual bound states were not searched for. 11 For the couplings calculated from the experimental values we vary both the mass and width by ±1σ and take the maximal positive and negative deviations as the uncertainties -25 -based on the following arguments: The mass and coupling are reasonably consistent with experiment, in particular, when considering the experimental errors and the systematic uncertainties in the lattice results. The mass is also close to the value obtained from the only previous lattice study of DD scattering [23] 12 , although the width and coupling are larger in the present work.
Narrow 0 ++ resonance χ DsDs c0 and its possible relation to χ c0 (3930), X(3915) We find a narrow 0 ++ resonance near the D sDs threshold. It has a large coupling to D sDs and a very small coupling to DD. The latter is responsible for its small decay rate to DD and the small total width. This state corresponds to the bound state in one-channel D sDs scattering discussed in Section 8.2. We compare the resulting resonance parameters The χ c0 (3930) with J P C = 0 ++ was very recently discovered in DD 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 DD was not observed [2]. They might represent the same state if X(3915) is a scalar. Both experimental states lie just below the D sDs threshold. One would naturally expect 0 ++ states with this mass to be broad, given the large phase space available to DD decay. Their narrow widths can only be explained if their decay to DD is suppressed by some mechanism. If the resonance found on the lattice is indeed related to X(3915)/χ c0 (3930), our results indicate that this state owes its existence to a large interaction in the D sDs channel near threshold, which naturally explains why its width is small and its decay to DD is suppressed. Note that a detailed quantitative comparison of lattice and experimental results in Eqs. (9.11) and (9.12) is not possible due to the unphysical masses of the quarks in the lattice study and due to the omission of decays to J/ψω and η c η, which may affect the determination of the width. The qualitative comparison, however, suggests the existence of a D sDs resonance with small coupling to DD. This could be further investigated experimentally by considering the D sDs invariant mass spectrum near threshold, where a peak (see Fig. 7) would be visible for a state just below threshold. The X(3915) was proposed to be a groundccss state within the diquark-antidiquark approach by Polosa and Lebed [9]. The identificationccss was considered also in phenomenological studies [58,59].

Conclusions
We presented a lattice study of coupled-channel DD-D sDs scattering in the J P C = 0 ++ and 2 ++ quantum channels with isospin 0. Using the generalized Lüscher method and a piecewise parametrization of the energy region from slightly below 2m D to 4.13 GeV, the coupled-channel scattering matrix S along the real energy axis was determined. The resulting S was then analytically continued to search for pole singularities in the complex energy 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 CLS ensembles with u/d and s quarks, spatial extents ∼2.07 fm and ∼2.76 fm and a single lattice spacing ∼0.086 fm.
In addition to χ c0 (1P ), the results suggest three charmonium-like states with J P C = 0 ++ below 4.13 GeV. One is a yet undiscovered DD bound state just below threshold. The second is a narrow resonance just below the D sDs threshold predominantly coupled to D sDs . This state is possibly related to the narrow resonance X(3915)/χ c0 (3930), which is also below the D sDs threshold in the experiments. The third feature is a DD resonance possibly related to the χ c0 (3860) observed by Belle, which is believed to be χ c0 (2P ). An overview of the resulting pole structure of the coupled-channel DD-D sDs scattering matrix in the complex energy plane is given in Fig. 9, and the possible implications of this singularity structure for experiments are illustrated in Figs. 7 and 8. The masses are compared to experiment in Fig. 10 and summarized in Section 9.
Turning to states with J P C = 2 ++ , the mass of the ground state χ c2 (1P ) was determined directly from the lattice energy and is compared with the experimental value in Eq. (9.3). We have assumed the 2 ++ resonance to be coupled only with the DD scattering channel 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 the conventional χ c2 (3930) in Eqs. (9.7) and (9.8). These are then fixed for the finite-volume coupled-channel analysis discussed in Section 8.4.2. We find the estimates for positions and residues for the poles with J P C = 0 ++ to be robust with the exclusion/inclusion of the l = 2 partial wave contribution to the analysis. The resulting pole positions and the residues 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 ensembles we used have heavier-than-physical light and charm quarks and a lighter-than-physical strange quark. This results in a smaller-than-physical splitting between the DD and the D sDs thresholds. In future studies, it will be important to systematically improve upon the current results by successively relaxing our assumptions, for example, by explicitly including η c η interpolating fields and adding this channel as well as J/ψω to the coupled-channel study. It is also essential to investigate additional parametrizations to test the model inde-pendence of our findings and this will require a larger set of ensembles with high statistics. With regard to the pole structure observed in this work, it would be particularly interesting to investigate how our observations, such as the shallow DD bound state and our X(3860) candidate evolve when simultaneously approaching the limit of physical quark masses and the continuum limit. funding No. P1-0035 and No. J1-8137). We are grateful to the Mainz Institute for Theoretical Physics (MITP) for its hospitality and its partial support during the course of this work.

A Error treatment
Central values for all quantitiesQ are obtained from the average of correlation matrices over the gauge ensemble, while the errors are based on N b = 999 bootstrap samples. The 1σ standard error formulae for a Gaussian-distributed quantity Q provides the range that captures the central 68% of the bootstrap samples, which is represented by the gray bands in various figures. 13 We present resonance masses, widths/couplings, pole positions, phase shifts and |t| with the asymmetric errorsQ Here cov is the modified correlation matrix defined as The central valueQ corresponds to the average over the gauge configurations and not to the average of the bootstrap samples, therefore it is possible thatQ is not within the range captured by the 68% of the bootstrap samples.
The sum indicated by a prime runs over the bootstrap samples b in which Q b i and Q b j are not among the 16% of the values excluded on either end. cov ii in Eq. (A.1) is equal to the standard covariance 1 for Gaussian-distributed quantities. cov ij also coincides with the standard covariance for completely correlated or uncorrelated Gaussiandistributed quantities i and j. We have verified that the standard covariance and cov render almost identical errors on the energies and most of the parameters of the scattering matrix. The advantage of the modified correlation matrix is the exclusion of outliers for non-Gaussian distributions with long tails. Such distributions might occur for the bootstrap samples of scattering-matrix parametrizations due to the highly non-linear nature of the box-functions B(E cm ) in Eq. (7.1).

B Eigen-energies
The Figure 11 compares the original eigen-energies E lat cm and the eigen-energies E calc cm obtained via Eq. (3.2). The energies E calc are taken as inputs to the scattering matrix according to our approach towards disretization errors, as outlined in Section 3.

C Fitting the parameters of the scattering matrix
The parameters of the matrixK are determined from the energies presented in Fig. 2 via the quantization condition (7.1) following the determinant residual method proposed in 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))    Figure 11: The lower figures present the original eigen-energies E lat cm in four irreducible representations considered. The upper figures present the eigen-energies E calc cm obtained via Eq. (3.2) and match those in Fig. 2. The E calc are inputs to our scattering analysis, while E lat are considered to be less reliable input according to the discussion in Section 3. irreps shown in Fig. 2: these are levels n = 2(3) from | P | = 0(1) on both volumes. The charmonium-like state obtained lies just below threshold, therefore the relative error on its binding energy given in Eq. (8.5) is large. We note that 6% of the bootstrap sam- ples do not render any poles on the real axes -this corresponds to the bootstraps for which p cot δ/E cm just fails to cross the orange line in Fig. 3a. An additional 6% of the bootstrap samples render a virtual bound state -this corresponds to the bootstraps for which p cot δ/E cm crosses ip/E cm = |p|/E cm rather than −|p|/E cm in Fig. 3a. Both of these scenarios happen in extreme cases and end up within the 32% of bootstrap samples that are excluded when computing the errors via (A.1). In Fig. 12, we present the pole positions along the (virtual) bound state constraints across the bootstrap samples showing the continuous distribution of the poles along the constraint curves and hence also across the Riemann sheets. The preferred fit with the scattering parameters (8.2) utilizes the 4 levels shown in violet in Fig. 3a. We also performed the fits using the 3 lowest levels, the 6 lowest levels and all 7 levels shown: the ensemble averages of the data lead to a bound state in all these fits and the binding energy is within the error given in Eq. (8.5).
This analysis of DD scattering near threshold includes only the eigen-states with energies close to the threshold and omits the eigen-state related to χ c0 (1P ), which is significantly below threshold. We are unable to constrain the DD scattering below the lowest violet point in Fig. 3a. Hence, the pole at around E cm 3.80 GeV, which would arise at the crossing of the orange and red curves (and would violate the consistency check, see Section VC of [25]), is below the region in which our analysis can reasonably be applied and also outside of the energy range of interest.
D.2 D sDs scattering with l = 0 near threshold in the one-channel approximation This analysis employs only those eigenstates whose overlaps are dominated by D sDs operators and that do not have significant overlap with DD operators. These are the four levels in Fig. 2 near the D sDs threshold in the A (+) 1 irreps: levels n = 3, 4 from | P | 2 = 0, 1 on the N L = 24 ensemble and levels n = 4, 7 from | P | 2 = 1, 2 on N L = 32. Here 97% of the bootstrap samples result in a bound-state pole, while 2.3% result in a virtual bound-state pole and 0.7% do not render any poles on the real axis -the latter two cases end up among the extremal 32% of bootstrap samples.
D.3 DD scattering with l = 2 DD scattering in partial wave l = 2 is the not the main focus of our study. It was considered in order to investigate and constrain its contribution to the A 1 irreps we have studied. Since this partial wave was initially not the goal of our study, we did not evaluate all irreps where it appears (for example E + and T + 2 for P = 0), instead we implemented only the B 1 irrep with | P | 2 = 1. The extraction of the phase shift in Eqs. (8.10,8.11) employs four lattice levels in the B 1 irrep with | P | 2 = 1: these are levels n = 3, 4 on both lattice volumes (levels n = 1, 2 correspond to the ground states with J P C = 2 ++ and 2 −+ , respectively). A fit using only three lattice levels (omitting the higher level on the smaller volume) renders the resonance pole position E p = (4.013 +0.013 −0.016 ) − i 2 (0.098 +0.044 −0.057 ) GeV. This is compatible with our main result (8.11) and has a larger central value for the width.
D.4 Coupled DD, D sDs scattering with l = 0 for E cm 3.9 − 4.13 GeV Fig. 13 shows an example of Ω(E cm ) (C.1) for the parameters of the coupled-channel scattering matrix given in Eq. (8.14). The values of E cm at which Ω crosses zero are indeed near the observed eigen-energies (indicated by the black circles). The number of crossings agrees with the number of observed levels in the relevant energy ranges.
In Figure 14, we present the pole distribution across the bootstrap samples for various poles we extract in the complex p Ds plane, where p Ds is the momentum of the D s meson in D sDs scattering in the CMF. The two islands of poles, one in sheet II and the other in sheet IV, lying close to Im(ap Ds ) are mutually exclusive and hence represent the same dynamics. The island in sheet II constitutes 70% of the samples, while a pole appears at a similar location on sheet IV in the remaining samples. Hence the results for the pole location related to χ DsDs  P 2 =2, L/a=32 Figure 13: The function Ω(E cm ) (defined in Eq. (C.1)) for the resulting scattering matrix of the coupled channels DD − D sDs (8.14) is given by the orange line. The observed eigen-energies are given by the circles: the black levels are employed to fit the parameters (8.14), while the blue circles are not.
from 2m D up to 4.13 GeV requires additional parameters. This is difficult as, with the statistics and the number of lattice QCD ensembles available to us, the fits become unstable. Instead, as a cross-check, we model the infinite-volume scattering matrix in the wider energy range using the parametrizations similar to those presented in Section 8. One of the aims is to verify that the resulting scattering matrix predicts the same number of finite-volume energy levels as observed in the actual simulation. The t-matrix elements are modeled in the energy range E cm 2m D − 4.13 GeV as shown forK −1 / √ s in Fig. 15. We require that they are continuous in energy and that they have continuous derivative. In the high energy region they asymptote to the linear dependence on s (7.2,8.13) and the parameters are fixed to the values (8.14) obtained from the coupled-channel analysis. Below we provide more details on each t-matrix element in  t 11 The energy region considered is divided into three intervals, as shown by the red line in the figure. t 11 asymptotes to the coupled and single channel results (of the main text) in the high and middle energy intervals, respectively. In order to ensure a smooth transition between two regions, we emply hyperbola-type shape forK −1 11 / √ s that smoothly asymptotes 14 to the linear dependencesK −1 11 / √ s = a 11 + b 11 s (8.14) andK −1 11 / √ s = a 11 + b 11 s (8.2). The four parameters a, b are fixed to the values in the main text. The value of the smoothing parameter c hyp is the only free parameter ofK −1 ij in this appendix. Its value c hyp = 0.00021 (2) is obtained from fitting the scattering matrix to all energy levels and the resulting fit has χ 2 /dof = 1.8. In the region below DD threshold, we choose a shape ofK −1 11 which prevents the occurrence of a second bound-state (this would correspond to the red and orange lines in the figure intersecting a second time). The exact form of this choice is not important as this is beyond the region of interest.
t 12 In the upper regionK −1 12 / √ s asymptotes to the constant value of the coupled channel analysis (8.14). In the region near DD threshold, where the effects from D sDs channel are expected to be negligible, it asymptotes to zero. The smooth transition between two constant values is ensured by using the sigmoid function.
t 22 This element is parametrized as in Eq. (8.14), for the entire energy region, see the black line in the figure. Note that we ignore any crossing of the D sDs bound state condition with the t 22 parametrization that occurs well below the DD threshold.
We find that the number of poles with the above-designed scattering matrix is the same as that obtained from the analysis of the separate energy regions. The pole locations on the various complex Riemann sheets and their residues are also almost unchanged. Following Lüscher's finite-volume analysis, we extracted the finite-volume spectrum from this scattering matrix. In Fig. 16, we present the Ω(E cm ) function defined in Eq. (C.1). The zeros of Ω(E cm ) are the predictions for the finite-volume spectrum derived from the scattering matrix. The points indicate the energy levels observed in the actual simulation. It is clear from the figure that the predicted spectrum agrees qualitatively with the lattice energy levels within the energy range of interest and that the same number of levels are obtained. P 2 =2, L/a=32