2023 Near-threshold structures in the D + s D − s mass distribution of the decay B + → D + s D − s K +

Two near-threshold peaking structures with spin-parities J PC = 0 ++ were recently discovered by the LHCb collaboration in the D + s D − s invariant mass distribution of the decay B + → D + s D − s K + . The ﬁrst of them is the resonance X (3960), whereas the second one X 0 (4140) is a structure with the mass around 4140 MeV. To explore their natures and model them, we study the hadronic molecule M = D + s D − s and calculate its mass, current coupling, and width. The mass and current coupling of the molecule are extracted from the QCD two-point sum rule analyses by taking into account vacuum condensates up to dimension 10. To evaluate its full width, we consider the processes M → D + s D − s and M → η c η ( ′ ) . Partial widths of these decays are determined by the strong couplings g i , i = 1 , 2 , 3 at vertices M D + s D − s and M η c η ( ′ ) . They are computed by means of the three-point sum rule method. Predictions for the mass m = (4117 ± 85) MeV and width Γ M = (59 ± 12) MeV of the molecule M are compared with the corresponding LHCb data, and also with our results for the diquark-antidiquark state X = [ cs ][ cs ]. We argue that the structure X 0 (4140) may be interpreted as the hadronic molecule D + s D − s , whereas the resonance X (3960) can be identiﬁed with the tetraquark X .


I. INTRODUCTION
Different X resonances discovered and studied during last years by the LHCb collaboration became unavoidable and important part of exotic hadron spectroscopy. Thus, the resonances X(4140), X(4274), X(4500) and X(4700) were seen in the process B + → J/ψφK + as peaks in the J/ψφ invariant mass distribution [1]. The structures X(4140) and X(4274) bear the quantum numbers J PC = 1 ++ , and in the tetraquark model are composed of the quarks ccss. The resonances X(4500) and X(4700) are scalar particles with spin-parities J PC = 0 ++ and the same quark content. It is worth noting, that X(4140) and X(4274) were observed previously in the decays B ± → J/ψφK ± by different collaborations [2][3][4], whereas the scalar resonances X(4500) and X(4700) were discovered for the first time by LHCb. The resonance X(4630) fixed in the J/ψφ invariant mass distribution of the decay B + → J/ψφK + is a vector member of the X tetraquarks' family [5].
The X 0 (4140) may be interpreted either as a new resonance with J PC = 0 ++ assignment or J/ψφ ↔ D + s D − s coupled-channel effect [6]. In the present work, we assume that the structure X 0 (4140) is a second resonance seen by LHCb in the D + s D − s mass distribution. The X resonances are interesting objects for theoretical investigations: Features of exotic mesons ccss were studied in numerous articles by employing different models and technical methods [7][8][9][10][11][12][13][14]. Some of these states was undergone to rather detailed exploration also in our publications. Thus, the axial-vector resonances X(4140) and X(4274) were analyzed in Ref. [12], in which they were modeled as states composed of scalar and axial-vector (anti)diquarks. In the case of X(4140) the constituent diquark (antidiquark) is the antitriplet (triplet) state of the color group SU c (3), whereas to model X(4274) we used (anti)diquarks from the sextet representation of SU c (3). Predictions for masses and full widths of these tetraquarks were compared with the LHCb data. It turned out, for masses of the resonances X(4140) and X(4274) these models led to nice agreements with the LHCb data. The model based on color triplet (anti)diquarks reproduced also the full width of X(4140), therefore it could be considered as a serious candidate to resonance X(4140). The width of the construction with color sextet constituents is wider than that of X(4274), which excludes it from a list of possible pretenders.
The discovery of structures X(3960) and X 0 (4140) activated investigations of hidden charm-strange resonances to account for features of these new states [15][16][17][18][19][20]. The X(3960) was considered as a coupled-channel effect [15], or as near the D + s D − s threshold enhancement by the conventional P -wave charmonium χ c0 (2P ) [20]. The hadronic D + s D − s molecule model was suggested in Ref. [17] to explain observed properties of the resonance X(3960).
In our paper [21], we examined the tetraquark X = [cs][cs] with quantum numbers J PC = 0 ++ and calculated its mass and full width. The spectroscopic parameters of X, i.e., its mass and current coupling, were found by means of the QCD two-point sum rule method. The full width of this state was evaluated using decay channels X → D + s D − s and X → η c η (′) . Predictions for the mass m = (3976 ± 85) MeV and width Γ X = (42.2 ± 8.3) MeV obtained in Ref. [21] allowed us to consider the diquark-antidiquark state X as an acceptable model for X(3960).
In the present article, we continue our studies of the resonance X(3960) and include into analysis also the structure X 0 (4140). We investigate the hadronic molecule M = D + s D − s by computing its mass, current coupling, and full width. The mass and current coupling of this state are calculated in the context of QCD two-point sum rule approach. The full width of M is estimated by considering the decay channels M → D + s D − s , M → η c η and M → η c η ′ . Partial widths of these processes, apart from parameters of the initial and final particles, depend also on strong couplings g i , i = 1, 2, 3 at vertices MD + s D − s , Mη c η ′ , and Mη c η, respectively. To extract numerical values of g i , we use QCD threepoint sum rule method. Predictions for the mass and width of the molecule M are compared with the LHCb data for the resonance X(3960) and structure X 0 (4140). They are also confronted with parameters of the diquarkantidiquark state X = [cs][cs]. This paper is organized in the following way: In Sec. II, we compute the mass and current coupling of the molecule M by employing the QCD two-point sum rule method. The processes M → D + s D − s and M → η c η (′) are considered in Sec. III, where we determine the couplings g i and partial width of these decays. The full width of M is evaluated also in this section. The Sec. IV is reserved for our conclusions. The mass and current coupling of the molecule M = D + s D − s can be evaluated using the QCD two-point sum rule method [22,23]. This approach works quite well not only for analysis of conventional particles, but leads to reliable predictions also in the case of multiquark hadrons.
The principal quantity in this method is an interpolating current for a hadron under analysis. In the case of the molecule D + s D − s , this current J(x) has a rather simple form where a, and b are color indices. This current belongs to [1 c ] sc ⊗[1 c ] cs representation of the color group SU c (3). It corresponds to a molecule structure with spin-parities J PC = 0 ++ , but may couple also to different diquarkantidiquark structures and other four-quark hadronic molecules [21,24,25].
To find sum rules for the mass m and coupling f of the molecule M, one has to start from calculation of the following correlation function At the first stage, we express Π(p) in terms of the spectral parameters of M. To this end, it is necessary to insert a complete set of states with J PC = 0 ++ into the correlation function Π(p) and carry out integration over x. These operations lead to the simple formula The expression derived by this way is a hadronic representation of the correlator Π(p), which forms the phenomenological side of sum rules. The term written down explicitly in Eq. (5) is a contribution of the ground-state particle M, whereas contributions coming from higher resonances and continuum states are denoted by dots. The function Π Phys (p) can be rewritten in a more convenient form using the matrix element Then, it is not difficult to find Π Phys (p) in terms of the parameters m and f The Lorentz structure of Π Phys (p) has a simple form, and consists of a term proportional to I. Then, the invariant amplitude Π Phys (p 2 ) corresponding to this structure is given by expression in the right hand-side of Eq. (7). The QCD side of the sum rules Π OPE (p) is determined by Eq. (4) calculated using the c and s-quarks propagators. To this end, we insert the explicit form of J(x) into Eq. (4), contract heavy and light quark fields, and write obtained expression using quark propagators. After these operations, we get Here S c (x) and S s (x) are the c and s-quark propagators, explicit expressions of which can be found in Ref. [26].
The correlation function Π OPE (p) is calculated by employing the quark propagators with some fixed accuracy of the operator product expansion (OPE). The Π OPE (p) has a trivial Lorentz structure ∼ I as well. Having denoted the corresponding invariant amplitude by Π OPE (p 2 ) and equated it to Π Phys (p 2 ), we get a sum rule equality, which can be undergone further processing. The ground-state term and ones due to higher resonances and continuum states contribute to this sum rule equality on equal footing. There is a necessity to suppress unwanted contributions of higher resonances and subtract them from this expression. For these purposes, we apply the Borel transformation to its both sides. This operation suppresses effects of higher resonances and continuum states, but at the same time generates dependence of obtained equality on the Borel parameter M 2 . Afterwards, using the assumption about quark-hadron duality, we perform continuum subtraction, which leads to additional parameter s 0 in formulas.
The Borel transformation of the main term in Π Phys (p 2 ) has a simple form For the Borel transformed and subtracted amplitude Π OPE (p 2 ), we find (10) The first term in Eq. (10) contains an essential part of contributions, and is expressed using the two-point spectral density ρ OPE (s), derived as an imaginary part of the correlation function. The second term Π(M 2 ) collects nonperturbative contributions extracted directly from Π OPE (p). Formulas of the functions ρ OPE (s) and Π OPE (p) are lengthy, therefore we refrain from writing down them explicitly.
The sum rules for m and f read and where Numerical calculations of m and f should be performed in accordance with Eqs. (11) and (12), but only after fixing different vacuum condensates and working windows for the parameters M 2 and s 0 . The quark, gluon and mixed condensates are universal and well known quantities [22,23,[27][28][29]. Their numerical values, extracted from numerous processes are listed below We have included the masses of the c and s-quarks into Eq. (13) as well. The correlation function Π(M 2 , s 0 ) is calculated by taking into account vacuum condensates up to dimension 10. The expression of Π(M 2 , s 0 ) is rather lengthy, therefore we do not provide it here explicitly. In numerical analysis we set m 2 s = 0, but take into account contributions proportional to m s . To carry out numerical analysis one also needs to choose working regions for the Borel and continuum subtraction parameters M 2 and s 0 . They should satisfy standard constraints of sum rule calculations. Thus, the parameters M 2 and s 0 employed in calculations have to guarantee the dominance of the pole contribution (PC) and convergence of OPE. The former can be defined by the expression whereas to make sure the operator product expansion converges, we utilize the ratio where Π DimN (M 2 , s 0 ) is a sum of a few last terms in OPE. The R(M 2 ) and PC are used to restrict the lower and upper bounds for the Borel parameter, respectively. In fact, at M 2 min the function R(M 2 min ) should be less than some fixed value, whereas at M 2 max the pole contribution PC has to overshoot minimally acceptable limit for this parameter. In our present analysis, we impose on PC and R(M 2 ) the following constraints The first criterion in Eq. (16) is usual for investigations of conventional hadrons and ensures the dominance of the pole contribution. It may be used in analysis of multiquark hadrons, as well. But in the case of multiquark particles this constraint reduces a window for M 2 . The second condition is required to enforce the convergence of OPE. Dominance of a perturbative contribution to Π(M 2 , s 0 ) over a nonperturbative term, as well as maximum stability of extracted physical quantities when varying M 2 are among constraints to determine working regions for the parameters M 2 and s 0 . Numerical tests carried out by taking into account these aspects of the sum rule analysis demonstrate that the windows for M 2 and s 0 , The region for s 0 together with M 2 has to provide the dominance of PC and convergence of the operator product expansion. The parameter √ s 0 carries also useful information on a mass m * of the first radial excitation of the molecule D + s D − s . Thus, in the "groundstate +continuum" scheme adopted in the present work, √ s 0 should be less than the mass m * of the first excited state. This fact allows us to estimate low limit for m * ≥ m + 480 MeV.  Our result for the mass of the molecule M overshoots the LHCb datum m 1exp , but nicely agrees with m 2exp .

III. FULL WIDTH OF
The mass and current coupling of the molecule M calculated in the previous section provide information to select its possible decay modes. Besides, one should take into account its quantum numbers J PC = 0 ++ . Because the structures X(3960) and X 0 (4140) were discovered in the invariant mass distribution of the D + s D − s mesons, we consider M → D + s D − s as a dominant decay mode of M. The two-meson threshold for this decay is ≈ 3937 MeV, which makes M → D + s D − s kinematically allowed channel for M. Two other processes which we study here are decays M → η c η ′ and M → η c η. The two-meson threshold 3941 MeV for these decays is below the mass m of the molecule M as well. It is easy to prove that these decays conserve the spin, P and C parities of the initial particle M.
The partial width of the decay M → D + s D − s is governed by a coupling g 1 that describes strong interaction at the vertex MD + s D − s . This partial width depends also on masses and decay constants of the molecule M and mesons D + s and D − s . The mass and coupling of M have been calculated in the present article, whereas physical parameters of the mesons D + s and D − s are known from independent sources. Therefore, only physical quantity to be found is the strong coupling g 1 .
To determine g 1 , we employ the QCD three-point sum rule method, and begin our exploration from the corre-lation function Here, J(x), J D + s (y) and J D − s (0) are the interpolating currents for M, and pseudoscalar mesons D + s and D − s , respectively. The currents J D + s and J D − s are given by the expressions where i and j are color indices. The four-momenta of M and D + s are labeled by p and p ′ : Then, the momentum of the meson D − s is equal to q = p − p ′ . To find g 1 , we apply standard recipes of the sum rule method and calculate the correlation function Π(p, p ′ ). For these purpose, we use the physical parameters of the molecule M and mesons involved into this decay. The correlator Π(p, p ′ ) obtained by this manner forms the physical side Π Phys (p, p ′ ) of the sum rule. It is easy to see, that with m Ds being the mass of the mesons D ± s . To get Eq. (21), we isolate contributions of the ground-state and higher resonances and continuum state particles from each other. In Eq. (21) the ground-state term is written down explicitly, whereas other contributions are denoted by ellipses.
The function Π Phys (p, p ′ ) can be rewritten in terms of the D ± s mesons matrix elements where f Ds is their decay constant. We model the vertex Using these expressions, it is not difficult to calculate the new expression of the correlation function Π Phys (p, p ′ ) The double Borel transformation of the function Π Phys (p, p ′ ) over variables p 2 and p ′2 is given by the formula The correlator Π Phys (p, p ′ ) and its Borel transformation has simple Lorentz structure ∼ I. Then the whole expression in Eq. (24) determines the invariant amplitude Π Phys (p 2 , p ′2 , q 2 ). To find the QCD side of the three-point sum rule, we express Π(p, p ′ ) in terms of quark propagators, and get The correlator Π OPE (p, p ′ ) is computed by taking into account nonperturbative contributions up to dimension 6, and as Π Phys (p, p ′ ), contains the Lorentz structure proportional to I. Denoting relevant invariant amplitude by Π OPE (p 2 , p ′2 , q 2 ) , equating its double Borel transformation BΠ OPE (p 2 , p ′2 , q 2 ) to BΠ Phys (p 2 , p ′2 , q 2 ), and performing continuum subtraction, we get the sum rule for the coupling g 1 (q 2 ).
After these manipulations, the amplitude Π OPE (p 2 , p ′2 , q 2 ) can be rewritten using the spectral density ρ(s, s ′ , q 2 ), which is extracted as an imaginary part of Π OPE (p, p ′ ) where M 2 = (M 2 1 , M 2 2 ) and s 0 = (s 0 , s ′ 0 ) are the Borel and continuum threshold parameters. The couples (M 2 1 , s 0 ) and s 0 = (M 2 2 , s ′ 0 ) correspond to the molecule M and D + s meson channels, respectively. The sum rule for g 1 (q 2 ) is determined by the formula Ds /M 2 2 Π(M 2 , s 0 , q 2 ). (28) Expression of g 1 (q 2 ) depends on the spectroscopic parameters of the molecule M, and masses and decay constants of the mesons D ± s : they are input parameters of numerical computations. Values of these parameters as well as masses and decay constants f Ds and f ηc which are necessary to study other decays are collected in Table I. The masses all of mesons are borrowed from Ref. [30]. For decay constant of the mesons D ± s , we employ information from the same source, whereas for f ηc use a prediction made in Ref. [31] on the basis of the sum rule method.
The partial width of the decay M → D + s D − s besides various input parameters, also is determined by the strong coupling g 1 (m 2 Ds ) at the mass shell q 2 = m 2 Ds of the meson D − s . At the same time, the sum rule computations of g 1 can be carried out in deep-Euclidean region q 2 < 0. For simplicity, we introduce a variable Q 2 = −q 2 and in what follows label the obtained function by g 1 (Q 2 ). An explored range of Q 2 covers the region Numerical calculations also require choosing of the working regions for the Borel and continuum subtraction parameters Regions for M 2 and s 0 are chosen in such a way that to minimize their effects of the coupling g 1 (Q 2 ). Results of computations are pictured in Fig. 4. It is seen that results for g 1 (Q 2 ) are extracted at the region Q 2 > 0, where sum rule gives reliable predictions. It has just been explained above, that we need g 1 (Q 2 ) at Ds . To this end, one has to introduce some fit function G 1 (Q 2 ), which at the momenta Q 2 > 0 gives the same values as the sum rule computations, but can be easily extrapolated to the region Q 2 < 0. There are different choices for such functions. In this article, we use G i (Q 2 ), i = 1, 2, 3 where G 0 i , c 1 i and c 2 i are parameters, which should be extracted from fitting procedures. Calculations demonstrate that G 0 1 = 1.29 GeV −1 , c 1 1 = 2.38, and c 2 1 = −1.84  lead a reasonable agreement with the sum rule's data (see, Fig. 4).
At the mass shell q 2 = m 2 Ds this function gives The width of the decay M → D + s D − s is calculated by means of the formula where λ = λ (m, m Ds , m Ds ) and The decays M → η c η ′ and M → η c η are studied by a way described above. But, here we take into account peculiarities of the η − η ′ mesons connected with mixing in this system due to U (1) anomaly [32]. This effect modifies a choice of interpolating currents and matrix elements for these particles. Although η − η ′ mixing can be considered in different approaches, we use the quarkflavor basis |η q = (uu + dd)/ √ 2 and |η s = ss, where the physical particles η and η ′ have simple decompositions [32][33][34] η = |η q cos ϕ − |η s sin ϕ, η ′ = |η q sin ϕ + |η s cos ϕ.
where ϕ is the mixing angle in the {|η q , |η s } basis. Such state-mixing implies that the same assumption is applicable to their currents, decay constants and matrix elements. Because in the decays M → η c η ′ , η c η participate only ss components of the mesons η and η ′ , relevant interpolating currents are given by the expressions where j is the color index. We start our analysis from the process M → η c η ′ . Then, we should consider the correlation function where J ηc (y) is the interpolating current of the meson η c The main contribution to the correlation function Π(p, p ′ ) has the following form where the ellipses stand for contributions of higher resonances and continuum states. The function Π Phys (p, p ′ ) can be rewritten using the the matrix elements where m ηc and f ηc are the mass and decay constant of the η c meson. The matrix element of local operator siγ 5 s sandwiched between the meson η ′ and vacuum states is denoted by h s η ′ [33]. The parameter h s η ′ follows the η − η ′ state-mixing pattern, and we get The parameter h s can be defined theoretically [33], but for our purposes it is enough to use values of h s and ϕ extracted from phenomenological analyses The vertex M η c η ′ has the following form where g 2 is the strong coupling at the vertex M η c η ′ . By employing these matrix elements, we obtain a new expression for Π Phys (p, p ′ ) The QCD side of the sum rule for g 2 (q 2 ) is given by the formula The sum rule for the coupling g 2 (q 2 ) is obtained using Borel transformations of invariant amplitudes Π Phys (p 2 , p ′2 , q 2 ) and Π OPE (p 2 , p ′2 , q 2 ) and is equal to Here, Π(M 2 , s 0 , q 2 ) is the Borel transformed and subtracted amplitude Π OPE (p 2 , p ′2 , q 2 ). The coupling g 1 (q 2 ) is calculated using the following Borel and continuum threshold parameters in the η c channel whereas as M 2 1 and s 0 for the M channel, we use Eq. (17). The strong coupling g 2 is defined at the mass shell of the η ′ meson. The fit function G 2 (Q 2 ) given by Eq. (30) has the parameters G 0 2 = 0.21 GeV −1 , c 1 2 = 5.08, and c 2 2 = −4.04. Computations yield The partial width of this decay can be evaluated using the formula Eq. (32), in which one should make substitutions g 1 → g 2 , m 2 Ds → m 2 ηc and λ (m, m Ds , m Ds ) → λ (m, m ηc , m η ′ ). Then, for the width of the decay M → η c η ′ , we find Γ [M → η c η ′ ] = (4.9 ± 1.1) MeV. (49) Analysis of the decay M → η c η can be performed in a similar manner. Avoiding further details, let us write down predictions obtained for key quantities. Thus, the strong coupling g 3 at the vertex Mη c η is given by the formula The width of the hadronic molecule M, within errors of calculations and measurements, agrees with the LHCb datum from Eq. (2).

IV. CONCLUSIONS
In this work, we have calculated the mass and width the scalar molecule M = D + s D − s in the framework of the QCD sum rule methods. The mass of M has been evaluated using the two-point sum rule approach. The full width of M has been computed by taking into account decay modes M → D + s D − s , and M → η c η (′) . Strong couplings g i that determine width of these decays have been found in the framework of the three-point sum rule method.
Our result for the mass m = (4117 ± 85) MeV of the molecule M exceeds considerably the corresponding LHCb datum m 1exp , but is comparable with m 2exp . The full width Γ M = (59 ± 12) MeV of M is consistent with the LHCb measurement Γ 2exp as well. It is evident, that the molecule M is considerably heavier than the resonance X(3960), which makes problematic its interpretation as X(3960).
The resonance X(3960) was examined in our article [21] as the tetraquark X = [cs][cs] with quantum numbers J PC = 0 ++ . We obtained the following predictions for the parameters of X: the mass m = (3976 ± 85) MeV and width Γ X = (42.2 ± 8.3) MeV. The parameters of the diquark-antidiquark state X are in nice agreements with the LHCb data given by Eq. (1), therefore in Ref. [21], X was identified with X(3960).
In the context of the sum rule approach the molecule D + s D − s was also studied in Ref. [17]. In accordance with this paper, the mass of such hadronic molecule is equal to (3980 ± 100) MeV and agrees with the LHCb data. It is worth noting that the authors did not analyze quantitatively the width of this state. Our results for parameters of M, even within existing errors of calculations, does not support molecule assignment for X(3960).
By taking into accout predictions for the mass m and full width Γ M of the molecule M = D + s D − s obtained in the present work, we argue that the molecule M may be a candidate to the structure X 0 (4140).