Towards a precise determination of the scattering amplitudes of the charmed and light-flavor pseudoscalar mesons

We study the scattering of the light-flavor pseudoscalar mesons ($\pi, K, \eta$) off the ground-state charmed mesons ($D,D_s$) within chiral effective field theory. The recent lattice simulation results on various scattering lengths and the finite-volume spectra both in the moving and center-of-mass frames, most of which are obtained at unphysical meson masses, are used to constrain the free parameters in our theory. Explicit formulas to include the $S$- and $P$-wave mixing to determine the finite-volume energy levels are provided. After a successful reproduction of the lattice data, we perform a chiral extrapolation to predict the quantities with physical meson masses, including phase shifts, inelasticities, resonance pole positions and the corresponding residues from the scattering of the light pseudoscalar and charmed mesons.


Introduction
The spectroscopy of the open charmed mesons is an active and interesting research topic in hadron physics. The discovery of the scalar charm-strange meson D * s0 (2317) [1][2][3] challenges the quark model description [4], which predicts a mass around 160 MeV heavier than the experimental value. Another puzzle is that the mass of D * s0 (2317) is almost the same as the mass of its non-strange partner D * 0 (2400). The scattering process of the ground-state charmed mesons (D, D s ) and the light pseudoscalar mesons (π, K, η) offers an excellent environment to explore the properties of the scalar charmed resonances D * s0 (2317), D * 0 (2400) and possible resonances with other quantum numbers as well.
the charmed and the light pseudoscalar mesons is carried out in Ref. [14]. It is found that the massive η ′ meson plays a minor role in the energy region considered, therefore we work in the conventional SU (3) chiral Lagrangian in this work.
We briefly introduce the relevant SU (3) chiral Lagrangians to set up our notations. The ground-state charmed-meson triplet P = (D 0 , D + , D + s ) is incorporated in the chiral Lagrangians as a matter field. The light pseudoscalar mesons π, K and η are treated as pseudo-Nambu-Goldstone bosons (pNGBs). The leading-order (LO) chiral Lagrangian describing the interactions between the pNGBs and the charmed mesons reads L (1) where M D denotes the mass of the charmed-meson triplet in the chiral limit. The covariant derivative D µ is given by where Here, F denotes the weak decay constant of the pNGBs in the chiral limit, with the normalization F π = 92.1 MeV. The NLO Lagrangian, with six additional low energy constants h i=0,...,5 , takes the form [7,9] L (2) Pφ with where s and p denote the scalar and pseudoscalar external sources, respectively. By taking (s + ip) = diag(m,m, m s ), withm the average of up-and down-quark mass and m s the strange quark mass, one can introduce the light-quark masses in the chiral Lagrangian. We do not consider any isospin violation effect in this work. At leading order, the quantity B in Eq. (5) is related to the light-quark condensate through 0|q i q j |0 = −F 2 Bδ ij . The LO squared masses of the pNGBs are then given by With different combinations of the strangeness (S) and isospin (I), the scattering amplitudes of the ground-state charmed mesons and the pNGBs are classified into seven different cases. See the first and second columns of Table 1 for the specific channels involved in each case. For the process D 1 (p 1 ) + φ 1 (p 2 ) → D 2 (p 3 ) + φ 2 (p 4 ) with definite strangeness and isospin, the general scattering amplitude takes the form V (S,I) where s = (p 1 + p 2 ) 2 = (p 3 + p 4 ) 2 , t = (p 1 − p 3 ) 2 = (p 4 − p 2 ) 2 , u = (p 1 − p 4 ) 2 = (p 3 − p 2 ) 2 correspond to the standard Mandelstam variables, and the functions H 24 (s, t, u) and H 35 (s, t, u) are given by H 24 (s, t, u) = 2h 2 (p 2 · p 4 ) + h 4 [(p 1 · p 2 )(p 3 · p 4 ) + (p 1 · p 4 )(p 2 · p 3 )] , The coefficients C i in Eq. (7) have been given in many works [9,21,22] and we show their expressions in Table 1 for the sake of completeness. The results from the generalization to the U (3) case with explicit η ′ meson have been given in Ref. [14].
In the present work we mainly focus on the S-wave scattering of the pNGBs and the charmed mesons. In order to obtain the partial-wave amplitudes, we need to perform the partial-wave projection of the full amplitudes in Eq.
where ϕ is the scattering angle of the incoming and outgoing states in the center-of-mass (CM) frame, and the Mandelstam variable t is related to ϕ through t(s, cos ϕ) = m 2 with λ(a, b, c) = a 2 + b 2 + c 2 − 2ab − 2bc − 2ac the Källén function. The S-wave amplitude can be obtained by taking J = 0 in Eq. (10). In later discussions the subscript J in the partial wave amplitude V (S,I) J, D 1 φ 1 →D 2 φ 2 (s) will be omitted for simplicity. The nonperturbative strong interactions of the pNGBs and the ground-state charmed mesons, which manifest themselves in the emergence of bound states or resonances, can be accounted for by restoring unitarity and the analytical properties associated with the unitarity cut of the perturbative partial-wave amplitudes in Eq. (10). In this work we use the unitarization approach that has been widely used to discuss the pNGBs and charmed mesons scattering in Refs. [8-10, 20, 22]. The unitarized amplitude for the two-body scattering process takes the form [36,37] T (s) = 1 − V(s) · G(s) −1 · V(s) , where V(s) denotes the partial-wave amplitude in Eq. (10) and for simplicity both the superscripts and subscripts are omitted. By construction, the G(s) function includes the two-body (S, I) . The quantum numbers of different channels are classified by strangeness (S) and isospin (I), as shown in the first column. unitarity/right-hand cut and it can be given by the loop function One can use a once-subtracted dispersion relation or dimensional regularization by replacing the divergence by a constant to calculate the explicit form of the G(s) function, which reads [36] G(s) DR = 1 16π 2 a(µ) + ln where and µ is the regularization scale. The superscript DR in Eq. (14) stresses that the G(s) function in this equation corresponds to the form obtained in dimensional regularization. The function G(s) DR does not depend on the regularization scale µ, since the explicit µ dependence in Eq. (14) is canceled by that from the subtraction constant a(µ). In later discussion we take µ = 1 GeV in order to allow for a comparison with the previous works [9, 13-15, 20, 22, 24]. The unitarized partial-wave amplitude in Eq. (12) can be easily extended to coupled-channel scattering, where one should promote V(s) and G(s) to n × n matrices in case of n channels. The matrix elements for V(s) are given by Eq. (10). G(s) becomes a diagonal matrix, with its diagonal elements given by Eq. (14) with the masses m 1 and m 2 in question. For easy comparison, we follow the previous works [14,24] for the convention of scattering the length. The S-wave scattering length is related to the unitarized chiral amplitude in Eq. (12) through where the superscripts for isospin and strangeness and the subscript for J = 0 are omitted for simplicity.

Chiral amplitudes in the finite volume
One of the main novelties in this work is to fully exploit the rich finite-volume spectra from the lattice simulations given in Ref. [28], in order to constrain the unitarized chiral amplitudes. In order to do so, we use the method proposed in Refs. [30,31] to introduce the finite-volume effects into the unitarized chiral amplitudes. As it was demonstrated in Ref. [38], this framework is quite efficient to fit the lattice finite-volume spectra for the coupled-channel scattering of πη, KK and πη ′ . In this work, we use the same method to study the coupled-channel Dπ, Dη and D sK system. Below, we briefly describe the method. The loop function G(s) in Eq. (13) is ultraviolet divergent and needs to be regularized. One way to do this is to perform the integral with the three-momentum cutoff q max . After integrating over the variable q 0 analytically, one gets where To obtain the above results when integrating out q 0 , it is convenient to choose the CM frame, by taking the total four-momentum P µ of the two-particle system as (P 0 , P = 0). Since the G(s) function in the infinite volume, i.e. Eqs. (13) or (17), is a Lorentz scalar, its final expression is the same in different frames. However due to the breaking of Lorentz invariance in the finite volume, one should distinguish the finite-volume quantities defined in different frames. The quantities in the CM frame will be denoted with an asterisk in the following. The finite-volume effects are introduced into the unitarized chiral amplitudes by discretizing the above three-momentum integral, defining the loop function. The allowed momenta q * in the cubic box of length L with periodic boundary conditions take the discrete values The three-momentum integral in Eq. (17) should be replaced by the sum of the allowed momenta. Hence, the finite-volume loop function reads Here we introduce a tilde on top of a symbol to distinguish it from the same quantity in the infinite volume. The finite-volume correction ∆G in the CM frame to the loop function G(s) is then given by It should be stressed that, as L → ∞, the quantity ∆G is independent of the three-momentum cutoff due to the cancellation of the q max -dependences of the two terms in this equation and, up to the terms that vanish exponentially at large L, can be related to the pertinent Lüscher zeta-function. In practice, for finite L, it was verified numerically (see Ref. [38]) that the cutoff dependence of ∆G is indeed rather weak. The final expression of the function G(s), used in our finite-volume analysis, takes the form where G DR and ∆G are explicitly given in Eqs. (14) and (21), respectively. As mentioned previously, although the loop function G DR (s) in the infinite volume is Lorentz invariant, the corresponding finite volume expression in Eq. (22) is not Lorentz invariant any more. As a result, one has to explicitly work out the different expressions for the loop functions in different frames, which are considered in Refs. [31,32,39,40]. Here we recapitulate the main results to set up the notation.
For the two-body system, moving with the four-momentum P µ = (P 0 , P ), the CM energy squared is s = E 2 = (P 0 ) 2 − | P | 2 and the three-momenta of the particles in the moving frame are q 1 and q 2 , respectively, with q 1 + q 2 = P . The corresponding three-momenta in the CM frame are denoted by q 1 * and q * 2 , respectively, with q * 1 = − q * 2 = q * . By performing the standard Lorentz boost, one obtains where the on-shell energies q * 0 1 and q * 0 2 take the form With these definitions, the finite-volume loop function in the moving frame reads [31] It is obvious that the expression in Eq. (25) in the moving frame recovers the formula of Eq. (20), defined in the CM frame with P = 0. In analogy with Eq. (22) in the CM frame, the final expression for the loop function used in the moving frame takes the form where with G DR , G MV and G cutoff given in Eqs. (14), (25) and (17), respectively 1 . In order to account for the higher partial waves in the determination of the finite-volume energy levels, the generalized G(s) functions are introduced [31] where | q on * | denotes the on-shell value for | q * |,q * = q * /| q * |, k = 0 (1) for ℓ + ℓ ′ = even (odd), and the Y ℓm denote the spherical harmonics functions with the normalization One can establish the relation of G MV ℓm,ℓ ′ m ′ in Eq. (29) with M ℓm,ℓ ′ m ′ (the linear combination of the Lüscher zeta functions) in Eq. (39) of Ref. [32]. See also Refs. [41,42] for further details on M ℓm,ℓ ′ m ′ .
Further, it is convenient to introduce the quantity which, up to the exponentially suppressed terms, is related to the quantity w ℓm , defined in Eq. (40) of Ref. [32], through In analogy to Eq. (27), we define where with G DR , G MV ℓm and G cutoff given in Eqs. (14), (31) and (17), respectively. It is easy to show that by taking ℓ = 0, m = 0 and P = 0, Eq. (33) reduces to the CM formula of Eq. (22), as it should be. In order to simplify the notation, we will denote G DR,MV ℓm in Eq. (33) by G ℓm in the following.
Due to the rotational invariance, different partial waves do not mix in the infinite volume. However, this feature is lost in a finite volume, and the different partial-wave amplitudes V J (s) in Eq. (10) get mixed. A more subtle issue is that the mixing patterns of the partial-wave amplitudes vary in different moving frames. In the following, we shall retain only the S-and Pwave amplitudes of the Dπ, Dη and D sK system, which should be a reasonable approximation up to the D sK threshold energy region [28].
The projection of the two-body quantization condition onto the irreducible representations of the different little groups of the octahedral group O h , corresponding to the different moving frames, has been carried out in all details in Ref. [32]. In this paper, we wish to adapt these results for the case of the unitarized ChPT in a finite volume. The pertinent formulas can be directly read off from Ref. [32], replacing w ℓm by the quantity G ℓm introduced above, and keeping track of the normalization factors. Of course, in the present work we consider the coupled-channel scattering, but this does not change the symmetry properties of the equation as, simply, in addition, the amplitudes become matrices in the channel space.
Below, we display the explicit equations in different frames. In the CM frame, there is no mixing between S-and P -wave scattering amplitudes. For the S-wave in the A + 1 irreducible representation, the finite-volume energy levels are given by the solutions of the equation [31,32] where I is the unit matrix and the matrix elements of V 0 (s) and G 00 can be calculated via Eqs. (10) and (33), respectively. Further, according to Ref. [28], there exists a bound state in the P -wave Dπ scattering. In order to consider the contribution of the P -wave to the energy levels, we use a simple ansatz to include the bound state where the superscripts (0, 1 2 ) of V 1 are omitted for later convenience, and g V and m D * shall be adjusted to reproduce the lattice energy levels. For the P wave in the T − 1 irreducible representation, the finite-volume energy levels are determined by [31] where V 1 (s) and G ℓm are given in Eqs. (36) and (33), respectively. We mention that the determinant in the above equation is in fact trivial, since the single-channel approximation is used for the P -wave scattering.
In the moving frame with the total three-momentum P = (2π/L) N , the S-and P -wave amplitudes will get mixed. For the moving frame with N = (0, 0, 1), the equation to determine the discrete energy levels in the irreducible representation A 1 is where Here, V 0 , V 1 and G ℓm should be understood as matrices in the scattering-channel space. To be more specific, the S-wave V 0 corresponds to a 3 × 3 matrix, spanned by the Dπ, Dη and D sK channels. For V 1 it is an ordinary function, since the single-channel approximation is taken for the P wave. As a result, the 4 × 4 matrix of V 0,1 · M A 1 0,1 in Eq. (38) is given by where i and j in the subscripts of V 0,ij and G ℓm,i are the channel indices. The Dπ, Dη and D sK channels are labeled by 1, 2, and 3, respectively. For other moving frames, the corresponding equations to determine the discrete energy levels for the irreducible representation A 1 can be obtained by replacing the M A 1 0,1 in Eq. (38) with the proper ones, which are given in Ref. [32]. We quote the explicit results below for completeness. For N = (1, 1, 0), it is For N = (1, 1, 1), it is The partial-wave scattering amplitudes V 0 and V 1 are the same as those in Eq. (38) for different moving frames and different irreducible representations. For the irreducible representations E when N = (0, 0, 1), B 1 and B 2 when N = (1, 1, 0) and E when N = (1, 1, 1), the S-wave amplitudes are decoupled and only the P wave enters. The general equation to determine the discrete energy levels is given by the solution of According to Ref. [32], M 1 takes different forms for different representations. For the irreducible presentation E when N = (0, 0, 1), it reads For the irreducible presentation B 1 when N = (1, 1, 0), it reads For the irreducible presentation B 2 when N = (1, 1, 0), it reads For the irreducible presentation E when N = (1, 1, 1), it reads The partial-wave amplitude V 1 (s) in different representations takes the same expression in Eq. (36). All formulas, which are relevant for further discussions, were listed above. The formulas for other irreducible representations will not be explicitly given here. We refer to Ref. [32] for further details.

Fits to the finite-volume spectra and scattering lengths from lattice calculations
In order to precisely determine the scattering amplitudes of the charmed and light pseudoscalar mesons, we perform global fits to the discrete finite-volume spectra and the scattering lengths from several lattice calculations [24][25][26]28]. To be more specific, we include the finitevolume spectra, which were used in Ref. [28] to study the S-and P -wave Dπ, Dη and D sK coupled-channel scattering with I = 1/2, and which amount to 47 data points in total (38 data points below the D sK threshold) 2 . In addition, the elastic scattering lengths obtained with m π < 600 MeV from Ref. [24], which amount to 15 data points, are incorporated in our fits.
The 2+1 flavor lattice calculation of the DK scattering length and the two energy levels well below the D s η threshold with (S, I) = (1, 0) from Refs. [25,26] are also considered in the global fits, which amount to 3 additional data points. In the previous reference, the DK scattering lengths with relatively small statistical uncertainties are obtained from the lowest two energy levels within the effective-range-expansion framework, which provides an efficient method to determine quantities at thresholds. We have also tried to only fit the lowest two energy levels, which turn out to be quite close to the present results.
Regarding the finite-volume spectra from the S-and P -wave coupled Dπ, Dη and D sK scattering with I = 1/2, in one fit strategy we use exactly the same 47 data points as those in Ref. [28] to determine the scattering amplitudes, which amounts to 65 data points in total. The present study is based on the chiral amplitudes with NLO local interactions and all the bound states or resonances are generated through the unitarization procedure. Therefore it is not expected that we can reliably describe the strong dynamics well above the scattering threshold. To make an estimate of the systematic error, in another fit strategy we only include the 38 points below the D sK threshold among the overall 47 data from Ref. [28] in our study, which amounts to 56 data points in total.
Before going to the details of the fits, we comment on the value of pion decay constant F π appearing in Eq. (7). One approach is to use the physical value F π = 92.1 MeV when fitting the lattice results calculated at unphysical pion masses, as done in Refs. [14,24]. Another approach is to use the unphysical F π values at the corresponding unphysical pion masses. For the lattice data in Ref. [24], the values of F π have been calculated in Ref. [43] and we take the values therein. The F π value for the lattice used in Ref. [28] has not been given from lattice calculation. We take the chiral extrapolated value F π = 105.9 MeV determined in our previous work [38]. The physical F π value is used to study the lattice data of Refs. [25,26], since the pion mass used in the lattice calculation is quite close to the physical value. From the chiral power counting point of view, there is no preference as to which approach to use up to the order considered in this work. The discrepancy resulting from the two approaches can be considered as a systematic uncertainty. One may also think of introducing different pNGB decay constants, such as F π and F K , to different channels. In practice we do not expect this effect as important as the differences between the physical and unphysical values of F π when performing the chiral extrapolation in next section. This has been explicitly verified in Ref. [38], where we have carried out the calculation to study the effects by using different pNGB decay constants in different channels, which indeed turn out to be small. One of the reasons is that the shift of F K when extrapolating m π from 391 MeV [28] to the physical value is very moderate, which is estimated to be from 115 MeV to 110 MeV in Ref. [38]. As a result, we do not introduce another type of fit to distinguish different pNGB decay constants in this work.
Four different types of fits are performed in our study by using different data sets and different F π values as discussed above. In the following we denote the four fits by Fit-1A, Fit-1B, Fit-2A and Fit-2B. In the notations, 1 and 2 stand for the two different data sets. 1 is for the 56 data points and 2 the 65 data points. A and B stand for the two different choices of the F π values. A means using the unphysical F π values for the lattice data at the unphysical pion masses, while B means using physical F π for all lattice data.
When only fitting the elastic scattering lengths from the lattice simulations in Refs. [24][25][26], it was found that for all the channels one common subtraction constant, defined in Eq. (14), is able to satisfactorily reproduce the lattice results [13,14,16,20,22,24]. The same value of the subtraction constant a(µ) determined from the elastic channels [24] was also used in the coupled-channel Dπ, Dη and D sK S-wave scattering to predict resonance poles of D * 0 (2400) in Ref. [15]. However, it is found that one common subtraction constant is not sufficient any more when simultaneously including the finite-volume spectra [28] and the elastic scattering lengths [24][25][26] in the global fits. For the S-wave coupled-channel Dπ, Dη and D sK scattering with I = 1/2, two subtraction constants a 0,1/2 Dπ and a 0,1/2 Dη are needed to reasonably describe the finite-volume spectra. We fix the subtraction constant in the D sK channel to be the same as in the Dη channel since, as seen a posteriori, the fit quality does not improve in general by introducing a free D sK subtraction constant. Further, in our study we use the single-channel formula to fit the two energy levels well below the D s η threshold in analogy to the lattice study of the I = 0 DK scattering [25,26], which did not consider the coupling of D s η channel. We find that a common subtraction constant a 1,0 DK for both the DK and D s η channels is able to well reproduce the lattice scattering length. For all other channels listed in Table 1, a common subtraction constant a EC is used and we find it is sufficient to describe the scattering lengths given by the lattice calculation [24].
The elastic P -wave scattering amplitude in Eq. (36) is incorporated in our study to describe the finite-volume spectra in Ref. [28]. According to the energy levels in Fig. 3 of that reference, clearly there is a bound state well below the Dπ threshold in the P -wave amplitude. Furthermore, the similarities of the lowest levels in Figs. 2 and 3 in Ref. [28], indicate that the S-and P -wave mixing effects are weak, which also justifies the elastic approximation of the P -wave amplitude. The lowest energy levels in Fig. 3 of Ref. [28] are dominated by the P -wave amplitude, which determines m D * = 2009 MeV in Eq. (36). For the coupling g V , we find that the fits are rather insensitive to its value. Therefore we fix m D * = 2009 MeV and g V = 3 in the following discussions. It is verified that the fits are barely affected by varying g V in a wide range from 0.5 to 5. The subtraction constant in the P -wave amplitude is fixed to be equal to the value in the S-wave case.
There are six LECs h i=0,...,5 in the NLO scattering amplitude. The values of h 0 and h 1 can be fixed to be h 0 = 0.033 and h 1 = 0.43 from the masses of D and D s [14], comparing with the slightly different values h 0 = 0.014 and h 1 = 0.42 used in Ref. [24]. We still have 8 parameters, i.e. the remaining four LECs h i=2,3,4,5 and the four subtraction constants a 0,1/2 Dπ , a 0,1/2 Dη , a 1,0 DK and a EC , which need to be determined from the fits to the lattice data. As has been done in Refs. [14,24], we redefine the LECs h i=2,3,4,5 as follows in order to reduce the correlations in the fits: whereM D ≡ (M phys D + M phys Ds )/2. Unlike the subtraction constants that each of them can only enter in a specific channel, every single chiral LEC could appear in all the scattering amplitudes. This is another reason that urges us to perform global fits by including the finite-volume energy levels of the coupled-channel Dπ, Dη and D sK scattering [28] as well as the scattering lengths of various channels given in Refs. [24][25][26]. The values of the parameters from the four types of fits are collected in Table 2. The results from Ref. [24] are also presented in the last column for comparison.
We would like to mention that the correlations between different energy levels within the same volume from Ref. [28] are included in our fits. If the correlations are neglected, the resulting χ 2 will be greatly reduced, which turn out to be 75.5, 87.6, 134.6 and 132.0 for Fit-1A, Fit-1B, Fit-2A and Fit-2B, respectively. Further, three sets of data from different lattice collaborations using rather different ensembles are included in our fits and they may introduce potentially large systematical uncertainties, which are difficult to estimate and hence are not  included in this work. This provides another explanation of the somewhat large χ 2 from our fits. In fact, we have tried to only fit the 47 data from the HSC for the coupled Dπ, Dη and D sK scattering with I = 1/2. By releasing all the six chiral low energy constants in Eq. (4), it is possible for us to obtain much smaller χ 2 values for the HSC data [28]. However the resulting values of h 0 and h 1 , are significantly different from the results of Refs. [14,24], which are determined by properly reproducing the masses of the grounds-state charmed mesons. Given the fact that the χ 2 /d.o.f. for the Fit-2A and Fit-2B are around 4 one could consider the possibility to double the relative errors of the data fitted in order to estimate the precision achieved by employing our parameterization based on unitarized ChPT. By taking into account that the relative errors for the data of Ref. [28] range in an interval of around 0.05 − 0.6%, this would imply that we are able to give a fair reproduction of the lattice QCD data at the level of a 0.1 − 1.2%, which indeed is a great achievement for a parameterization based on unitarized SU (3) NLO ChPT. The latter is expected to be affected by errors from higher-order corrections at the level of [(m π ∼ m K )/1 GeV] 3 , i.e. around 6 ∼ 15% for the unphysically large meson masses used here. One should take into account that by unitarizing ChPT the resulting parameterization is expected to be more precise, particularly if the data reflect the presence of resonances that are properly reproduced with the nonperturbative approach. One possible way to improve the discussions is to generalize the present study to next-to-next-to-leading order [13,16], which is clearly beyond the scope of this work. As a result, we shall focus on the more constrained fits shown in Table 2 in the following discussions. As shown in Table 2, the χ 2 values resulting from Fit-2A and Fit-2B, which include the finite-volume energy levels above the D sK threshold from Ref. [28], are clearly larger than those from Fit-1A and Fit-1B that only include the finite-volume spectra below the D sK threshold. We further verify that the extra amounts of the χ 2 from Fit-2A and Fit-2B are mainly contributed by the finite-volume energy levels from the Dπ, Dη and D sK coupledchannel scattering. Comparing the four types of fits with the results in the last column in Table 2, we observe that the parameters from Fit-2B are the closest to the values given in Ref. [24]. According to the large N C arguments, h 2 is expected to be 1/N C suppressed comparing with h 3 . The same expectation is also applied to h 4 and h 5 . In Ref. [45], it provides another useful theoretical criteria to discriminate different parameter sets, which relies on the positive constraints of the scattering amplitudes. If one considers the N C argument and the positivity bound, it is plausible that Fit-2B is the preferred one comparing with the other three fits in Table 2. We also find an additional solution for Fit-2A, which gives similar total χ 2 . However, the additional solution gives a worse description of the elastic scattering lengths in Ref. [24] and the energy levels of the DK scattering in Refs. [25,26] than the other fits in Table 2. In this respect, the other additional solution of Fit-2A is considered to be disfavored and we refrain from discussing the results from that solution. Red squares denote the results from Fit-2B and the gray shaded areas denote the corresponding one-sigma error bands. The lattice data are taken from Ref. [28].
With all the parameters determined from the fits, we can reproduce the finite-volume energies of the scattering channels considered in this study. The reproduced energy levels as a function of the box size L in various channels together with the lattice data are presented in Figs. 1, 2 and 3. Fig. 1 is for the I = 1/2 coupled-channel scattering of Dπ, Dη and D sK . Fig. 2 is for the I = 1/2 P -wave Dπ scattering and the I = 0 S-wave DK scattering. Fig. 3 is for the I = 3/2 S-wave Dπ scattering, which is not included in the fits and is a prediction of our study. One can see that our theoretical formalism can well reproduce the lattice results. Similarly, we can also reproduce the scattering lengths given by the lattice calculations [24,26]. This is shown in Fig. 4. Notice that only the data points with m π < 600 MeV from Ref. [24] are included in the fits. We explicitly show the fit results from Fit-1B and Fit-2B, where F π is fixed at its physical value. Both the central values and one-sigma statistical error bands from Fit-2B are explicitly given. In order to not overload the figures, we only show the central-value curves for Fit-1B. It is clear that our theoretical formalism can also well describe the various scattering lengths from the lattice calculations.
Having determined all the unknown parameters and verified the reliability of our fits, we proceed to discuss the resonance structures in the scattering amplitudes in the next section.

Scattering amplitudes and resonances at unphysical meson masses
In this part, we study the infinite-volume amplitudes of the Dπ, Dη and D sK coupledchannel scattering obtained at the unphysically large meson masses used in Ref. [28]. The phase shifts (δ) and inelasticity parameters (ε) are related to the S matrix, which is given by the unitarized scattering amplitude T in Eq. (12) through with To be more specific, the phase shifts δ kk and δ kl and the inelasticity parameters ε kk and ε kl , with k = l, are related to the matrix elements S kk and S kl through For the inelasticity parameters ε kk , one has 0 ≤ ε kk ≤ 1.
The phase shifts and inelasticities of the S-wave coupled-channel Dπ, Dη and D sK scattering with (S, I) = (0, 1/2) are given in the left and right panels in Fig. 5, respectively. We show the representative results with both central values and the statistical uncertainties at the one-sigma level from Fit-2B. Within uncertainties, we observe that there are two branches of the Dπ phase shifts in the energy region around E > 2530 MeV. The two branches of phase shifts in fact correspond to similar physical dynamics, since they differ by 180 degrees. Although the Dπ phase shifts around 2530 MeV show large uncertainties, the inelasticities in the same energy region almost vanish, indicating that the underlying dynamics of the S matrix in this region shows a unique feature. The central-value plots and uncertainties for the phase shifts and inelasticities of the Dη and D sK channels from Fit-2B are also shown in Fig. 5. In order to not overload the figure, we only give the central-value phase shifts and inelasticities for the Dπ channel from Fit-1A, Fit-1B and Fit-2A. Roughly speaking, the Dπ phase shifts below the first inelastic channel from different fits are quite compatible. In the region above 2450 MeV, the Dπ phase shifts and inelasticities start to show different behaviors from different fits. The results from Fit-1A and Fit-1B, which include the same lattice simulation data up to the D sK threshold, show somewhat similar behaviors. The results from Fit-2A and Fit-2B, which includes the lattice data above the D sK threshold from Ref. [28], are clearly different from those from Fit-1A and Fit-1B in the energy region above 2450 MeV. The resulting plots from Fit-2A and Fit-2B, which include the same data sets, are compatible with each other within the uncertainties. The scattering amplitudes in the inelastic energy region are clearly affected by the different data included in the fits. In contrast, the amplitudes in the elastic energy region show quite consistent behaviors.
In order to study the resonance poles, we need to perform the analytical continuation of the scattering amplitudes to the complex energy plane. In our formalism, this can be easily achieved by modifying the G(s) function in Eq. (14). For each channel, one can define two Riemann sheets (RS's) for the G(s) function. The formula in Eq. (14) is the corresponding expression on the first RS. The expression on the second RS is given by [46] with G(s) DR and σ(s) defined in Eqs. (14) and (15), respectively. In this convention σ(s) has to be calculated with Im σ(s) > 0 in the complex s plane. This implies that the signs of Im G(s) along the real s axis above the threshold on the first and second RS's are opposite. For a n-channel problem, one can then define 2 n RS's. The first RS will be denoted as (+, +, +, · · · , +). The second, third and fourth RS's are labeled as (−, +, +, · · · , +), (−, −, +, · · · , +) and (+, −, +, · · · , +), respectively. The plus and minus signs correspond to the G(s) function of this channel evaluated on the first and second RS's, respectively. Apart from the pole position s P itself, one can also calculate the residues γ at s P , which are given by with γ an n row vector and its transpose γ T = (γ 1 , γ 2 , · · · , γ n ). The residues correspond to the coupling strengths of the resonance pole to the interacting channels and encode important information of the resonance. With the pole position and its residues, one can then further discuss the composition of the resonance. In Ref. [47], the calculation for the compositeness coefficient is generalized to the resonances by extending Weinberg's bound-state compositeness relation [48]. Then one can interpret the values of the compositeness X as the probabilities to find the twobody components inside the resonances and bound states. The prescription to calculate the compositeness coefficient X i contributed by the ith channel in Ref. [47] is where the function G(s) is given in Eq. (14) or (53), depending on the location of the pole. The total compositeness X is then given by X = i X i , with the sum only spanning on the open channels for the resonance in question, which are the channels below the real part of the pole. We mention that the prescription in Eq. (55) only applies to the canonical resonance, in the sense that the resonance pole should reside in the RS that can be directly accessed from the physical RS. For a near-threshold bound state, Eq. (55) recovers the Weinberg compositeness [48]. We refer to Ref. [47] for further details and also to Ref. [49], where the general framework for the calculation of the compositeness for poles is developed.
In Table 3, we give both the resonance pole positions and their residues for the coupledchannel Dπ, Dη and D sK S-wave scattering with (S, I) = (0, 1/2) obtained at the unphysically large meson masses used in Ref. [28]. The most robust conclusion from Table 3 is that there is a bound state just below the Dπ threshold. Within uncertainties the four fits lead to compatible results for the pole positions and also the corresponding residues. Our determinations of the bound state pole are close to the value in Ref. [28], which is 2275.9 ± 0.9 MeV. Furthermore, the results from Fit-2A and Fit-2B, which include exactly the same lattice data of the Dπ coupled-channel scattering as Ref. [28], are perfectly compatible with the value given in the former reference. This presents a nice crosscheck with our chiral amplitudes assisted finitevolume study, comparing with the K-matrix assisted Lüscher formula in Ref. [28]. While in Ref. [15], the bound state pole is predicted to be 2264 +8 −14 MeV, with much larger error bars than ours and those in Ref. [28]. The possible reason is that the values of the chiral LECs and the subtraction constant used in Ref. [15] are taken from Ref. [24], which were determined by only including the elastic scattering lengths of the latter reference. Regarding the coupling strengths of the bound state, our study shows that this pole is more strongly coupled to the D sK channel than to the Dη one, which is also the case of Ref. [15] but differs from the results of Ref. [28]. By applying Eq.  Table 3. Our results for the resonance poles on the third RS and their residues are consistent with the determinations in Ref. [15]. However the resonance poles, either on the second or the third RS's, are not reported in Ref. [28]. Apart from the poles shown in Table 3, we find that there are also other heavier poles in the region around or above 2600 MeV. The heavier poles appear in different RS's depending on the different fits and the uncertainties of their masses are usually large. Since these heavier resonances are much less constrained in our fits and show large model dependences, we refrain from discussing further about their properties.
According to Ref. [47], the prescription in Eq. (55) can be applied to the third-sheet poles from Fit-2A and Fit-2B in Table 3 For the P -wave Dπ scattering, a bound state pole with the mass in the range 2008.2 ∼ 2009.8 MeV is found by combining the results from the four fits. Our determination is in good agreement with the value 2009 ± 2 MeV in Ref. [28].
Having shown the scattering amplitudes and the resonance structures at the unphysical meson masses, we proceed the study for the physical meson masses in the following section.

Chiral extrapolation to the physical meson masses
By assuming that the free parameters in Table 2 are independent on the light-flavor meson masses, it is straightforward to perform the chiral extrapolation to the physical meson masses in our study. The phase shifts and inelasticities of the S-wave coupled-channel Dπ, Dη and D sK scattering obtained at physical meson masses are given in Fig. 6. The central-value plots and the statistical uncertainties at the one-sigma level from Fit-2B are shown. As in Fig. 5, we give the central-value curves for Fit-1A, Fit-1B and Fit-2A in Fig. 6. In order to clearly demonstrate the resonance structures, the magnitudes of the scattering T matrices are provided in Fig. 7. One can clearly see the discrepancies from different fits.
Another subtlety issue on the chiral extrapolation is about the pion mass dependences of the parameters in Table 2. Note that the strange-quark mass is basically kept fixed to its physical value in the lattice QCD here considered. Clearly the chiral LECs by definition are independent of the pion mass. The subtraction constant a, on the contrary, could possibly vary with different pion masses, although many previous works simply assume the constant behavior of a when performing the chiral extrapolation [13][14][15][16]19,20,22,24]. One possible way to estimate the pion mass dependences is by comparison of the function G DR (s) in Eq. (14) to the three-momentum-cutoff version of G(s) [50].
First, the on-shell three-momentum is denoted by q(s) with We introduce the function δG(s) by rewriting G DR (s) in Eq. (14) as Next, we denote by G C (s) the function that results by evaluating the divergent integral of Eq. (13) with a three-momentum cutoff q max . An algebraic expression for G C (s) can be found in Ref. [50] Now, let us consider the possible pion mass dependence of the subtraction constants. One can work out explicitly the (nonrelativistic) limit of the functions G DR (s), Eq. (14), and G C (s) [50] for s → (m 1 + m 2 ) 2 , which implies |q(s)| ≪ m 1 , m 2 , q max . In this limit these functions are simply a constant plus −iq(s)/8π(m 1 + m 2 ) plus quadratic and higher-order terms in threemomentum. Therefore, we can write where α is a constant. For the case q(s) = 0, the following expression for the subtraction constant a is obtained The constant α does not depend on µ, since a(µ)−log µ 2 is µ independent, cf. Eq. (14). Eq. (59) reflects the splitting in the mechanisms underlying the contributions to the subtraction constant a. On the one hand, the last two contributions in this equation stem from rescattering effects (by taking q max around 1 GeV), associated with the unitarity cut. On the other hand, the former contribution (α) is associated with properties of contact terms (short-range physics). The crossed-channel contributions involving the explicit degrees of freedom in the effective field theory will be accounted for order by order in V(s) (12). Thus, the variation of α with the masses is at least quadratic in the pion mass. There is some remnant cut-off dependence in the splitting of Eq. (59) that could be ascertained by varying q max around q max ≃ 1 GeV, the typical scale for hadronic rescattering.
The fact that the leading correction to Eq. (59) is linear in m 2 implies a linear change in the pion mass for the Dπ subtraction constants, as m π /m D . Differently, for the other channels involving a K or an η the change will be just quadratic in m π and, therefore, much less important. The linearized version of Eq. (59) with respect to the smaller mass m 2 for two sets of values of the masses m 1 and m 2 (a ′ and a for m ′ i and m i , in order) gives The subtraction constant a is said to have its natural value [37] when the constant α has an absolute value much smaller than 1 for q max ≃ 1 GeV. One then has [37] a(q max ) = −2 log 1 + 1 + m 2 where the ellipses indicates higher order terms in the nonrelativistic expansion and in m 2 /m 1 with m 2 ≪ m 1 , as it follows directly from Eq. (59).
We take the Fit-2B as a concrete example to check the shifts of a by varying pion masses. When Eq. (60) is applied to a Dπ subtraction constant from a pion mass of 391 MeV to its physical value of 138 MeV we have a variation in the subtraction constant of −0.12. Compared to the value reported in Table 2 for a 0,1/2 Dπ , we obtain a mild effect of around a 10%. For the subtraction constants in the Dη and D sK channels, their values are kept fixed. In Fig. 8, we explicitly show the results by taking the extrapolated value a 0,1/2 Dπ = −1.64 obtained from Eq. (60), together with the figures by assuming pion mass independence of the subtraction constants from Fit-2B and also the results from Ref. [24]. The three different sets of plots in Fig. 8 reveal qualitatively similar resonant behaviors, although the heights of the peaks around the resonances at 2.1 GeV and 2.45 GeV are different. Comparing the curves from Fit-1A, Fit-1B and Fit-2A in Fig. 7 with those in Fig. 8, we observe that the discrepancies among the three different types of plots in Fig. 8, which include the additional uncertainties of the chiral extrapolation and the fitting results from the previous work [24], are clearly smaller than the differences of the four types of fits in the present work. Therefore in the following discussions, we shall concentrate on the results from the four types of fits in Table 2 without introducing the pion mass corrections to the subtraction constants. In Fig. 9, we show the phase shifts and inelasticities of the DK channel obtained at physical meson masses. In this case, the central-value results from Fit-1A and Fit-1B are almost identical, which lead to a sharp rise of the phase shifts near the DK threshold. The central-value behavior from Fit-1A and Fit-1B indicates a virtual pole near the DK threshold. In contrast, the phase shifts from Fit-2A and Fit-2B fall rapidly, which implies a bound state pole near the threshold. 3 The resonance poles and the corresponding residues from the S-wave coupled-channel Dπ, Dη and D sK scattering with (S, I) = (0, 1/2) obtained at physical meson masses are collected in Table 4. The first lesson we learn is that all the four fits give robust resonance poles on the second RS with the mass around 2100 MeV and the half width lying between 100 and 200 MeV. Furthermore, for all the four fits we also find their shadow poles on the third RS, whose masses and widths are quite close to the values on the second RS. Heavier resonance poles lying around 2.4 GeV are found as well. In Table 4, the relevant poles that are mostly responsible for the resonant behaviors on the physical sheet are given. For Fit-1A, the relevant pole is located on the second RS, and its shadow pole lies on the third RS. While for the Dπ are labeled as Fit-2B and a ′ , respectively. See the text for details. The results using the parameters from Ref. [24] are labeled as LOGHM.  [15]. We stress that the poles and their residues in Table 4 are only slightly affected by the pion mass dependences of the subtraction constants. We also take Fit-2B to demonstrate this point. Taking the chiral extrapolated value a 0,1/2 Dπ = −1.64 as explained previously, the pole around 2.1 GeV on the second RS is found at (2112.5 − i127.0) MeV, with the residues |γ 1 | = 9.9 GeV, |γ 2 /γ 1 | = 0.08, |γ 3 /γ 1 | = 0.58. The pole around 2.4 GeV on the third RS is at (2475.7 − i108.9) MeV, with the residues |γ 1 | = 6.7 GeV, |γ 2 /γ 1 | = 1.15, |γ 3 /γ 1 | = 2.07. These results are consistent with the values from Fit-2B that assumes the pion mass independence of the subtraction constants given in Table 4.
Next we discuss the compositeness for the resonance poles obtained at physical meson masses. According to Ref. [47], the prescription in Eq. (55) can be applied to the poles around 2.1 GeV on the second RS in Table 4.  Table 4, the poles on the third RS from Fit-1B, Fit-2A and Fit-2B are not far away from the region of validity according to Ref. [47]. As a rough estimate, we also   Table 3 for details.
For the D * s0 (2317) in the S-wave coupled-channel DK and D s η scattering, its poles and residues obtained at physical meson masses are summarized in Table 5. First we stress that each of the parameter configurations from all of the fits only gives one pole for the D * s0 (2317), either a bound state or a virtual state. For the parameters from Fit-2A and Fit-2B within onesigma uncertainties, all the parameter configurations only give the bound state poles on the first RS. While for Fit-1A and Fit-1B, within one-sigma uncertainties parts of the parameter configurations give the bound state poles on the first RS and others give the virtual poles on the second RS. E.g., with the central values of the parameters from Fit-1A and Fit-1B we only obtain the virtual poles. This tells us that the interactions are strong enough to produce a prominent enhancement below the DK threshold. However at the present stage we can not definitely conclude that the enhancement is caused by a bound or a virtual state. The findings of the bound and virtual states are consistent with the behaviors of the phase shifts shown in Fig. 9.
For the bound state poles, we can use Eq. (55) to calculate the compositeness coefficients for the D * s0 (2317). The compositeness coefficients contributed by the DK and D s η for the Fit-2A case are 0.72 +0.14 −0.13 and 0.16 +0.04 −0.07 , respectively. The corresponding values from Fit-2B are 0.77 +0.11 −0.13 and 0.11 +0.03 −0.04 , which are compatible with those found in Ref. [47]. The robust conclusion from these numbers is that the DK component is the dominant one inside the D * s0 (2317). Before ending the phenomenological discussion, we give the predictions for the S-wave scattering lengths of various channels at the physical meson masses in Table 6 Table 6, we provide the values with statistical uncertainties. This is because within the one-sigma fitted parameter configurations the scattering lengths for these four channels vary from huge negative values to huge positive values. The reason behind is that for these channels parts of the parameter configurations could lead to bound state poles near threshold, which correspond to large negative scattering lengths, and others could give virtual poles near threshold, which correspond to large positive scattering lengths. These findings are consistent with the pole contents discussed in Table 5 for the D * s0 (2317). We also verify that similar situations happen for the S-wave DK scattering with (S, I) = (−1, 0). The results from Fit-2B in Table 6, which gives the closest values of the parameters to Refs. [14,24], are qualitatively compatible with the numbers of the former references within uncertainties.  Table 5: Poles and their residues obtained at physical meson masses from the S-wave coupledchannel DK and D s η scattering amplitudes with (S, I) = (1, 0). For Fit-1A and Fit-1B, parts of the parameter configurations within one-sigma uncertainties give bound state poles for the D * s0 (2317) and others give virtual poles. In these cases, we simply show the ranges for the bound-and virtual-state poles obtained for the parameter configurations from Fit-1A and Fit-1B within one-sigma uncertainties. For Fit-2A and Fit-2B, all the parameter configurations within one-sigma uncertainties give bound state poles. The indices of the residues γ i=1,2 correspond to the channels DK and D s η in order. The physical thresholds of the DK and D s η channels are 2362.8 and 2516.2 MeV, respectively. For the definition of the different RS's, see the text.

Conclusions
In this work we simultaneously analyzed the lattice finite-volume energy levels and the scattering lengths for the scattering of the charmed and light pseudoscalar mesons from Refs. [24][25][26]28] within the chiral effective field theory. Several different fit strategies, by using different values for the pion decay constant and including different data sets in the fits, are explored in our study. Through the fits of the lattice data, we fix the values the chiral low-energy constants up to the next-to-leading order and the subtraction constants introduced in the unitarization procedure. The updated values for the low energy constants and subtraction constants are collected in Table 2, which provide a useful starting point for future study on the charmed resonance dynamics in various processes.
The scattering amplitudes and the resonance poles obtained at the physical meson masses provide the most important outputs of this work. Regarding the resonance poles in the coupledchannel Dπ, Dη and D sK S-wave scattering with (S, I) = (0, 1/2), a robust conclusion in  our study is that there is pole with the mass around 2100 MeV and the width more than 200 MeV, see Table 4. Another type of heavier poles lying between 2300 MeV and 2500 MeV, depending on different fits, are also found, with their widths varying from 70 MeV to 200 MeV. According to Fig. 6, the physical S-wave Dπ phase shifts and inelasticities with I = 1/2 show somewhat different behaviors from different fits, specially in the energy region above 2350 MeV. To implement the scattering amplitudes obtained here in the semileptonic charmed meson decays [51] or the phenomenological study of B decays [44] may offer us another way to further discriminate the different fits. For the phase shifts and inelasticities of the S-wave coupled-channel DK and D s η scattering with (S, I) = (1, 0), we find two different types of solutions. In one solution, i.e. the lower branch of the phase shifts in Fig. 9, the D * s0 (2317) corresponds to a bound state. While the upper branch of the phase shifts in the former figure implies a virtual state nature of the D * s0 (2317). Future lattice simulations with more energy levels in this channel may enable us to discriminate the two different solutions.