Doubly-hidden scalar heavy molecules and tetraquarks states from QCD at NLO

Alerted by the recent LHCb discovery of exotic hadrons in the range (6.2 -- 6.9) GeV, we present new results for the doubly-hidden scalar heavy $(\bar QQ) (Q\bar Q)$ charm and beauty molecules using the inverse Laplace transform sum rule (LSR) within stability criteria and including the Next-to-Leading Order (NLO) factorized perturbative and $\langle G^3\rangle$ gluon condensate corrections. We also critically revisit and improve existing Lowest Order (LO) QCD spectral sum rules (QSSR) estimates of the $({ \bar Q \bar Q})(QQ)$ tetraquarks analogous states. In the example of the anti-scalar-scalar molecule, we separate explicitly the contributions of the factorized and non-factorized contributions to LO of perturbative QCD and to the $\langle\alpha_sG^2\rangle$ gluon condensate contributions in order to disprove some criticisms on the (mis)uses of the sum rules for four-quark currents. We also re-emphasize the importance to include PT radiative corrections 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. Our LSR results for tetraquark masses summarized in Table II are compared with the ones from ratio of moments (MOM) at NLO and results from LSR and ratios of MOM at LO (Table IV). The LHCb broad structure around (6.2 --6.7) GeV can be described by the $\overline{\eta}_{c}{\eta}_{c}$, $\overline{J/\psi}{J/\psi}$ and $\overline{\chi}_{c1}{\chi}_{c1}$ molecules or/and their analogue tetraquark scalar-scalar, axial-axial and vector-vector lowest mass ground states. The peak at (6.8--6.9) GeV can be likely due to a $\overline{\chi}_{c0}{\chi}_{c0}$ molecule or/and a pseudoscalar-pseudoscalar tetraquark state. Similar analysis is done for the scalar beauty states whose masses are found to be above the $\overline\eta_b\eta_b$ and $\overline\Upsilon(1S)\Upsilon(1S)$ thresholds.


I. INTRODUCTION
QCD spectral sum rules (QSSR)à la SVZ [1][2][3] have been applied since 41 years 1 to study successfully the hadron properties (masses, couplings and widths) and to extract some fundamental QCD parameters (α s , quark masses, quark and gluon condensates,...). In previous series of papers [18][19][20][21][22][23], we have used the inverse Laplace transform (LSR) [25][26][27][28] of QSSR to predict the couplings and masses of different heavy-light molecules and tetraquarks states by including next-to-next nonleading order (N2LO) factorized perturbative (PT) corrections where we have emphasized the importance of these corrections for giving a meaning of the input heavy quark mass which plays an important role in the analy-sis though these corrections are small in the M S-scheme. However, this feature (a posteriori) can justify the uses of the M S running masses at LO in some channels [24] if the α n s -corrections are small, especially in the ratios of moments used to extract the hadron masses where these corrections tend to compensate [4,5].
In this paper, we pursue the analysis for the fully / doubly-hidden heavy quarks (QQ)(QQ) molecules and (QQ)(QQ) tetraquarks states, where the effect of the quark mass value and its definition are (a priori) important as we have four heavy quarks which bound these states.
We separate explicitly the factorized and nonfactorized contributions to the four-quark correlators at LO of PT QCD and for the lowest dimension gluon condensate ⟨α s G 2 ⟩ contributions. We add the contribution of the NLO perturbative corrections from the factorized part of the diagrams which as we shall see is a good approximation. We also include the triple gluon condensate ⟨G 3 ⟩ contributions in the Operator Product Expansion (OPE).
We use these QCD results using the LSR sum rules within different stability criteria used successfully in some arXiv:2008.01569v4 [hep-ph] 6 Aug 2023 other channels to extract the masses and couplings of the previous molecules and tetraquarks states assumed to be resonances.
Our results as improved estimates of the LO ones given in the QSSR literature.
We expect that these will be an useful guide for further experimental searches of these exotic states and for identifying the different new states found by LHCb [29,30].

II. THE INVERSE LAPLACE SUM RULES
A. The QCD molecule local interpolating currents We shall be concerned with the following QCD local interpolating currents of dimension-six: where f H M is the meson decay constant; J H M (x) is the lowest dimension bilinear quark currents and H ≡ S, P, V, A.
For the scalar (0 ++ ) molecule states, these currents are : Interpolating currents constructed from bilinear (pseudo)scalar currents are not renormalization group invariants such that the corresponding decay constants possess anomalous dimension: where :f (S,P ) M is the renormalization group invariant coupling and −β 1 = (1/2)(11 − 2n f /3) is the first coefficient of the QCD β-function for n f flavours. a s ≡ (α s /π) is the QCD coupling. k f = 2.028(2.352) for n f = 4(5) flavours.

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 : where m Q is the heavy quark mass, τ is the LSR variable, n = 0, 1 is the degree of moments, t c is the threshold of the "QCD continuum" which parametrizes, from the discontinuity of the Feynman diagrams, the spectral function Im Π H M (t, m 2 Q , µ 2 ) where Π H M (t, m 2 Q , µ 2 ) is the scalar correlator defined as :
We have cross-checked that, using our calculation method, we recover the results for charmonium where the heavy quark condensate contribution is already included into the gluon condensate one through the relation [1,2,[31][32][33][34]: The LO ⊕ ⟨G 2 ⟩ expressions of the other molecules spectral functions are given in Appendix B. The one of the ⟨G 3 ⟩ condensates which are lengthy are not quoted.
We note that the inclusion of higher dimension condensate contributions (d ≥ 8) , as abusively done in the current literature, does not help, except in some few cases, because the OPE is often convergent at the optimization scale while the size of higher dimension condensates are not under control due to the violation of factorization for the four-quark [35][36][37][38] and to the inaccuracy of the dilute gas instanton estimate of higher dimensions gluon [39][40][41] condensates.

B. NLO PT corrections to the Spectral functions
We extract the next-to-leading (NLO) perturbative (PT) corrections by considering that the molecule /tetraquark two-point spectral function is the convolution of the two ones built from two quark bilinear currents (factorization) as illustrated in Fig. 4. This is a good approximation because we have seen for the LO that the non-factorized part of the QCD diagrams gives a small contribution and behaves like 1/N c where N c is the number of colours.
where : with the phase space factor: and m Q is the on-shell / pole perturbative heavy quark mass.
The NLO perturbative expressions of the bilinear equal masses pseudoscalar spectral functions are known in the literature [4,5,10,44].
We estimate the N2LO contributions assuming a geometric growth of the numerical coefficients [45]. We consider this contribution as an estimate of the error due to the truncation of the PT series.

C. From the On-shell to the M S-scheme
We transform the pole masses m Q to the running masses m Q (µ) using the known relation in the M Sscheme to order α 2 s [46][47][48][49][50][51][52][53][54]: for n l = 3 : u, d, s light flavours. In the following, we shall use n f =4 or 5 total number of flavours for the numerical value of α s respectively for the charm and bottom quarks.

IV. QCD INPUT PARAMETERS
The QCD parameters which shall appear in the following analysis will be the QCD coupling α s the charm and bottom quark masses m c,b , the gluon condensates ⟨α s G 2 ⟩. Their values are given in Table I.

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

B. c and b quark masses
For the c and b quarks, we shall use the recent determinations [56,57] of the running masses and the corresponding value of α s evaluated at the scale µ obtained using the same sum rule approach from charmonium and bottomium systems. We use the recent estimate obtained from a correlation with the values of the heavy quark masses and α s which can be compared with the QSSR average from different channels [55].

V. THE SPECTRAL FUNCTION
In the present case, where no complete data on the spectral function are available, we use the duality ansatz: (17) for parametrizing the molecule spectral function. M H and f H are the lowest ground state mass and coupling analogue to f π . The "Continuum" or "QCD continuum" is the imaginary part of the QCD correlator from the threshold t c . Within a such parametrization, one obtains: indicating that the ratio of moments appears to be a useful tool for extracting the mass of the hadron ground state [4][5][6][7]16].
This simple model has been tested in different channels where complete data are available (charmonium, bottomium and e + e − → I = 1 hadrons) [4,5,12]. It was shown that, within the model, the sum rule reproduces well the one using the complete data, while the masses of the lowest ground state mesons (J/ψ, Υ and ρ) have been predicted with a good accuracy. In the extreme case of the Goldstone pion, the sum rule using the spectral function parametrized by this simple model [4,5] and the more complete one by ChPT [59] lead to similar values of the sum of light quark masses (m u +m d ) indicating the efficiency of this simple parametrization.
An eventual violation of the quark-hadron duality (DV) [60][61][62] has been frequently tested in the accurate determination of α s (τ ) from hadronic τ -decay data [35,61,63], where its quantitative effect in the spectral function was found to be less than 1%. Typically, the DV behaves as: where κ, α, β are model-dependent fitted parameters but not based from first principles. Within this model, where the contribution is doubly exponential suppressed in the Laplace sum rule analysis, we expect that in the stability regions where the QCD continuum contribution to the sum rule is minimal and where the optimal results in this paper will be extracted, such duality violations can be safely neglected. Therefore, we (a priori) expect that one can extract with a good accuracy the masses and decay constants of the mesons within the approach. An eventual improvement of the results can be done after a more complete measurement of the corresponding spectral function which is not an easy experimental task.
In the following, in order to minimize the effects of unkown higher radial excitations smeared by the QCD continuum and some eventual quark-duality violations, we shall work with the lowest ratio of moments R c 0 for extracting the meson masses and with the lowest moment L c 0 for estimating the decay constant f H . Moment with negative n will not be considered due to their sensitivity on the non-perturbative contributions at zero momentum.

VI. OPTIMIZATION CRITERIA
For extracting the optimal results from the analysis, we have used in previous works the optimization criteria (minimum sensitivity) of the observables versus the variation of the external variables namely the τ sum rule parameter, the QCD continuum threshold t c and the subtraction point µ.
Stability on the number n of heavy quark moments have also been used [39][40][41]57].
One should notice in the previous works that these criteria have lead to more solid theoretical basis and noticeable improvement of the sum rule results. The quoted errors in the results are conservative as the range covered by t c from the beginning of τ -stability to the one of t c -stability is quite large. However, such large errors induce less accurate predictions compared with some other approaches (potential models, lattice calculations) especially for the masses of the hadrons. This is due to the fact that, in most cases, there are no available data for the radial excitations which can be used to restrict the range of t c -values. However, the value of t c used in the "QCD continuum" model does not necessarily coincide with the 1st radial excitation mass as the "QCD continuum" is expected to smear all higher states contributions to the spectral function. This feature has been explicitly verified by [38] in the ρ-meson channel.
VII. THE SCALAR χ c0 χc0 MOLECULE Using the previous QCD expression given in Eq. 6 and adding the PT NLO contribution, we study the dependence of the coupling and mass on the LSR parameter τ , the continuum threshold t c and the subtraction scale µ. We shall also study the relative contribution of the continuum versus the ground state one.

B. µ-stability
Fixing t c = 70 GeV 2 and τ = (0.35 − 0.38) GeV −2 , we show in Fig. 6 the µ behaviour of the mass and coupling where we note an inflexion point at : in agreement with the one quoted in [56,75,76] using different ways and/or from different channels.

C. QCD continuum versus lowest resonance
To have more insights on the QCD continuum contribution, we study the ratio of the continuum over the lowest ground state contribution as predicted by QCD : We found that for t c ≥ 55 GeV 2 , the continuum contribution is less than 60% of the ground state one and decreases quickly for increasing t c indicating a complete dominance of the ground state contribution in the sum rule.  Table I.

D. PT series and higher order terms
We compare in Fig. 7 the LO and NLO perturbative contributions. As the input definition of the quark mass is ambiguous at LO, we use the running mass evaluated at µ = 4.5 GeV and the corresponding on-shell / pole mass M (µ = M ) = 1.53 GeV. We see that, for the coupling, the two mass definitions lead to about the same predictions but there is a difference about 400 MeV for the mass prediction. This systematic error is never considered in the literature where a running mass is often used ad hoc with not any justification. This ambiguity is avoided when the PT corrections are added.
Comparing the predictions for the running mass at given τ ≈ 0.17 GeV −2 , t c ≃ 70 GeV 2 and µ = 4.5 GeV, one can parametrize numerically the result as : where the PT corrections tend to compensate in the ratio of moments used to determine the mass of the meson. We have estimated the N2LO contributions from a geometric growth of the PT coefficients [45] which we consider as an estimate of the uncalculated higher order terms of the PT series.
One can notice, like in the case of the two-point functions of the scalar quark bilinear currents, that the coefficients of radiative corrections are large for the decay constant [4,5,27]. However, the PT series converge numerically at µ = 4.5 GeV but induce a relatively large systematic error when the higher order terms of the PT series are estimated using a geometric growth of the numerical coefficients.
VIII. THE η c ηc, J/ψJ/ψ, χ 1c χ1c MOLECULES The τ , t c and µ behaviours of the coupling and mass of these molecules are very similar to the one of χ 0c χ 0c and will not be repeated here. The values τ -and t c at the stability regions are shown in Table XII where one can notice that, for the η c η c , the stabilities are reached at earlier values of t c which is dual to the lower value of the η c η c molecule mass.
In all cases, the inclusion of the ⟨G 3 ⟩ condensate shift the τ -stabilty to smaller values. In the case of the η c η c , it becomes 0.36 GeV −2 for the coupling (minimum) and 0.34 GeV −2 for the mass (inflexion point).
The main difference with the χ 0c χ 0c as shown in Figs.7 is the almost equal position of the τ minima for the LO and LO ⊕ NLO contributions as shown in Fig. 8, which can be attributed to the different reorganisation of the terms in each channel.
Our results also emphasize the importance to add radiative PT corrections for a proper heavy quark input (pole or M S running) mass definition. In the M S scheme, the α s correction is small as can be seen explic-  itly in this numerical parametrization : The µ-stability is reached at µ = 4.5 GeV. The results of the analysis are shown in Table II.
The extension of the analysis to the b quark channel is straigthforward. We show in this example the details of the analysis.
The µ-stabilty is shown in Fig. 11.  Table I.

B. µ-stability
Fixing t c = 460 GeV 2 and τ = 0.17 GeV −2 , we show in Fig. 11 the µ behaviour of the mass and coupling, where we find a clear inflexion point for the coupling but a slight for the mass at : in agreement with the one quoted in [56,75,76] using different ways and/or from different channels.

C. LO versus NLO contributions
We compare in Fig 12 the LO and LO ⊕ NLO contributions. We note (as expected) that the radiative corrections is smaller for b than for c as the coupling and mass are evaluated at higher µ-values. Using this result, we can numerically parametrize the previous observables  Table I. as: where the PT corrections tend to compensate in the ratio of moments while, compared to the c-quark channel, the PT corrections are relatively small. As in the previous cases, we have estimated the N2LO contributions from a geometric growth of the PT coefficients [45] which we consider as an estimate of the uncalculated higher order terms of the PT series.

X. ⟨G 3 ⟩ AND TRUNCATION OF THE OPE
We have included the ⟨G 3 ⟩ condensate contribution into the sum rule. We have cross-checked that with our method of calculation we reproduce the results of [33] for charmonium sum rules.
We have noticed that in the χ c0 χ c0 channel, the contribution of the ⟨G 3 ⟩ condensate is relatively small and does not modify the shape of the mass and coupling curves versus the variation of τ and for different values of t c . It only decreases the decay constant by 0.4 keV and increases the mass by 14 MeV.
However, this is not the case of some other channels which will be analyzed later on where the ⟨G 3 ⟩ contribution can be large and modify the minimum of the mass found for ⟨α s G 2 ⟩ into an inflexion point (see Fig. 13) and vice-versa for the coupling. This feature renders the mass result quite sensitive to the localisation of this inflexion point. An analogous effect of ⟨G 3 ⟩ has been also observed e.g in the analysis of charmonium sum rules [39][40][41] and the inclusion of the ⟨G 4 ⟩ condensates which act with an opposite sign restores the stability of these sum rules.
To circumvent this problem and due to the difficulty for evaluating the ⟨G 4 ⟩ contribution, we consider the optimal result at the value of τ where the coupling presents a minimum. Then we consider as a final result (here and in the following), the mean obtained with and without the ⟨G 3 ⟩ contribution. The error induced in this way will be included as the systematics due to the truncation of the OPE as quoted in TableII.
The analysis of these scalar molecules is very similar to the analysis presented above. The QCD expressions of their corresponding two-point functions are given in Appendix A. One should mention that in these channels the PT radiative corrections and the contribution of the ⟨G 3 ⟩ condensate are small indicating a good convergence of the PT series and of the OPE at the optimization scale. The results are quoted in Table II where the LSR parameters used to get them are shown in Table XII.

XII. THE SCALAR TETRAQUARK STATES
We repeat the previous analysis for the case of tetraquark states with same choice of diquark currents as in [77] : in order to make a direct comparison with their LO results. We do not consider the current associated to σ µν which corresponds to a two-point correlator of higher dimension. We shall also consider the four-quark operator : in order to make a direct comparison with [78]. One should notice that due to the epsilon-tensor, most of the currents used by [77] are not present in [78].
The QCD expressions of their corresponding two-point functions are given in Appendix B.
The behaviours of different curves are very similar with the ones of the corresponding molecule case.
We quote the results in Table II and the optimal LSR parameters used to get them in Table XII. These results are compared with the ones in [77,78] in Table IV.   Table I. ∆µ are given in Eqs. 20 and 24. We take |∆τ | = 0.02 GeV −2 . In the case of asymetric errors, we take the mean value. The inclusion of the ⟨G 3 ⟩ contribution and the way to estimate the systematics induced by the truncation of the OPE are explained in Section X.

A. The quest of factorization and Landau singularities
We have shown explictily in Eq. 6 that the contributions from the non-factorized diagrams appear already to LO of perturbative series and for the lowest dimension ⟨α s G 2 ⟩ gluon condensate contributions. This result does not support the claims of [79,80] that non-factorized contributions start to order α 2 s . However, this effect shown in Fig. 14 is numerically small(about 3%≈ 1/(10N c )) of the sum of factorized ⊕ non-factorized contributions as expected from large N c limit and Fierz transformations. This feature has been already observed explicitly in our previous work [19,23]. This small effect of the nonfactorized contribution justifies the accuracy of our ap-proximation by only using the factorized diagrams in the NLO perturbative contributions. We do not also see the relevance / appearance of the Landau singularities mentioned by [79,80] in the analysis using the OPE in the Euclidian region. However, the two-point function analyzed in [79,80] has nothing to do with the one analyzed in our paper as it corresponds to a four-point function compacted into a twopoint function but with four legs i.e with two incoming and two-outgoing momenta. This four-point function is more relevant for the analysis of hadron-hadron scatterings (see the example of ππ and γγ in [81,82]), while in this case, a two-point function enters differently via a gluonium intermediate state [83].
From the analysis of Eq. 21, we have shown that the postulated lowest mass ground state dominates the spectral function. This feature indicates that the nonresonant states do not play a crucial role in the analysis. This conclusion may go in line with the answer of [84,85] on some of the comments of [79,80].

B. Systematic errors
As mentioned in Section V, one expects that at the optimization region, an eventual duality violation is expected to be negligible and the QCD continuum contribution which parametrizes non-resonant states is dominated by the lowest resonance as can be checked from Eq. 21. Therefore, the high-energy tail of the spectral function cannot bring a sizeable systematic error.
The error due to the truncation of the PT series cannot be quantified with a good accuracy as the LO contributions are quite sensitive to the quark mass definition (pole or running) in some other channels. Using an approach similar to the one leading to Eq. 22 where a geometric growth of the a s -coefficients has been assumed, we deduce the error estimate in Table II. We have estimated the unknown higher dimension condensates contributions in the OPE quoted in Table II as discussed in Section VII.

N ew compared with available QSSR results
Compared to previous QSSR LO results given in the literature (see Table IV): We have included (for the first time) the NLO corrections which is mandatory for giving a sense on the definition and numerical values of the input heavy quark mass which plays a crucial role in the analysis.
We have added the contributions of the dimensionsix ⟨G 3 ⟩ condensates, which are quite large for the η q η q and J/ψJ/ψ, ΥΥ molecules and for theV q V q andP q P q tetraquark states.
Our results are shown in Table II where systematic analysis of some possible configurations of the 0 ++ molecule and four-quark states have been done.

C. LSR versus the ratio of MOM results
Taking, the example of the χ c0 χ c0 molecule andS q S qtetraquark, we use the ratio of moments as in [77]: where M T ,M is the molecule or tetraquark mass. We take e.g Q 2 0 = 4m 2 Q . Then, we find that the LO and LO ⊕ NLO results are about the same as from the LSR obtained in the previous sections. To NLO and including ⟨α s G 2 ⟩, one obtains in units of GeV : 19.29 , (29) compared to the ones from LSR in Table II, indicating that the two methods give (within the errors) the same results.
D. On the ratio of MOM results of Ref. [77] Using the QCD expression of theS q S q tetraquark twopoint function given in Appendix A, we have also compared our LO ⊕ ⟨α s G 2 ⟩ MOM results : with our LO LSR results given in Table IV where we find (within the errors) a good agreement. However, by comparing these LO MOM results with the ones from [77] quoted in Table IV, one can see that the results of [77] are about 0.34 GeV (resp. 1.08 GeV) for the charm (resp. beauty) case lower than the ones   [77] are from Moments at LO and of Ref. [78] from LSR at LO. As already mentioned earlier, we notice that the choice of the numerical values of the M S running quark masses used at LO is not justified due to the ambiguous quark mass definition to that order. One may also have equally used a pole / on-shell mass which naturally appears in the expression of the spectral function evaluated using on-shell quark mass.
in Eq. 30. More generally, compared to our LO ones, the LO results of [77] have the tendancy to underestimate the mass results.
With the inclusion of the NLO QCD corrections, our predictions agree (within the errors) with the LO results of [77] for the charm and for theP q P q beauty channels. For theS q S q ,Ā q A q andV q V q beauty ones, the disagreements persist and range from 0.77 to 1.41 GeV. We cannot trace back the origin of a such discrepancy as the comparison of the QCD expressions of the full correlator given in [77] with the one using the spectral function is not easy due to the choice of variables used by the authors.
Therefore, unlike Ref. [77], we expect, like in the charm case, that the tetraquark beauty states are above the η b η b and Υ(1S)Υ(1S) thresholds. The future experimental findings of these beauty states may select among these theoretical predictions.
E. Comparison with the LSR results of Ref. [78] We have also compared our results for theĀ q A q scalar tetraquark with the LO ones of [78] using the current in Eq. 27.
The PT QCD expressions agree each other at LO. There is a slight difference for the ⟨α s G 2 ⟩ contribution for higher values of t to all values of the heavy quark mass but this difference affects only slightly the predictions.
At LO and including the ⟨α s G 2 ⟩ contribution, our values of theĀ q A q couplings of about 287 (resp. 78) keV for the charm (resp. bottom) are comparable with the ones of [78] (289 (resp. 54) keV) if the (unjustified) choice of M S-mass is used.
For the charm, theĀ q A q mass of [78] is (460 550) MeV lower than the one of [77] and our LO result, while for the bottom it is 670 MeV lower than our LO result but 380 MeV higher than that of [77] (see Table IV). However, the origin of this discrepancy does not come from the QCD input parameters as we use about the same values. This example puts a question mark on the unusual treatment of the sum rules by the author in [78].
His choice of the subtraction scale µ ≃ (1.2 ∼ 2.2) GeV for the charm (resp. (2.3 ∼ 3.3) GeV for the bottom) based, for instance, on the identification of the sum of the PT running mass (m c + m b )(µ) with the value of the B c -mass [86] is difficult to justify in the absence of NP-contributions (binding energy). However, such low values of µ are quite dangerous as, at this low scale, the PT radiative corrections are expected to be large and can strongly affect the final result. This is indeed the case for the coupling where, at the µ-stability (4.5 GeV for the charm and 7.25 GeV for the bottom) the NLO corrections increase it by 59% for the charm and 83% for the bottom. This effect is obviously larger for smaller values of µ.
Moreover, using only the µ-dependence of the running values of α s and m Q into the PT LO expression of the sum rule is also inconsistent while the identification of the QCD continuum threshold with the mass of the first radial excitation can be inaccurate as the QCD continuum is expected to smear all higher state contributions.
It is also remarkable to notice from Tables II and IV the (almost) independence of our results on the form of the current for theĀ q A q tetraquark.
For a consistency check of our results, we compare our result for theĀ q A q tetraquark mass MĀ q Aq ≃ 6.47 (resp. 19.72) GeV from the current of [78] within a3 c ⊗3 c color representation with the one from the combination of molecule currents 2(S q S q +P q P q ) +V q V q −Ā q A q given there. Using a quadratic mass relation, we deduce at NLO ⊕G3: MĀ q Aq ≃ 6.38 (resp. 19.49) GeV in agreement (within the errors) with our predicted tetraquark masses.

F. Some phenomenological implications
One can notice from Tables II and IV that : Our different QSSR predictions cannot disentangle (within the errors) the mass of a molecule from a tetraquark state as already found in some of our previous works [18][19][20][21][22].
Our results do not favour the ones from some potential models where the exotic states are below the η c η c meson thresholds. Instead, our results may explain the existence of a 0 ++ broad structure around (6.2 6.7) GeV which can be due to η c η c , χ c1 χ c1 and J/ψJ/ψ molecules or /and to scalar-scalar, vector-vector and axial-axial scalar tetraquark states.
If the new LHCb peak candidate [29,30] around (6.8 6.9) GeV is a 0 ++ state, the value of its mass suggests that it is likely a χ c0 χ c0 molecule or a pseudoscalarpseudoscalar tetraquark states. Its signature from a J/ψJ/ψ invariant mass may come from the di-χ c0 decaying to di-γJ/ψ.
In the case of a χ c1 χ c1 molecule, the predicted mass is below the χ c1 χ c1 threshold while our NLO predictions for the beauty states indicate that all of them are above the η b η b and Υ(1S)Υ(1S) thresholds.
We plan to calculate the spectra of some other 0 − , 1 ± and 2 ++ channels and eventually their widths in a future work.
Appendix A: ⟨G 3 ⟩ contributions to the χ 0q χ0q spectral function The ⟨G 3 ⟩ contributions to theχ 0q χ 0q spectral function are given by the Feynman diagrams drawn in Figs. 15 and 16. As the expression is quite lengthy, we shall only present the one for χ 0q χ 0q but not for some other molecules.