Charmonium-like resonances with JPC = 0++, 2++ in coupled DD¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathrm{D}\overline{\mathrm{D}} $$\end{document}, DsD¯s\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\mathrm{D}}_{\mathrm{s}}{\overline{\mathrm{D}}}_{\mathrm{s}} $$\end{document} scattering on the lattice

We present the first lattice investigation of coupled-channel DD¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ D\overline{D} $$\end{document} and DsD¯s\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {D}_s{\overline{D}}_s $$\end{document} scattering in the JPC = 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üscher’s formalism. Lattice QCD ensembles from the CLS consortium with mπ ≃ 280 MeV, a ≃ 0.09 fm and L/a = 24, 32 are utilized. The resulting scattering matrix suggests the existence of three charmonium-like states with JPC = 0++ in the energy region ranging from slightly below 2mD up to 4.13 GeV. We find a so far unobserved DD¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ D\overline{D} $$\end{document} bound state just below threshold and a DD¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ D\overline{D} $$\end{document} resonance likely related to χc0(3860), which is believed to be χc0(2P). In addition, there is an indication for a narrow 0++ resonance just below the DsD¯s\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {D}_s{\overline{D}}_s $$\end{document} threshold with a large coupling to DsD¯s\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {D}_s{\overline{D}}_s $$\end{document} and a very small coupling to DD¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ D\overline{D} $$\end{document}. This resonance is possibly related to the narrow X(3915)/χc0(3930) observed in experiment also just below DsD¯s\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {D}_s{\overline{D}}_s $$\end{document}. The partial wave l = 2 features a resonance likely related to χc2(3930). We work with several assumptions, such as the omission of J/ψω, ηcη 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.

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 ), S = 1 + 2 i ρt = e 2iδ , (2.1) 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 , 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 JHEP06(2021)035 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 figure 1, along which the experimental measurements are made. Figure 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. denoted H105 and U101, have an inverse gauge coupling β = 6/g 2 = 3.4, corresponding to a lattice spacing a = 0.08636(98) (40) fm and volumes N T × N 3 L = 128 × 24 3 and 96 × 32 3 , respectively [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).

JHEP06(2021)035
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 section 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 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 in the moving lab frames also receive contributions from states with J P [λ] = 2 + [0] and JHEP06(2021)035 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 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.
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 twomeson operators involving spin 1 mesons, such as J/ψω and D * D * (see table 2). For nonzero 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 figure 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 section 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 -9 -JHEP06(2021)035 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 figure 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 figure 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 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 figure 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 ± .

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 finite- 5 Here, the energy in the lattice frame En stands for E calc n obtained according to eq. (3.2) or eq. (17) 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 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. volume 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 2 × 2 scattering matrix S(E cm ) for these coupled channels via the generalization of this formalism [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 functionK(E cm ) for one-channel scattering (2.1) and a real symmetric 2 × 2 matrixK(E cm ) for two coupled channels (2.3), as detailed 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]

JHEP06(2021)035
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 eq. (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 -12 -

JHEP06(2021)035
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.

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 figure 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  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 figure 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. Figure 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 -14 -

JHEP06(2021)035
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.

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 ( The resulting fit is shown in figure 4a. The scattering matrix has a bound state pole at the energy E cm = m where condition (8.3) is satisfied, see figure 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]). We find a shallow D sDs bound state at , 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 figure 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.

DD 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 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)  figure 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 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.

3.93-4.13 GeV
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 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 5 parameters of the scattering matrix (8.13)   -18 -

JHEP06(2021)035
The lattice energy levels and the levels predicted by this parametrization are compared in appendix D.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 and their locations (E p cm ) in the complex energy plane are shown in the left pane of figure 6. The couplings c i to both channels are extracted from the behavior of t ij near the poles using eq. (2.5) and are given in the same figure.

Analysis including l = 0 and l = 2
In the following we present our main result for the coupled channel scattering in the energy region 3.93 − 4.13 GeV. We fix the scattering amplitude for DD scattering in partial wave l = 2 to the values in eq.
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    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 figure 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 figure 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 figure 6 for the pole presented in red. The resonance parameters are 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. This resonance is related to the bound state in the analysis of D sDs -scattering in the one-channel 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 figure 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 figure 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.

Summary of the resulting hadrons and their relation to experiment
Below we summarize the properties of the charmonium-like hadrons found in this simulation. These are denoted by the conventional names χ c0 (1P ), χ c2 (1P ), χ c2 (3930) when 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 figure 9, while the corresponding masses are compared to experiment in figure 10. Due to the unphysical quark masses employed in the simulation, the hadron masses obtained are not compared to experiment directly. Instead  [2], where χ c0 (3930) [3] and X(3915) [2] may be the same state.

JHEP06(2021)035
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.

χ c0 (1P ) and χ c2 (1P )
These states lie significantly below the DD threshold and their masses are extracted from the ground state energies m = E 1 (P = 0) in irreps A ++ 1 and E ++ , with J P C = 0 ++ and 2 ++ , respectively. We obtain 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 figure 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 Experimental searches for this state in inclusive final states are challenging since the region above the DD threshold can be populated by DD from X(3872) → DD * → DDπ (see, for example, ref. [52]). Various strategies for the experimental search of such a state in exclusive decays were proposed: B 0+ → D 0D0 K 0+ [53], ψ(3770) → D 0D0 γ [54], γγ → DD [55], ψ(3770, 4040) → ηη γ and e + e − → ηη J/ψ [56].
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 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. The experimental couplings g of χ c0 (3860) and χ c2 (3930), quoted in eqs. (9.10) and (9.8), respectively, are similar, and their widths differ mainly because of the different phase space. On the lattice side, the χ c0 coupling in eq. (9.9) is smaller than the χ c2 (3930) coupling in eq. (9.7), however, both have large uncertainties.
Narrow 0 ++ resonance χ D sDs 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 figure 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 -26 -

JHEP06(2021)035
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 figure 9, and the possible implications of this singularity structure for experiments are illustrated in figures 7 and 8. The masses are compared to experiment in figure 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 figure 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 independence 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. 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 (A.1) 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. 13 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.   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 figure 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.

C Fitting the parameters of the scattering matrix
The parameters of the matrixK are determined from the energies presented in figure 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)) irreps shown in figure 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 samples 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 figure 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 figure 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 figure 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 figure 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 figure 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 figure 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)  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. 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. 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 c0 given in eq. (8.18) are based on the samples for which this pole appears on sheet II, whereas the errors on the rates and phase shifts presented in figures 7 and 8 are computed from the entire bootstrap samples.

E Coupled DD, D sDs scattering in a wider energy region
In section 8, we discussed the analysis of rather narrow energy ranges with a linear and/or constant parametrization in s for the elementsK −1 ij (s)/ √ s (see eq. (8.13)). A single description of coupled DD − D sDs scattering encompassing the whole of the energy range 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 n the actual simulation.   P 2 =2, L/a=32