Coupled-Channel $D\pi$, $D\eta$ and $D_{s}\bar{K}$ Scattering from Lattice QCD

We present the first lattice QCD study of coupled-channel $D\pi$, $D\eta$ and $D_{s}\bar{K}$ scattering in isospin-1/2 in three partial waves. Using distillation, we compute matrices of correlation functions with bases of operators capable of resolving both meson and meson-meson contributions to the spectrum. These correlation matrices are analysed using a variational approach to extract the finite-volume energy eigenstates. Utilising L\"uscher's method and its extensions, we constrain scattering amplitudes in $S$, $P$ and $D$-wave as a function of energy. By analytically continuing the scattering amplitudes to complex energies, we investigate the $S$-matrix singularities. Working at $m_\pi \approx 391$ MeV, we find a pole corresponding to a $J^{P} = 0^{+}$ near-threshold bound state with a large coupling to $D\pi$. We also find a deeply bound $J^{P} = 1^{-}$ state, and evidence for a $J^{P} = 2^{+}$ narrow resonance coupled predominantly to $D\pi$. Elastic $D\pi$ scattering in the isospin-$3/2$ channel is studied and we find a weakly repulsive interaction in $S$-wave.


Introduction
Over the last few years, a number of new states have been observed in both the charm-light (isospin-1/2, strangeness-0) D meson and the charm-strange (isospin-0, strangeness-1) D s meson systems and experiments continue to investigate their properties and find additional states [1,2]. The intermediate mass scale of the charm quark means that these systems provide a window on heavy-light dynamics away from the heavy-quark limit. The low-lying excitations are generally in agreement with expectations from quark-potential models [3] with some notable exceptions: the lightest scalar, D s0 (2317), and axial vector, D s1 (2460), charm-strange mesons were expected to be broad and above the relevant strong-decay threshold (DK and D K respectively), but they were both observed to be narrow and below threshold. A number of possible explanations have been put forward [2,4]. On the other hand, the corresponding charm-light mesons, D 0 (2400) and D 1 (2430), were both observed, as expected, to be broad resonances. The dynamics in the charm-light and charm-strange sectors are expected to be similar but the different masses that enter and the relative position of thresholds appear to be playing an important role. In this study, we investigate charm-light states as a step towards understanding these differences.
Lattice QCD provides a method for performing first-principles non-perturbative computations of the masses and other properties of hadrons within QCD. Correlation functions are computed numerically by Monte-Carlo sampling gauge configurations on a discretised Euclidean spacetime of finite volume, yielding a discrete spectrum of energy eigenstates. One virtue of lattice QCD is that it is systematically-improvable, permitting increasingly precise spectra to be obtained through efficient use of computational resources.
Recent lattice QCD investigations of charm-light mesons beyond the lightest pseudoscalar and vector include Refs. [34][35][36][37], but these calculations were not sensitive to meson-meson energy levels and so could not robustly determine states close to threshold or properly take into account the resonant nature of states above threshold. Elastic Dπ scattering was investigated to a limited extent in Ref. [38], and in Ref. [39] the isospin-3/2 Dπ scattering length was computed and used as an input to a chiral unitary approach to indirectly calculate the isospin-1/2 Dπ scattering length.
Here we present a lattice QCD investigation of isospin-1/2 coupled-channel Dπ, Dη, D sK scattering relevant for charm-light mesons, the first coupled-channel calculation using ab-initio methods in the charm sector: the Dπ channel opens first with Dη and D sK found (L/a s ) 3 × T /a t N cfgs N tsrcs N vecs 16 3 × 128  479  4  64  20 3 × 128  603  3  128  24 3 × 128  553  2-3  162   Table 1. The gauge field ensembles used in this study. The volume is given by, (L/a s ) 3 × (T /a t ), where L and T are respectively the spatial and temporal extents of the lattice. The number of gauge field configurations used, N cfgs , and the number of time-sources used per configuration, N tsrcs , are shown. N vecs refers to the number of eigenvectors used in the distillation framework.
close together a little higher in energy. We compute precise finite-volume spectra in many different symmetry channels for various momenta on multiple lattice volumes. From these spectra, we use extensions of the Lüscher method to determine infinite-volume scattering amplitudes. Considering coupled channels enables us to constrain the amplitudes over a larger range in energy than would be possible with elastic scattering; the extensive spectra allow us to determine these amplitudes robustly and to constrain the effect of higher partial waves. We also study elastic Dπ scattering in the exotic-flavour isospin-3/2 channel for which preliminary results have already appeared [40]. The remainder of this paper is laid out as follows: in Section 2 we give a brief description of the lattice ensembles used in this work, along with an overview of our methodology for extracting finite-volume spectra from two-point correlation functions. We then discuss how we obtain scattering amplitudes from finite-volume spectra. In Section 3 results for isospin-1/2 coupled-channel Dπ, Dη, D sK scattering are presented and in Section 4 we show our isospin-3/2 Dπ results. We summarise in Section 5.

Calculation Details
We use an anisotropic lattice formulation where the temporal lattice spacing, a t , is finer than the spatial lattice spacing, a s ≈ 0.12 fm, with ξ ≡ a s /a t ≈ 3.5. The finer temporal resolution is crucial in allowing us to accurately resolve the time dependence of two-point correlation functions enabling a precise determination of finite-volume energies. The gauge sector is described by a tree-level Symanzik-improved anisotropic action while in the fermionic sector a tadpole-improved anisotropic Sheikholeslami-Wohlert (clover) action, with stout-smeared gauge fields and N f = 2 + 1 flavours of dynamical quarks, is used. For these ensembles, m π ≈ 391 MeV, while the heavier dynamical quark is tuned to approximate the physical strange quark. The three different spatial volumes used are summarised in Table 1. Full details of the formulation are given in Refs. [41,42].
The same action is used for valence charm quarks as for the light and strange quarks (with tadpole-improved tree-level clover coefficients), where the charm quark mass parameter has previously been tuned using the physical η c mass [43]. By fitting to a relativistic dispersion relation, the anisotropy measured from the D-meson has been determined to be ξ D = 3.454(6) [34], which agrees with the value measured from the pion ξ π = 3.444(6) [44]. In this work, we use ξ π as the anisotropy and present our determination of scattering amplitudes incorporating its statistical uncertainty.
When we quote values in physical units, we set the scale by comparing the Ω-baryon mass determined on these ensembles, a t m Ω = 0.2951 [45], to the physical mass, m phys Ω , via atm Ω , leading to a −1 t = 5.667 GeV.

Finite-Volume Spectra
To determine the discrete spectrum of finite-volume energies we compute Euclidean twopoint correlation functions, where the interpolating operators, O † and O, are chosen to have the quantum numbers of the states of interest. In order to robustly extract many energy levels in each channel we follow our well established procedure [46]: a matrix of two-point correlation functions, C ij (t), is computed using a basis of operators with the relevant quantum numbers. A variational procedure [47] is employed, which amounts to solving a generalised eigenvalue problem, The energies then follow from analysing the time dependence of the eigenvalues (known as principal correlators), λ n (t, t 0 ). We fit each principal correlator to the form , where the fit parameters are E n , A n and E n ; the second exponential proves useful in stabilising the fit by "mopping up" excited state contamination. The eigenvectors, v n j , are related to the operator-state overlaps, Z n i ≡ n|O † i |0 , and also give weights for constructing a variationally-optimal operator to interpolate state n, Working in a finite cubic volume with periodic boundary conditions quantises the allowed momenta, P = 2π L (n x , n y , n z ), where (n x , n y , n z ) is a triplet of integers. We will use a shorthand notation when labelling momenta in which we omit the 2π L factor, e.g. P = [n x , n y , n z ] or [n x n y n z ]. The finite lattice volume also breaks the rotational symmetry of an infinite-volume continuum: for mesons at rest the relevant symmetry is that of a cube, the octahedral group with parity O h , whereas for mesons at non-zero momentum, P , the symmetry is reduced further to that of the little group, LG( P ) [48]. As a result, the continuum spin, J, is no longer a good quantum number and instead states must be labelled by the irreducible representations (irreps) of O h or LG( P ). The consequences of this for scattering will be discussed below in Section 2.2.
To reliably extract the many energy levels required to map out scattering amplitudes, we compute C ij (t) for large bases of interpolating operators with various structures. These include fermion-bilinearqq operators,ψΓ ← → D . . . ψ [46,49], as well as those resembling the combination of two-mesons, p 1 , [25,44], where Ω M i ( p i ) is a variationally-optimal linear combination of fermion-bilinear operators to interpolate meson M i with momentum p i , and C is a generalised Clebsch Gordan coefficient. For isospin-1/2 scattering we useqq operators along with Dπ, Dη and D sK "two-meson" operators. For  Figure 1. Principal correlators, labelled by the extracted energy, determined on the 20 3 volume in the [000]A + 1 irrep with isospin-1/2. Points show λ n (t, t 0 = 13) and error bars correspond to the one-sigma statistical uncertainty. In each plot the dominant time-dependence, e −En(t−t0) , has been divided out so that a horizontal line is observed when a single exponential dominates. Curves show fits to the form described in the text; the red curves show the fit range and blue points are not included in the fit.
isospin-3/2 we only use Dπ operators -there are noqq operators with this isospin. The operator bases we use are listed in Tables 5, 6, 7, 8 and 9 in Appendix B.
We use the distillation framework [50] which enables us to efficiently compute correlation functions involving operators with various structures where each operator is projected onto a definite momentum. In Table 1 we give the number of distillation vectors, N vecs , used on each lattice volume, along with the number of time-sources used per configuration, N tsrcs .
As an example of the quality of the signals extracted, in Fig. 1 we show the principal correlators from the [ P ]Λ P = [000]A + 1 irrep with isospin-1/2 on the 20 3 volume; the resulting spectrum is shown in the leftmost panel of Fig. 2. In each plot, the dominant time-dependence, e −En(t−t 0 ) , has been divided out and we observe a horizontal line when a single exponential dominates. The principal correlators shown here are representative of all those determined within our calculation.

Scattering Amplitudes from Finite-Volume Spectra
Having determined the finite-volume spectra, we relate these to infinite-volume scattering amplitudes using the Lüscher method [5,6] and its extensions to moving frames [7,10,13,51] and coupled-channels [16][17][18]52]. In this approach, the dependence of the spectra on finite volume is used as a tool but exponentially-suppressed corrections in the volume are neglected -typically the leading such corrections fall off as e −mπL and, since our volumes have m π L ∼ 4 to 6, we can safely neglect these.
For lattice irrep Λ and overall momentum P = 2π L d, the relation between the finitevolume spectra and the infinite-volume scattering t-matrix is given by where i and j label the scattering channels, ρ i = 2k i /E cm is the phase-space factor for channel i, k i denotes the the momentum in the centre of momentum frame, and t ( ) ij is the infinite-volume t-matrix for partial wave . M(q 2 i ) is a matrix of known functions of q 2 i = (k i L/2π) 2 where and are partial waves that can subduce into irrep Λ and the index n indicates the n'th subduction of partial wave (similarly for n and ) -we show the pattern of subductions in Table 4 in Appendix A. The mixing between different partial waves, encoded in M, is a consequence of the reduced symmetry of the finite cubic volume. From the table it is clear that even and odd partial waves mix when the overall momentum is non-zero (this is a consequence of the unequal mass of the scattering mesons).
Given an infinite-volume t-matrix, Eq. 2.3 can be solved to find the finite-volume spectrum, {E cm }. However, the reverse is in general an under-constrained problem: for N coupled channels there is one equation for each energy level but (N 2 + N )/2 energydependent parameters in the t-matrix. Hence, even neglecting the mixing between different partial waves, there is insufficient information to determine the t-matrix directly. In order to circumvent this difficulty, we follow Ref. [32] and parametrise the energy-dependence of the t-matrix with a relatively small number of parameters. Using Eq. 2.3, this parametrised tmatrix gives a spectrum {E par cm } and we vary the parameters to fit {E par cm } to our computed spectrum {E cm }, minimising the χ 2 function defined in Eq. 8 in [32]. By analytically continuing the resulting t-matrix into the complex s = E 2 cm plane, we can determine the pole and residue content of the scattering amplitude, which are arguably the least methoddependent quantities that can be compared between analyses. We consider a wide range of different parametrisations to ensure the final answer does not depend on a particular form used.
When considering elastic scattering, the t-matrix can be described by a single energydependent parameter, the scattering phase shift, δ (E cm ), where t ( ) = 1 ρ e iδ sin δ . Two commonly used parametrisations in this case are the effective range expansion and the relativistic Breit-Wigner. The first is given by where the constants a and r are known as the scattering length and the effective range respectively. The second, which is commonly used to describe a resonance, is given by where m R is known as the Breit-Wigner mass. Γ (s) is the energy-dependent width which can be parametrised in terms of the coupling g R , Γ (s) = For coupled-channel scattering, the relationship between the t-matrix, phases and inelasticities becomes more complicated. For each channel i, the phases, δ i , can be extracted from the diagonal elements of the t-matrix for each partial wave , which also provides a convention for determining the inelasticities η i . When parametrising the t-matrix in the coupled-channel case, we make use of the K-matrix formalism, where the inverse of the t-matrix for a partial wave is given by The factors (2k i ) − ensure correct behaviour in the proximity of kinematic thresholds [53], while we parametrise the symmetric matrix K −1 ij (s). There is of course some freedom in this parametrisation and in this work we will make use a variety of forms which can be written generally as ij are real free parameters. Unitarity of the t-matrix is ensured when K −1 ij (s) is real for real values of s, Im[I ij (s)] = −δ ij ρ i (s) for energies above the kinematic threshold of channel i and Im[I ij (s)] = 0 for energies below that same threshold. Since unitarity does not directly constrain Re[I ij (s)], there is some freedom in its choice, with the simplest option being being Re[I ij (s)] = 0, i.e. I ij (s) = −iρ ij (s). A different choice is the Chew-Mandelstam prescription [54], which uses the known imaginary part of I ij (s) to determine the real contribution through a dispersion relation. In this scheme, which captures many of the correct analytic properties of scattering amplitudes, the dispersion integral is made finite by subtraction at an arbitrary point. For elastic S-wave Dπ scattering we subtract at threshold and in all other cases we subtract at the K-matrix pole parameter (s = m 2 ). 1 Details of our implementation are given in Ref. [32]. In this work, we will only consider energies far from the left-hand cut (which arises from t-channel exchanges). As a consequence, we do not consider the effects of this cut.

Results: Isospin-1/2
We now present the results of our calculations in the isospin-1/2 channel. First we discuss the finite-volume spectra obtained from the variational procedure, before moving on to discuss both elastic Dπ and coupled-channel Dπ, Dη and D sK scattering. We end the section with an interpretation of our results in terms of poles of the infinite-volume scattering matrix.

Finite-Volume Spectra
Following the procedure described in Section 2.1, the large bases of interpolating operators listed in Tables 5, 6 and 7 in Appendix B are used to determine finite-volume spectra in a meson a t m π 0.06906(13) Table 2. Relevant stable meson masses and kinematic thresholds on our ensembles [33,34,43,44]. Those shown in italics do not contribute to pseudoscalar-pseudoscalar scattering in S-wave.
number of lattice irreps, [ P ]Λ (P ) , where parity P is only a good quantum number when the system has zero overall momentum. These are shown in Figs. 2, 3 and 4, where the black and grey points correspond to the extracted finite-volume energy levels; only black points are used in the subsequent scattering analyses. The red, green and blue curves (dashing) are the non-interacting Dπ, Dη and D sK energies (thresholds) respectively. The grey dotted lines show the opening of channels for which we have not included interpolating operators in our variational procedure. Relevant meson masses and multi-meson thresholds are given in Table 2.
The leftmost panel of Fig. 2 shows the spectrum obtained in the [000]A + 1 irrep, whose lowest partial wave contribution comes from = 0. The lowest energy level, which has a large overlap with our D 000 π 000 operator (defined in Table 5), shows a volume dependent shift away from the Dπ threshold. Furthermore, there appears to be an "extra" energy level compared to the number of non-interacting multi-meson levels expected in the energy region below the Dη threshold. These features may point to a non-trivial meson-meson interaction in S-wave within the energy region around Dπ threshold. In Fig. 3, we show spectra extracted in irreps whose lowest partial wave contribution comes from = 1. In this case we observe a level far below the Dπ threshold at a t E cm ≈ 0.35, suggesting that a vector bound state is present. In Fig. 2, we also show the [ P = 0]A 1 irreps where both S and P -waves can contribute. Here we find further evidence for a non-trivial S-wave interaction near the Dπ threshold and the deeply bound vector state; we observe volume dependent shifts near the Dπ threshold along with the appearance of an "extra" energy level, while also finding an energy level far below the Dπ threshold at a t E cm ≈ 0.35.
The spectra shown in Fig. 4 have = 2 as their lowest contributing partial wave. Within the energy range 0.44 a t E cm 0.46, we observe significant shifts of the energy levels away from non-interacting energies along with the presence of an "extra" energy level, indicative of non-trivial interactions in D-wave.  Table 5 (with error bars showing the statistical uncertainty); those coloured black are used in our subsequent scattering analysis whereas those coloured grey are not. Solid curves represent noninteracting energies while dashed lines correspond to channel thresholds. The colour coding is as follows: red corresponds to Dπ, green to Dη and blue to D sK . The grey dashed and dotted lines show the opening of channels for which we have not included operators in our variational procedure.

Elastic Dπ Scattering
To begin our scattering analysis we consider the region where only the Dπ channel is kinematically open, that is below the Dη threshold when the system has zero overall momentum and below the D π threshold when the system has overall non-zero momentum. In the near-threshold region, higher partial waves are suppressed in proportion to k 2 +1 in the absence of resonances. This is important because partial waves can mix; for Dπ scattering in a finite volume, the allowed partial waves are P = 0 + , 1 − , 2 + , ..., where parity, P , is only a good quantum number when the system has overall zero momentum.

S-wave
By following the procedure described in Section 2.2, we determine the S-wave scattering amplitude using only the [000]A + 1 irrep, neglecting ≥ 4 contributions. Taking the lowest two levels in each volume and parametrising the t-matrix using a K-matrix containing a pole term and a constant, the χ 2 function is minimised, obtaining   Table 7.
For each parameter, our convention is that the first uncertainty reflects the χ 2 minimisation and the second uncertainty is obtained by varying the scattered meson masses and the anisotropy within their statistical uncertainties. The matrix shows the correlation between each parameter. Although this form does not demand the existence of a nearby pole in the t-matrix, it permits one to arise with parametric efficiency and the well-determined K-matrix pole parameter might suggest the presence of a t-matrix pole. Furthermore, parametrisations that do not permit poles were unable to successfully reproduce our finitevolume spectra. In Section 3.2.4 we explore further forms of the t-matrix and interpret our results in terms of poles in Section 3.5.

P -wave
In each of the irreps with an = 1 contribution, we find an energy level at a t E cm ≈ 0.35 well below Dπ threshold. The first source of possible inelasticity in = 1 comes from D π contributions, so we exclude extracted levels above D π threshold. We consider the ten energy levels below the D π (and also below the Dπ) threshold from irreps that have = 1 as their lowest contributing partial wave, these are the black points in Fig.  3. Assuming that contributions coming from ≥ 2 are negligible, the t-matrix can be parametrised yielding a P -wave scattering amplitude around a t E cm ≈ 0.35. Using a K-matrix description that includes only a pole term, we obtain 23 10−2 = 1.40.
All of the energy levels found at a t E cm ≈ 0.35 have a large overlap withqq operators subduced from those with J = 1 in the infinite-volume continuum. Along with the welldetermined pole parameter in Eq. (3.2), this may indicate the presence of a deeply bound vector state consistent with what was previously found in Ref. [34] (which did not include multi-meson operators). We find no other levels in the elastic scattering region meaning that we can not constrain the P -wave scattering amplitude further without including irreps that have both S and P -wave contributions.

S and P -waves
We now determine S and P -wave scattering amplitudes simultaneously.  Separately for each of the = 0 and = 1 parts of the t-matrix, we use a K-matrix containing a pole term and a constant. This leads to where the parameters with a subscript 1 denote the P -wave. In Fig. 5, we show the phase shifts, δ Dπ 0 and δ Dπ 1 , along with the finite-volume energy levels (black points) and those given by the parametrisation in Eq. (3.3) (orange). Although the points clustered around a t E cm ≈ 0.35 are not shown, they are included in our fit. Just above the Dπ threshold the S-wave phase shift shows a rapid variation, a feature indicative of a nearby pole. On the other hand, the P -wave phase shift varies slowly throughout the energy range shown. In the upper (lower) panel of Fig. 6, we show the quantity k Dπ cot δ Dπ 0 (k 3 Dπ cot δ Dπ 1 ) determined from the parametrisation in Eq. (3.3). The dotted curves correspond to the quantity ik 2 +1 Dπ , which should intersect the bands at the location of a subthreshold t-matrix pole on the physical sheet. The intersection with k Dπ cot δ Dπ 0 provides evidence that the possible pole near the Dπ threshold actually lies just below it, while the intersection with k 3 Dπ cot δ Dπ 1 suggests a bound state around a t E cm ≈ 0.35. We defer further discussion to Section 3.5 where we investigate the singularity content of these scattering amplitudes.

Parametrisation Variation
To assess the extent to which the scattering amplitudes depend upon the choice of parametrisation, we repeat the procedure above for a variety of forms of both the S and P -wave parts of the t-matrix. A selection of these forms is summarised in Table 10 in Appendix C, where we also show the χ 2 /N dof obtained when minimising our χ 2 function using the same 33 energy levels as above. Several other parametrisations have also been tested, such as higher orders of the forms shown or the inverse polynomial K-matrix used in Ref. [32]; these are found to be highly correlated or contained parameters consistent with zero so we do not show them.
In Fig. 7, we show the quantity ρ 2 Dπ |t DπDπ | 2 (left), which is proportional to the Dπ → Dπ cross-section, and the phase shift δ Dπ (right). The size of the bands encompasses the variation and uncertainty coming from all parametrisations with χ 2 /N dof < 1.9. From a  Figure 7. The quantity ρ 2 Dπ |t DπDπ | 2 (left) and the phase shift δ Dπ (right). The red and orange bands correspond to S and P -wave respectively -they encompass the various parametrisations with χ 2 /N dof < 1.9 in Table 10 as well as the uncertainty from the anisotropy and the scattered meson masses.
comparison of the phase shifts with Fig. 5, it is clear that the parametrisation dependence is almost negligible. In the quantity ρ 2 Dπ |t DπDπ | 2 , we observe further evidence of an S-wave pole near the Dπ threshold, as the large "peak" almost saturates the unitarity bound, which is unity in our normalisation.

Coupled-Channel Scattering in S-wave
We now go beyond the the elastic Dπ energy region and consider the case of coupledchannel Dπ, Dη and D sK scattering in S-wave in combination with elastic Dπ scattering in P -wave; the P -wave constraint comes only from energy levels below the coupled-channel region. The energy levels in the coupled-channel region are all from [000]A + 1 where the first contamination comes from = 4 and is expected to be highly suppressed. The irreps with P = 0 can have contributions from = 2 but these are later shown to be negligible in the elastic region. Using the 47 energy levels coloured black in Figs. 2 and 3, we parametrise the t-matrix by a coupled-channel K-matrix in S-wave and an elastic K-matrix in P -wave, giving  Note that the P -wave parameters are consistent with what is obtained in Eq. (3.3); as one might expect, the inclusion of the coupled-channel region in S-wave appears to have a negligible effect on the P -wave amplitudes. In Fig. 8 we show a comparison between our finite-volume spectrum in the [000]A + 1 irrep (black points) and the spectrum coming from the parameters in Eq. (3.4) (orange points).
In the upper left panel of Fig. 9, we show the S-wave phase shifts δ Dπ 0 (red), δ Dη 0 (green), and δ DsK 0 (blue) corresponding to the parametrisation in Eq (3.4). By comparing to the elastic case in the right panel of Fig. 7, we see that our results for the elastic Dπ region are largely unaffected when we allow for the Dη and the D sK channels. However, at the opening of the Dη threshold we do observe a noticeable "kink" in the Dπ phase shift suggesting a non-zero coupling between the two channels. We see a much smaller  Table 11 with χ 2 /N dof < 1.9. The black points between the upper and lower panels show the location of the finite-volume energy levels used to constrain the parametrisations.
effect at the opening of the D sK threshold. The non-zero coupling between channels is further demonstrated in the lower left panel of Fig. 9, which shows a clear deviation of the inelasticities from unity. The upper (lower) left panel of Fig. 10 shows ρ i ρ j |t ij | 2 for i = j (i = j) determined from the parametrisation in Eq (3.4); this quantity is proportional to the cross section for scattering of channel i → i (i → j). We see that just above Dπ threshold, as in the elastic case, the unitarity bound is almost saturated for Dπ → Dπ.

Parametrisation Variation
We now assess the extent to which our results depend upon our choice of parametrisation of the t-matrix. Table 11 in Appendix B shows a selection of parametrisations of the tmatrix we considered with the χ 2 /N dof obtained in each case. Note that, we have also attempted several other parametrisations, such as ij s k and those with higher order terms of the forms shown. However, these were found to either have insufficient freedom to describe our finite-volume spectra or to give highly correlated parameters.
In the upper (lower) right panel of Fig. 9, we show the Dπ, Dη and D sK S-wave phase shifts (inelasticities) where the size of the bands include the one-sigma statistical uncertainty coming from all parametrisations with a χ 2 /N dof < 1.9 as well as the statistical uncertainty coming from the scattered meson masses and the anisotropy. It appears that there is almost no parametrisation dependence up to around a t E cm ≈ 0. 46 Table 11 with χ 2 /N dof < 1.9. The black points show the location of the finite-volume energy levels used to constrain the parametrisations. many energy levels to constrain the scattering amplitude and hence we see a dramatic reduction on the constraint we can place on the phases 2 .
In the right panels of Fig. 10, we show the quantity ρ i ρ j |t ij | 2 , where the size of the bands include the parametrisations from Table 11 with χ 2 /N dof < 1.9. There appears to be very little parametrisation dependence in this quantity with all features, including the large "peak" in the Dπ → Dπ channel just above the Dπ threshold, remaining intact.

Coupled-Channel Scattering in D-wave
We now turn our attention to coupled-channel scattering in D-wave. The spectra shown in Fig. 4, namely [000]E + , [000]T + 2 , [001]B 1 and [001]B 2 , have = 2 as the lowest contributing partial wave. Although we have included operators for many of the open channels within the energy region shown, there are some for which we have not included operators, notably the relatively low-lying D π channel. Nevertheless, D-wave channels are known to open slowly since they are suppressed in proportion to k 2 +1 , suggesting that it may still be possible to apply the Lüscher formalism; Ref. [55] showed that, at least in S-wave, the breakdown of the formalism above an inelastic threshold results in values of the phase shift clearly incompatible with those below the threshold. We have thoroughly checked for the presence of such effects but do not find any evidence for them. Encouraged by this, we cautiously proceed with our scattering analysis.
Utilising the 28 energy levels up to around the D ππ threshold in Fig. 4, we parametrise the D-wave part of the t-matrix using a three-channel K-matrix form. We find a reasonable  (3.5) In Fig. 11, we compare the finite-volume spectra in the [000]E + and [000]T + 2 irreps to those determined from the parametrisation in Eq. (3.5).
As before, we assess the extent to which our results depend upon a given parametrisation; Table 12 in Appendix B shows a selection of parametrisations of the K-matrix used to determine the D-wave scattering amplitude. The upper left panel of Fig. 12 shows the phase shift, δ Dπ 2 , where the size of the band incorporates all of the parametrisations shown in Table 12. We observe that the phase shifts are small and well determined in the elastic region, justifying our neglect of D-wave contributions when determining the S and P -wave amplitudes above. The lower left panel shows the corresponding inelasticities, where a clear decoupling of the channels is observed. This enables a one-to-one correspondence between the finite-volume energy levels and the phases. In the right panel of Fig. 12, we show the phase shift points determined for each energy level superimposed onto the phase shifts from the upper left panel 3 . The agreement between the two approaches further indicates the lack of coupling of the Dη and D sK channels to the resonance.  Table 12. The right panel shows the same phase-shift bands superimposed with points determined for each energy level as described in the text.
The narrowness of the phase shift and the apparent decoupling of the channels suggest that a Breit-Wigner parametrisation may also be capable of describing the resonance. By selecting only those levels identified as belonging to the Dπ channel, we obtain a reasonable description with the parameters,

Poles and Interpretation
Finite-volume energy levels determined from Euclidean two-point correlation functions are real and in the sections above we have used these energies to constrain the t-matrix at real values of energy. However, the t-matrix can also be considered a function of complex energies, where bound states and resonances can be associated with poles in the complex plane. In the proximity of a pole, s pole , the t-matrix is dominated by the term where the factorised residues, c i , are complex numbers that give a measure of the "coupling" of the pole to channel i.
In terms of complex energies, branch cuts appear in the t-matrix for each scattering threshold leading to 2 N Riemann sheets for N coupled-channels. Sheets can be labelled by the sign of the imaginary part of the momenta, k i , for each channel i. Poles that correspond to a resonance occur in complex conjugate pairs on "unphysical sheets", where Im[k i ] < 0 in at least one channel. The only poles permitted to occur on the "physical sheet", where all Im[k i ] > 0, are those corresponding to bound states. Bound states far below threshold are unlikely to influence physical scattering, but one sufficiently close to threshold can cause noticeable effects. We now proceed to interpret our results in terms of poles we find in our parametrised t-matrices.

S-wave
To begin, we investigate the pole structure of the S-wave parametrisations used to describe elastic Dπ scattering in Section 3.2.4. In all of the parametrisations we considered, we consistently find a bound-state pole on the real axis of the physical sheet extremely close to the Dπ threshold. In the central panel of Fig. 13, we show the location of this pole for each parametrisation along with the χ 2 /N dof in each case. By averaging over the pole positions from parametrisations with χ 2 dof < 1.9, we find a t √ s pole = 0.40155 ± 0.00015 , (3.8) where the quoted uncertainty encompasses the uncertainties from the individual parametrisations. Although the central value lies below the Dπ threshold, which is located at a t E cm = 0.40171 ± 0.00015, they overlap within uncertainties. The effect of the pole is seen in our S-wave amplitude as shown in Fig. 7, where we observe a rapid variation coincident with the Dπ threshold. We also extract the residue of the pole, measuring the strength of its coupling to the Dπ channel; this is shown in the right panel of Fig. 13. Since the central value of the pole is below threshold, the residue has no imaginary part. We find that, like the pole, the residue is very stable across parametrisations, and averaging over all parametrisations shown we obtain a t c Dπ = 0.110 ± 0.025 . (3.9) We do not find any further poles in the region where we have constrained the amplitudes. In all of the coupled-channel parametrisations used in Section 3.2.4, we find a pole consistent with that of the elastic case described above. Averaging over the coupledchannel parametrisations with a χ 2 /N dof < 1.9, we find our final location of the pole to be a t √ s pole = 0.40161 ± 0.00015 .

P -wave
In the P -wave part of the t-matrices determined in Section 3.3.1, we consistently find a bound-state pole on the real axis of the physical sheet well below threshold. Averaging over parametrisations with χ 2 /N dof < 1.9, we obtain our final position for the pole to be a t √ s pole = 0.35440 ± 0.00023 . (3.12) As discussed in Section 3.2.2, this pole can be associated with the stable J = 1 − state found at almost exactly the same energy in Ref. [34]. Because this state is far below the Dπ threshold it is not expected to strongly influence Dπ scattering; this is consistent with the small P -wave amplitudes we find.

D-wave
In Section 3.4, we determined D-wave scattering amplitudes by considering the Dπ, Dη and D sK channels 4 . In all of the parametrisations of the t-matrix we considered, we found what appeared to be an "extra level" in close proximity to both the Dη and D sK thresholds along with a rapid phase shift through 180 • in the Dπ channel. We find an isolated resonance pole on each of the Riemann sheets with Im[k Dπ ] < 0. As shown in Fig. 14 Table 12.  Figure 15. The coupling, c i , of the Dπ (red), Dη (green) and D sK (blue) channels to the D-wave resonance pole shown in Fig. 14. The inset shows the coupling to the Dπ channel for the various parametrisations listed in Table 12 with the black point corresponding to the parametrisation in Eq. (3.5). We also determine the couplings of each channel to the pole. As shown in Fig. 15, we find a non-zero coupling only in the Dπ channel; as expected from Section 3.4, the Dη and D sK channels are decoupled from the resonance. Averaging over all of our parametrisations, we obtain our final value for the coupling a t c Dπ = (0.0431 ± 0.0015) · exp iπ(−0.0106 ± 0.0013) . (3.14) 4 Results: Isospin-3/2 We now change focus and present the results of elastic Dπ scattering in the isospin-3/2 channel. In our calculation we include only Dπ interpolating operators: there are noqq operators with this isospin.

Finite-Volume Spectra
As before, we determine energy levels from a variational procedure applied to a matrix of two-point correlation functions. We construct these correlation functions using the Dπ operators listed in Tables 8 and 9 in Appendix B. In Fig. 16  In all four spectra, we observe small positive shifts of our energy levels from the non-interacting energies, which is usually indicative of a weakly repulsive interaction.

Elastic Dπ Scattering
We begin by using the energy levels coloured black in Fig. 16 to map out the S-wave phase shift δ 0 as a function of energy. To do so, we ignore higher partial waves and solve Eq. (2.3) for each considered energy level, resulting in a corresponding value for δ 0 at that energy. We show the outcome of this procedure in Fig. 18 irreps. The form of the phase shift is consistent with that of a weakly repulsive interaction, which may be expected in the isospin-3/2 channel. As mentioned above, we have neglected contributions of higher partial waves; in determining the black points we have neglected contributions coming from ≥ 4, which is justified due to the large angular momentum suppression, and in determining the grey points we have neglected contributions coming from ≥ 1. To justify the use of the grey points, we assume negligible inelasticity into D π and compute the magnitude of the P -wave phase shift, δ 1 , within the energy range 0.43 ≤ a t E cm ≤ 0.45 using the energy levels in Fig. 17. We find that |δ 1 | ≤ 5 • at a t E cm ≈ 0.45 and |δ 1 | ≤ 3 • at a t E cm ≈ 0.43. We expect this trend to continue down to a t E cm ≈ 0.41, justifying our earlier assumption.
We can also follow the approach taken in the isospin-1/2 section and use the determined spectra to constrain the scattering amplitude as a function of energy. By again ignoring the negligible contribution coming from ≥ 4, we use the [000]A + 1 energy levels in Fig. 16 and parametrise the t-matrix using a scattering length and effective range. We find that the parameters are sufficient to describe our determined spectrum. In Fig. 19 we show the comparison between our determined spectrum and the spectrum resulting from the parameters in Eq. (4.1). We repeat this process for various parametrisations, and in Table 3 show a selection that were found to obtain a reasonable χ 2 /N dof . Our results suggest that the best description of our finite-volume spectrum is achieved using two free parameters. Fig. 20 shows the S-wave phase shift, δ 0 , where the size of the band includes parametrisations which were found to result in a χ 2 /N dof < 1.9. The points in Fig. 20 are taken from Fig. 18 and superimposed to demonstrate the consistency between this approach and the approach discussed above.
To get a handle on the P -wave amplitudes we can include the energy levels coloured black from the [ P = 0]A 1 irreps in Fig. 16. However, we find that for all our forms of the t-matrix, the P -wave parameters are always consistent with zero. Furthermore, by setting the P -wave parameters to zero, we we see no significant variation in any of our S-wave parameters when we include energy levels coming from the [ P = 0]A 1 irreps. This agrees Effective range expansion k cot δ = 1 a 1 3.66 k cot δ = 1 a + 1 2 r 2 k 2 2 1.29 Table 3. A selection of S-wave parametrisations for elastic isospin-3/2 Dπ scattering. Only those with χ 2 /N dof < 1.9 are used in Fig. 20. Figure 20. The S-wave phase shift, δ 0 . The extent of the band encompasses the phase shifts resulting from the parametrisations with χ 2 /N dof < 1.9 in Table 3 as well as the systematic uncertainty coming from the anisotropy and the scattered meson masses. The points from Fig. 18  with what we found using the approach discussed above; the contribution coming from ≥ 1 below the D π threshold is negligible.
In physical units our result for the scattering length is a 0 = −0.19 ± 0.05 fm. This value is in agreement with the results of Ref. [39] where a chiral unitary approach is used to interpolate between different values of the pion mass.

Summary and Outlook
We have presented the first lattice QCD study of coupled-channel Dπ, Dη and D sK scattering. Utilising the distillation framework and large bases ofqq, Dπ, Dη and D sK -like interpolating operators, we have determined finite-volume energy levels in many lattice symmetry channels. These were used to constrain scattering amplitudes as a function of energy which we analytically continue to complex energies. We determined the pole structure of these amplitudes and, from their residues, the coupling of poles to channels. Figure 21 shows the resulting S, P and D-wave isospin-1/2 poles.
In the J P = 0 + channel, we find a pole just below the Dπ threshold signalling the presence of a bound state. Although this pole appears to couple predominantly to Dπ, we obtain significant couplings to the Dη and D sK channels. As shown in Fig. 21, we find the pole to be at (2275.9 ± 0.9) MeV, which is statistically indistinguishable from the Dπ threshold which is at (2276.4 ± 0.9) MeV on our ensembles. As a consequence, the effect of the pole can be seen above threshold; as shown in Fig. 22, we observe a large "peak" in ρ 2 Dπ |t Dπ,Dπ | 2 almost saturating the unitarity bound. Since m π = 391 MeV in this calculation, we only make a qualitative comparison with experiment. Although we find a near-threshold bound state, it shares similarities with the experimental D 0 (2400) resonance; both states couple dominantly to Dπ and influence a similar broad energy range [1]. Noting that the light quark mass in this calculation lies between the physical light and strange quark masses, and ignoring the differences due to flavour, the relative position of the pole to the threshold lies between what is observed experimentally for the resonant D 0 (2400) and the bound D s0 (2317). Qualitatively, such behaviour is anticipated from unitarised chiral perturbation theory amplitudes [39,56,57].
As shown in Fig. 21, we find a pole at (2009 ± 2) MeV in the 1 − channel corresponding to a deeply-bound state, consistent with what was found in Ref. [34]. Experimentally, the near-threshold D (2007) resonance, which is narrow and decays predominantly to Dπ, has a mass of (2006.97 ± 0.08) MeV [1].
In the 2 + channel we find a narrow resonance coupled to Dπ. As shown in Fig. 21, we determine its pole mass and width to be (2527 ± 3) MeV and (8.2 ± 0.7) MeV respectively. Experiment finds a relatively narrow tensor resonance, the D 2 (2460), coupled to Dπ. However, this also couples to D π, a kinematically open channel which we have neglected in the determination of our 2 + state.
We have also performed a study of elastic Dπ scattering in the isospin-3/2 channel. We find that the weakly-repulsive S-wave interaction can be successfully described using a scattering length and effective range parametrisation with a 0 = −0.19 ± 0.05 fm and r 0 = −0.9 ± 0.4 fm.
This work, which is the first ab initio coupled-channel study including charm quarks, has taken a significant step towards understanding the striking differences between the D 0 (2400) and D s0 (2317). A complementary study of DK scattering is already underway and in the near future we will perform calculations with lighter pion masses -these will enable a more direct comparison with experiment and allow us to study the dependence of the pole position on the light-quark mass. These calculations will also include additional  Figure 22. ρ i ρ j |t ij | 2 for S-wave scattering in the isospin-1/2 channel. The bands encompass all the parametrisations with χ 2 /N dof < 1.9 in Table 11 along with the uncertainties coming from the scattered meson masses and the anisotropy. Black points show the location of the finite-volume energy levels used to constrain the scattering amplitudes. channels, such as D π, whose role may become increasingly important as the light quark mass is reduced.

Appendices A Lattice Irreps and Partial Waves
In Table 4 we show the pattern of contributing partial waves, , for various overall momentum types, P , and lattice irreps, Λ, relevant for the scattering of two unequal-mass pseudoscalars. Note that, because the mesons have different masses, even and odd partial waves mix for non-zero overall momentum.

B Operator Lists
In Tables 5, 6 and 7 we list the interpolating operators used to determine the finite-volume energy levels in the isospin-1/2 channel shown in Figs. 2, 3 and 4 respectively. Tables 8  and 9 show the operators used in the isospin-3/2 channel to determine the finite-volume spectra shown in Figs. 16 and 17 respectively; note that there are noqq operators with this isospin.

C Parametrisation Variations
In Table 10 we show the elastic S and P -wave parametrisations used in Section 3.2.4. In Tables 11 and 12 we show the coupled-channel S and D-wave t-matrix parametrisations used in Sections 3.3.1 and 3.4 respectively.
LG( P ) is the double-cover little group and the corresponding single-cover little group relevant for only integer spin is given in parentheses. Also shown are the various J ≤ 4 or |λ| ≤ 4 that appear in each of the relevant irreps. The J P values and |λ|η = 0 − in italics are in the "unnatural parity" [P = (−1) J+1 ] series and do not contribute to pseudoscalar-pseudoscalar scattering.

Parameters
ReI(s) Npars χ 2 /N dof m g  Table 12. As Table 11 but for the D-wave parametrisations used in Section 3.4.