On chiral extrapolations of charmed meson masses and coupled-channel reaction dynamics

We perform an analysis of QCD lattice data on charmed meson masses. The quark-mass dependence of the data set is used to gain information on the size of counter terms of the chiral Lagrangian formulated with open-charm states with J^P= 0^- and J^P =1^- quantum numbers. Of particular interest are those counter terms that are active in the exotic flavour sextet channel. A chiral expansion scheme where physical masses enter the extrapolation formulae is developed and applied to the lattice data set. Good convergence properties are demonstrated and an accurate reproduction of the lattice data based on ensembles of PACS-CS, MILC, ETMC and HSC with pion and kaon masses smaller than 600 MeV is achieved. It is argued that a unique set of low-energy parameters is obtainable only if additional information from HSC on some scattering observables is included in our global fits. The elastic and inelastic s-wave pi D and eta D scattering as considered by HSC is reproduced faithfully. Based on such low-energy parameters we predict 15 phase shifts and in-elasticities at physical quark masses but also for an additional HSC ensemble at smaller pion mass. In addition we find a clear signal for a member of the exotic flavour sextet states in the eta D channel, below the bar K D_s threshold. For the isospin violating strong decay width of the D_{s0}(2317) we obtain the range (104-116) keV.

In this work we wish to study the size of such chiral counter terms. First rough studies [5,6] suffer from limited empirical constraints. Additional information from first QCD lattice simulation on a set of s-wave scattering lengths was used in a series of later works [9][10][11][12]. Results are obtained that in part show unnaturally large counter terms and/or illustrate some residual dependence on how to set up the coupled-channel computation. Here we follow a different path and try to use the recent data set on the quark-mass dependence of the D meson ground-state masses [13][14][15][16][17][18][19]. This dynamics is driven in part by the counter terms that also have significant impact on the open-charm coupled-channel systems as discussed above. One may hope to obtain results that are less model dependent in this case.
However, it is well known that chiral perturbation theory formulated with three light flavours does not always show a convincing convergence pattern [20][21][22][23][24][25][26][27]. How is this for the case at hand? Only few studies are available in which this issue is addressed for open-charm meson systems. In a recent work the authors presented a novel chiral extrapolation scheme for the quark-mass dependence of the baryon octet and decuplet states that is formulated in terms of physical masses [28][29][30][31]. It is the purpose of our study to adapt this scheme to the open-charm sector of QCD and apply it to the available lattice data set. This requires in particular to consider the D mesons with J P = 0 − and J P = 1 − quantum numbers on equal footing. For a given set of low-energy constants each set of the four D meson masses has to be determined numerically as a solution of a non-linear system.
The work is organized as follows. In section 2 the part of the chiral Lagrangian that is relevant here is recalled. It follows a section where the one-loop contributions to the D meson masses are derived in a finite box. In section 3 and 4 power counting in the presence of physical masses is discussed. The application to available lattice data sets is presented in sections 5 and 6. Lattice data taken on ensembles of PACS-CS, MILC, HPQCD, ETMC and HSC are considered. In section 7 we present our predictions for phase shifts and inelasticities based on a parameter set obtained form the considered lattice data. In addition the fate of possible exotic states but also the isospin violating strong decay width of the D s0 (2317) is discussed. With a summary and outlook the paper is closed.

II. THE CHIRAL LAGRANGIAN WITH OPEN-CHARM MESON FIELDS
We recall the chiral Lagrangian formulated in the presence of two anti-triplets of D mesons with J P = 0 − and J P = 1 − quantum numbers [1,32]. In the relativistic version the Lagrangian was developed in [4][5][6]. The kinetic terms read Following [6] we represent the 1 − field in terms of an antisymmetric tensor field D µν . The covariant derivative · ∂ µ involves the chiral connection Γ µ , the quark masses enter via the symmetry breaking fields χ ± and the octet of the Goldstone boson fields is encoded into the 3 × 3 matrix Φ. The parameter f is the chiral limit value of the pion-decay constant.
Finally the parameters M and M + ∆ give the masses of the D and D * mesons in that limit with m u = m d = m s = 0.
We continue with first order interaction terms the parameterg P in (3) can not be extracted from empirical data directly. The size ofg P can be estimated using the heavy-quark symmetry of QCD [1,32]. At leading order one expectsg P = g P .
Second order terms of the chiral Lagrangian were first studied in [5,6], where the focus was on counter terms relevant for s-wave scattering of Goldstone bosons with the D mesons.
A list of eight terms with dimension less parameters c i andc i was identified. This list was extended by further terms relevant for p-wave scattering in [33]. A complete collection of relevant terms is where the parameter M and M + ∆ are the D and D * meson masses as evaluated at m u = m d = m s = 0. In the limit of a very large charm-quark mass M → M 0 ← M + ∆ a common mass M 0 arises. All parameters c i andc i are expected to scale linearly in the parameter M 0 . As illustrated in [6] it holdsc i = c i in the heavy-quark mass limit.
A first estimate of some parameters can be found in [6] based on large-N c arguments.
Since at leading order in a 1/N c expansion single-flavour trace interactions are dominant, the corresponding couplings should go to zero in the N c → ∞ limit, suggesting to first QCD lattice computations for s-wave scattering lengths of the Goldstone bosons with the D mesons. It is remarkable that their range for c 3 + c 5 1 is quite consistent with the earlier estimates [5,6] based on the empirical πD invariant mass spectrum. The c 3 parameter is of crucial importance for the physics of two exotic sextets of J P = 0 + and J P = 1 + resonances. Such multiplets are predicted by the leading order chiral Lagrangian (1), which entails in particular the Tomozawa-Weinberg coupled-channel interactions of the Goldstone bosons with the D mesons [4]. The latter predicts weak attraction in the flavour sextet channel. If used as the driving term in a coupled-channel unitarization exotic signals appear. A reliable estimate of the correction terms proportional to c 3 andc 3 is important in order to arrive at a detailed picture of this exotic sector of QCD [5,6].
We close this section with a first construction of the symmetry breaking counter terms proportional to the product of two quark masses: Such terms are relevant in the chiral extrapolation of the D meson masses. For the pseudoscalar mesons we provide the tree-level contributions to the polarization Π where we consider the isospin limit with m u = m d = m. Analogous expressions hold for the vector mesons polarization Π H∈ [1 − ] where the replacements c i →c i and d i →d i are to be applied to (9). Withc i = c i andd i = d i and ∆ → 0 the heavy-quark spin symmetry is recovered exactly.
We need to mention a technical issue. The propagator S αβ µν (p) of our 1 − fields involves four Lorentz indices, which are pairwise antisymmetric. Either interchanging α ↔ β or µ ↔ ν generates a change in sign. A mass renormalization from a loop contribution arises from a particular projection Π(p 2 ) of the polarization tensor Π µν αβ (p) with where d is the space-time dimension. This is the part which is used in (9) and will be used also in the following.

III. ONE-LOOP MASS CORRECTIONS IN A FINITE BOX
The chiral Lagrangian of section I is used to compute the D-meson masses at the one-loop level. In order to prepare for a comparison of QCD lattice data this computation is done in a finite box of volume V . A direct application of the relativistic chiral Lagrangian in the conventional M S scheme does lead to a plethora of power-counting violating contributions.
There are various ways to arrive at results that are consistent with the expectations of power counting rules [28,[34][35][36].
We follow here the χ-MS approach developed previously for the chiral dynamics of baryons [28,37,38], which is based on the Passarino-Veltman reduction scheme [39]. Recently this scheme was generalized for computations in a finite box [30]. Our results will be expressed in HQ and a set of generic loop functions. While the index H or R runs either over the triplet of pseudo-scalar or vector D mesons the index Q runs over the octet of Goldstone bosons (see Tab. I and Tab. II). In our case there will be two tadpole integralsĪ Q from the Goldstone bosons and the scalar bubble-loop integralĪ QR . In addition there may be tadpole contributionsĪ • any tadpole integral involving a heavy particle is dropped • the scalar bubble-loop integral requires a single subtraction.
The required loop functions have been used and detailed in a previous work [28,30,40] for finite box computations. For the readers' convenience we recall the loop functions in the infinite box limit [28,40] with where we note that in the infinite volume limit the two tadpole integralsĪ Q turn dependent and can no longer be discriminated in that case.
We point at the presence of the additional subtraction term γ H as suggested recently in [40] in the analogous case of a baryon self-energy computation. The subtraction term depends on the chiral limit values M and M + ∆ of the D and D * meson masses only. It was not yet imposed in earlier computations [28][29][30]. As was discussed in [40] the request of such a term comes from a study of the chiral regime with Within a counting scheme with m Q ∼ ∆ ∼ Q there is no need for any additional subtractions beyond the ones enforced by the χ-MS approach. However, to arrive at consistent results for m Q ∆ this subtraction is instrumental. While for ∆ ∼ m Q ∼ Q and γ H R = 0 the scalar bubble scales withĪ QR ∼ Q as expected from dimensional analysis, in the chiral regime with m Q ∆ and m Q ∼ Q one would expectĪ QR ∼ Q 2 ∼ m 2 Q . This expectation turns true only, for γ H R = 0 as chosen in (11).
We are now well prepared to collect all contributions to the D-meson self energies at the one-loop level. Consider the bubble loop and tadpole contributions. The Passarino-Veltman Π tadpole where the loop functions are expressed in terms of physical meson masses. The sums in (14,16) extend over intermediate Goldstone bosons (Q) and pseudo-scalar or vector D mesons QR are specified in Tab. I. In the contributions from the tadpole diagrams the sums in (15,17) Tab. II. The results (14,16) deserve a detailed discussion. First let us emphasize that a chiral expansion of the loop function as they are given confirms the leading chiral power as expected from dimensional counting rules. All power-counting violating contributions are subtracted owing to the χ-MS approach. Here we adopted the conventional counting rules which is expected to be effective for ∆ ∼ m Q . Our results (14,16) are model dependent, as there are various subtraction schemes available to obtain loop expression that are compatible with dimensional counting rules. Most prominently there is the infrared regularization of Becher and Leutwyler [35] and the minimal subtraction scheme proposed by Gegelia and collaborators [36]. Following our previous work on the chiral extrapolation of the baryon masses we will attempt to extract a model independent part of such loop expressions. This goes in a few consecutive steps. The driving strategy behind this attempt is to keep the physical masses inside the loop function.
Consider first the terms that are proportional to the tadpole loop functionĪ Q . There are two distinct classes of terms. The coefficient in front of anyĪ Q is either proportional to The terms proportional to m 2 QĪ Q or also toĪ (2) Q in (14,16) have the same form as the corresponding structures in (15,17) and therefore renormalize the low-energy parameters c n andc n with We conclude that the terms proportional to m 2 QĪ Q orĪ (2) Q in (14,16) may be dropped if we use the renormalized low-energy parameters c r n andc r n in the tadpole contributions (15, 17) but also in (24). Note, however, that by doing so some higher order terms proportional to with n ≥ 1 are neglected in Π H∈[1 − ] . We argue that the latter terms would cause a renormalization scale dependence that can not be absorbed into the available counter terms at the considered accuracy level. In order to avoid a model dependence such terms should be dropped.
We are left with the terms proportional to (M 2 R − M 2 H )Ī Q . If the charm meson masses are decomposed into their chiral moments the leading renormalization scale dependence of such terms can be absorbed into the Q 2 counter terms c 0,1 andc 0,1 . Similarly the components of order Q 4 can be matched with counter terms d n andd n . Most troublesome, however, are the subleading contributions proportional to m 5 Q log µ in such a strict chiral expansion of the vector D meson masses. There is no counter term available to remove such a scale dependence. In fact, only within a two-loop computation this issue is resolved in a conventional approach. Instead we keep the charm meson masses unexpanded in the terms (M 2 R − M 2 H )Ī Q and follow the strategy proposed in [40]. For those terms we provide the following decomposition where the second term depending on the heavy-meson tadpoleĪ R can be systematically dropped without harming the chiral Ward identities. We end up with the following renor- which will be the basis for our following studies. Note yet the additional subtraction terms α H QR in (22). Such terms were suggested in [40] for the analogous case of a baryon selfenergy computation. In order to arrive at consistent results for m Q ∆ the terms α H QR are instrumental: where the functions α i ,α i and γ i ,γ i depend on the ratio ∆/M only. They are listed in Appendix A and Appendix B. While the rational functions α i andα i all approach the numerical value one in the limit ∆ → 0, the functions γ i andγ i show a logarithmic divergence in that limit. We summarize the convenient implications of our subtraction scheme • the chiral limit values of the D meson masses are not renormalized • the low-energy parameters c 0,1 andc 0,1 are not renormalized • the wave-function factor of the D mesons are not renormalized in the chiral limit We close this section with a brief discussion on the role of the renormalization scale µ.
Given our scheme a scale dependence arises from the tadpole terms only. Such terms need to be considered in combination with the tree-level contribution Π where identical results hold for thec i andd i coupling constants. However, it is evident that scale invariant results follow with (24) only if the meson masses in the tadpole contributions are approximated by the leading order Gell-Mann Oakes Renner relations with m 2 π = 2 B 0 m and m 2 K = B 0 (m+m s ) for instance. This is unfortunate since we wish to use physical masses inside all loop contributions. Recalling our previous work [31] there may be an efficient remedy of this issue. Indeed the counter term contributions can be rewritten in terms of physical masses such that scale invariance follows without insisting on the Gell-Mann Oakes Renner relations for the meson masses. Such a rewrite is most economically achieved in terms of suitable linear combinations of the low-energy constants With Tab. III our rewrite is specified in detail. We assure that replacing the meson masses in the table by their leading order expressions the original expressions as given in (9) are recovered identically. We note a particularity: at leading order the effects of c 0 in G

IV. SELF CONSISTENT SUMMATION APPROACH
The physical D meson masses M H are determined from a set of coupled equations. This is so since the renormalized loop functions depend themselves on the physical masses of the D mesons. In a conventional chiral expansion scheme the meson masses inside the loop would be expanded to a given order so that a self-consistency issue does not arise. This is fine as long as the expansion is rapidly converging. For a slowly converging system such a summation scheme is of advantage even though this may bring in some model dependence [28][29][30]40].
Let us be specific on how the summation scheme is set up in detail. There is a subtle point emphasized recently in [40] which needs some discussions. The coupling constant g P was determined in [6] from the pion-decay width of the D * meson using a tree-level decay amplitude. Alternatively the decay width can be extracted from the D * meson propagator in the presence of the one-loop polarization Π bubble D * . The latter has imaginary contributions proportional to same coupling constant g 2 P that reflect the considered decay process. In the absence of wave-function renormalization effects one would identify a Breit-Wigner width by where the loop function is evaluated at the D * meson mass M D * . Both determinations would provide identical results. However, in the presence of a wave-function renormalization effect from the loop function this would no longer be the case. Following [40] we will therefore use the following form of the Dyson equation where we takeΠ 2 for the pseudo-scalar and vector D mesons respectively. The second order termsΠ (2) H are the tree-level contributions (9) proportional to the quark masses as written in terms of the parameters c 0 , c 1 andc 0 ,c 1 . The fourth order termsΠ (4−χ) H are the tree-level contributions (9) proportional to the product of two quark masses. Here the parameters d i andd i are probed. We recall that the wave-function renormalization Z H has a quark-mass dependence which cannot be fully moved into the counter terms of the chiral Lagrangian.
We provide a first numerical estimate of the importance of the various terms in (28).
We putΠ = 0 since the associated counter terms are not known reliably.
Insisting on the large-N c relations for the ratios of the Goldstone boson masses over the D meson masses. The mass differences of either pseudo-scalar or vector mesons  (22) are evaluated with the coupling constants g P =g P 0.57 and the physical isospin averaged meson masses. A decomposition according to (30,31) and (32) is performed.
can also be counted without controversy. In (31) we use a notation H R requesting Less obvious is how to treat the mass differences of a pseudoscalar and a vector D meson.
There are different schemes possible. Technically most straight forward is the extreme which can be motivated in the limit of a large charm quark mass where ∆ → 0 and therewith ∆ m π . In (32) we use a notation H ⊥ R implying that either H ∈ accurate to order Q 5 . The coefficients α H QR and γ H R were given already in (23) and (11). In Tab  How to improve on the counting rule (32). Before presenting a universal approach we consider yet two further interim power counting scenarios. First we work out the extreme chiral region where all Goldstone boson masses are significantly smaller than ∆ ∼ 200 MeV.
In this case the counting rules are used. Since the extreme chiral region is not realized in nature, such an assumption is not expected to provide any significant results for quantities measurable in experimental laboratories.
Since at some stage lattice QCD simulations may be feasible at such low strange quark masses we provide the corresponding expressions for the loop function nevertheless. Here we decompose all meson masses into their chiral moments in application of a strict chiral expansion. At third order The dimension less coefficients γ are expressed in terms of the basic coefficients α n , γ n andα n ,γ n in Appendix A and B. Again they depend on the ratio ∆/M only. We note that the rational functions α n andα n approach one in the limit ∆/M → 0. In contrast the γ n andγ n have contributions proportional to log ∆/M and do not approach one in the heavy-quark mass limit. All terms in (36) that are proportional to γ d can be viewed as a renormalization of the low-energy parameters d n andd n . This is illustrated in Appendix A and B, where explicit expressions are provided. We note that the fifth order terms can also be readily constructed. For the vector D mesons we derivē where the dots stand for additional terms extracted from (36) with the replacement Π we conclude that for pion masses smaller than ∆ the successive orders (dashed, dotted and dash-dotted lines) approach the exact solid line convincingly. Unlike the consequences of the power-counting ansatz (34) as illustrated in the previous Fig. 1 this is clearly not the case for (36) in the large pion mass domain with m π > ∆.
Neither of the extreme counting assumptions (32) nor (34) generates an expansion scheme that converges for physical up, down and strange quark masses. A step forward may be provided by the following conventional ansatz suggested originallay by Banerjee and collaborators [41,42] for the chiral expansion of baryon masses. Even though the authors demonstrated in a recent work [31] that such an expansion is not suitable to arrive at a meaningful expansion for the baryon octet and decuplet masses at physical values of the up, down and strange quark masses, it deserves a closer study whether it may prove significant for a chiral expansion of the D meson masses. The counting rules (38) lead to somewhat more complicated expressions. Again we derive the third, fourth and fifth order terms. We find andΠ bubble−4  (22,40) are evaluated with the coupling constants g P =g P 0.57 and the physical isospin averaged meson masses. A decomposition according to (30,31) and (38) is performed.
with ∆ Q of (38  (38) is not suitable for a chiral extrapolation of the D meson masses. We note that (38) neither reproduces the results of (32) nor those of (34). We further demonstrate our claim by a plot of the loop functionΠ H in the flavour limit with m π = m K = m η as was done in the previous figures 1 and 2. Fig. 3 demonstrates that for m π > ∆ no quantitative reproduction of the solid line is obtained.
We finally present our counting ansatz that is expected to be applicable from small to medium size quark masses uniformly. It is an adaptation of the framework developed recently for the chiral extrapolation of the baryon octet and decuplet masses [31] and implements the driving idea to formulate the expansion coefficients in terms of physical masses. It is supposed to interpolate the two extreme counting rules (32) and (34). The counting rules are where the sign ± is chosen such that the last ratio in (41) vanishes in the chiral limit. The implications of (41) are more difficult to work out. The counting rules (41) as they are necessarily imply which is at odds with the assumption in (34). Therefore we supplement (41) by the request that the implications of (41) are recovered in the chiral regime. This requires a particular summation of terms proportional to (∆/M ) n with n = 1, 2, 3, · · · .
There is yet another issue pointed out in [31]. The chiral expansion of the scalar bubble function is characterized by an alternating feature. We recall from [31] the following approximation hierarchy where we denoted x = m Q /M H and M R = M H . As was discussed in [40] the terms with even and odd powers in x have opposite signs always. This implies a systematic cancellation effect amongst terms proportional to x n and x 1+n , where the effect is most striking for n = 1.
Therefore it is useful to always group such terms together. Even though the need of such an reorganization is not very strong for the D meson systems under consideration we adapt this strategy in the following. Note that the convergence domain of (43) was proven to be limited by |x| < 2 only, a surprisingly large convergence circle. Given this scheme accurate results can be obtained by a few leading order terms. We construct the third order contributions from the one-loop diagrams.
with ∆ Q and ∆ H as introduced in (41 (22) are evaluated with the coupling constants g P =g P 0.57 and the physical isospin averaged meson masses. A decomposition according to (30,31) and (41) is performed. This leads to (44,45,58,63).
with ∆ Q and ∆ H already introduced in (44).
In Tab  to the previous figures 1, 2 and 3 and demonstrate that for any reasonable pion mass, say 0 ≤ m π < 600 MeV, a quantitative reproduction of the solid line is obtained. We conclude that it is justified to identify the full loop expressions as the loop function to be used at chiral order Q 4 without any significant error from the incomplete 5th order terms.

VI. FIT TO QCD LATTICE DATA
In this section we will determine the low-energy constants c i and d i of the chiral Lagrangian from lattice QCD simulations of the D meson masses. Open-charm mesons have been extensively studied on different QCD lattices [16,[43][44][45][46][47][48][49]. For a recent review we refer to [50]). There exits a significant data set for D-meson masses at various unphysical quark masses. We consider data sets where the pion and kaon masses are smaller than about 600 MeV only. Once we determined the LECs in our mass formula, the D-meson masses can be computed at any values for the up, down and strange quark masses, sufficiently small as to justify the application of the chiral extrapolation.
While for instance in [14,16]  A comprehensive published data set is from Mohler and Woloshyn [14] based on the PACS-CS ensembles [13]. The Fermilab approach is employed in implementing the valence charm-quark [51,52]. In this approach, heavy-quark mass dependent counter terms are added in the heavy-quark action to systematically reduce discretization effects. The valence charm-quark mass dependence is parameterized by a hopping parameter κ c , which is tuned to match the average of the physical kinematic D-meson masses. In Tab and [13]. Statistical errors are given only. The results are based on ensembles from PACS for which their estimate of the lattice spacing is a = 0.0907 (13) fm.
available unpublished results, which allow us to independently extrapolate their lattice data to the physical charm quark mass. For each ensemble, the four D-meson masses but also the η c and J/Ψ masses are computed at two different values of the charm valence-quark mass µ c . As a consequence of the discretization procedure there are corresponding pairs of meson masses that turn degenerate in the continuum limit. We use the notation (±, ∓) and (±, ±) from [16,17]. In this work we focus on the (±, ∓) states and use the masses of the partner states (±, ±) only as a rough estimate for the size of the discretization error. In the vicinity of the physical charm quark mass a linear behavior is expected to hold for all hadron masses. Since the chosen charm quark masses are close to the physical one the ansatz (46) should be justified to sufficient accuracy. The parameters as recalled in Tab. IX is assumed. A typical example for this procedure is shown in Fig. 5 where a sizable uncertainty for the extracted value of µ c is observed.
How such an uncertainty propagates into the masses of the D mesons is shown in Tab. splittings of the two modes, (±, ∓) and (±, ±), for the latter.
It is immediate from Tab. XI and Tab. X that the D meson masses are quite sensitive to the precise charm quark mass used but also to the lattice scale a assumed. We note that, for instance, there exist two distinct values for the lattice spacing for the coarsest ensembles: The value a = 0.0885(36) fm obtained from the pion decay constant [54] and a = 0.0920 (21) fm obtained from the nucleon mass [53]. We conclude that it may be of advantage to determine the lattice scale and the charm-quark mass from the D meson masses directly.
Such a procedure is expected to minimize the discretization errors for the D meson masses. This is what we will do in the following. All information required for such a strategy is provided with Tab. XI and Tab. X, from which the parameters α H and β H in (46) can be read off.
There are yet three further sources of QCD lattice data on the D meson masses, which we will discuss briefly [12,15,18]. The two data sources [12,15] [55,58]. The results are recalled from [12] in units of the lattice spacing a. The lattice spacing is a = 0.1241 (25) fm.
results of [12] rely on previous studies by the LHP collaboration [55], who use a mixed action framework with domain-wall valence quarks but staggered sea-quark ensembles generated by MILC [56][57][58][59]. For the charm quark they use a relativistic heavy-quark action motivated by the Fermilab approach [51,52]. In Tab. XII we summarize the relevant masses that are considered in our study.
The results of the HPQCD Collaboration [15] are based on a highly improved staggered quark (HISQ) action. The HISQ action has since been used very successfully in simulations involving the charm quark such as for charmonium and for D and D s meson decay constants. In Tab. XIII we collect the relevant masses in units of the lattice spacing for the configurations on three coarse and two fine lattices.

Most recently the Hadron Spectrum Collaboration (HSC) computed the excited open-
charm meson spectrum in a finite QCD box [18,19]. Results for the for the D mesons masses based on an ensemble with a pion mass of about 390 MeV are published in [19] and recalled in Tab. XIV. For an additional ensemble at smaller pion masses studies are on going [18].
We note that the charm-quark mass in [12], [15] and [19] was not adjusted to the D meson masses. While in [12] the spin average of the physical J/Ψ and η c meson mass was used, in [15] the charm quark mass was tuned to the physical η c mass. In both cases we cannot exclude uncertainties significant to our analysis. In order to minimize any bias from a possibly imprecise charm-quark mass determination we consider only mass differences from Tab. XII, Tab. XIII and Tab. XIV in our fits. In addition we fine tune the lattice scales.
As we have seen in case of the ETMC results such a procedure reduces any possible bias significantly. taken from [15] and [60]. The lattice spacing is a = 0.119(2) fm and a = 0.0846(7) fm for the two sets of data respectively.
We introduce a universal parameter ∆ c of the form which is supposed to fine tune the choice of the charm quark mass. In principle the values of H depends on the type of D meson considered but also the β QCD value of the ensemble considered. argue that a precise determination of ∆ c and therewith the physical charm quark mass for a given ensemble requires the quantitative control of the chiral extrapolation formulae for the D meson masses.
We do not implement discretization effects in our chiral extrapolation approach since this would introduce a significant number of further unknown parameters into the game. For each lattice group such effects have to be worked out in the context of our chiral extrapolation scheme. As a consequence a fully systematic error analysis is not possible yet in our present study. Here we follow the strategy suggested in [30,31] where the statistical error given by the lattice groups is supplemented by a systematic error in mean quadrature. We perform fits at different ad-hoc values for the systematic error. Once this error is sufficiently large the χ 2 per data point should be close to one. In our current studies we arrive at the estimate  of 5-10 MeV. In anticipation of our analysis of the lattice data set we collect the result of four representative fits. Their characteristics and defining assumptions will be discussed in more detail in the next sections.
Our fit procedure goes as follows. For a given lattice ensemble we take the pion and kaon masses as given in lattice units and then determine from the one-loop expressions (28) in [31] the quark masses for that ensemble. They depend on the three particular linear combinations of the low-energy constants of Gasser and Leutwyler [61]. One combination can be fixed by the request that the η meson mass is reproduced at physical quark masses. The other two are determined by our fit to lattice data. This is analogous to [31] where those low-energy constants are determined from a fit the lattice data on baryon masses. In Tab. XV we show our results for four distinct fit scenarios, which are reasonably close to the results of [31]. In Tab  shown with either green, blue or red filled symbols. In case that for a considered lattice ensemble there is no lattice result for the considered D meson mass available our theory prediction is presented with a yellow filled symbol.
In Fig. 6 we scrutinize the lattice results of [12,14]. XV is considered. From Fig. 6 we conclude that all masses from [12,14] are recovered well is currently computing various scattering observables. We will return to this issue below. are available yet. Note that we show the HPQCD data in units of their spatial lattice spacing but the HSC data in units of 3.5 times their temporal lattice spacing.

VII. LOW-ENERGY CONSTANTS FROM QCD LATTICE DATA
We report on our efforts to adjust the low-energy parameters to the D meson masses as evaluated by the various lattice groups. Our first observation is that the available data set is not able to determine a unique parameter set without additional constraints. Therefore it would be highly desirable to evaluate the D meson masses with J P = 0 − and J P = 1 − quantum numbers on further QCD lattice ensembles with unphysical pion and kaon masses.
Typically solutions can be found with similar quality in the lattice data reproduction but quite different values for the low-energy parameters. This problem is amplified by the unknown size of the underlying systematic error from discretization effects. Almost always the size of the statistical errors given by the lattice groups is negligible, and it is expected that the systematic error is dominating the total error budget. In turn it is unclear whether a parameter set with a better χ 2 value is more realistic than a solution with a worse χ 2 .
The D meson masses may be over fitted.
To actually perform the fits is a computational challenge. For any set of the low-energy parameters four coupled non-linear equations are to be solved on each lattice ensemble considered. We apply the evolutionary algorithm of GENEVA 1.9.0-GSI [62] with runs of a population size 4000 on 100 parallel CPU cores.
In Tab. XVI we collect four distinct fit scenarios which are constrained by additional input from first lattice results on some scattering observable. All sets reproduce the D meson masses with a χ 2 /N close to one given an estimate for the systematic error in the range 5-10 MeV. In all fit scenarios the four low-energy constants c 0,1 andc 0,1 are adjusted to recover the isospin averaged physical D meson masses with J P = 0 − and J P = 1 − quantum numbers from the PDG [63]. This implies that deviations from leading order large-N c or heavy-quark symmetry sum rules are considered for c 0,1 andc 0,1 . In turn we must not impose the heavy quark-symmetry relations d n =d n for all n = 1, ...
the remaining scenarios Fit 2 and Fit 4 keep those parameters unrelated.   The quality with which the four scenarios reproduce the D meson masses from the lattice ensembles is summarized in Tab. XVII. All low-energy parameters are in qualitative agreement with the first rough estimates in (7). On the other hand we find significant tension with the low-energy parameters as obtained in [9,33,64,65]. The parameters of Fit 2 are reasonably close to the two sets claimed in [64] with the notable exception of c 1 which differs  by about a factor 2. Despite the considerable variations in the low-energy constants we deem all four parameter sets acceptable from the perspective of describing the D meson masses.
We repeat that it is unclear whether Fit 1 should be trusted more, only because it would be compatible with a discretization error slightly smaller than the one for Fit 4. After all a 5 MeV systematic error would be an astonishingly small value.
We take up the additional constraints considered. In [12] a set of s-wave pion and kaon scattering lengths was computed on 4 different lattice ensembles as recalled in Tab. XII. Since only for the first two ensembles the kaon mass is smaller than our cutoff choice of 600 MeV, we include into our χ 2 function only the scattering lengths from the first two ensembles of that table. The scattering lengths are computed in the infinite volume limit based on the parameter sets collected in Tab. XVI. We apply the coupled-channel framework pioneered in [4][5][6] which is based on the flavour SU(3) chiral Lagrangian. It relies on the on-shell reduction scheme developed in [38] which can be justified if the interaction is of short range nature or the long-range part is negligible [66,67]. Fortunately this appears to be the case for the s-wave interactions of the Goldstone bosons off any of the D mesons. An alternative chain of works based on a somewhat different treatment of the coupled-channel effects are [9-12, 33, 65]. We did a careful comparison of the three available sources for the flavour structure of the coupled-channel interaction [5,9,12]. We find two discrepancies amongst the original work [5] and [12] where we do take into account the different phase conventions used in the two works for the isospin states. The two discrepancies are in the (I, S) = (1/2, 0) sector. One is traced as a misprint, in C W T of Tab. 2 of [5], in which the two entries 13 and 22 need to be interchanged (see [4]). The second one we attribute to a misprint in [12]. Unfortunately, we were not able to relate to the flavour coefficients shown in [9]. As compared to [5] and [12] there are more than 10 unresolved contradictions.
In Tab. XVIII we collect the χ 2 /N values that characterize how well we reproduce the s-wave scattering length of [12] in our four fit scenarios. Note that we use here our estimates for the lattice scales a LHPC as shown in Tab. XV. The table is complemented by Fig. 9 where a direct comparison of our results with the lattice data is provided for prediction. Our systematic error estimate is implied by a variation of the matching scale µ M around its natural value [4][5][6]. The error bars are implied by ∆µ M = ±100 MeV with It is important to recall that dialing the matching scale slightly off its natural value does not affect our self consistent determination of the D meson masses. The latter is a convenient tool to estimate the uncertainties of the unitarization process.
In the upper panels of Fig. 9 we show the channels that are dominated by a repulsive Tomozawa-Weinberg interaction term [4]. In terms of a flavour SU(3) multiplet classification they belong to a flavour 15plet, that can not be reached within the traditional quark-model picture. A minimal four quark state configuration is required. In contrast in the lower panels, channels are presented that belong to the exotic flavour sixtet sector in which the leading Tomozawa-Weinberg interaction shows a weak attraction [4]. As pointed out in [2,[4][5][6] depending on the size of chiral correction terms exotic resonance states may be formed by the chiral dynamics. Final state interactions distort the driving leading oder term and ultimately generate the more complicated quark mass dependence as seen in the figure. We discriminate results based on ensembles with a kaon mass larger or smaller than 600 MeV by distinct colored symbols. With red symbols we indicate that the kaon mass is larger than our cutoff value, and therefore chiral dynamics is not expected to be reliable.
A fair reproduction of all relevant scattering lengths is seen in (9). Our predictions for the scattering lengths at the physical point is also included by the additional yellow filled points farthest to the left.
We would conclude that with the constraints set by scattering lengths of [12] we cannot rule out any of our four fit scenarios in Tab. XVI.

VIII. SCATTERING PHASE SHIFTS FROM QCD LATTICE DATA
In this section we finally present an additional constraint on the low-energy parameters that provide a clear criterion which of the four fit scenarios is most reliable and should be used in applications. Recently HSC computed πD phase shifts in both isospin channels.
The results are based on the ensemble recalled in Tab. XIV. Given our four parameter sets we can compute those observable at the given unphysical pion and kaon masses. We do this for all four parameter sets.
It is necessary to explain how we compare with those lattice results. Ultimately one should compute the various discrete levels the collaboration computed and than apply the Lüscher method [68,69] to extract the coupled-channel scattering amplitudes. This requires an ansatz for the form of the reaction amplitudes. In the case of a single channel problem this can be analyzed in a model independent manner. In turn for πD scattering in the I = 3/2 channel we can compare our results with the single energy phase shifts as taken from Fig. 20 of [19]. They are to be compared with the four lines from our four fit scenarios.
In the figure of Tab. XIX we see that the two red lines are significantly off the lattice data points, where with those Fit 1 and 2 are presented. This is the case even though in Fit 2 an attempt was made to reproduce the πD phase shifts from [19]. Note that in Fit 1 we ignored any of the latter. We assure that our conclusions are stable against a reasonable variation of the matching scale in this sector.
Based on this observation we made our ansatz for the scattering amplitudes more quantitative by the consideration of an additional set of low-energy constants relevant at chiral order three. Such terms were constructed in [70,71] to take the form Our motivation to consider such terms is slightly distinct to the one followed in [70,71].
From the previous work [6] we expect the light vector meson degrees of freedom to play a crucial role for the considered physics. Ultimately we would like to consider them as active degrees of freedom. This is beyond the scope of the current work. Here we consider the low-energy constants as a phenomenological tool to more accurately integrate out the light vector meson degrees of freedom. In scenario Fit 3 and Fit 4 the contributions of the g n are  worked into the coupled-channel interaction. Their values are displayed in Tab. XIX, which consecutively lead to a significantly improved reproduction of the scattering phase shift.
We proceed by the coupled-channel πD system with I = 1/2 for which its determination of the three phase shifts and in-elasticities is more involved. Some model dependence may enter the analysis. In [19] an estimate of the latter was accessed by allowing a quite large set of different forms of the ansatz for the coupled-channel amplitudes. That then leads to two error bands in their plotted phase shifts and in-elasticity parameters. The smaller one shows the statistical uncertainty, the larger one includes also the systematic error. In Fig. 9 and Fig. 10 of [19] it is shown in addition, on how many levels their results are based on in a given energy bin. Above the πD and below the KD s threshold there are three clusters of levels. We take their center and translate those into single energy phase shifts and in-elasticities with error bars taken from the estimated uncertainties. In Fig. 10 those 'lattice data' points are shown and confronted with our results from the four fit scenarios.
In addition a fourth lattice data point at energies above the KD s threshold is also included in the figure, but shown in red symbols. We do have some reservation towards those points, since the number of close-by energy levels is quite scarce. This is particularly troublesome since here it is a true three channel system that would need more rather than a fewer number of levels to unambiguously determine the scattering amplitude. In turn, the particular choice of ansatz is expected to play a much more significant role in the determination of the red lattice data points. We conclude that the error bars must be significantly underestimated for those points. There is a further piece of information provided by HSC in the given ensemble. The mass M B of a bound state just below the πD threshold is predicted. It is a member of the conventional flavour anti-triplet, which formation was predicted by chiral dynamics unambiguously [4,5]. Within the given error it is not distinguishable from the πD threshold value. The following bound is derived from data published by HSC at the one sigma level. We compute this value in the four fit scenarios with where we find discrepancies for the bound state mass of the order of our resolution of 5-10 MeV. As a consistency check we exploit the uncertainties in the unitarization process, by tuning the matching scale to meet the condition (50)  In the following we take our best fit scenario Fit 4 and provide a thorough documentation of its consequences. In Fig. 11 and Fig. 12 all phase shift and in-elasticity parameters are shown for all possible combinations of (I, S). In Fig. 11 we present the channels in which no exotic signals are expected. Indeed, the evolution from the two HSC ensembles of Tab passes through 90 o in between the ηD andKD s thresholds. We checked that the amplitudes ηD → ηD but alsoKD s →KD s show a well defined resonance structure, with a width significantly smaller than the 300 − 400 MeV of the flavour anti-triplet partner at lower masses. We find this to be a spectacular confirmation of the leading order prediction of this state advocated since 15 years ago by one of the authors (see [2]). It is amusing to see that the clear signature of this state at the physical point may not be seen at the studied HSC ensemble with unphysically large pion masses. Most exciting is the most recent claim in [72] that this state can be seen in data from LHCb [73,74].

IX. ISOSPIN VIOLATING DECAY OF D s0 (2317) FROM QCD LATTICE DATA
A most striking prediction of chiral dynamics is the formation of the D s0 (2317) as a coupled-channel hadronic molecule with significant components in theKD and ηD s twobody states [4]. At leading order in a chiral expansion the coupled-channel interaction is predicted by the Tomozwa-Weinberg term that is parameterized only by the pion-decay or kaon-decay constants, f π or f K , driven into their chiral flavour SU (3) limit with f π,K → f .
This term dominates the s-wave coupled-channel force of the Goldstone bosons with the pseudo-scalar and vector D mesons. The force is short ranged: it may be visualized in terms of a vector meson t-channel exchange process with properly adjusted coupling constants. In contrast to a widespread confusion in the field there are hadronic molecular states that are not driven by a long-range force as provided by an exchange process involving the pion. The challenge is to control and predict such short range forces.
The original work [4] was taken up by many authors [5-7, 9-11, 65, 72, 75, 76] who confirm this universal picture. The challenge is to make this approach more quantitative by controlling chiral correction terms. A first attempt was made in [5,6] based on rough assumptions on the πD invariant mass distributions. A more sophisticated approach was pursued in [10,12] where first QCD lattice data on some s-wave scattering lengths were used. With the significantly improved and extended lattice data set the determination of the low-energy constants, as achieved in our work, is expected to be more controlled and reliable.
In this section we focus on a particular property of the D s0 (2317), its isospin violating hadronic decay width. Since its mass is below the KD threshold and it carries isospin zero it can decay into the πD channel only via isospin violating processes. Estimates of that width within typical quark-model approaches predict such a width of less than 10 keV [77]. This is contrasted by estimates from chiral-coupled channel approaches. Here, already the leading order Tomozawa-Weinberg predicts a width of about 75 keV as demonstrated first in [6]. A corresponding computation with similar physics input but less stringent framework arrived at a similar value [75]. This is to be compared to the significantly larger values of about 140 keV in [6] and later with even an error estimate of (133 ± 22 ) keV [12]. The latter two works implemented chiral correction terms, where the more sophisticated approach [12] was based on additional constraints from some early lattice data. The results of our study for the decay width is collected in Tab. XX for all four fit scenario. Since the mass of the D s0 (2317) was not tuned in any of our fits we again use the uncertainty in the unitarization and adjust the matching scale as to recover the precise mass of the D s0 (2317). This is achieved with 50 MeV < ∆µ M < 90 MeV in the four scenarios. Beside the low-energy constants determined in our work the computation of the width parameter depends crucially on the mixing angle of the π 0 − η system. According to [61] it is determined by the quark masses as follows sin (2 ) cos ( While in [6] the value = 0.010(1) was taken from [61] an updated estimate = 0.0129(7) was used in [12]. Here we consider the impact of a recent and more precise lattice determination of the quark masses by ETMC [54]. This leads to a significantly lower estimate = 0.0061 (8) which our faithful results in Tab. XX are based on.
Since we argued that the lattice data of HSC rule out Fit 1 and Fit 2, we estimate the isospin violating hadronic width of the D s0 (2317) with (45 − 49) keV, significantly lower than the previous claimed value of (133 ± 22) keV [12].

X. SUMMARY AND CONCLUSIONS
We studied the chiral extrapolation of charmed meson masses based on the three-flavour chiral Lagrangian formulated with pseudo-scalar and vector charmed fields. Here the recent approach by the authors constructed for the chiral extrapolation of the baryon ground state masses was adapted to the charm sector successfully, where good convergence properties for the chiral extrapolation are observed. Within the framework the chiral expansion is formulated in terms of physical masses. All D meson masses arise in a manifest scale invariant manner. The framework was applied to lattice data such that an almost unique set of low-energy constants was established.
The low-energy parameters were adjusted to QCD lattice data at N 3 LO, where large-N c sum rules or relations that follow in the heavy charm-quark mass limit were used systematically. We considered lattice data based on ensembles of PACS-CS, MILC, HPQCD, ETMC and HSC with pion and kaon masses smaller than 600 MeV. Besides taking into account constraints from the D meson masses from the various lattice groups, we also considered first results on scattering observables in particular from HSC. Only with the latter, in particular their estimate of the ηD phase shift, we arrive at a rather well defined parameter set, in terms of which we make predictions. The data set on the D meson masses together with constraints from s-wave scattering lengths is not sufficient to nail down the set of low-energy constants.
We computed 15 phase shifts and in-elasticities at physical quark masses but also for an additional HSC ensemble. Such results can be scrutinized by lattice QCD with available computing resources and technology. In addition we predict the isospin violating strong decay width of the D s0 (2317) to be (45 − 49) keV. Given our favorite set of low-energy parameters we find a clear signal for a member of the exotic flavour sextet states in the ηD channel, below theKD s threshold.
To further substantiate the claimed chiral low-energy parameters it is necessary to take additional data on QCD lattices in particular at unphysical quark masses. Our predictions are relevant for the PANDA experiment at FAIR, where the width of the D s0 (2317) may be accessible by a scan experiment [78]. Also the invariant ηD mass distribution, in which we expect a signal from an exotic flavour sextet state, may be accessed by the efficient detection of neutral particles with the available calorimeter.  table are provided to us by the authors of [16].