Lowest lying even-parity $\bar B_s$ mesons: heavy quark spin-flavor symmetry, chiral dynamics, and constituent quark model bare masses

The discovery of the $D^\ast_{s0}(2317)$ and $D_{s1}(2460)$ resonances in the charmed-strange meson spectra revealed that formerly successful constituent quark models lose predictability in the vicinity of two-meson thresholds. The emergence of non-negligible effects due to meson loops requires an explicit evaluation of the interplay between $Q\bar q$ and $(Q\bar q)(q\bar q)$ Fock components. In contrast to the $c\bar s$ sector, there is no experimental evidence of $J^P=0^+,1^+$ bottom-strange states yet. Motivated by recent lattice studies, in this work the heavy-quark partners of the $D_{s0}^\ast(2317)$ and $D_{s1}(2460)$ states are analyzed within a heavy meson chiral unitary scheme. As a novelty, the coupling between the constituent quark model P-wave $\bar B_s$ scalar and axial mesons and the $\bar B^{(\ast)}K$ channels is incorporated employing an effective interaction, consistent with heavy quark spin symmetry, constrained by the lattice energy levels.


I. INTRODUCTION
The low-lying positive parity charmed-strange spectrum moved into spotlight in 2003 after the experimental observation of the D * s0 (2317) and D s1 (2460) states by the BaBar [1] and Cleo collaborations [2] in the D ( * )+ s π 0 invariant mass spectrum with J P = 0 + and 1 + quantum numbers, respectively. The interest in such states arose from their low masses and widths, contrary to constituent quark model (CQM) [3][4][5][6] and lattice QCD (LQCD) [7][8][9][10] expectations. Consequently, the nature of those resonances became a matter of intense research, mainly being interpreted as naive cs [11][12][13] or meson-meson and four-quark states [14][15][16][17][18][19][20][21]. An indirect hint for the nonperturbative nature of the DK chiral amplitudes near threshold, and the possible existence of a bound state, was obtained in Ref. [22] from an Omnès analysis of the LQCD data for the scalar form factor in the D →Kℓν ℓ semileptonic decay. More recently, lattice QCD simulations [23][24][25] and the quark model calculations [26,27] emphasized the importance of including the D ( * ) K dynamics when describing the P-wave cs mesons. These meson-meson channels produce large mass shifts which improve the description of the experimental masses.
The combination of a heavy quark and a light antiquark in the D s orB s spectra is a great advantage when it comes to describe the system. In such mesons, heavy quark spin symmetry (HQSS) [28][29][30][31] (see also the book [32]) is in good approximation fulfilled by QCD, being thus the parity and the total angular momentum of the light antiquark jq approximately conserved. The decoupling of the spin of the heavy quark (s Q ) and jq generates the arrangement of states in doublets labelled by their parity 1 and jq values, so members within a doublet are governed by the same dynamics and become degenerate in mass, up to Λ QCD /m Q corrections, with m Q the heavy quark mass and Λ QCD , a typical scale accounting for the dynamics of the light degrees of freedom. The D * s0 (2317) and D s1 (2460) are members of a positive-parity jq = 1 2 doublet, which will strongly couple to S-wave D ( * ) K pairs, being the dynamics of these latter meson pairs in turn governed by chiral symmetry. Within a CQM scheme, these resonances will correspond to P-wave states, where the spin and angular momentum of the light antiquark are coupled to a total jq = 1 2 . There will be another HQSS doublet with j P q = 3 2 + .
Besides, in the m Q → ∞ limit the dynamics of systems containing a single heavy quark becomes also independent of the flavor of the heavy quark [32]. 2 Hence the bottom-strange sector is expected to present a pattern similar to that of the cs one, and in particular there should exist heavy-flavor partners of the D * s0 (2317) and D s1 (2460) resonances. Furthermore, since the b quark is heavier than the c quark, the O(m −1 Q ) corrections are expected to be smaller, and thus the HQSS relations should be more accurate. In these circumstances, the bs j P q = 1/2 + doublet is a perfect scenario to discuss the interplay between CQM states and meson-meson channels with thresholds located close to the former. This study is relevant to unveil the nature of the D * s0 (2317) and D s1 (2460), where such interplay turns out to be essential to understand the dynamics of these states. Hence, we expect also a strong influence of the close continuum twomeson channels on the properties of the bottom partners of these even-parity resonances.
Unfortunately, unlike the cs spectrum, the lowest-lying positive parity bs states have not been fully discovered. While experimental searches have successfully observed the j P q = 3 2 + doubletB s1 (5830) andB * s2 (5840) [33,34], the lower mass j P q = 1 2 + doublet states still wait to be observed. Note that the dynamics of the j P q = 3 2 + doublet is not governed by chiral symmetry, since theB ( * ) K * channel, involving the light vector meson K * , should be taken into account. Under this lack of experimental data, many theoretical predictions have been produced within a wide variety of techniques (quark models [6,35,36], effective field theory approaches (EFTs) [37][38][39][40][41][42][43][44], and LQCD [45,46]). A special attention deserves the recent LQCD study of the even-parity isoscalar bs energy-levels carried out in Ref. [46]. There, clear signatures for theB s1 (5830) and B * s2 (5840) are found above theB ( * ) K thresholds. Below these thresholds, two QCD bound states are identified using a combination of quark-antiquark andB ( * ) K interpolating fields, and assuming that the mixing with B ( * ) s η and the isospin-violating decaysB ( * ) s π are negligible. A J P = 0 + bound state, with mass 5.711 ± 0.023 GeV (adding errors in quadrature) is predicted and with some further assumptions, it is also found a 1 + state with a mass of 5.750 ± 0.025 GeV [46].
In this work, we will pay attention to the 0 + and 1 + isoscalar bottom-strange sector. We will use a heavy meson chiral unitary scheme to describe the isoscalar Swave elasticB ( * ) K →B ( * ) K T −matrix. The scattering amplitudes will be obtained by solving a renormalized Bethe-Salpeter equation (BSE) with an interaction kernel determined from leading order (LO) heavy meson chiral perturbation theory (HMχPT) [47][48][49][50]. We will couple the two-meson channels with the CQM P-wavē B s scalar and axial mesons using an effective interaction consistent with HQSS. To that end, we will follow the scheme detailed in Ref. [51], where the (DD * + h.c.) two-meson channel was coupled to the χ c1 (2P ) charmonium state, and the consequences for the X(3872) and its spin-flavor partners, were examined. Finally, we will use the LQCD energy levels reported in Ref. [46] to constrain the undetermined low-energy-constants (LECs) of the present approach. As a final outcome, we will present robust predictions for the lowest-lying bs J P = 0 + and 1 + states, that can serve as an important guidance for experimental searches and to shed light into the situation in the analog charm sector.

A. HQSS fields
To study theB ( * ) K interactions, we first introduce the matrix field H which combines the isospin doublet and singlet of pseudoscalar heavy-mesons P (Q) a = Qū, Qd, Qs fields and their vector HQSS partners P * (Q) a . We use the isospin conventionū = |1/2, −1/2 andd = −|1/2, +1/2 , which induces a minus sign between the D + (andB 0 ) particle and isospin states. Besides v is the four-velocity of the mesons, and the vector field satisfies v · P * (Q) a = 0. Note that in our normalization the heavy-light meson field, H (Q) , has dimensions of E 3/2 (see Ref. [32] 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 Qq states, associated to the j P q = 1 2 + HQSS doublet, are described by the matrix field J (Q) a [52], with v µ Y * (Q) a µ = 0. The Y a and Y * a fields respectively annihilate the 0 + and 1 + meson states belonging to the Under a parity transformation we have The field H (Q) a transforms as a (2,3) under the heavy spin ⊗ SU(3) V symmetry [47], this is to say: The hermitian conjugate field is defined by: and transforms as [47]: The matrix field J (Q) a satisfies transformation relations identical to those in Eqs. (3)-(6).

B. Interactions
In S-wave, the spin-parity quantum numbers of thē BK B * K meson pair are 0 + (1 + ), while the light degrees of freedom are coupled to spin-parity 1/2 + . As in Ref. [46], we will neglect theB ( * ) s η channels, and thus the (elastic) isoscalarB ( * ) K →B ( * ) K interaction potential V (s) consists of a chiral contact term [V c (s)] plus the exchange [V ex (s)] of a bare bs state, with s, the center of mass (c.m.) energy squared. At LO in the chiral counting, V c (s) is given by the Weinberg-Tomozawa Lagrangian that reads [47][48][49][50] (omitting from now on the (Q) superscript), where the ξ field contains the pseudoscalar light mesons, and can be written as ξ = exp( iM √ 2f ), with f ∼ 93 MeV, and the M matrix defined in Ref. [48]. The Lagrangian density in Eq. (8) is invariant under SU(3) L ×SU(3) R chiral symmetry, Lorentz transformations, HQSS and parity. From the previous Lagrangian, the isoscalar contact term contribution V c (s) can be easily derived, and the result after projecting into S-wave reads: 3 (9) Neglecting theB * −B mass difference, the interactions in the J = 0 and 1 sectors are identical, as expected from HQSS because they correspond to the same configuration (1/2 + ) of the light degrees of freedom.
The exchange term in Eq. (7) is determined by the coupling between theB ( * ) K meson pairs and the P-wave bare quark model states described by the matrix field J (Q) a introduced in Eq. (2). At LO in the heavy quark expansion, there exists only one term invariant under Lorentz, parity, chiral and heavy quark spin transformations, 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 both the heavy-quark flavor and 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, the same LEC governs the interplay between two-meson and quark model degrees of freedom in all isospin and strangeness channels. Paying attention to the isoscalar bottom-strange sector, we find a bs →B ( * ) K coupling in S-wave, where • m bs is the mass of the bs meson without the effect of theB ( * ) K meson loops. This mass is the same, up to small HQSS breaking corrections, for both J = 0 and J = 1 sectors, and it can be in principle obtained from CQMs. In what follows, we will denote it as the bare mass of the state. 4 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 theB ( * ) K 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 EFT approaches. We will come back to this point later.
The Lagrangian density in Eq. (10) allows to compute the V ex (s) term contribution to the full potential, Eq. (7), that accounts forB ( * ) K scattering via the exchange of intermediate even-parity bottom-strange mesons [51]: The HQSS consistent potential V (s) detailed above is used to obtain theB ( * ) K elastic unitary amplitude, T (s). This is done by solving the BSE within the so-called onshell approximation [53]. We use [22]: where G(s) is the two-meson loop integral, regularized with a Gaussian cut-off, Above, P(· · · ) stands for the principal value of the integral and with ω( q ) and ω ′ ( q ) the energies of theB ( * ) and K mesons, respectively. Finally, k 2 is the square of the c.m. on-shell momentum, C. Bound, resonant states, couplings and the compositeness condition for bound states The dynamically-generated meson states appear as poles of the scattering amplitudes on the complex s−plane. The poles of the scattering amplitude on the first Riemann sheet (FRS) that appear on the real axis below threshold, s + , are interpreted as bound states. The poles that are found on the second Riemann sheet (SRS) below the real axis and above threshold are identified with resonances. The SRS is simply obtained by analytical continuation of the amplitude in the physical FRS across the unitarity cut, where Note that the cuts for k(s)/ √ s go along the real axis for −∞ < s < s − and s + < s < ∞. The function is chosen to be real and positive on the upper lip of the second cut, s + < s < ∞, and hence it satisfies: The mass and the width of the state can be found from the position of the pole on the complex energy plane. Close to the pole, the scattering amplitude behaves as The mass M R and width Γ R of the state result from √ s R = M R − i Γ R /2, while g (complex in general) is the coupling of the state to theB ( * ) K channel.
Information on the compositeness of the bound states can be obtained from the derivative of the meson loop function and the residue at the pole position. Indeed, it can be shown [54,55], inspired by the Weinberg compositeness condition [56][57][58], that the probability of finding theB ( * ) K molecular component in the bound state wave function is given by where M b is the mass of the bound state and g the coupling of the state to theB ( * ) K meson pair. The above probability deviates from one because of the energy dependence of the potential [Eq. (12)] [55,59], which is enhanced by the exchange of intermediate quark model (bare) bottom-strange mesons [51]. We do not extend this discussion to resonances, restricting it to the case of bound states. For poles located in the complex plane the strict probabilistic interpretation is lost, since the definition in Eq. (22) gives rise to complex numbers (see for instance the discussion in Ref. [51]).

D. Finite volume
To compare with LQCD simulations, we consider the T -matrix [Eq. (13)] in a finite box of size L. The boundaries of the box impose a quantization condition for the momentum, q = 2π L n, with n ∈ Z 3 . The loop function G(s) is thus replaced with its finite volume counterpart, G(s, L) [60,61], Up to the order we are considering in this work, the potential does not receive finite volume corrections, and hence the finite volume T -matrix, denoted as T (s, L), reads: The energy levels obtained in LQCD simulations can be computed within our approach from the poles of T (s, L).
To better describe the energy levels reported in the LQCD simulation carried out in Ref. [46], we use the masses and the energy-momentum dispersion relations given in that work. In particular, we will employ a modified energy-momentum dispersion relation for theB ( * ) mesons, where the parameters appearing in the above equation can be found in Table 1 of Ref. [46]. The lattice size and spacing in that simulation are 32 3 × 64 and a = 0.0907 ± 0.0013 fm, respectively, while the pion mass is 156 ± 7 MeV. For the kaon, the ordinary relativistic dispersion relation is used with an unphysical mass of m K a = 0.2317 ± 0.0006 (m K = 504 ± 7 MeV) [24]. To compute the potentials (chiral+exchange) at finite volume, we set theB ( * ) mass to m 1 , introduced in the modified dispersion relation of Eq. (25). Finally k 2 , that appears in the Gaussian regulator needed to render G(s, L) finite is obtained from Eq. (17) using the lattice masses.

III. RESULTS
We will fit our model to the levels 1,2,3 and 1,3,4 given in Table 3 of Ref. [46] for the 0 + and 1 + sectors, respectively. 5 We will consider the energy-levels in lattice units. Hence, and to properly take into account the uncertainty on this scale, we introduce it as an additional best-fit 5 Note that the 1 + level 2 is interpreted in [46] to be the j P q = 3 2 + state with J P = 1 + , which does not couple to B * K in S-wave in the heavy quark limit.  [46], together with masses and properties of the low-lying j P q = 1 2 +B s meson doublets deduced from these parameters in the infinite volume case. In this latter case, physicalB ( * ) and K masses have been used, and we have searched for poles in the FRS (bound states) and SRS (resonances) of the isoscalar S-waveB ( * ) K amplitudes. Besides, we find a th = 0.0909 ± 0.0013 fm and 0.0910 ± 0.0013 fm for sets (a) and (b), respectively, which compare rather well with the lattice spacing (a = 0.0907 ± 0.0013 fm) quoted in Ref. [46]. The isoscalar 0 + and 1 + scattering lengths (a0) and the isoscalar S-waveB ( * ) K →B ( * ) K amplitudes at threshold are related by T (s+) = −8πa0 √ s+, with s+ = (mB( * ) + mK ) 2 . TheB ( * ) K molecular probability PB( * ) K is calculated using Eq. (22) and it is given only for bound states. The coupling g, defined in Eq. (21), is also calculated only for the bound state. Errors on the fitted parameters show 68% confidence levels (CLs), which are obtained from distributions computed by performing a large number of best fits to Monte Carlo synthetic datasets. The synthetic datasets are generated from the original energy levels of Ref. [46] and the lattice spacing assuming that the data points are Gaussian distributed. The 68% CLs are obtained by discarding the higher and the lower 16% tails of the appropriate distributions. These parameter distributions are used to estimate the uncertainties on the derived quantities in the infinite volume case.

Parameters
Infinite volume predictions a0  We compare our predictions for the difference E −m with the results of Ref. [46], given also in Table 3 of that work. We have usedma = 1.62897 (43), as in Ref. [46]. We have evaluated energy levels ( T (s, L) presents poles at √ s = E for a given L, when the lattice masses are employed) using both sets of parameters compiled in Table I, which are shown by dashed blue and solid red lines for sets (a) and (b), respectively. Black lines stand for the first and second non-interactingB ( * ) K levels, while the data points show the lattice levels reported in Ref. [46]. The uncertainty bands in the predicted energy levels mark 68% CLs obtained from a Monte Carlo simulation considering the parameter distributions of Table I. parameter, a th , constrained by the central value and error quoted above and taken from Ref. [46]. Thus, we minimize the following χ 2 , The sum runs over the six 0 + and 1 + energy-levels (given in lattice units) specified above. We could instead have fitted directly to the energy levels in physical units, but in that case, the errors on the levels inherited from the lattice spacing need to be treated as correlated ones, since variations in the lattice spacing affect to all the energy levels in the same manner. In addition to the use of correlated errors, one would have also to consider that the lattice meson masses, appearing in the theoretical T (s, L), will change with a as well because of their conversion into physical units. The χ 2 introduced in Eq. (26) accounts for all these correlations. Indeed, the latter make the uncertainties on (Ea) lat i to become the most relevant ones. Note that if these correlations induced by the lattice spacing are not taken into account, one will end up with large and unrealistic errors. On the other hand, one might treat m K a, (m 1,2,4 )a also as best-fit parameters, supplementing appropriately the χ 2 . Results do not change significantly and for simplicity we have fixed all these masses in lattice units to the central values reported in Ref. [46]. This little sensitivity can be expected since the errors on m K a and m 1 a, which determine the threshold and the chiral potential, are much smaller than those on a. Indeed, the largest uncertainty in the magnitude of these quantities is induced from the error on the lattice spacing.
Besides the lattice spacing, the parameters of the model are the bare masses of the CQM bs 0 + and 1 + mesons, • m bs , the LEC c that gives us the strength of the coupling of the latter states with the two-meson B ( * ) K channels, and the UV Gaussian regulator Λ. We expect to get reasonable estimates for the bare masses from the predictions of CQMs [6,35,36]. These kind of models find masses in the ranges 5800 − 5850 MeV and 5840 − 5880 MeV for the J P = 0 + and J P = 1 + sectors, respectively. TheBK andB * K thresholds are located approximately at 5775 MeV and 5820 MeV, respectively. Thus, in principle, we expect the quark model states to be relatively close to, but above, the respectiveB ( * ) K thresholds, which would produce attractiveB ( * ) K interactions for energies below the bare masses [51].
We will explore the ranges of bare masses mentioned above, and we will perform two independent fits using • m bs values close to their respective lower and upper limits. To maintain a consistent 0 + − 1 + mass splitting, we will use the predictions of a widely used non-relativistic CQM [62][63][64]. This quark model was already employed to study the low-lying P-wave charmed-strange mesons [27]. In that reference, since the D * s0 1 3 P 0 (n 2S+1 L J ) bare state was found significantly above the experimental level (2511 versus 2317.7 MeV), an additional one-loop correction to the one-gluon exchange (OGE) potential was introduced. This extra term was motivated from the studies of Refs. [13,65], where a spin-dependent term was added to the quark-antiquark potential affecting only mesons in the case of unequal quark masses. Such correction is in general rather small, except for the 0 + sector, where a large mass shift was found (around 128 MeV in the case of the D * s0 ). Hence, as commented above, we will consider two sets of bare masses • m bs . In the set (a), we will use bare masses of 5851 MeV and 5883 MeV for the 0 + and 1 + states, respectively, as deduced from the CQM of Refs. [62][63][64] without including any correction to the OGE potential. For the second set, (b), we will fix the 0 + and 1 + bare masses to 5801 MeV and 5858 MeV, as obtained when the latter CQM is supplemented with the OGE one-loop terms discussed in Ref. [13,65,66]. Since the LQCD simulation carried out in Ref. [46] uses nonphysical meson masses, the CQM bare masses have been corrected using the difference between the experimental and the prediction of Ref. [46] for the spin-average mass m = 1 4 (mB s + 3mB * s ).
For each set of bare masses, the values of the other two parameters, c and Λ, are obtained from a simultaneous fit to the first three J P = 0 + and 1 + energy levels reported in the LQCD study of Ref. [46]. In Fig. 1, we present the predicted energy levels as a function of the lattice size L, and the values of the fitted parameters are compiled in Table I. As can be seen, we find an excellent description of the LQCD levels of Ref. [46] in both J P sectors, despite the large deviations from the free levels [67]. The LEC c is rather insensitive to the used bare masses, taking a value about 0.75 with an error of around 6%. The central values of the UV regulator present however a mild dependence, and we find Λ = 730 ± 40 MeV for set (a), while for set (b) the fitted value is Λ = 650 ± 30 MeV. We should remind that the CQM bare masses depend on the renormalization scheme, in particular on the UV regulator, or equivalently the UV regulator is expected to depend on the bare masses. Nevertheless, set (a) and (b) UV regulators turn out to be almost compatible within errors.
When the loop function is renormalized by a suitable subtraction, instead of using a Gaussian regulator, the physical results showed in Table I and Fig. 1 do not appreciably change, besides some variation of the renormalization-scheme-dependent low energy constant c. Thus, a similar good reproduction of the LQCD energy levels is achieved. Note that the finite volume loop function, in both renormalization schemes, is related to the Lüscher function [68,69], as it is shown in Refs. [60,61]. Hence, the continuous volume dependent curves in Fig. 1 are essentially the Lüscher curves obtained from the phase shift by solving δ(q) + φ(q) = nπ (27) withq = qL/2π and φ(q) determined by the Lüscher function (see Eq. (6.13) of Ref. [69]). Next, and once the parameters have been fixed, we search for poles in the FRS (bound states) and SRS (resonances) of the isoscalar S-waveB ( * ) K amplitudes for the infinite volume case. Pole positions are also compiled in Table I, together with the 0 + and 1 + isoscalar scattering lengths and the probabilities of the molecular B ( * ) K component in the bound states. For both sets of parameters, and for both J P = 0 + and 1 + sectors, we find a bound state (FRS) and a resonance (SRS).
In all cases, the mass of the bound state is rather independent of the UV regulator, or equivalently of the bare quark model mass, and it is located more than 100 MeV below the corresponding bare pole, consequence of the strong attraction produced by the chiral potential. This is a first hint of the importance of the meson loops in the dynamics of the bound state, which can be also inferred from its large (∼ 50%) molecular component. From the results of the Table I, we predict masses of 5709 ± 8 MeV and 5755 ± 8 MeV for theB * s0 andB s1 states, respectively. These states are the heavy flavor partners of the charmed-strange D * s0 (2317) and D s1 (2460) resonances, and are clear candidates for future experimental searches.
The masses obtained in this work are in excellent agreement with the estimations given in Ref. [46], and mentioned in the introduction. 6 They are also quite compatible within errors with other HMχPT predictions, where the explicit coupling (LEC c) of the two-meson channel and the bare quark model state is not considered 7 [39-41, 43, 44]. In all cases a similar binding energy around 60 − 70 MeV is obtained, which favors a molecular interpretation of such states, where one would expect a (B s1 −B * s0 ) mass splitting similar to the (B * −B) one. The latter is around 45 MeV, while in our calculation we find (mB s1 − mB * s0 ) ∼ 41 MeV for set (a) and ∼ 49 MeV for set (b), around 4 MeV below and above the pseudoscalar-vector mass splitting, respectively. This is a clear indication of the non-canonical quark model nature of theB s1 andB * s0 states. It is interesting, though, to note that the molecular proportion in the 0 + state (∼ 50%) is below the EFT estimations for the D * s0 (2317), predicted to be around 70% of D ( * ) K [21,25,72].
In Ref. [66] 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. TheB ( * ) K interactions are derived from the same CQM used to compute the bare states, instead of using HMχPT. The 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. The 0 + and 1 + states reported in Ref. [66] are around 30 − 40 MeV less bound than those found here and in the LQCD study of Ref. [46]. Presumably, this is because theB ( * ) K interactions derived in the CQM of Ref. [66] are weaker than those used here. Molecular probabilities are reported in Ref. [66] to be around 30 − 40%, which are smaller than those found in the present approach.
Regarding the isoscalar scattering lengths, we predict (combining the results of both sets) a 0 = −0.89 ± 0.07 fm for both J P = 0 + and 1 + sectors, which compares well with the results aB K 0 = −0.85 ± 0.10 fm and aB * K 0 = −0.97 ± 0.16 fm, obtained in the analysis carried out in the LQCD study of Ref. [46]. In the approach of Ref. [66], the 0 + and 1 + scattering lengths turn out to be 6 Note that the uncertainties obtained here are smaller than those quoted in Ref. [46] because we go beyond the effective range approximation and determine a potential (see the discussion in Ref. [61]). 7 In these schemes, such effect is encoded either in the renormalization subtraction constants or in higher order LECs, appearing at next-to-leading order (NLO) in the HMχPT expansion.
(Note that in the present work, we obtain reasonable UV cutoff values ∼ 700 MeV, which do not hide large higher order contributions [70,71].) Despite the bare quark model pole was located above, and relatively close to, theB ( * ) K threshold, we find the bound state significantly below • m bs . Hence, the bare pole induces a mild energy dependence in the vicinity of the physical bound states, which can be accounted for by means of local terms in the potential [51]. The bare pole, however, should produce a relevant energy dependence in the amplitudes above threshold and close to its position. ∼ −1.18 fm and −1.35 fm, respectively, which are larger (in absolute value) than those found here and in Ref. [46]. This is expected, since the bound states in Ref. [66] lie closer to the respective thresholds.
We now pay attention to the extra poles found in the SRS, located well above (∼ 400-500 MeV) their respective thresholds. From the very beginning we should take these results with some caution, since most likely they should be affected by sizable NLO and higher thresholdchannels corrections. In other words, they are not as theoretically robust as those concerning the lowest-lyingB s1 andB * s0 states. These resonances, likely, are originated from the bare bs-quark-model poles that are dressed by theB ( * ) K meson loops. In that case, the bare pole has been highly renormalized, moving to significant higher masses (∼ 6.2 − 6.3 GeV) and acquiring a significant B ( * ) K width (∼ 70 − 80 MeV). We should also bear in mind that radial excitations of the CQM states orB ( * ) K * two-meson loops, neither of them taken into account in this study, might lie in this region of energies. Further theoretical and experimental studies will be helpful in shedding light on the possible existence and properties of these resonant states.

IV. CONCLUSIONS
We have adopted a chiral unitary approach, based on leading-order HMχPTB ( * ) K interactions, and for the very first time in this context, the two-meson channels have been coupled to the corresponding CQM P-wavē B s scalar and axial mesons using an effective interaction consistent with HQSS. We have examined two different sets of masses for the bare quark model poles, and in each case, successfully fitted the rest of parameters to the recently reported LQCD isoscalar bs 0 + and 1 + energylevels [46]. Results turned out to be rather independent of the bare masses, showing that the changes can be easily re-absorbed by means of reasonable variations of the UV regulator.
We have focused on the scalar and axialB * s0 andB s1 states, which form a HQSS j P q = 1/2 + meson doublet, heavy-flavor partner of that in the charmed-strange sector integrated by the D * s0 (2317) and D s1 (2460) resonances. We have searched for bound states (poles in the FRS) of the isoscalar S-waveB ( * ) K amplitudes and found masses of 5709 ± 8 MeV (0 + ) and 5755 ± 8 MeV (1 + ) for these states. Therefore, theB * s0 andB s1 appear well below theBK andB * K thresholds, being in this way the lowest-lying mesons with these quantum numbers and stable through strong interactions. These states are clear candidates for experimental search in the LHCb experiment, B−factories or future high-luminosity proton-antiproton experiments.
We have also predicted the isoscalar elastic S-waveBK andB * K scattering lengths to be similar and approximately equal to −0.89 ± 0.07 fm, in good agreement with the findings of Ref. [46].
We have obtained sensible UV cutoff values ∼ 700 MeV, which do not hide large higher order contributions. In addition, and within the renormalization scheme adopted in this work, we have determined the dimensionless LEC c that controls the strength of the coupling between theB ( * ) K meson pairs and the P-wave bare quark model states. This LEC, though it depends on the orbital angular momentum and radial quantum numbers of the CQM state, is in principle independent of both the heavy-quark flavor and the light SU(3) flavor structure of the vertex. Thus, up to Λ QCD /m Q corrections, it can for example be also used to address the interplay between meson-loops and CQM degrees of freedom in the case of the D * s0 (2317) and D s1 (2460) resonances. Moreover, the same LEC governs the interplay between two-meson and quark model degrees of freedom in all isospin and strangeness channels. Next, we have looked at the Weinberg compositeness condition. Thanks to this admixture between CQM and two-meson degrees of freedom, we could realistically estimate the molecular component (B ( * ) K) of theB * s0 andB s1 , which turned out to be of the order of 50%. This is a clear indication of the non-canonical quark model nature of these states. Finally, we have further predicted the volume dependence of the isoscalar bs 0 + and 1 + energy-levels, which could be useful for future LQCD simulations.