Z_c -like spectra from QCD Laplace sum rules at NLO

We present a global analysis of the observed Z_c, Z_cs and future Z_css-like spectra using the inverse Laplace transform (LSR) version of QCD spectral sum rules (QSSR) within stability criteria. Integrated compact QCD expressions of the LO spectral functions up to dimension-six condensates are given. Next-to-Leading Order (NLO) factorized perturbative contributions are included. We re-emphasize the importance to include PT radiative corrections (though numerically small) for heavy quark sum rules in order to justify the (ad hoc) definition and value of the heavy quark mass used frequently at LO in the literature. We also demonstrate that, contrary to a na\"ive qualitative 1/N_c counting, the two-meson scattering contributions to the four-quark spectral functions are numerically negligible confirming the reliability of the LSR predictions. Our results are summarized in Tables III to VI. The Z_c(3900) and Z_cs(3983) spectra are well reproduced by the T_c(3900) and T_cs(3973) tetramoles (superposition of quasi-degenerated molecules and tetraquark states having the same quantum numbers and with almost equal couplings to the currents). The Z_c(4025) or Z_c(4040) state can be fitted with the D*_0D_1 molecule having a mass 4023(130) MeV while the Z_cs bump around 4.1 GeV can be likely due to the (D^*_s0D_1+ D^*_0D_s1) molecules. The Z_c(4430) can be a radial excitation of the Z_c(3900) weakly coupled to the current, while all strongly coupled ones are in the region (5634-6527) MeV. The double strange tetramole state T_css which one may identify with the future Z_css is predicted to be at 4064(46) MeV. It is remarkable to notice the regular mass-spliitings of the tetramoles due to SU(3) breakings M_{T_cs}-M_{T_c}= M_{T_css}-M_{T_cs= (73- 91) MeV.


I. INTRODUCTION
Beyond the successful quark model of Gell-Mann [1] and Zweig [2], Rossi and Veneziano have introduced the four-quark states within the string model [3] in order to describe baryon-antibaryon scattering, while Jaffe [4] has introduced them within the bag models for an attempt to explain the complex structure of the I = 1, 0 light scalar mesons (see also [5][6][7]).
In earlier papers, QSSR has been used to estimate the I = 0 light scalar mesons (σ, f 0 ,) masses and widths [8] assumed to be four-quark states. However, the true nature of these states remains still an open question as they can be well interpreted as glueballs/gluonia [9][10][11][12].
After the recent discovery of many exotic XYZ states beyond the quark model found in different accelerator experiments 1 , there was a renewed interest on the fourquarks and molecule states for attempting to explain the properties of these new exotic states 2 .
In previous works [22][23][24], we have systematically extracted the couplings and masses of the XY Z states using QCD spectral sum rules (QSSR)à la Shifman-Vainshtein-Zakharov (SVZ) [25,26] 3 where the N2LO factorized perturbative and QCD condensates up to dimension-six in the Operator Product Expansion (OPE) corrections have been included 4 . In so doing, we have used the inverse Laplace transform (LSR) [38][39][40] version of QSSR within stability criteria where we have emphasized the importance of the PT corrections for giving a meaning on the input heavy quark mass value which plays an important role in the analysis, though these corrections are numerically small, within the M S-scheme.
More recently, we have applied the LSR for interpreting the new states around (6.2-6.9) GeV found by the LHCb-group [41] to be a doubly/fully hiddencharm molecules (QQ)(QQ) and (QQ)(QQ) tetraquarks states [42], while the new states found by the same group from the DK invariant mass [43] have been interpreted by a 0 + and 1 − tetramoles (superposition of almost degenerate molecules and tetraquark states having the same quantum numbers and couplings) slightly mixed with their radial excitations [44]. In this paper, we pursue the analysis using LSR by studying the recent data from BESIII where the K + recoil or invariant D * s D ⊕ D * D s mass (see Fig. 1) [45] is a good Z cs (1 + )-like state candidate. A narrow peak is experimentally found at (in units of MeV): We shall be concerned with the QCD local currents O µ H (x) of dimension-six given in Table I for 1 + axialvector molecules and tetraquarks where q ≡ u, s. b is a free mixing parameter where its optimal value was found to be zero [23,24]. The appropriate 1 + hadron H couples to the current as : where f H is the hadron decay constant analogue to f π and ϵ µ is the axial-vector polarisation. In general, the four-quark operators mix under renormalization and acquire anomalous dimensions [46]. In the present case where the interpolating currents are constructed from bilinear (pseudo)scalar currents, the anomalous dimension can be transfered to the decay constants as : where :f H is the renormalization group invariant coupling and −β 1 = (1/2)(11 − 2n f /3) is the first coefficient

B. Form of the sum rules
We shall work with the Finite Energy version of the QCD Inverse Laplace sum rules (LSR) and their ratios : dt t n e −tτ 1 π Im Π (1) where m c is the charm quark mass, τ is the LSR variable, n = 0, ... is the degree of moments, t 0 is the quark/hadronic threshold. t c is the threshold of the "QCD continuum" which parametrizes, from the discontinuity of the Feynman diagrams, the spectral function Im Π (1) H (t, m 2 c , µ 2 ) where Π (1) H (t, m 2 c , µ 2 ) is the transverse scalar correlator corresponding to a spin one hadron :

III. QCD TWO-POINT FUNCTION
Using the SVZ [25] Operator Product Expansion (OPE), we give in the Appendix the QCD expressions to lowest order (LO) of the two-point correlators associated to the currents given in Table I up to dimension d = 6 condensate contributions.
A. LO perturbative (PT) and 1/Nc counting Motivated by the criticisms raised in Ref. [47] based on the large N c -limit which state that the non-factorized contribution of the two-point correlator starts at order α 2 s and that the non-resonant (scattering states) dominate the sum rules, we check explicitly these statements for finite N c here and in the following subsection C, which we complete in Section X by comparing the D * D molecule resonance and the non-resonating D * and D channels contributions to the sum rule.
The LO perturbative contributions are given by the diagrams in Fig. 2. Explicit evaluations of the Trace appearing in the two-point function indicates that the lowest order (LO) PT contribution behaves like N 2 c (or 2N 2 c (1 − 1/N c ) if the current has an ϵ-tensor like the one in Table I where the 1/N c term arises from the ϵcontraction) as expected from large N c .
However, non-factorised contribution appears at LO both from PT and condensate contributions when one has two or more identical quark flavours because one has more possibilities to do the Wick's contraction 5 . Moreover, some care has to be taken when applying the large N c analysis to the case of baryons and tetraquark states with string junctions [3]  The QCD condensates entering in the analysis are the light quark condensate ⟨qq⟩ and the SU (3)breaking parameter κ ≡ ⟨ss⟩/⟨qq⟩, the gluon condensates ⟨α s G 2 ⟩ ≡ ⟨α s G a µν G µν a ⟩ and ⟨g 3 G 3 ⟩ ≡ ⟨g 3 f abc G a µν G b νρ G c ρµ ⟩, the mixed quark-gluon condensate g⟨qGq⟩ ≡ ⟨qgσ µν (λ a /2)G a µν q⟩ = M 2 0 ⟨qq⟩ and the fourquark condensate ρ⟨qq⟩ 2 , where ρ ≃ (3 ∼ 4) indicates the deviation from the four-quark vacuum saturation. Their different contributions within the SVZ expansion to LO are shown in Figs. 4 to 8.
Unlike often used in the literature, we have not included higher dimension d ≥ 8 condensates contributions 7 due to our poor knowledge of their size. Indeed, a violation of the vacuum saturation for the four-quark condensates [48][49][50][51][52] has been already noticed in different light quark channels. In addition, the mixing of different four-quark condensates under renormalization [46] does not also favour the vacuum saturation estimate. On the other, the inaccuracy of a simple dilute gas instanton estimate [25,53] has been also observed from the phenomenological estimate of high-dimension gluon condensates [54][55][56].

C. Convolution representation and Matching
We have explicitely proved in our previous papers [22-24, 42, 44] that the non-factorized contribution to the four-quark correlator which appears at lowest order α 0 s of PT QCD to order in 1/N c (but not to order α 2 s as claimed by [47]), gives numerically negligible contribution to the sum rule. Therefore, we can consider that the molecule /tetraquark two-point spectral function is well approximated by the convolution of the two ones built from two quark bilinear currents (factorization) as illustrated in Fig. 9 where the diagrams in the two sides of Fig. 9 are of the order N 2 c . In order to fix the matching factor k 2 , we consider the example of the D * D molecule current where the QCD expression is given in the Appendix. The bilinear currents and the corresponding spectral functions entering in the RHS of Fig. 9 are : in the limit where m 2 c ≪ t. In this way, we obtain, for a spin 1 state, the convolution integral [57][58][59]: with the phase space factor: (8) √ t 10 and √ t 20 are the quark / hadronic thresholds and m c is the on-shell / pole perturbative charm quark mass.
The appropriate k-factor which matches this convolution representation with the direct perturbative calculation of the molecule spectral function given in the appendix is 8 : which comes from the dynamics of the Feynman diagram calculations and which is missed in a standard large N cresults of [60]. Related phenomenology will be discussed in Section X.

D. NLO PT corrections to the Spectral functions
We extract the next-to-leading (NLO) perturbative (PT) corrections by approximating the molecule /tetraquark two-point spectral function with the convolution of the two ones built from two quark bilinear currents (factorization) illustrated in    The NLO perturbative expressions of the bilinear unequal masses (pseudo)scalar and (axial)-vector spectral functions are known in the literature [27,28,31,[61][62][63].
In the following, we shall use n f =4 total number of flavours for the numerical value of a s ≡ α s /π.

A. QCD coupling αs
We shall use from the M χ0c − M ηc mass-splitting sum rule [64]: which is more precise than the one from M χ 0b − M η b [64] : These lead to the mean value quoted in Table II, which is in complete agreement with the world average [65] :

B. Quark masses
We shall use the recent determinations of the running masses m s (µ) and m c (m c ) quoted in Table II and the corresponding value of α s evaluated at the scale µ obtained using the same sum rule approach.

C. QCD condensates
Their values are quoted in Table II. One should notice that taking into account the anomalous dimension of the ⟨qq⟩ condensate, α s ⟨qq⟩ 2 has a very smooth 1/ log (τ Λ 2 ) 1/25 behaviour for n f = 4 flavours such that, contrary to simple minded, its contribution is not suppressed by 1/ log (τ Λ 2 ) for light mesons [48,49] and τdecays [50] permitting its reliable phenomenological estimate while its extraction from light baryons [36,51,52] occurs without the α s factor like in the case of the fourquark currents discussed here.

V. THE SPECTRAL FUNCTION
A. The minimal duality ansatz In the present case, with no complete data on the spectral function, we use the minimal duality ansatz: for parametrizing the molecule spectral function. M H and f H are the lowest ground state mass and coupling analogue to f π . The "QCD continuum" is the imaginary part of the QCD correlator (as mentioned after Eq. 4) from the threshold t c which is assumed to smear all higher states contribution. Then, it insures that both sides of the sum rules have the same large t asymptotic behaviour which leads to the Finite Energy Sum Rule (FESR) in Eq. 4. Within a such parametrization, one obtains: indicating that the ratio of moments is a useful tool for extracting the mass of the hadron ground state [27][28][29].
The corresponding value of t c approximately corresponds to the mass of the 1st radial excitation. However, one should bear in mind that a such parametrization cannot distinguish two nearby resonances but instead will consider them as one "effective resonance".
This simple model has been tested successfully in different channels where complete data are available (charmonium, bottomium and e + e − → I = 1 hadrons) [27,28,33]. It was shown that, within the model, the sum rule reproduces quite well the integrated data while the masses of the lowest ground state mesons (J/ψ, Υ and ρ) have been predicted within a good accuracy.
In the extreme case of the pseudoscalar Goldstone pion, the sum rule using the spectral function parametrized by this simple model and the more complete one by Chiral Perturbation Theory (ChPT) [76] lead to similar values of the sum of light quark masses (m u + m d ) indicating the efficiency of this simple parametrization [27,28].
An eventual violation of the quark-hadron duality (DV) [77,78] tested from hadronic τ -decay data [50,78,79] is negligible here thanks to the double exponential suppression of this contribution in the Laplace sum rule (see e.g [42] for details).
The uses of this model for the tetraquarks and molecule states are also quite successful compared to the recent data (see e.g [42,44] and references therein). Then, we (a priori) expect to extract with a good accuracy the masses and couplings of the mesons within the approach.
In order to minimize the effects of radial excitations smeared by the QCD continuum, we shall work with the lowest moment L c 0 and ratio of moments R c 0 for extracting the meson masses and couplings f H . Moments with n < 0 will not be considered due to their sensitivity on the non-perturbative contributions at zero momentum.
However, once we have fixed the ground state parameters, we attempt to extract the mass and coupling of the first radial excitation by using a :"two resonance" +(Θ(t − t c1 ) "QCD continuum " parametrization where t c1 is above the t c -value obtained for the ground state.

B. Optimization Criteria
As τ, t c and µ are free external parameters , we shall use stability criteria (minimum sensitivity on the variation of these parameters) to extract the hadron masses and couplings. Results based on these stability criteria have lead to successful predictions in the current literature (see [27][28][29] and original papers).

VI. REVISITING fD * D AND MD * D
We start by revisiting and checking the results obtained in [24] where they have included the factorized contributions to N2LO of perturbative series. In this example, we show explictly our strategy for extracting the mass and coupling. The same strategy will be used in some other channels discussed later in this paper. This example is also a test of the efficiency of the method by confronting the prediction with the : Z c (3900) [80,81], Z c (4020) [82], Z c (4025) [83], Z c (4050) [84], Z c (4226), Z c (4257) by BESIII [85] and Z c (4430) [86,87] found earlier, where one can notice that only the Z c (3900) and Z c (4430) have been retained as established in the Meson Summary Table of Particle Data Group (PDG) [65].  Table II.

A. τ and tc-stabilities
We show in Fig. 11 the τ -behaviour of the coupling and of the mass for different values of t c and fixing the value of the subtraction constant µ at 4.65 GeV (see next subsection) where a µ stability (4.5 ± 0.5) GeV has been found in [24]. From this figure, one can see that the coupling presents minimas in τ and the mass inflexion points.
-In a first step, we use the experimental mass M Zc = 3900 MeV for extracting the value of the coupling f D * D shown in Fig. 11a).
-In a 2nd step, we take the value of τ at the minimum of the coupling and use it for extracting the value of M D * D from Fig. 11b).
-In a 3rd step, we take the common range of t c where both curves present stabilities in τ . In the present case, this value ranges from t c = 22 GeV 2 (beginning of τstability) to t c = 38 GeV 2 (beginning of t c -stability) where the range is given in Table IV. For the mean t c = 30 GeV 2 , it is : τ R ≃ 0.38(resp. 0.34) GeV −2 for LO (resp. NLO) QCD expression.
-The errors given in Table III Table VII are the mean of the ones from the previous two extremal values of t c .

B. µ-stability
For doing the analysis, we shall fix t c = 30 GeV 2 which is the mean of the two extremal values delimiting the stability region and take M D * D = 3900 MeV for fixing the coupling. The analysis is shown in Fig. 12. One can see a µ-stability for : at which we shall evaluate the results quoted in Table III.
This value of µ from a more refined analysis is more precise than the conservative one (4.5 ± 0.5) GeV quoted in [24]. The results of the analysis are given in Tables III  and VII. One should note that the resummed QCD expression of the sum rule which obeys an "Homogeneous" Renormalization Group Equation (RGE) is obtained by putting µ 2 = 1/τ in the QCD expression of the sum rule and where the parameters having anomalous dimension γ run as 1/(log τ Λ 2 ) γ/β (see e.g. [40]). We have often used this choice in the past (see [27][28][29]) which corresponds here to the value: at the τ -minimum for t c =30 GeV 2 . However, this value is outside the µ-stability region obtained previously and then does not correspond to the optimal choice of µ. This is the reason why we have abandoned this choice µ 2 = 1/τ . This result does not support the argument of Ref. [88] based on the observation that, when the term of the type (α s (τ )/α s (µ)) γ/β1 disappears for µ = 1/ √ τ , one would obtain the best choice of µ. Indeed, the PT series behaves obviously much better in the µ-stability region where the value of µ is about 3 times higher at which the radiative corrections are more suppressed.
Moreover, the physical meaning of the relation between µ with the so-called bound energy or virtuality µ 2 = M 2 Z − 4M 2 c used in his different papers (see e.g. [88,89]) remains unclear to us where M c is the PT constituent or pole charm quark mass taken by the author to be about 1.84 GeV which corresponds to µ ≃ 1.3 GeV.
Indeed, one may expect from this formula that the difference between the resonance and PT quark constituent masses has a non-perturbative origin which (a priori) has nothing to do with the scale µ where the PT series and the Wilson coefficients of the condensates are evaluated.
However, one may also consider µ as a scale separating the calculable PT Wilson coefficients and the NPT noncalculable condensates in the OPE [25] though one has to bear in mind that the d ≤ 6 condensates appearing in the present analysis are renormalization group invariant (µ independent like m c ⟨qq⟩, m 2 c ⟨qq⟩ 2 , ) or have a weak dependence on µ (⟨α s G 2 ⟩, ⟨gqGq⟩) [27] (Part VII page 285) such that the truncation of the PT series does not affect much their values. This feature indicates that the separation of the condensates from the PT Wilson coefficients are not ambiguous while the size of the non-perturbative condensate is almost independent on the scale at which the PT series is truncated.

C. NLO and Truncation of the perturbative series
We have mentioned in previous works that the inclusion of the NLO perturbative corrections is important for justifying the choice of heavy quark mass definition used in the analysis where an ad hoc value of the M S running mass is frequently used in the literature while the spectral function has been evaluated using an on-shell renormalization where the pole (on-shell) quark mass naturally enters into the LO expression.
We show the analysis for the mass and coupling in Fig. 13 at fixed value of the continuum threshold t c and subtraction constant µ, where one can find that the use of the pole mass at LO decreases by 30% at the minimum (resp. increases by 0.5% at the inflexion point) the value of the coupling (resp. mass) obtained using the M S running mass at LO while the NLO correction is relatively small within the M S-scheme.
The smallness of radiative corrections for the ratio of moments demonstrates (a posteriori) why the use of the M S running mass at LO leads to a surprisingly good prediction for the mass.
Assuming that the PT series grow geometrically [90][91][92][93], one can deduce from the previous analysis an estimate of higher order PT contributions given in Table III which can be compared with the ones in [23,24].  Table II. obtains for t c =30 GeV 2 : One can notice the important role of the dimension -4 and -6 condensates which are dominated by the chiral condensates ⟨qq⟩ and ρ⟨qq⟩ 2 in the extraction of the coupling and mass where their strengths are more pronounced for the coupling.
We estimate the sytematic errors due to the truncation of the OPE from the size of the dimension-6 condensate contributions rescaled by the factor m 2 c τ /3 where 1/3 comes from the LSR exponential form of the sum rule. It can be compared with the contributions of the known d = 8 ⟨qq⟩⟨qGq⟩ condensate obtained in [23,24] but bearing in mind that this is only a part of the complete d=8 condensate ones where the validity of the vacuum saturation used for its estimate is questionable.
The quoted errors in Table III are the mean from the two extremal values of t c delimiting the (τ, t c ) stability region. We expect that the systematics for the m s ̸ = 0 cases are similar as confirmed explicitly in Table III.
The values of mc(mc) and the condensates are given in Table II. FAC means that we use the vacuum saturation assumption for the estimate of the four-quark condensates.

E. Vacuum saturation of the four-quark condensate
The four-quark condensates have been demonstrated to mix under renormalization [46] and [27] (Part VII page 285) at finite N c while its vacuum saturation estimate is only valid in the large N c -limit. A such estimate from light mesons [48,49], light baryons [36,51,52] and τdecays [50] has been shown to underestimate the actual value of the four-quark condensates by a factor 3-4 (see also the comments in Subsection IV C). For completeness, we also show in Fig. 14 the effect of this estimate, where we see that the position of the minimum and inflexion point are shifted at τ ≈ 0.4 GeV −2 . At this value and for t c = 30 GeV 2 , one obtains an effect of +38% for the coupling and -7% for the mass leading to : (18) which can be compared with the one in Eq. 17. One can notice that the results are in perfect agreement with the previous ones in [24]. The slight difference in the error calculation is due to the fact that, in Ref. [24], we have estimated the error by choosing t c = 38 GeV 2 but not considering the one due to t c = 22 GeV 2 . The result for the mass is in a very good agreement with the BELLE [80] and BESIII [81] data Z c (3900) MeV which we shall discuss later on.
Here, we also revisit the estimate of the D * 0 D 1 , A cd masses and couplings done in [24]. We repeat exactly the same procedure as in the previous section.
The τ and t c -behaviours of the coupling and mass are similar to the previous case and will not be shown here. The τ -stability region ranges from τ = 0.21 GeV −2 for t c = 28 GeV 2 (beginning of τ -stability) to 0.30 GeV −2 for t c = 40 GeV 2 (beginning of t c -stability). One can notice that the τ -stability starts at a larger value of t c than the one t c = 22 GeV 2 for the case of D * D which will imply a larger value of M D * 0 D1 than of M D * D . The µ-stability is shown in Fig. 15 where the optimal value is the same as in Eq. 15. The result : (20) differs with the one given in [23,24] (see Table VII) originated from the unprecise value of the τ used there for extracting the optimal value which affect in a sensible way the value of the mass in this channel. The τ and t c -behaviours of the A cd -mass and coupling are similar to the one in Fig 11. The τ -stability starts  TABLE III. Sources of errors and predictions from LSR at NLO for the couplings and masses of the molecules and tetraquark ground states. The errors from the QCD input parameters are from Table II. We take |∆τ | = 0.01 GeV −2 and ∆µ = 0.05 GeV. The quoted errors are the mean from the two extremal values of tc delimiting the stability region quoted in Table IV.  for t c =22 GeV 2 at τ = 0.25 GeV −2 while the t c -stability starts at t c = 38 GeV 2 where τ = 0.36 GeV −2 .
The µ-behaviours are the same as the ones of D * 0 D 1 where the µ-stability is also at 4.65 GeV as in Eq. 15.
The estimates of the errors and the results for the coupling and mass : (17) keV M A cd = 3889(58) MeV (21) are collected in Tables III and VII. The mass value is in good agreement with the one 3888(130) MeV obtained in Ref. [24] but with a smaller error due to a better localisation of the τ and µ stability points. The A cd mass also coincides with the observed Z c (3900). The different sets of (t c , τ ) used to get the previous optimal results are summarized in Table IV. The τ and t c -behaviours of the coupling and mass are also similar to the previous cases as shown in Fig. 16. The τ -stability region ranges from τ = 0.22 GeV −2 for t c = 22 GeV 2 (beginning of τ -stability) until 0.36 GeV −2 for 38 GeV 2 (beginning of t c -stability).
The µ-stability is shown in Fig. 17 where the optimal value is the same as in Eq. 15.
The sources of the errors and the results are quoted in Table III Tables III and VII. C. The D * s0 D1 and D * 0 Ds1 molecules Similar analysis leads to (see Tables III and VII) : (24) where one can notice that the two molecule states are almost degenerated and have the same couplings to the currents.

D. The Acs tetraquark
We pursue the previous analysis for the A cs tetraquark. The behaviours of the different curves are similar to the previous ones and will not be shown. The result quoted in Tables III and VII is: where one can notice that it is almost degenerated to the D * D s and D * s D and has almost the same couplings.

IX. THE (1 + ) (cs)(cs) AND (cs)(cs) STATES
In this section, we revisit and improve our previous estimate of the masses and couplings of these abovementioned states and give a new estimate of the D * s D s radial excitation mass and coupling.
A. The D * s Ds molecule The analysis of the τ and t c behaviours is shown in Fig. 19. The µ-behaviour is shown in Fig. 20. These behaviours are similar to the previous ones. The (τ, t c ) stabilities are obtained for t c inside the range 24 to 40 GeV 2 . We deduce (see Tables III and VII) : In this section, we revise and check the results obtained in [23]. The behaviours of the τ, t c and µ-behaviours of the masses and couplings are similar to the previous cases and will not be shown.
For the D * s0 D s1 molecule, the set of (τ, t c )-stabilities are obtained from (0.20,28) to (0. 30,44) in units of (GeV −2 , GeV 2 ), where one notice the sensitivity of the results in the change of τ which is quantified by the large error induced by the variation of τ as shown in Table III.  One obtains from Tables III and VII : which agree within the errors with the previous results in [23].
For the A css tetraquark, τ -stability at 0.27 GeV −2 starts for t c = 24 GeV 2 while the t c -stability for τ = 0.38 GeV −2 starts for t c = 40 GeV 2 . These values are about the same as the ones for A cd and A cs indicating like in the case of molecule states a good SU (3) symmetry for the values of the LSR parameters. The optimal results are given in Tables III and VII which read: f Acss = 114 (16) keV , M Acss = 4014 (86) MeV . (28) Compared to the previous non-strange case where the molecule D * D is quasi-degenerated with the tetraquark A cd , which leads us to conclude that the observed state is a tetramole T c , in this channel, we find the tetraquark A css is almost degenerated (within the errors) with the molecule D * s D s and has the same couplings implying that it can also be a tetramole T css state with a coupling and mass: The D * s0 D s1 molecule is about 100 MeV slightly higher but has a weaker coupling to the corresponding current.

X. TWO-MESON SCATTERING STATES
To study this contribution, we take the example of the D * D molecule 9 .
We saturate the RHS of Fig. 9 by the two non-resonant (scattering) states D * and D which we shall compare with the one due to the molecule D * D (LHS) appropriately matched with the k 2 -factor.
Then, one obtains the two-resonance scattering LSR moment : where I(t, M 2 D , M 2 D * ) is the integral (see Eq. 6): to be compared with the molecule sum rule result :  [28,94]. We have neglected 9 Our conclusion remains valid for some other molecules and tetraquarks states.
the equal and small contributions to the two sum rules from the QCD continuum above the continuum threshold taken to be : t c ≈ 30 GeV 2 which is the mean of the two extremal values delimiting the stability region. These results indicate that : -For finite N c , the non-resonant contribution is about one order of magnitude smaller than the one of the resonance molecule (a similar conclusion using an alternative approach has been reached in [88]) and disprove the claim of Ref. [60] based on large N c once an appropriate matching of the two correlators via the k 2 factor is done.
-A posteriori, the existence of the stability region or "sum rule window" in the LSR analysis where one can extract the (postulated) resonance mass and coupling using the spectral function parametrization within the duality ansatz "one resonance" + QCD continuum is a strong indication of the duality between the QCD-OPE and the phenomenological side of the sum rule, where the resonace contribution is dominant over the one of the nonresonant states. This fact does not also support the claim of Ref. [60].

XI. RADIAL EXCITATIONS
We present in this section a new estimate of the couplings and masses of the first radial excitations using the lowest moments L c 0 and R c 0 10 . In so doing, we use a "two resonance" parametrization of the spectral function.
A. The (D * D)1 first radial excitation We insert the previously obtained values of the lowest ground state D * D mass and coupling and study the τand t c -stability of the sum rules by fixing the subtraction constant µ as in Eq. 15. The analysis is shown in Fig. 21 where the stability region is delimited by t c = 39 and 50 GeV 2 to which corresponds respectively the τ -stability of (0.14-0. 18 Table V. The set of (t c , τ ) is compiled in Table VI.
One can notice that the mass value is roughly about the (expected) one of √ t c ≃ (4.7 − 6.2) GeV inside the stability region where the lowest ground state mass has been extracted. 10 We note that the moment R c 1 which is expected to be more sensitive to the radial excitation contribution has a behaviour similar to R c 0 and will not presented here.  TABLE V. Sources of errors and predictions from LSR at NLO and for the decay constants and masses of the molecules and tetraquark radial excitation states. The errors from the QCD input parameters are from Table II. We take |∆τ | = 0.01 GeV −2 and ∆µ = 0.05 GeV. The quoted errors have been estimated fixing tc at the mean of the two extremal values delimiting the stability region quoted in Table VI. We notice that the relative large values of some individual errors in the case of D * s Ds and (D * 0 D1)1 are mainly induced by the shift of the position of the minima compared to the one of the central values.   We show the analysis in Fig. 22. The curves have similar behaviour as in the case of (D * D) 1 but the stabilities are reached for higher values of t c which imply a higher value of the (D * 0 D 1 ) 1 radial excitation mass. The results are quoted in Table V and the set of (t c , τ ) is compiled in Table VI. C. The (A cd )1 radial excitation The behaviours of the coupling and mass versus τ and t c are similar to the case of the one of (D * D) 1 and will not be repeated here. The results are quoted in Table V. In this subsection, we study the first radial excitation (D * s D s ) 1 . The τ and t c behaviours of its mass and coupling are shown in Fig. 23 for µ = 4.65 GeV. The τ and t c -stabilities are reached for the set (τ, t c ) =(0.28,42) to (0.30,54) in units of (GeV −2 , GeV 2 ). We deduce (see Table V One can notice that the coupling of the 1st radial excitation is similar to the previous radial excited states which are relatively large compared to the ones of lowest ground states. This feature differs from the case of ordinary mesons built from bilinear currents. The mass is also found to be relatively high.  Table II.

F. Comments on the radial excitations
We have shown previously that the couplings of the excited states to the corresponding currents are as large as the one of the ground states (see Table V) which is a new feature compared to the case of ordinary hadrons.
We have also shown that the mass-splittings between the first radial excitation and the lowest ground state are (see Table V) : which is much bigger than the one (500 ∼ 600) MeV for ordinary mesons but comparable with the one obtained for the DK state in [44]. The authors in Ref. [89] have also noticed a such anomalously large value of the mass of the radial excitation M R1 ≈ √ t c ≥ 4.5 GeV in the analysis of the Z c (4020, 4025). Then, they have concluded that the Z c (4020, 4025) cannot be a tetraquark state as the value of t c used to extract these masses is much larger than the one expected from the empirical relation : where the value 0.5 GeV has been inspired from ordinary mesons and from the results of [95]. From our result, one can already conclude that the Z c (4020) to Z c (4430) are too low to be the radial excitations of the Z c (3900) unless they couple weakly to the interpolating current such that their effect is tiny in the LSR analysis.

G. Can there be a weakly coupled radial excitation
If one literally extrapolates the phenomenological observation from ordinary hadrons, one would expect a radial excitation with a mass: which is relatively low compared to the previous prediction in Eq. 33 but seems to fit the Z c (4430). To understand why it can have been eventually missed in the previous analysis, we shall determine the Z c (4430) coupling to the current and use as input the experimental mass value 11  In so doing, we reconsider the L c 0 moment. We use a two-resonance parametrization of the spectral function and introduce as inputs the previous values of the D * D molecule mass and coupling. We use the experimental mass of the Z c (4430). We include the high-mass (D * D) 1 radial excitation obtained previously into the QCD continuum contribution.
We show the τ -behaviour of the coupling for different t c -values and for fixed µ = 4.65 GeV in Fig. 24. At the (τ, t c ) stability regions (0. 30,27) to (0.34, 46) (GeV −2 , GeV 2 ), we deduce : where the different sources of the errors are given in Table V. The coupling is indeed relatively small compared to that of f D * D = 140 keV. It can even be consistent with zero due to the large errors mainly induced by the coupling of the ground state and of the dimension-6 condensates. This result may support an eventual radial excitation interpretation of the Z(4340) which couples very weakly to the current and having a mass much lower than the strongly coupled [f (D * D)1 = 197 MeV] radial excitation with a mass 5709 MeV (see Table V). This weak coupling disagrees with the one obtained in Ref [96] and may originate from the fact that the latter analysis has been done at a smaller value of t c and at the scale µ = 1.5 GeV which is too low compared to the optimal choice obtained in the present work. However, as we have already mentioned in Subsection VI B, we do not find any convincing theoretical basis for justifying this low choice of µ. We expect that similar results can be obtained in some other channels. 11 No stability region is obtained for an the attempt to determine this mass from the LSR R c 0 .

XII. VERSUS OUR PREVIOUS RESULTS
We compare our results with our previous ones from Refs. [23,24] in Table VII. Notice that in [23], the double ratio of moments has been also used to extract directly the SU (3) breaking contributions to the masses and couplings. One can notice a good agreement between the different results. The exception is the central value of M D * 0 D1 which slightly moves in the 3 papers though the results are in agreement within the errors. This is due to the difficult localization of the inflexion point which we have identified in the present work with the minimum of the coupling like in some other channels.

XIII. CONFRONTATION WITH THE DATA
As proposed in [44], we consider that the physical states are superposition of quasi-degenerated hypothetical molecules and tetraquark states having the same quantum numbers J P C and having almost the same coupling strength to the currents. We have denoted these observed states as Tetramoles T (tetraquarks ⊕ molecules).

A. Zc(3900) as a tetramole state
One may define the tetramole T c as a superposition of the D * D molecule and A cd tetraquark state with the parameters : M Tc = 3900(42) MeV, f Tc = 155(11) keV, (40) which are the mean of the two previous couplings and masses. We identify this tetramole state with the Z c (3900) found by BELLE [80] and BESIII [81].  Tables VII and III might be identified with the Z c (4025) found by BES [83] or with the Z c (4050) found by BELLE [84]. However, one should mention that these states are not well established and have not yet been retained in the PDG Summary Table [65].
Pursuing our confrontation with the data, we note that the Z c (4430) [86,87] is too low to be the strongly coupled (f (D * D)1 = 196(41) MeV) first radial excitation of the D * D expected to be in the range (5.4 ∼ 5.8) GeV (see Table III). However, it can be fitted by the weakly coupled (f (D * D)0 = 46(56) MeV) low mass one obtained in Subsection XI G.

States
Couplings Masses Ref. [24] Ref. [23] New Ref. [24] Ref. [  We consider the tetramole T s as the combination of the D * s D, D * D s molecules and A cs tetraquark states have almost the same mass and same coupling to the currents (see Table III). Taking the mean values of these parameters, we obtain: f Tcs = 136(9) keV, (42) which we can identify with the recent Z cs (3983) found recently by BESIII [45].
E. The Zcs(4100) bump as a D * s0 D1 ⊕ D * 0 Ds1 molecule Inspecting our predictions for the D * s0 D 1 and D * 0 D s1 masses in Tables III and VII, we find that they are almost degenerated and have almost the same couplings. Taking their combination having a mean mass : we are tempted to identify this state with the Z cs (4100) bump observed by BESIII [45]. Its relatively small coupling to the current compared to the one of the ground state Z cs (3983) may indicate that it can be large using a Golberger-Treiman-type relation argument where the hadronic coupling behaves as the inverse of coupling (see e.g: [97]).
F. The future Zcss states From our predictions in Table III, one can also define a tetramole T css which is a superposition of the D * s D s molecule and A css tetraquark states with : M Tcss = 4068(48) MeV, f Tcss = 114(10) keV. (44) One also expects to have a D * s0 D s1 molecule at a higher mass value: where its coupling to the current is relatively small indicating that it can be relatively large using a Golberger-Treiman-type relation argument where the hadronic coupling behaves as 1/f D * s0 Ds1 (see e.g: [97]). These predicted states are expected to be seen in the near future experiments and can be considered as a test of the predictions given in this paper.

XIV. ON SOME OTHER AND LSR RESULTS
After the publication of the recent BESIII results on the observation of the Z cs (3982) candidate [45], many papers using different models appear in the literature [98] for attempts to explain the true nature of this state.
Besides the pioneer QSSR estimate of the D * D s molecule and tetraquark states [99] 12 , some recent papers using LSR come to our attention (see e.g [100])where we notice some common caveats : -All analysis is done at lowest order (LO) of perturbation theory where the choice of the value of the M S running mass in favour of the pole mass is unjustified because the definitions of the two masses are undistinguishable at this order. Moreover, the calculation of the spectral functions using on-shell renormalization would (a priori) favour the choice of the on-shell mass. The difference on the effect of this choice is explicitly shown in Fig. 13. Hopefully, the effects of NLO corrections in the M S-scheme are tiny (see Fig. 13) in the LSR analysis confirming (a posteriori) the intuitive choice of the M S running mass at LO.
-The value of t c ≡ s 0 used to determine the mass and coupling of about (18 ∼ 23) GeV 2 is relatively low as it corresponds to the beginning of the τ -stability region (see previous figures) where one also notice that the predictions increase until the t c -stability value. As a consequence, the absolute value and the error in the extraction of the mass and coupling have been underestimated.
-In general, the way how the errors from different sources has been estimated are not explained in details, raises some doubts on the real size of the quoted errors having in mind that the extraction of the different errors quoted in Tables III and V require some painful works.
-The value of the four-quark condensates from the vacuum saturation assumption is often used. However, though this estimate is correct in the large N c limit, it has been found phenomenologically from different lightquarks and τ -decay channels [36,[48][49][50][51][52] that this estimate is largely violated at the realistic case N c = 3. Moreover, it is also known that the dimension-six quark condensates which mix under renormalization [46] (Part VII page 285) does not support the previous assumption.
-The previous papers extend the OPE to highdimension up to d = 10 vacuum condensates but only for some classes of high-dimension condensate contributions and by assuming the validity of the factorization assumption for estimating their sizes. However, it has been shown in [46] [27] (Part VII page 285) that the structure of these high-dimension condensates are quite complex due to their mixing under renormalization such that their inclusion in the QSSR analysis should deserve more care.
In addition to these caveats, we also notice that : -Differentiating the neutral from the charged Z cs is purely academic from the approach as the two states are almost degenerated within the errors.
-In Ref. [88,89], the author introduces a relation between the PT subtraction scale µ and the so-called energy bound which is quite obscure to us as discussed in Subsection VI B.

XV. SUMMARY AND CONCLUSIONS
We have : -Systematically studied the spectra and couplings of the (cq)(q ′ c) molecules and (cq)(q ′ c) tetraquark states where q, q ′ ≡ d, s are light quarks.
-Improved our previous predictions obtained in Refs. [23,24] by a much better localization of the t c , τ and µ-stability points and by using updated values of some QCD input parameters.
-Emphasized that the localization of the inflexion point for extracting the values of the masses can be fixed more precisely in most channels by identifying it with the value of τ corresponding to the minimum of the curve where optimal value of the coupling is extracted.
-Introduced the tetramole states as a superposition of quasi-degenerated molecules and tetraquark states having the same quantum numbers J P C with almost the same couplings strengths to the interpolating currents. It is remarkable to notice that the mass-spliitings of the tetramoles T c , T cs and T css due to SU (3) breakings are successively about (73 ∼ 91) MeV.
-Completed the analysis with new predictions of some first radial excitation masses and couplings. One can notice that the mass gap of about (1.7 ∼ 2.35) GeV, between the lowest mass ground state and the 1st radial excitation strongly coupled to the current, is quite large compared to (0.5 ∼ 0.7) GeV for ordinaryqq mesons. Similar results have been obtained for the DKlike states [44]. However, a weakly coupled radial excitation (see Section XI G) having a lower mass of about 4.4 GeV is not excluded from the approach. These features may signal some new dynamics of these exotic states which can be found in these high-mass regions.
-Compared successfully our predictions with the observed Z c and Z cs -spectra.
-Given new predictions for the future Z css states.