Contribution of constituent quark model $c\bar{s}$ states to the dynamics of the $D^\ast_{s0}(2317)$ and $D_{s1}(2460)$ resonances

The masses of the $D^\ast_{s0}(2317)$ and $D_{s1}(2460)$ resonances lie below the $DK$ and $D^\ast K$ thresholds respectively, which contradicts the predictions of naive quark models and points out to non-negligible effects of the $D^{(\ast)}K$ loops in the dynamics of the even-parity scalar ($J^\pi=0^+$) and axial-vector ($J^\pi=1^+$) $c\bar s$ systems. Recent lattice QCD studies, incorporating the effects of the $D^{(\ast)}K$ channels, analyzed these spin-parity sectors and correctly described the $D^\ast_{s0}(2317)-D_{s1}(2460)$ mass splitting. Motivated by such works, we study the structure of the $D_{s0}^\ast(2317)$ and $D_{s1}(2460)$ resonances in the framework of an effective field theory consistent with heavy quark spin symmetry, and that incorporates the interplay between $D^{(\ast)}K$ meson-meson degrees of freedom and bare P-wave $c\bar s$ states predicted by constituent quark models. We extend the scheme to finite volumes and fit the strength of the coupling between both types of degrees of freedom to the available lattice levels, which we successfully describe. We finally estimate the size of the $D^{(\ast)}K$ two-meson components in the $D^\ast_{s0}(2317)$ and $D_{s1}(2460)$ resonances, and we conclude that these states have a predominantly hadronic-molecular structure, and that it should not be tried to accommodate these mesons within $c\bar{s}$ constituent quark model patterns.

The presence of the heavy charm quark in the D s states implies the validity, up to Λ QCD /m Q corrections, of Heavy Quark Spin Symmetry (HQSS) [16][17][18][19][20], with m Q the heavy quark mass (the charm quark mass in this case), and Λ QCD a typical scale related to the dynamics of the light degrees of freedom. Thus, in good approximation, the spin of the heavy quark s Q is decoupled from the total angular momentum of the light degrees of freedom jq, and hence they are separately conserved. This gives rise to the arrangement of charmedstrange mesons in doublets, classified by the total angular momentum and parity, j π q , of their light degrees of freedom content, and with total spin J = jq ± 1/2 and parity π. For the P-wave D s mesons, the expected HQSS doublets are, 1 on the one hand, j π q = 1 2 + with J π = 0 + , 1 + mesons, which in S-wave couple to DK and D * K, respectively; and, on the other hand, j + q = 3 2 + with J π = 1 + , 2 + mesons, where the 1 + (2 + ) meson can couple to D * K (DK and D * K) in D-wave. Axial 1 + states will also couple to D * K states in D-wave. In addition, D ( * ) K * pairs in S and D waves might also couple to the latter D s states, but however the dynamics of these pairs will not be governed by chiral symmetry since light-vector mesons are involved. Experimentally, the positive parity jq = 1 2 doublet would be composed of the D * s0 (2317) and D s1 (2460) mesons, which were expected to be almost degenerated and broad, decaying to D ( * ) K through S-wave. However, neither of these two properties are empirically ob-served. This caused an intense debate on the nature of the resonances, producing a wide variety of interpretations among which we highlight their assignment to cs states [21][22][23] and two-meson or four-quark systems [24][25][26][27][28][29][30][31][32][33][34][35].
The latest lattice QCD (LQCD) simulations [36][37][38] have achieved a good description of these charmedstrange resonances when DK interpolators are included in the set of used operators. Notably, the mass of the D * s0 (2317) was found to be overestimated if the DK interpolators were omitted, which gives further support to the idea of a necessary interplay between constituent quark model (CQM) configurations and nearby D ( * ) K thresholds.
On the other hand, if one accepts the predictions of generally successful CQMs, one should expect the charmed-strange J π = 0 + ground state to lie much closer to the DK threshold than the physical D * s0 (2317), so the latter meson pair could act as an essential dynamical agent to reduce the mass of the bare meson state closer to the experimental value, as suggested by some authors [39]. Hence, in this picture, the physical D * s0 (2317) resonance would be the result of a strong renormalization of a bare cs component, rather than a new dynamical state generated from a strongly attractive DK interaction. Nevertheless, since the required renormalization would be quite significant, one would expect, even in this context, that the D * s0 (2317) resonance will acquire a sizable two-meson molecular probability. Indeed, the low-lying P-wave charmed-strange mesons were studied in Ref. [40] employing a widely used CQM [41][42][43], where the coupling between the quark-antiquark and meson-meson degrees of freedom was modeled with the 3 P 0 transition operator [44]. In that work, where all the parameters were constrained from previous studies on hadron phenomenology, the bare 1 3 P 0 cs state 2 developed a large mass-shift as a consequence of its coupling with the DK−meson pair, becoming its mass closer to that of the physical D * s0 (2317) resonance. On the other hand, the dressed state contained a large molecular component that gave rise to a DK−meson pair probability of around 33% [40] in the final configuration of the meson.
In sharp contrast, the LQCD energy-levels reported in Refs. [37,38] were analyzed in Ref. [45], employing an auxiliary potential method, where DK molecular probabilities for the D * s0 (2317) much higher, of the order of 70%, were found. This result was consistent with the previous values obtained in Ref. [46]. The authors of this work performed a LQCD calculation of several heavy-light meson-Golstone boson scattering lengths, that they used to fit the LECs entering in the unitarized next-to-leading (NLO) heavy meson chiral perturbation theory (HMChPT) coupled-channel T −matrix derived in Ref. [47]. The latter amplitudes were employed to estimate the D * s0 (2317) molecular component. These high values for the DK probabilities 3 are similar to those obtained in Ref. [49] from the analysis of the experimental DK invariant mass spectra of the reactions B + →D 0 D 0 K + , B 0 → D − D 0 K + [51] and B s → π +D0 K − [52] measured by the BaBar and LHCb Collaborations, respectively. In all cases an enhancement right above the threshold is seen and it is related in Ref. [49] to the presence of the D * s0 (2317) The latter is dynamically generated when the leading order (LO) HMChPT amplitudes are used as the kernel of a Bethe-Salpeter equation (BSE), whose renormalized solutions fulfill exact elastic unitarity in coupled-channels.
The predominantly molecular structure of the D * s0 (2317) and D s1 (2460) resonances has recently received an indirect robust theoretical support in Refs. [53,54] (see also the review of Ref. [55]). In the first of these two references, the heavy-light pseudoscalar meson J π = 0 + scattering in the strangeness-isospin (S, I) = (0, 1/2) sector was studied, and a strong case for the existence of two poles in the D * 0 (2400) energy-region was presented. The affirmative evidence came from a remarkably good agreement between a parameter-free predicted low-lying levels calculated using NLO HMChPT unitarized amplitudes [47] and the LQCD results reported in Ref. [56]. The dynamical origin of this two-pole structure was elucidated from the light-flavor SU(3) structure of the interaction, and it was found that the lower pole would be the SU(3) partner of the D * s0 (2317). Thus, this latter state will have also clear hadron-molecular origin. A similar pattern was found for J π = 1 + and in the bottom sector. This in fact might solve a long-standing puzzle in charm-meson spectroscopy, since it would provide an explanation of why the masses, quoted in the PDG [57], of the non-strange mesons D * 0 (2400) and D 1 (2430) are almost equal to or even higher than their strange siblings. In the second reference [54], it is shown that the wellconstrained amplitudes for Goldstone bosons scattering off charm mesons used in Ref. [53] are fully consistent with recent high quality data on the B − → D + π − π − and B 0 s →D 0 K − π + final states provided by the LHCb experiment in Refs. [58] and [52], respectively. Indeed, all these results suggest a new paradigm for heavy-light meson spectroscopy that questions their traditional qq CQM interpretation [54].
Most of these latter works [46,49,53,54] do not incorporate explicitly the bare cs degrees of freedom, whose effects are, in principle, encoded in the low energy constants (LECs) that appear beyond LO in the chiral expansion, and in the non-perturbative resummation employed to restore elastic coupled-channels unitarity. However, and depending on the proximity of the CQM states to the energy-region under study, this approximation might not be sufficiently accurate.
Such radically different pictures of the inner structure of the D * s0 (2317) makes timely a re-analysis of this resonance, paying special attention to the interplay between meson molecular and CQM degrees of freedom. We will employ here the scheme designed in Ref. [59] for its bottom heavy-flavor partner, and we will try to describe the charm 0 + and 1 + LQCD energylevels obtained in Ref. [36]. Such comparison will also serve to further 4 constrain/determine the LECs that appear in the approach. The scheme of Ref. [59] started from a unitary ansazt for the heavy-light-meson-Goldstone-boson 0 + and 1 + amplitudes, based on LO HMChPTB ( * ) K interactions, and computed for finite volumes. In addition, and for the very first time in this context, the two-meson channels were coupled to the corresponding CQM P-waveB s scalar and axial mesons using an effective interaction consistent with HQSS. In the m Q → ∞ limit and besides HQSS, Quantum Chromodynamics (QCD) acquires also an approximate heavy flavour symmetry that ensures that the dynamics of the system containing a single heavy quark is, up to O(Λ QCD /m Q ) corrections, independent of the flavor of the heavy quark [20]. As a consequence, the bottomstrange and charm-strange systems are expected to have a similar behavior. Hence, the analysis of the B s lowlying spectrum carried out in Ref. [59], can be readily used here to study the charmed-strange j π q = 1 2 + HQSS doublet.
The recent LQCD simulation of Ref. [36] reported finite volume energy-levels from a high statistics study of the J π = 0 + and 1 + charmed-strange mesons, D * s0 (2317) and D s1 (2460), respectively, where the effects of the nearby DK and D * K thresholds were taken into account by employing the corresponding four-quark operators. As we will discuss below, the work of Ref. [36] represents a clear improvement on the pioneering ones of Refs. [37,38]. Some of the energy-levels reported in Ref. [36] lie in the proximity of, when not above, the expected CQM bare masses of the ground scalar and axial charmedstrange states, being thus interesting to include explicitly these degrees of freedom in the scheme, since their effects might not be properly taken into account by simply including LECs. Moreover, in the present study, we will also include the next (S = 1, I = 0) higher thresholds, D ( * ) s η, since they appear at energies below some of the finite-volume levels computed in Ref. [36].
This work is structured as follows. After this Introduction, the theoretical formalism is described in Sec. II, 4 As it will be discussed below, LECs in the charm and bottom sectors might be related by heavy-quark flavor symmetry.
where details about the D ( * ) K and D ( * ) s η scattering amplitudes, the coupling of the meson-pair degrees of freedom with the bare CQM cs spectrum, the restoration of unitarity and the extension of the scheme to finite volume are discussed. In Sec. III our results for the finite volume 0 + and 1 + energy levels are presented and compared with the ones reported in Ref. [36]. The properties of the D * s0 (2317) and D s1 (2460) mesons, in particular their molecular content, are discussed in this section. We also compute the energy levels obtained using the unitarized NLO HMChPT amplitudes derived in Refs. [46,47], and extensively compare the predictions of this latter scheme with those deduced by including a bare CQM pole. The section concludes with predictions for DK S-wave scattering phase-shifts, and a discussion about a few aspects of the renormalization dependence of the results obtained in this work. Finally, our conclusions and a summary of results are presented in Sec. IV. As mentioned in the Introduction, we will extend the formalism of Ref. [59] for theB ( * ) K and bs system to the charm sector. We will briefly review here the most relevant aspects, paying attention to the inclusion of the D ( * ) s η channels, whose counterparts were not considered in Ref. [59].

A. Interactions
We are interested in the S-wave Goldstone boson (K + , K 0 and η) scattering off charm mesons P ( * ) a ≡ (D 0( * ) , D +( * ) , D +( * ) s ). We will generically refer as φ to the former mesons and P ( * ) to the latter ones. The heavy-light charm mesons are described in terms of HQSS matrix field H a , with v µ the four-velocity of the heavy mesons, and a a light-flavor SU(3) index. The HQSS field combines the isospin doublet and singlet of pseudoscalar heavy-mesons P (c) a = cū, cd, cs fields and their vector HQSS partners P * (c) a . Transformations of these fields under parity and heavy spin and chiral rotations can be found for instance in Ref. [59]. The Weinberg-Tomozawa Lagrangian (WTL) describes the S-wave P ( * ) φ chiral interaction at LO, and it reads [60][61][62][63], withH a = γ 0 H † a γ 0 the hermitian conjugate of H a , and ξ the Goldstone-boson matrix field given by [61], with normalization f ∼ 93 MeV and the matrix M reads, where we explicitly consider the η 1 and η 8 unflavored SU(3) states. Note that in our normalization the heavylight meson field, H, has dimensions of E 3/2 (see Ref. [20] for details). This is because we use a non-relativistic normalization for the heavy mesons, which differs from the traditional relativistic one by a factor √ M H . On the other hand, within the HQSS formalism, the even parity CQM bare cq states, associated to the j π q = 1 2 + HQSS doublet, are described by the matrix field J a [64], with v µ Y * a µ = 0. The Y a and Y * a fields respectively annihilate the 0 + and 1 + bare states belonging to the 1 2 + doublet. Since in this work we will be interested in the (S, I) = (1, 0) sector, in what follows and for simplicity we will denote Y ( * ) cs as Y ( * ) . The mass of CQM bare states, • m cs , is a renormalization-dependent parameter of the scheme [65] that will be discussed below. At LO in the heavy quark expansion, there exists only one term invariant under Lorentz, parity, chiral and heavy quark spin transformations [59], where c is a dimensionless undetermined LEC that controls the strength of the vertex. This LEC, though it depends on the orbital angular momentum and radial quantum numbers of the CQM state, is in principle independent of the spin of the quark-model heavy-light meson, and of the light SU(3) flavor structure of the vertex. Thus, up to Λ QCD /m Q corrections, it can be used both for J = 0 and J = 1 in the charm and bottom sectors. Moreover, in the SU(3) limit, the same LEC governs the interplay between two-meson and quark model degrees of freedom in all isospin and strangeness channels. This LEC was found to be c = 0.74 ± 0.05 in Ref. [59] from a fit to the LQCD isoscalar bs 0 + and 1 + energy-levels computed in Ref. [66]. Let us consider the transitions involving pseudoscalar (P ) heavy-light and Goldstone (φ) mesons, for which the tree level isoscalar 5 amplitude (V c ) deduced from the WTL of Eq.
(2) reads [67] V c (s, where channels 1 and 2 are DK and D s η, respectively, s and u are the usual Mandelstam variables and we have considered the η − η ideal mixing [68], After projecting into J = 0, we replace where M I(F ) and m I(F ) are the masses of the initial (final) heavy-light charm and Goldstone mesons, respectively.
The WTL leads to similar isoscalar amplitudes for the transitions involving vector (P * ) heavy-light mesons, Indeed, one gets in this case a potential 6 like that of Eq. (8), supplemented by a term − I · F , where I(F ) is the polarization four-vector of the initial (final) heavylight meson. This latter factor reduces, at LO, to 1 after projecting into J = 1 (S-wave, i.e., zero orbital angular momentum). In summary, the amplitudes given in Eq. (8), together with the projection implicit in Eq. (10) provide the coupled-channel contact potential, V c (s), both in the 0 + and 1 + sectors. Next, we consider transitions between CQM bare states [Y ( * ) ] and P ( * ) φ meson pairs, Y ( * ) → P ( * ) φ. From Eq. (6) we find with M an m the masses of the P ( * ) and φ mesons, respectively. 5 The phase convention for isospin states |I, I 3 used in this work isū = |1/2, 1/2 andd = −|1/2, 1/2 , which implies |D + = − 1 2 , + 1 2 , while the other meson states are defined with a positive sign in front of the |I, I 3 state. Moreover, we use the order DK, as in Ref. [67], to construct the isoscalar amplitudes. 6 Now channels 1 and 2 are D * K and D * s η, respectively.
Note that, here, by bare mass, we mean the mass of the CQM states when the LEC c is set to zero, and thus it is not a physical observable. In the sector studied in this work, the coupling to the P ( * ) φ meson pairs renormalizes this bare mass, as we will discuss below. Since, in the effective theory, the ultraviolet (UV) regulator is finite, the difference between the bare and the physical resonance masses is a finite renormalization. This shift depends on the UV regulator since the bare mass itself depends on the renormalization scheme. The value of the bare mass, which is thus a free parameter, can either be indirectly fitted to experimental observations, or obtained from schemes that ignore the coupling to the mesons, such as some CQMs. In this latter case, the issue certainly would be to set the UV regulator to match the quark model and the HMChPT approaches [65].
The vertices in Eq. (12) can be used to compute the contribution [V ex (s)] to P ( * ) φ scattering via the exchange of intermediate even-parity charmed-strange mesons. It is given by [59,65] Finally, the full effective potential, V (s), consistent with HQSS is given by that incorporates the interplay between meson-pairs and CQM degrees of freedom in the P ( * ) φ dynamics [65]. Note that V c (s) is obtained from V c (s, u) in Eq. (8) using the projection implicit in Eq. (10). The LO potentials, V c (s) and V ex (s), are identical in the heavyquark limit in both the 0 + and 1 + sectors. Here, we will include some HQSS breaking corrections stemming from the mass difference between pseudoscalar and vector heavy-light mesons and the masses of scalar and axial CQM bare states. The scheme derived here is similar to that followed in Ref. [59], though here the D ( * ) s η channel has been also included, while the counterpart of this channel in the bottom-strange sector was not considered in the latter work.

B. Unitarity in coupled-channels
Elastic unitarity in coupled channels is restored by solving for each J π sector a BSE, using as kernel the HQSS effective potential of Eq. (14). The BSE is solved within the so-called on-shell approximation [69] and using a Gaussian cut-off, Λ, to regularize/renormalize its UV behaviour. Thus, the unitary scattering amplitude T (s) in coupled channels is obtained from the matrix equation where G(s) is the regularized (Λ) two-meson loop function. Normalizations are fixed thanks to the relation between the S and T −matrices, and the relation of the latter matrix with phase-shifts [δ(s)] and inelasticities [η(s)]. Namely, we use with Σ(s) = diag (σ 1 (s), σ 2 (s)) and the function σ i (s) is defined as: where λ(x, y, z) = x 2 + y 2 + z 2 − 2xy − 2xz − 2yz and (a, b) = (D ( * ) , K) and (D ( * ) s , η) for i = 1 and 2, respectively. Thus, unitarity in coupled-channels can be expressed as The expression for G was given in Eqs. (14) and (15) of Ref. [59] for the B ( * ) K pair, and here in coupled channels it is trivially modified to compute any of the elements of the diagonal matrix G(s) = diag G D ( * ) K (s), G D ( * ) s η (s) .

C. CQM states
We have seen that through the coupling of the CQM cs and the P ( * ) φ degrees of freedom, the effective interaction incorporates a term (V ex (s)) driven by the exchange of bare CQM states (Y ( * ) ). Such a term introduces a pole in the two-meson tree-level amplitudes, Eq. (14), located at the bare mass value, √ s = • m cs . As already mentioned, it should be interpreted as the mass of the CQM state in the limit of vanishing coupling to the P ( * ) φ mesonpairs (c → 0), and therefore it is not an observable. The interaction with the meson cloud dresses the CQM state through loops [Eq. (15)], renormalizing its mass, and the dressed state might also acquire a finite width, when it is located above threshold. At energies far enough • m cs , the contribution of V ex can be regarded as small contact interaction that can be accounted for by means of a LEC. However, the exchange contribution becomes more important for higher energies approaching • m cs , and its energy dependence might then not be safely ignored.
A priori, the value of • m cs is a free parameter of the present approach, and moreover it should depend on the renormalization scheme [65]. We will take predictions from quenched CQMs, which in principle do not include couplings with nearby meson-meson channels. In the J π = 0 + sector, quark models predict, in general, cs bare masses well above M D + m K , which lead to attractive V ex exchange interactions at DK threshold, which might help in forming the D * s0 (2317) resonance.
Most quark models predict masses in the range of 2.45 − 2.51 GeV for the 1 3 P 0 cs state [3,[6][7][8], significantly far from the D * s0 (2317) experimental mass. With the aim of improving these predictions, some other models [5,23] incorporated a one-loop correction diagram to the one-gluon exchange (OGE) potential, adding a spindependent term to the quark-antiquark potential which affects mesons with different flavor quarks, such as the cs mesons. In these works, it is shown that this correction is rather small except for the 0 + sector, where large shifts are found and the mass of the CQM 1 3 P 0 state is significantly lowered (∼ 100 MeV). Indeed, it is found in the 2.35-2.38 GeV region-closer to but still above the experimental D * s0 (2317) mass. Nevertheless, these predicted states would be still above the DK threshold and their width would be large due to the decay into final DK, and difficult to reconcile with the experiment that currently provides an upper bound of a few MeV for its total width, 7 as mentioned in the Introduction. Moreover, quark models including these modified OGE potential will still face difficulties to describe the J π = 1 + sector, where these corrections are quite small, and the experimental mass pattern is similar to that found in the isoscalar-scalar sector.
One should bear also in mind, as the one-loop OGE correction brings the bare state closer to the DK threshold, the interplay between the two-meson channel and the CQM degrees of freedom might have a major impact on the description of the resonance properties and LQCD energy-levels.
Thus, in this study we will explore both types of CQMs. On one side, we will take the • m cs values for the j π q = 1 2 + charmed-strange meson doublet predicted in the CQM calculations of Refs. [43,72] (Set A in Table I). Such CQM is based on the assumption that the light constituent quark mass and the exchange of pseudo-Goldstone bosons arise as a consequence of the spontaneous breaking of the chiral symmetry in QCD. Besides, the dynamics is completed with a perturbative OGE potential and a non-perturbative screened confining interaction [41,43]. On the other hand, another set of bare masses (Set B in Table I) will be employed, predicted within the same CQM but supplemented by the oneloop OGE corrections derived in Ref. [5] (see Ref. [72] for further details). In contrast to the bottom sector studied in Ref. [59], we see in Table I  s η averaged threshold energies (in MeV units) used in this work, taken from Ref. [72]. In addition, we take an average kaon mass of mK = 495.6 MeV.

D. NLO corrections
We will also use the O(p 2 ) (NLO) HMChPT amplitudes derived in Ref. [47], with LECs determined from the lattice calculation in Ref. [46] of the S-wave scattering lengths in several (S, I) sectors. After unitarization, the scheme provides an accurate description of the P φ interactions in coupled channels. For instance, as it is shown in Ref. [53] and already mentioned in the Introduction, the finite volume energy-levels in the (S, I) = (0, 1/2) channel calculated with the unitarized O(p 2 ) amplitudes, without adjusting any parameter, are in an excellent agreement with those recently reported by the Hadron Spectrum Collaboration [56]. These chiral amplitudes predict the existence of two scalar broad resonances, instead of only one, with masses around 2.1 and 2.45 GeV, respectively [53]. The lower pole would form, together with the D * s0 (2317) (also correctly predicted in [53]), a SU(3) 3 multiplet. This scheme provides also an excellent description of the recent LHCb data [58] on the B − → D + π − π − final states and give solution of a number of puzzles [54]. Moreover, these NLO unitarized amplitudes have been used in Ref. [73], as input of coupledchannel Muskhelishvili-Omnès integral equations, whose solutions produced scalar form factors of the semileptonic heavy meson decays D → π¯ ν , D →K¯ ν ,B → π ν andB s → K ν in good agreement with LQCD and lightcone sum rule predictions.
The consideration of the contributions from explicit exchanges of CQM states, in addition to the LO HM-ChPT amplitudes, contemplated in the previous subsections provides a different perspective to the physiognomy of the P ( * ) φ interactions. Thus, it is worth discussing if there exists an energy-regime where the LO & CQM scheme might mimic the NLO amplitudes, and when the CQM effects cannot be properly accounted for by the LECs that appear beyond LO in the chiral expansion. Such study could provide further insights on the energy range of applicability of the unitarized effective theory. Hence, in addition to results obtained with a LO & CQM interaction, we will also show results from the (S, I) = (1, 0) NLO input derived in Refs. [46,47], and will compare both sets of predictions with the LQCD energy-levels reported in Ref. [36].

E. Poles, couplings and the compositeness condition for bound states
The interplay between meson-meson and CQM cs states might dynamically generate new states that arise as poles of the scattering amplitudes on the complex s−plane. There exist two thresholds s where, for future purposes we have also defined s + . In our two-channel problem, Riemann sheets (RSs) are denoted as (ξ 1 ξ 2 ), ξ i = 0, 1, and are defined in the whole complex plane through analytical continuations of the loop functions [74]: where the cuts for σ i (s) go along the real axis for −∞ < s < s The function is chosen to be real and positive on the upper lip of the second cut, s With all these definitions, (00) is the physical RS, while the SRS is defined by requiring continuity across the unitarity cut between its fourth quadrant and the first one of the FRS. Therefore, the definition of the SRS of the amplitudes varies below and above the highest threshold (branch point of the T −matrix) [74], and it corresponds to (10) or (11) when the real part of s is above s (1) + , but below s (2) + , or above both thresholds, respectively (see Ref. [74] for some more details).
The mass and the width of the bound state/resonance can be found from the position of the pole on the complex energy plane. Close to the pole, the T -matrix behaves as The quantity √ s R = M R − i Γ R /2 provides the mass (M R ) and the width (Γ R ) of the state, and g i is the complex coupling of the resonance to the channel i that it is obtained from the residue. The residues can be used to get information on the compositeness of the bound states. Motivated by the Weinberg compositeness condition [48,75,76], the probability of finding the D ( * ) K or D ( * ) s η molecular component in the bound state wave function is given by [50,77], where M b is the bound state mass. The energy dependence of the potential produces probabilities in Eq. (23) that deviate from one. We will restrict the discussion to the case of bound states. The evaluation of Eq. (23) for resonances gives rise to complex values of P j , loosing then a straightforward probabilistic interpretation. (Further details can be found, for instance, in section 4.2 of Ref. [65]. See also Refs. [78,79].) F. Finite volume and details of the simulation of Ref. [36] As mentioned in the Introduction, the simulation of Ref. [36] reported finite volume energy-levels from a high statistics study of the D * s0 (2317) and D s1 (2460) resonances, taking into account effects of the nearby DK and D * K thresholds by employing appropriate fourquark operators. Six ensembles with N f = 2 nonperturbatively O(a) improved clover Wilson sea quarks at lattice-spacing a = 0.071 fm were employed in Ref. [36], covering different spatial volumes and pion masses: linear lattice size (L) of 1.7 fm to 4.5 fm were realized for m π = 290 MeV and 3.4 fm and 4.5 fm for an almost physical pion mass of 150 MeV. Thus, the work of Ref. [36] represents a clear improvement on the pioneering ones of Refs. [37,38], where an ensemble with m π = 156 MeV, at a fairly coarse lattice spacing of a = 0.09 fm and a small spatial lattice extent of L = 2.9 fm (Lm π = 2.29) were analyzed using the effective range approximation to extract infinite volume results. Thus, the use of a finer lattice spacing in Ref. [36] is important since discretization effects can be substantial for observables involving charm quarks, while exploring the dependence on the spatial volume is needed, since contributions which are exponentially suppressed in Lm π (that are ignored in the Lüscher formalism) may not be small for Lm π = 2.29.
To compare with the energy-levels reported in Ref. [36], we consider our scheme, based on unitarized HMChPT and the contribution of CQM states, in a cubic box of side L, and periodic boundary conditions for the fields. The three momentum is quantized q = 2π L n , where ω i ( q ) = M 2 i + q 2 1 2 , with M i the mass of the heavy-light meson or the Goldstone boson for i = P ( * ) or φ, respectively. (We have adopted relativistic dispersion relations as in Ref. [36].) Up to the order considered in this work, there are no finite volume corrections to the potential, so the full volume dependence is carried by the loop function G defined in a finite box. The T -matrix in finite volume, T (s, L), is given by, Thus, the energy-levels E(L) (s = E 2 , E ∈ R) are computed from the poles of T (s, L) for each size of the box. The spectrum becomes discrete, with levels that, in principle, can be associated to two-meson (P ( * ) φ) scattering states. In the non-interacting case, the free energies, E free Hence, the continuous volume dependent curves that will be presented below are essentially the Lüscher curves obtained from the phase shift by solving withq = qL/2π and φ(q) determined by the Lüscher function (see Eq. (6.13) of Ref. [82]). On the other hand, for a proper comparison with the results of Ref. [36], it is necessary to use the lattice meson masses obtained in that simulation. That work reported two different sets of results that correspond to two different pion masses used to compute the energylevels. We will label here the two sets of LQCD levels by Ensembles I and II, for m π = 290 MeV and 150 MeV, respectively. We take the π, K, D and D * masses, for the different box-sizes considered in each of the ensembles, from Table I of Ref. [36]. Besides, the 0 − and 1 − D ( * ) s masses for the Ensemble II are taken from Table VII (L/a = 64) of Ref. [36], while for Ensemble I the masses of the charmed-strange heavy-light mesons have been obtained using the values of (m 0 + − m 0 − ) and Fig. 13 of the latter reference, and taking for the m 0 + and m 1 + masses the L/a = 64 values reported in Table III of the same work (Ref. [36]). Hence, we neglect any dependence on L in the masses of 0 − and 1 − cs ground states. 8 The masses used in this work are compiled in Table II.
The mass of the η meson is not reported in Ref. [36] either, so we estimated its value, as function of the volume and pion mass, using the Okubo mass formula [83] 8 It is not straightforward to extract the masses of the 0 − and 1 − D ( * ) s mesons, for each of the pion masses and volumes studied in Ref. [36], from the results reported in that work. Note, however, that volume effects in the mases of these mesons are expected to be even smaller than for the D−meson, and that these masses enter only through the small effects originated from the coupledchannels dynamics, when the D  as follows with m 2 8 = 4(m 2 K − m 2 π )/3. Finally for unphysical pion masses, the CQM bare masses are also corrected using the difference between the experimental and lattice spinaveraged D s masses, Uncertainties in all input lattice mesons masses (Table  I of Ref. [36] and Table II of this work), as well as in • m lat cs , are taken into account in the error budget of our final results, as we will detail in the next section.
To end this subsection, we would like to stress that considering the LQCD meson masses, for finite volumes and unphysical pions, it is important to set up correctly the thresholds and to properly compute the loop function in a finite box. However, the current approach will still suffer from some systematic errors, mostly because we have not considered the dependence of the Goldstone decay constant, that appear in the WTL interaction, on the volume and the unphysical pion mass (such information is not given in Ref. [36]), and we still use its value in the infinite-volume chiral limit. 9 Nevertheless, some of this dependence might be partially reabsorbed in the parameters fitted to the LQCD energy-levels, and we certainly benefit from the fact that the pions simulated in Ref. [36] are quite light and close to the physical one.

III. RESULTS AND DISCUSSION
A. LO HMChPT+CQM analysis

Fit details
The S-wave D ( * ) K interaction in the LO+CQM interaction scheme depends on three, a priori, free parameters: the Gaussian cut-off Λ, the LEC c (Eq. (12)) and the masses of the bare CQM cs state, • m cs . Nevertheless, we will consider different CQM meson masses in the 0 + and 1 + sectors (Sets A and B in Table I), as discussed in Sec. II C.
To determine the values of the LEC c, that controls the interplay between CQM and meson-pairs degrees of freedom, and the Gaussian cut-off, we perform, for each set (A and B) of bare CQM masses and lattice ensembles (I and II) a combined fit to the 0 + and 1 + energy-levels, using an uncorrelated merit function defined as, where the sum spans over the 0 + and 1 + energy-levels compiled in Table III of Ref. [36]. For both 0 + and 1 + sectors, only two energy-levels for each volume are fitted.
We have not considered in the fits the third 1 + level given in the last column of that table. It is rather insensitive to the spatial volume, suggesting only a small coupling to the D * K threshold and, in Ref. [36], it is identified with the D s1 (2536) resonance, that would presumably have a large overlap with the j π q = 3 2 + HQSS state.
We show in Figs. 1 and 2 the 0 + and 1 + energylevels obtained using sets A and B of CQM bare masses, respectively. Fitted parameters and best fit χ 2 /dof values are collected in Table III.
We see that the Set A of bare CQM masses provides a fairly good description of the volume dependence of the LQCD energy-levels in both the J π = 0 + and 1 + sectors, despite the large deviations from the free levels. There exists a very mild dependence of the UV cutoff and LEC c on the pion mass, which however is not statistically significant. The D ( * ) s η coupled channel effects are negligible, except perhaps for the highest levels calculated with the heaviest pion mass ensemble at the smallest of the volumes, since that threshold is located sufficiently more higher than the measured levels in Ref. [36]. (Note that D ( * ) s η correlators are not considered in the LQCD study of the latter reference.) We find c = 0.61 (9) and Λ = 710(70) MeV from the fit to the lightest pion mass ensemble. This parameter was also determined in Ref. [59] from a similar analysis of the LQCD low-lying 0 + and 1 + B s −energy-levels calculated in Ref. [66]. There, it was found c = 0.75 (6), with an UV cutoff in the range 620-770 MeV, which points out to some small dependence of the LEC c on the heavy-flavor mass.
On the other hand, when the Set B of 0 + and 1 + CQM bare masses are used, we find unacceptable fits, with χ 2 /dof values above 18 for both pion mass ensembles. Indeed, as can be seen in Fig. 2, the set B leads to a really poor description of the LQCD data. For the latter, the HQSS breaking corrections thus look i) compatible with those encoded in the current scheme when the Set A of bare masses is used, but ii) much smaller than those implemented by the Set B of bare masses. The one-loop corrections [5] to the OGE potential implemented in Ref. [72] produce a 1 + − 0 + shift of the bare CQM masses of around 190 MeV, while it is around only 80 MeV when these corrections are neglected (see Table I). This is because, as we already mentioned, this correction mostly affects the 0 + sector [5]. Actually, because of the denominator in Eq. (13) and for fixed c and Λ, the decrease of the 0 + bare mass produces an enhancement of the attraction close to the DK threshold from the exchange of the CQM state. This effect is much less important in the 1 + sector, and thus the current scheme using Set B of bare masses produces a visible tension between the predicted 0 + and 1 + levels and the LQCD data (Fig. 2). The scheme tends to overestimate (underestimate) the attraction in the former (latter) energy-levels by around 10-20 MeV, amount significantly larger than the errors of the LQCD data. Therefore, this study strongly disfavors the bare masses found in Ref. [72], where the one-loop corrections derived in Ref. [5] are taken into account. These oneloop corrections to the OGE potential produced a much smaller HQSS breaking of the bare masses in the bottom sector, and thus nothing statistically meaningful could be concluded about this issue in the analysis carried out in Ref. [59] of the B s −energy-levels reported in Ref. [66]. Due to the former discussion, we will always make reference to the results obtained from the Set A of CQM bare masses in the rest of the work.  Table IV  FIG. 1. Black points: LQCD J π = 0 + (bottom panels) and J π = 1 + (top panels) energy-levels, taken from Ref. [36], for different volumes and pion masses. LQCD data for Ensembles I and II are depicted in the left and right panels, respectively. Solid lines: Energy-levels obtained from 0 + and 1 + combined fits to Ensembles I and II using the Set A of CQM bare masses, as a function of the box-size, L (some volume interpolated meson masses are used to compute the continuous energy-levels for values of L different than those employed in Ref. [36]). Details of the fits and some derived quantities are collected in Table III. Statistical 68%-confident level (CL) bands are also shown. They are calculated from the distributions obtained from a sufficiently large number of fits to synthetic sets of LQCD data, randomly generated assuming that each of the LQCD energy-levels is Gaussian distributed. (Note that possible correlations between the different energy-levels are not considered, since these are not provided in Ref. [36].) In addition, in all synthetic fits, the LQCD meson masses for each volume are randomly chosen as well. For comparison, the free energies E free j , j = D ( * ) K, D ( * ) s η, for each volume and set of LQCD meson masses are also shown (dashed lines).

Properties of the lowest-lying states: masses, molecular probabilities and couplings
We first pay attention to the mesons of the lowest-lying j π q = 1 2 + D s doublet. We see that the predicted mass of the 0 + bound state, using Set A of CQM bare masses, is only around 15 MeV higher than the experimental mass of the D * s0 (2317) and 75-80 MeV smaller than the CQM bare mass, while that of the 1 + nicely agrees, taking into account the errors, with the experimental mass of the D s1 (2460) state. Nevertheless, this level of discrepancy can be well attributed to the presence of discretization errors, or some uncertainties in the determination of the mass of the charm quark in the lattice simulation. Actually, the results compiled in Table I of Ref. [36] show discrepancies of the order of 10 MeV between the experimental masses of the D and D * mesons and the LQCD ones, determined from the lightest quark mass, which provides almost physical pion masses. We should also note some differences (1-2 σ's) between the 0 + and 1 + masses found in this work and those reported in Table   VII of Ref. [36], taken from the m π = 150 MeV and L = 64a ensemble. These might be due to the use here of physical meson masses, and also because the LQCD ones are accessed via the Lüscher's relation [Eq. (27)] using the effective range approximation. The approach followed here, where the two meson loop function is computed in a finite volume, the unknown LECs are determined from fits to the LQCD data, and finally poles are searched for in the infinite volume unitarized chiral amplitudes, provides a theoretically sound tool to analyze the LQCD energy-levels. A good example of the latter affirmation can be found in Ref. [53], where such approach led to the existence of two D * 0 (2400)−poles, instead of only one reported in the original lattice work of the Hadron Spectrum Collaboration [56], where the LQCD energylevels were calculated (see also the discussion in Ref. [54], where it is emphasized how the two-pole pattern of the D * 0 (2400), together with their SU(3) structure, provides a natural solution to a number of puzzles).
Interestingly, we appreciate in Fig. 1 a quite significant dependence of the lowest-lying LQCD energy-levels on the pion mass (differences between left and right plots), as it is also evident in the results of Table IV of Ref. [36]. Thus, the LQCD 0 + and 1 + masses reported in that table vary between 30 to 50 MeV, when the pion mass is reduced from 290 MeV down to 150 MeV. These changes are likely related to the modifications of the DK and D * K thresholds. All of this clearly indicates that the D * s0 (2317) and D s1 (2460) states should have a sizable molecular component, and that any CQM cs component in their dynamics cannot be dominant, because it could not accommodate such visible dependence of their masses on the light quark mass, as exhibited by the LQCD data. These findings are corroborated by the molecular probabilities collected in Table III. Using the modified Weinberg compositeness condition of Eq. (23), we derive the molecular S-wave D ( * ) K probabilities for the D * s0 (2317) and D s1 (2460), which turn out to be around 65% and 56% for the scalar and axial states, respectively. On the other hand, D ( * ) s η components are small for both mesons, of the order of 2%. The LQCD studies of Ref. [36][37][38] showed a non-zero overlap of the energy-levels related to the D * s0 (2317) and D s1 (2460) and meson-meson lattice interpolating fields, but no precise information was provided on the percentage of mesonmeson channels in the wave function of these states. Only in the latest work of Ref. [36], the compositeness probability is studied, and found to be 1 within errors (around 20-30%) for both states.
The large molecular probabilities found in this work confirm those reported in previous works [45,46,49] that, employing also unitarized meson-meson amplitudes, had already obtained molecular components of around 70% for the D * s0 (2317), as mentioned in the Introduction. In what respects the D s1 (2460), the authors of Ref. [45] found a molecular probability of 0.57 ± 0.22 also in good agreement with our findings, although with a much larger error. The interesting and novel aspect of the present calculation is that the LO HMChPT interactions have been supplemented by those driven by the exchange of even-parity charmed-strange CQM mesons, and thus the couplings of CQM cs and P ( * ) φ degrees of freedom have been explicitly taken into account. 10 In Ref. [40] two-meson loops and CQM bare poles are also coupled. For the latter, the values of the bare masses are the same as those used here. The D ( * ) K interactions are derived from the same CQM used to compute the bare states, instead of using HMχPT. The 10 Other studies have done something similar (e.g. Ref. [45]) by including in the interactions a Castillejo-Dalitz-Dyson pole [85] to visualize a genuine (CQM) state that couples weakly to a meson-meson component. However, those studies do not make use of the HQSS to relate the interplay between both types of degrees of freedom in the 0 + and 1 + sectors, which will be fundamental to disfavor the Set B of CQM bare masses used in the study of Ref. [40], that will be discussed in the next paragraph, and that claimed a much smaller (∼ 30%) molecular components for the D * s0 (2317). Moreover, in some of these studies chiral symmetry is not fully used to constraint the P ( * ) φ interactions, and different solutions were obtained with many sets of parameters, obviously correlated, though the claim in Ref. [45] was that the particular values of the parameters did not have a special significance, and all of them led to similar hadronic-molecular probabilities [45]. III. Best fit LO+CQM parameters, together with infinite volume properties (masses, D ( * ) K and D ( * ) s η molecular components and couplings) of the lowest-lying j π q = 1 2 + Ds charm-strange meson doublet, determined from the fits to the lattice energy-levels obtained for each of the two lattice pion mass ensembles (I and II) calculated in Ref. [36], and using either Set A or B of bare CQM masses (see discussion in Sec. II C). S-wave isoscalar D ( * ) K scattering lengths (a) are also given, which are related to the amplitudes at threshold by T s = (M D ( * ) + mK ) 2 = −8πa (M D ( * ) + mK ), as in Ref. [36]. All these infinite volume quantities have been computed using physical meson masses. LQCD energy-levels and those determined in this work are shown in Figs. 1 and 2. Statistical 68%-CL errors on the best fit parameters and derived quantities are calculated from the distributions obtained after performing a sufficiently large number of fits to synthetic sets of LQCD data, as explained in the caption of Fig. 1. In addition, the c − Λ correlation coefficients are −0.81, −0.93, 0.08 and −0.80 for fits AI, AII, BI and BII, respectively. Besides, in the Set C rows, we give the results obtained from a one-parameter (UV cutoff Λ)-fit that corresponds to an scheme where the LQCD energy-levels are described using finite-volume untarized LO HMChPT amplitudes. This is to say, the LEC c is fixed to zero, and therefore the contributions to the amplitudes of the exchange of even-parity charmed-strange CQM mesons are neglected. The volume dependence of the 0 + and 1 + energy-levels determined within this latter scheme are shown in Fig. 3 for the two lattice pion mass ensembles (I and II). 3 P 0 model is employed to couple both types of degrees of freedom, and the quark model wave functions provide form-factors that regularize the meson loops. Thus, all the inputs in this approach are constrained/determined from previous studies. The masses of the 0 + and 1 + states found in Ref. [40] are about 10 and 25 MeV higher than the experimental ones, respectively. The coupling of the CQM mesons to the DK and D * K thresholds is crucial to simultaneously lower the masses of the corresponding D * s0 (2317) and D s1 (2460) states predicted by the naive quark model. Such effects are of the order of 60 and 85 MeV in the 0 + and 1 + sectors, respectively. However, in the study carried out in Ref. [40], the one-loop corrections derived in Ref. [5] are taken into account, and they lower the predicted mass of the D * s0 (2317) by more than 100 MeV. At this respect, we should repeat once more that the simultaneous analysis of the 0 and 1 + LQCD energy-levels of Ref. [36] carried out in this work strongly disfavors such corrections. Molecular probabilities are reported in Ref. [40] to be around 33% and 54% for the D * s0 (2317) and D s1 (2460), respectively. Though the latter one turns out to be in a nice agreement with our results, the former one is around twice smaller than that found here and in previous works [45,46,49], and it would contradict a dominant molecular picture for the D s0 (2317). Moreover, this disparity between the molecular components in the scalar and axial states might also question that the D * s0 (2317) and D s1 (2460) could form a HQSS doublet. Within our scheme, it is however natural to assign these two states to the j π q = 1/2 + HQSS doublet, assignation that gets further support from the observation that the experimental mass splitting between these two resonances is remarkably similar to that between the D and D * mesons. Furthermore, interpreting the D * s0 (2317) and D s1 (2460) as DK and D * K bound states, the binding energies of both states will be very similar (approximately 46 MeV versus 54 MeV).
The couplings of the D * s0 (2317) and D s1 (2460) to D ( * ) K and D ( * ) s η are also compiled in Table III. We see that the coupling of both states to the latter channel, though around a factor two smaller than to D ( * ) K, is not negligible. 11 The analysis adopted in the original LQCD work of Ref. [36] led to g D ( * ) K couplings of 11.0±1.3 GeV and 13.8 ± 1.3 GeV for the D * s0 (2317) and D s1 (2460), respectively. These values are in good agreement with the values found in this work. We would like to stress that the clear similarities between the couplings of both resonances reinforces our conclusion that they form a HQSS doublet. Moreover, the D * s0 (2317) and D s1 (2460) should be heavy-flavor partners of theB s scalar and axial mesons found in Ref. [59] at 5709 ± 8 and 5755 ± 8 MeV, respectively. Note that the mass shift, due to the breaking of HQSS, is much smaller in the bottom sector, and it turns out to be quite similar to (MB * −MB), as expected.

D ( * ) K scattering lengths
The D ( * ) K scattering lengths are negative (see Table III), compatible with the interpretation of the D * s0 (2317) and D s1 (2460) as bound states. Indeed, the zero-range approximation, a 0 = −1/(2µB) 1 2 [B > 0 and µ are the D ( * ) K binding energy and reduced mass, respectively], provides already the first significant digit (−1 fm). This simple formula also anticipates that |a DK | > |a D * K |. Our predictions for the scattering lengths are consistent, within errors, with previous lattice determinations [36][37][38], and the bulk of the small differences existing among central values can be explained in terms of the differences between binding energies. The uncertainties on our estimates are significantly smaller than those affecting the LQCD ones. This is in line with the previous discussion about the errors on the masses of the D * s0 (2317) and D s1 (2460) resonances, and it shows, once more, that the present analysis of the LQCD energylevels leads to more precise results than those based on 11 The much larger differences found for the molecular probabilities are due, not only because it appears the square of the couplings, but also because the large distance of the D the Lüscher's relation using the effective range approximation.

Second pole: resonances
Within the current scheme, the amplitudes include an explicit pole. It is therefore reasonable to assume that the CQM bare state does not disappear, but it gets dressed by the meson-meson interaction moving into the complex plane. In addition to a mass shift, the new state acquires a width since it can decay into S-wave D ( * ) K and D ( * ) s η meson-pairs. The position and couplings of these extra poles, located in the SRS of the amplitudes, are collected in Table IV. We see that both in the 0 + and 1 + sectors, the resonances are relatively broad (85(4) and 98(5) MeV), respectively) and have similar couplings 12 to D ( * ) K and D ( * ) s η. On the other hand, the couplings of these resonances to D ( * ) K are around a factor of three smaller than those of the D * s0 (2317) and D s1 (2460) states.
The masses of these higher states are 175 (330) and 110 (270) MeV above the D ( * ) s η (D ( * ) K) threshold. Thus, we should take these results with some caution, since most likely they should be affected by sizable higher order chiral corrections and higher threshold-channels corrections. In other words, they are not as theoretically robust as those concerning the lowest-lying D * s0 (2317) and D s1 (2460) states. As mentioned, these additional resonances are likely originated from the bare cs-quarkmodel poles that are dressed by the D ( * ) K and D ( * ) s η meson loops. In that case, the bare poles have been highly renormalized, moving to significant higher masses and acquiring a significant width. We should also bear in mind that radial excitations (2 3 P 0 ) of the CQM states [8] or D ( * ) K * two-meson loops, neither of them taken into account in this study, might lie in this region of energies, then having a certain impact in the dynamics of these resonances.

LO HMChPT energy-levels
In addition to the results shown in the previous sections, where a CQM pole was added to the LO D ( * ) K interaction, it is enlightening to discuss whether we are able to describe the lattice data without any CQM cs contribution. To explore this scenario, we have performed an additional one-parameter (UV cutoff Λ)-fit, where the LQCD energy-levels are described using finite-volume  Fig. 1, but in this case the predictions have been obtained after neglecting the contribution to the amplitudes of the exchange of intermediate even-parity charmed-strange CQM mesons, i.e., setting c = 0 and evaluating the energy-levels using finite-volume unitarized LO HMChPT amplitudes. The besfit UV cutoff in this scenario and some derived quantities are given in the Set C rows of Table III. untarized LO HMChPT amplitudes. Thus the LEC c is set to zero, and consequently the contribution to the amplitudes of the exchange of the CQM mesons vanish as well. This fit is labeled as Set C in Table III, and the obtained 0 + and 1 + energy-levels, as a function of the box-size, are shown in Fig. 3. We find a quite reasonable description of the LQCD data, and the infinite volume properties of the lowestlying states agree well with those deduced using the Set A of CQM bare masses, though molecular probabilities and couplings of the D s1 (2460) and D * s0 (2317) are now much more similar. 13 Note that the CQM exchange potential induces some HQSS breaking corrections, driven by the 0 + and 1 + cs bare mass-shift, and the fact that these contributions have been eliminated might explain the found pattern of probabilities and couplings. The more distinctive difference, however, is that the UV cutoff is around 1100 MeV. This is to say, the new UV cutoff is around 400 MeV higher than the values needed when the contribution of the CQM meson exchanges are kept. That reveals that higher order chiral corrections, 13 We find P D ( * ) K ∼ 65% and P D ( * ) s η ∼ 6% for both states, and adding the probabilities, we obtain molecular components above 70% in the wave-functions of both mesons. On the other hand, the higher D ( * ) s η channel becomes also more important in their dynamics. previously effectively accounted for by means of the CQM pole, are not negligible. Conversely, taking into account explicitly the exchange of (bare) CQM mesons is not crucial to describe the D * s0 (2317) and D s1 (2460) states, because such contributions can be accommodated by appropriately modifying the finite contributions derived from short-distance physics. This is expected since the CQM bare states lie far from the latter physical states, for which the unitarity meson loops play a fundamental role.
The UV cutoff Λ is expected to be larger than the wave number of the states (∼ 200 MeV) and, at the same time, small enough to prevent it from inducing large HQSS breaking corrections. From this perspective, one might think that values in the range of 1.1 GeV, comparable to the mass of the charm quark, do not seem appropriate within the spirit of an EFT based on HQSS. However, one should bear in mind that here we are using a Gaussian UV regulator, which will approximately correspond to a sharp-cutoff, q max , smaller by a factor 14 π/8 [50]. Thus, in terms of a sharp-cutoff, we are dealing with values of the order of 700 MeV that are more acceptable from the point of view of a HQSS EFT.
The Set A of Gaussian regulators found in Subsec. III A would correspond to sharp cutoffs of the order of 400-450 MeV, which are still larger than the wave number of the D * s0 (2317) and D s1 (2460) states. Nevertheless, we should remind here that the CQM bare masses are not observables, and they depend on the UV renormalization. Here we have fixed the CQM masses, and thus the fitted cutoffs should be understood as those needed to effectively account for higher order chiral corrections, when these bare poles are ncorporated [65]. This also hints to a certain scale for which the CQM of Ref. [40] might match the chiral EFT.

NLO HMChPT energy-levels
Within this context, it seems appropriate to calculate the energy-levels obtained from the unitarized HMChPT NLO amplitudes [47] described in Subsec. II D. As in the previous subsection, the exchanges of CQM bare poles are not included. Indeed, as we have discussed above, considering such contributions, together with the NLO corrections, might lead to some double-counting problems. The UV divergences are renormalized in Ref. [47] by using subtracted meson loop functions instead of a sharp-cutoff. However, both schemes can be related (see for instance the discussion in Appendix A of Ref. [73]), and the subtraction constants determined in [46] correspond, in good approximation, to a sharp-cutoff q max = 0.72 +0.07 −0.06 GeV, similar to the values mentioned in the previous discussion.
Results are shown in Fig. 4, where we see that this scheme provides a more than acceptable description of the 0 + and 1 + LQCD data, for both pion mass ensembles, despite having set all LECs to the values determined in Ref. [46]. We emphasize that, therefore, the energy levels shown in the different panels of the figure are predictions and do not imply any fine adjustment of any type of parameter. 15 We find this remarkable, and together with the similar good description found in Ref. [53] of the (S, I) = (0, 1/2) LQCD low-lying energy levels calculated in Ref. [56], provides a great support for the finite-volume amplitudes obtained after unitarizing the NLO HMChPT amplitudes derived in Refs. [46,47]. Hence, the D * 0 (2400) two-pole structure and the SU(3) pattern 16 of the 0 + 15 We should mention that the untarized NLO HMChPT description of the LQCD energy levels shown in Fig. 4 might be improved by allowing for a slight variation of the LECs, which could, for example, effectively account for some discretization/finite volume errors etc. Note that, in addition, these systematic errors could be also different to those affecting the lattice study of scattering lengths carried out in Ref. [46]. Nevertheless, the second levels are quite far from the thresholds, and one might need to explicitly include higher order chiral corrections. Otherwise, the re-fitted NLO LECs would be biased, since they would effectively account for those contributions. All of this is beyond the scope of the present work. 16 Among other predictions, we point out that the lower D * 0 (2400) and 1 + heavy-light sectors claimed in Ref. [53] seem rather robust from the theoretical point of view. All these results reinforce the new paradigm to study the spectrum of heavy-light mesons [54], and that questions its traditional interpretation in terms of constituent Qq degrees of freedom.

C. DK S-wave Phase shifts
Here we will show predictions for DK S-wave scattering phase-shifts [see Eq. (16)], and will take the opportunity to discuss few aspects of the renormalization dependence of the scheme. For the sake of brevity, we will not address the rest of channels. In this subsection, we will always show results obtained for infinite volume, using physical meson masses, and the Set A of bare masses to incorporate the CQM degrees of freedom.
Phase shifts at threshold depend strongly on the position of the S-wave DK pole. Thus, and to make the discussion meaningful, we will only consider approaches leading to the same D * s0 (2317) mass value, 2315 +18 −28 , obtained in the unitarized NLO HMChPT approach [46,47], and whose central value is quite close to the experimental one. Details can be found in Table V, while the related phase shifts are shown in the left panels of Fig. 5.
Looking at the results of Table V, we first stress the dependence of the UV cutoff on the LEC c, or viceversa the dependence of c on Λ, for a fixed CQM bare mass. The latter should be also understood as a renormalization scheme dependent quantity. All these three parameters (c, Λ, • m cs ) should effectively incorporate higher order chiral corrections beyond LO, and not accounted by the unitarity loops. One expects that these further corrections should be rather independent of the short distance physics at energies moderately far from threshold. However predictions might significantly differ, lets us say above 2450 or 2500 MeV. Indeed, we have an indication from the masses and widths of the possible second (higher resonance) state compiled in the table. We see that for c = 0.3, a narrow resonance (28 MeV) is generated at around 2557 MeV, that would produce clear signatures, not seen, in the second energy level calculated in Ref. [36]. The other schemes either do not generate any resonance (c = 0 and NLO) or it is located close to 2700 GeV (c = 0.61) and it is sufficiently broad to become unnoticed for energies below 2600 MeV. The conclusion is that the c = 0.3 predictions above 2500 MeV are unreliable, and that the exact location of the resonance generated for c = 0.61 is not well constrained by the data of Ref. [36].
resonance (located at 2105(8)−i102(12) MeV) and the D * s0 (2317) state would be siblings forming a3 SU(3) multiplet. The same can be said of the new D * 1 resonance (located at 2240(6) − i93(9) MeV) and the D s1 (2460) state in the 1 + sector. All these features have counterparts in the bottom sector as well.  [46,47] compared to the lattice results of Ref. [36]. Exchanges of intermediate CQM mesons are not included, and the distribution of panels is the same as in Fig 1. There are no fitted parameters involved in these predictions since all LECs that appear in the definition of the chiral amplitudes were determined from the lattice calculation [46] of the S-wave scattering lengths in several (S, I) sectors. The 68% CL uncertainty bands depicted in the plots are inherited from the errors on the LECs quoted in Ref. [46]. ) obtained from the unitarized LO HMChPT+CQM approach. The UV cutoff Λ is fitted to reproduce the above mass, that it is deduced from the unitarized NLO HMChPT T −matrix [46]. The uncertainties in the D * s0 (2317) mass, inherited from the errors on the LECs, are taken into account and lead to the errors on Λ. We first generate synthetic sets of LECs, according to the correlated error distributions given in [46]. Next, we find the position of the D * s0 (2317) pole for each set of LECs, and finally, for each parameter c, we tune the UV cutoff to reproduce this mass. In this way, we determine the error distributions of the cutoffs that are also used to estimate the uncertainties on the derived quantities (S-wave DK scattering length, mass and width of the higher -dressedresonance, D * s0 (2317) molecular probabilities and couplings) given in the table. Predictions for the latter quantities from the unitarized NLO HMChPT approach [46,47] are compiled in the last row. Finally note that the choice c = 0.61 is motivated by the Set A of results in Table III On the other hand, we appreciate a variation pattern for c consistent with the physical interpretation of this LEC, and thus the total molecular probability decreases from (67 ± 10)% down to (60 ± 9)%, when c varies from 0 to 0.61. Scattering lengths are mostly determined by the common D * s0 (2317) mass and do not show statistically significant variations, as it also occurs for the couplings.
Paying attention now to the NLO results, we find that they coincide reasonably well with those found using LO HMChPT amplitudes, thanks to the freedom in the latter to re-adjust the cutoff. Finally, it could be surprising that the DK molecular probability, obtained within the NLO scheme, is only 54 ± 4%, when a value of ∼ 70% was claimed in the original work of Ref. [46]. This dis- Note that in the case c = 0, the approach reduces to the unitarized LO HMChPT approach followed in [59] for the bottom-strange sector, except that in this latter case theBsη channel was not considered. In the left (right) panels the UV divergences, that appear in the unitarization of the LO HMChPT+CQM amplitudes, have been renormalized using a Gaussian cutoff (subtraction constant). For details see Tables V and VI. The UV behaviour of the NLO unitarized amplitudes is always renormalized using subtraction constants [46,47]. Statistical 68%-CL error-bands are generated from the uncertainties of the LECs that enter in each scheme.
crepancy is due to the use in the latter work of the Weinberg compositeness rule [48], that provides the molecular probability in terms of the scattering length and the DK wave number γ B , which leads to one, in the limiting case when the scattering length is approximated by −γ −1 B (very loosely bound states), neglecting finite effective range corrections. Indeed, using the relation of Eq. (31), one obtains P DK ∼ 70 ± 4 %, as obtained in Ref. [46]. Note, however, that Eq. (31) has corrections due to the D * s0 (2317) binding energy [45], which is not too small (∼ 46 MeV), and possible inelastic effects [76], which should bring the 70% down to the more accurate estimate found in the present work.
Coming back to the phase shifts shown in the left panels of Fig. 5, we see that different schemes produce compatible phase shifts close to threshold, while the differences become larger as the energy increases. Thus for instance, the c = 0.3 phase shifts suddenly change curvature above 2500 MeV, but as we discussed above such behaviour, produced for a narrow resonance at 2557 MeV (see Table V), is not compatible with the higher 0 + LQCD energy level reported in Ref. [36]. Above 2450 MeV the rest of schemes lead to differences in the phase shifts of few degrees, and at most of around 15 • at 2550 MeV. However, taking into account the uncertainties of the different predictions, phase shifts are almost compatible up to this latter energy.
We now discuss about the dependence of the analysis on the regulator scheme. To this end, we consider the loop matrix function G renormalized by one subtraction, as in Refs. [46,47]. Suppressing the indices, the loop function is written for each channel as,  (33) where ν is the scale of the dimensional regularization. Changes in the scale are, in principle, reabsorbed in the subtraction constant a(ν), so that the results remain scale independent. Here we have taken ν = 1 GeV and a common subtraction constant a(1 GeV) = α for both DK and D s η channels, as in Refs. [46,47]. We have now constrained the LEC α to obtain the same D * s0 (2317) mass (2315 +18 −28 ), as in Table V, for each value of the parameter c. The results are presented in Table VI. Comparing these latter results with those in Table V, we see that the predictions, within uncertainties, are consistent in the two renormalization schemes. There is a slight dependence, and the D * s0 (2317) coupling to DK and the modulus of the scattering length are smaller and closer to those deduced from the NLO approach. Molecular probabilities are somewhat larger, specially P Dsη that becomes almost twice bigger. As a consequence, the total hadronic molecular component is now roughly 70 ± 10%. The dependence on c follows a similar pattern as in Table V, and the importance of the D s η channel decreases as the value of the LEC c increases. Finally, the results concerning the higher-dressed-resonance pole position are similar in the two renormalization schemes. Thus from the previous discussion, the c = 0.3 predictions for energies above 2500 MeV turn out to be little reliable also in this scheme.
The phase shifts deduced from the various possibilities discussed in Table VI are shown in the right panels of Fig. 5. First, the c = 0.3 phase shift above 2500 MeV presents the same pathologies as in the left top panel, where an UV Gaussian cutoff is used. It is interesting to note that the c = 0 and c = 0.61 phase shifts have smaller errors, the two sets of phase shifts are still statistically compatible, but in addition, they now agree quite well with the NLO predictions. Thus, the renormalization scheme dependence, while it is not much relevant for the D * s0 (2317) properties, turns out more important for the phase shift at energies above 2450 MeV, as one could reasonably expect from the previous discussions. This clearly could be regarded as a source of systematic uncertainty, though up to 2550 MeV remains smaller/comparable to the other uncertainties of the predictions accounted for the error bands depicted in Fig. 5.

IV. SUMMARY AND CONCLUDING REMARKS
In this work we have first carried out a coupled channel study of the 0 + and 1 + charm-strange meson sectors employing a chiral unitary approach based on LO HMChPT P ( * ) φ interactions, and that incorporates, consistently with HQSS, the interplay between intermediate CQM bare cs and P ( * ) φ degrees of freedom. We have extended the scheme to finite volumes and fixed the strength of the coupling between both types of degrees of freedom to the available 0 + and 1 + LQCD energy-levels [36], which we have successfully described. On the other hand and at variance to the situation in the bottom-sector reported in Ref. [59], we have found that the 0 + and 1 + CQM bare masses (denoted as Set B in this work) obtained in Ref. [72] using the one-loop corrections to the CQM OGE potential proposed in Ref. [5], lead to a really poor description of the LQCD data. This is because the HQSS breaking corrections induced by this modification of the OGE potential are inconsistent with the LQCD energy levels calculated in Ref. [36].
We have estimated the size of the the D ( * ) K twomeson components in the D * s0 (2317) and D s1 (2460), and conclude that these states have a predominantly hadronic-molecular structure. Furthermore, we have observed a quite significant dependence of the lowest-lying LQCD energy-levels of Ref. [36] on the pion mass, which is difficult to accommodate by a dominant CQM cs component. This is, however, consistent with having a large influence of the P ( * ) φ loops in the D * s0 (2317) and D s1 (2460) structure.
In addition, we have found one extra resonance, in both the 0 + and 1 + sectors, arising from the dressed CQM states. Our predictions for these states are not as robust as those for the low lying D * s0 (2317) and D s1 (2460), and moreover they are relatively broad, which might complicate their discovery. Some experimental efforts are needed to clarify their possible existence.
The LEC c depends on the radial quantum number, but not on the heavy flavor, up to Λ QCD /m Q corrections. Thus, the value determined here for this parameter should be similar to that found in the bottom-strange sector in Ref. [59]. There, it was obtained c = 0.75 (6), which is quite compatible with the values in the range 0.52-0.70 found in this work for the Set A of CQM bare masses. Note that in addition to heavy flavor symmetry breaking corrections, there might be also some discretization errors. Nevertheless, we have shown that taking into account explicitly the exchange of (bare) CQM mesons is not fundamental to describe the D * s0 (2317) and D s1 (2460) states, since such contributions can be accommodated by appropriately modifying the finite contributions derived from short-distance physics. This is natural because the CQM bare states lie far from the latter physical states, for which the unitarity meson loops play a fundamental role.
We have discussed how the approach followed here, where the two meson loop function is computed in a finite volume, the unknown LECs are determined from fits to the LQCD data, and finally poles are searched for in the infinite volume unitarized amplitudes using physical meson masses, provides a theoretically sound tool to analyze the LQCD energy-levels. We have shown that such procedure leads to more precise predictions that those obtained via the Lüscher's relation using the effective range approximation.
We have also calculated the energy-levels obtained from the unitarized HMChPT NLO amplitudes derived in Ref. [47], without including any contribution from the exchanges of CQM bare poles. We have shown ( Fig. 4) that this scheme provides a more than acceptable description of the 0 + and 1 + LQCD energy levels of Ref. [36], despite having fixed all LECs to the values previously determined in Ref. [46] (not fitted to the energy levels). These findings, together with the similar good description found in Ref. [53] of the (S, I) = (0, 1/2) LQCD low-lying energy levels calculated in Ref. [56], provide a great support for the amplitudes obtained after unitarizing the NLO HMChPT amplitudes derived in Refs. [46,47]. Hence, the D * 0 (2400) two-pole structure and the SU(3) pattern of the 0 + and 1 + heavy-light sectors claimed in Ref. [53] seem rather robust from the theoretical point of view. All these results reinforce a new paradigm to study the spectrum of heavy-light mesons [54], that questions its traditional interpretation in terms of constituent Qq degrees of freedom.
Finally, we have predicted S-wave DK phase shifts and discussed few aspects of the renormalization dependence of our results.