NLO results with operator mixing for fully heavy tetraquarks in QCD sum rules

We study the mass spectra of $\bar{Q}Q\bar{Q}Q\ (Q=c,b)$ systems in QCD sum rules with the complete next-to-leading order (NLO) contribution to the perturbative QCD part of the correlation functions. Instead of meson-meson or diquark-antidiquark currents, we use diagonalized currents under operator renormalization. We find that differing from conventional mesons $\bar qq$ and baryons $qqq$, a unique feature of the multiquark systems like $\bar{Q}Q\bar{Q}Q$ is the operator mixing or color configuration mixing induced by NLO corrections, which is crucial to understand the color structure of the states. Our numerical results show that the NLO corrections are very important for the $\bar{Q}Q\bar{Q}Q$ system, because they not only give significant contributions but also reduce the scheme and scale dependence and make Borel platform more distinct, especially for the $\bar{b}b\bar{b}b$ in the $\overline{\rm{MS}}$ scheme. We use currents that have good perturbation convergence in our phenomenological analysis. With the $\overline{\rm{MS}}$ scheme, we get three $J^{PC}=0^{++}$ states, with masses $6.35^{+0.20}_{-0.17}$ GeV, $6.56^{+0.18}_{-0.20}$ GeV and $6.95^{+0.21}_{-0.31}$ GeV, respectively. The first two seem to agree with the broad structure around $6.2\sim6.8$ GeV measured by the LHCb collaboration in the $J/\psi J/\psi$ spectrum, and the third seems to agree with the narrow resonance $X(6900)$. For the $2^{++}$ states we find one with mass $7.03^{+0.22}_{-0.26}$ GeV, which is also close to that of $X(6900)$, and another one around $7.25^{+0.21}_{-0.35}$ GeV, which has good scale dependence but slightly large scheme dependence.


QCD Sum Rule
In this section, we briefly review the framework of the QCD sum rules used to calculate the mass of the tetraquark ground state. See Ref. [55] for more details. We start with a two-point correlation function where D denotes the spacetime dimension, Ω denotes the QCD vacuum and J is the (pseudo-)scalar tetraquark current to be defined later.
On the one hand, the correlation function Π(q 2 ) can be related to the phenomenological spectrum by the Källén-Lehmann representation [55], where ρ(s) denotes the physical spectrum density. Taking the narrow resonance approximation for the physical ground state, one can parametrize the spectrum density as a pole plus a continuum part where M H and λ H denote the mass of the ground state and pole residue, respectively. ρ cont (s) denotes the continuum spectrum density, which could also contain information of higher resonances. s h is the threshold of the continuum spectrum.
On the other hand, in the region where −q 2 = Q 2 Λ 2 QCD , one can calculate correlation function Π(q 2 ) using the operator product expansion (OPE), which reads where C 1 and C i are perturbatively calculable Wilson coefficients, and O i is a shorthand of the vacuum condensate Ω|O i |Ω , which is a nonperturbative but universal quantity. The relative importance of the vacuum condensate is power suppressed by the dimension of the operator O i . In our calculations, we will only keep the relevant vacuum condensates up to dimension four, which gives the approximated expression of the OPE as Π(q 2 ) = C 1 (q 2 ) + C GG (q 2 ) g 2 sĜĜ , (2.5) where g 2 sĜĜ denotes the gluon-gluon (GG) condensate Ω|g 2 sĜĜ |Ω . According to Eq. (2.2), one can relate the physical spectrum density to the imaginary part of Π(q 2 ) in Eq. (2.5) using the dispersion relation, which gives Π(q 2 ) = ds ρ(s) s − q 2 − i = 1 π ∞ s th ds ImC 1 (s) + ImC GG (s) g 2 sĜĜ s − q 2 − i , (2.6) where s th = 16m 2 Q is the QCD threshold for theQQQQ system, and the integral in the second line has been assumed to be convergent. Then by employing the quark-hadron duality and Borel transformation [55], we obtain a sum rule for Π(q 2 ), where s 0 is the threshold parameter and M B is the Borel parameter. They are introduced into the formula due to the qurak-hadron duality and Borel transformation, respectively. By differentiating both sides of Eq. (2.7) with respect to − 1 where ρ 1 = 1 π ImC 1 and ρ GG = 1 π ImC GG . Similar to Eq. (2.1), for the (axial-)vector and tensor tetraquark currents J µ and J µν (to be defined later), one can introduce two-point correlation functions as For J P = 1 − vector particle and J P = 1 + axial vector particle, the correlation function Π V µν and Π A µν can be decomposed as where θ µν = g µν − qµqν q 2 and ω µν = qµqν q 2 . In this paper we use Π V (A) 1 and Π T 1 to construct sum rules, as they project out the spin-1 and spin-2 degrees of freedom we are interested in. The calculation of the corresponding ground state masses is similar to that in Eq. (2.9).
3 Calculation of C 1 and C GG In QCD Sum Rules, there are two kinds of expansions: the OPE and the perturbative expansion in α s . For the OPE, we only consider the most important contributions, the purely perturbative term C 1 and the GG condensate term C GG g 2 sĜĜ , because other higher dimensional operators are power suppressed in the OPE. According to Eq. (2.9), we need to calculate the imaginary parts of C 1 and C GG perturbatively. We can expect that the LO contribution of C 1 is the dominant one, and the next important contribution can be the NLO corrections for C 1 or the LO contribution of C GG . Therefore, the NLO corrections to C 1 need to be considered in the calculation in order to reduce theoretical uncertainties. For convenience, we will call the sum of the LO of C 1 and C GG as the LO contribution and the NLO corrections to C 1 as the NLO contribution in the following.
We use FeynArts [64,65] to generate Feynman diagrams and Feynman amplitudes of C 1 and C GG . Some representative Feynman diagrams at the LO and the NLO are shown in Fig. 1 and Fig. 2, respectively. The calculation procedure for C 1 and C GG are summarized below: where we use the column vector J to represent the basis in Eq. (4.1) or (4.2). A physical state can well be a mixture of all possible currents that share the same quantum numbers. The operator mixing has significant effects on the QCD sum rules calculations, say, for the heavy baryon spectrum [62]. However, there are no natural standards to pin down the mixing scheme only based on the LO calculation of C 1 and C GG . Thanks to the NLO calculations, the currents are mixed with each other naturally under the renormalization. If one choose the basis which diagonalizes the anomalous dimension matrix, then the operators in the basis have universal anomalous dimensions, separately. Thus, inserting these operators into the calculations in QCD sum rules, the dependence on the renormalization scale µ tends to be cancelled out in the righthand side of Eq. (2.9), which is desirable since the left-hand side M 2 H is a physical quantity. couple to the state with J P C = 0 −+ . The two types bases can be associated with each other by the Fierz transformation in 4 dimension, which is given as (4.10) Choosing the currents in Eq. (4.8) as the operator basis, one can get the anomalous dimension matrix (4.11) To diagonalize the matrix in Eq. (4.11), one needs the transformation matrix And one can get a unique set of the diagonalized currents The anomalous dimension matrix of J Dia P is diagonal, which is given by (4.14) 4.3 J P = 1 + For the J P = 1 + axial vector system, there are four independent interpolating currents. The operator basis, in the color singlet meson-meson type currents, can be chosen as where J M-M A,1 and J M-M A,2 couple to states with J P C = 1 ++ , while J M-M A,3 and J M-M A,4 couple to states with J P C = 1 +− . Alternatively, one can choose the diquark-antidiquark type currents [23] as the basis, which are given by (4.17) Choosing the currents in Eq. (4.15) as the operator basis, one can get the anomalous dimension matrix To diagonalize the matrix in Eq. (4.18), one needs the transformation matrix And one can get a unique set of the diagonalized currents (4.20) The anomalous dimension matrix of J Dia A is diagonal, which is given by For the J P = 1 − vector system, there are four independent interpolating currents. The operator basis, in the color singlet meson-meson type currents, can be chosen as couple to states with J P C = 1 −+ . Of course, one can choose the diquark-antidiquark type currents [23] as the basis, which are given by couple to states with J P C = 1 −+ .
In the calculation, J Di-Di V can be associated with J M-M V by Fierz Transformation in 4 dimension, which is given by If we choose the currents in Eq.
which make the anomalous dimension matrix A M-M V diagonal.

J P = 2 +
For the J P C = 2 ++ tensor system, there are three independent interpolating currents. The operator basis, in the color-singlet meson-meson type currents, can be chosen as (4.26) One could also construct the following operators but they won't contribute to Π T 1 . This suggests that these operators can not correspond to a tensor particle and we discard them from our analysis.
Of course, one can choose the diquark-antidiquark type currents [23] as the basis, which are given by (4.28) In the calculation, J Di-Di T can associate with J M-M T by the Fierz Transformation in 4 dimension, which is given by Choosing the currents in Eq. (4.26) as the operator basis, one can get the anomalous dimension matrix To diagonalize the matrix in Eq. (4.30), one needs the transformation matrix And one can get a unique set of the diagonalized currents (4.32) The anomalous dimension matrix of J Dia T is diagonal, which is given by Thus, all diagonalized currents can be determined uniquely.

Phenomenology
In our numerical analysis, we choose the following parameters [62,[76][77][78][79],  It is worth emphasizing that α s (µ) and the heavy quark mass m MS Q (µ) are obtained through two-loop running. Note that we don't need to consider the running of g 2 sĜĜ for the LO GG condensate contribution, as its anomalous dimension vanishes up to this order. As a typical choice, we set µ = M B in our phenomenological analysis [53,80], but the renormalization scale dependence will also be discussed. On-Shell (OS) masses m OS c and m OS b are extracted from the QCD sum rules analysis of the J/ψ and Υ(1S) spectrum, respectively, in which the mass renormalization scheme and truncation order of α s are the same as this paper.
According to Eq. (2.9), numerical result M H also depends on other two parameters: s 0 and M B . However, the physical value of M H should be independent of any artificial parameters. So a credible result should be obtained from an appropriate region where the dependence of s 0 and M B is weak. On the other hand, the choice of M B and s 0 should ensure the validity of the OPE and ground-state contribution dominance, which constrain the two parameters to be the so-called "Borel window". Within the Borel window, one should find the region, the so-called "Borel platform", in which M H depends on s 0 and M B weakly.
To search for the Borel window, we define the relative contributions of the condensate and continuum as and impose the following constraints: The two constraints guarantee the validity of OPE and the ground-state contribution dominance, respectively. In addition to the conditions given in Eq. (5.3), we also impose the following constrain on s 0 : around the point (x 0 , y 0 ) up to 10% in magnitude. It should be emphasized that the central point (x 0 , y 0 ) may lies on the margin of the Borel window in some cases. Therefore, the parameter space used to estimate errors of M H may exceed the Borel window, and also, the upper and the lower errors are usually asymmetric.

Numerical results and discussions for thecccc system
Our main results are shown in Fig. 3, where the 19 diagonalized currents, which should be more reasonable to be used in the QCD sum rules, are clustered by different quantum numbers. We set µ = M B and choose the MS renormalization scheme, and errors of M H include only that originated from uncertainties of s 0 and M 2 B . In the plot we also indicate the mass of X(6900) and the double J/ψ threshold. Let us first emphasize the importance of the NLO corrections. On one hand, NLO corrections to hadron masses are significant, which are larger than 0.5 GeV in both MS and on-shell schemes for almost all the currents involved in Tabs. 9-18. On the other hand, with the NLO corrections, the quark mass scheme dependence of M H tends to be reduced, especially for some diagonalized currents. To see this, we examine the difference of the predicted hadron masses between the two schemes, From Tabs. 9-18, for almost all the currents, one can find that the mass difference at LO is about ∆M LO H ≈ 1.2 GeV, which implies a roughly linear dependence of ∆M LO H on the quark mass difference between the two schemes. One can also find that the NLO corrections to M MS H are positive while those to M OS H are negative, therefore, the scheme dependence of M H tends to be reduced with the NLO corrections. Taking the J P = 0 + system as an example (see Tab. 9 and 10), for the three diagonalized currents J Dia S,2,3,4 , the NLO mass difference |∆M NLO H | < 0.4 GeV, which is explicitly smaller than that at LO. As for the currents J Dia S,1 and J Dia S,5 , the NLO corrections to M MS H are larger than 1 GeV, which implies that there are genuine large corrections other than the quark mass renormalization effects and the perturbation convergence may be bad for these currents.
The convergence of perturbation can also be explored by the µ dependence of the NLO results. One would expect that the µ dependence of the NLO result will be significantly reduced in comparison with the LO one for the current for which the perturbation convergence is good, since the truncation of the perturbation series up to NLO has weak effect on the result. On the other hand, the large µ dependence of the NLO result may imply that the perturbation convergence is bad, say, the next-to-next-to-leading-order (NNLO) corrections should be important in this case. In Tabs. 9-18, we have chose µ = M B in the MS scheme. To study the µ dependence, we vary µ = k M B with k ∈ (0.8, 2.0), where the range is chosen with the requirement that the Borel platform can be achieved and the perturbative expansion is under good control. We investigate the µ dependence for the diagonalized operators with J P C = 0 ++ , which are shown in Fig. 26 As for states other than J P C = 0 ++ , there are only three diagonalized operators that satisfy ∆M NLO the µ dependence of NLO results are improved significantly, compared with that of the LO one.
In cases where convergence of perturbation is bad, uncertainties from higher order corrections should be large. As higher order corrections, such as the NNLO ones for C 1 , are beyond the scope of this paper, we will only choose diagonalized operators that have good perturbation convergence in the following analysis. For the results of the diagonalized operators J Dia S,2,3,4 with J P C = 0 ++ , J Dia P,2 with J P C = 0 −+ , J Dia A,4 with J P C = 1 +− and J Dia T,1 with J P C = 2 ++ , we also estimate the uncertainties coming from errors of the quark masses in Eq. (5.1), and the uncertainties are shown in Tabs. 1-6.

Current
Order

Current
Order  Table 3.

Current
Order

Current
Order

Current
Order  Table 6. Phenomenologically, it is interesting to compare our calculations with the LHCb measurements of the possiblecccc tetraquark states in the J/ψJ/ψ spectrum [5]. The most likely quantum numbers J P C for the tetraquark states are 0 ++ and 2 ++ , since they can couple to J/ψJ/ψ in S-wave. The predicted NLO MS masses for the two operators J Dia S,3 and J Dia S,4 are 6.35 +0.20 −0.17 GeV and 6.56 +0.18 −0.20 GeV, respectively, which might account for the broad structure around 6.2 ∼ 6.8 GeV measured by the LHCb collaboration [5]. As for the narrow resonance X(6900) [5], the central value of the mass is consistent with the NLO MS mass for the operator J Dia S,2 , which gives 6.95 +0.21 −0.31 GeV. Moreover, the predicted NLO MS mass, 7.03 +0. 22 −0.26 GeV, for the operator J Dia T,1 with J P C = 2 ++ is also close to that of X(6900), so we can not assert that the quantum number of X(6900) is 0 ++ , while 2 ++ may also be possible. Since the quality of the Borel platform in On-Shell scheme is worse than that in MS one, which can be seen from Figs. 7-25, we only use the corresponding MS masses in the above analysis. As for the NLO On-Shell masses of J Dia S,2,3,4 and J Dia T,1 (see Table 1-3 and 6), they all lie on the the broad structure around 6.2 ∼ 6.8 GeV measured by the LHCb collaboration [5].
Since the abovecccc states are pure heavy-quark systems, their non-relativistic (NR) attributes should be important to understand them. In our calculations, although the amplitude are calculated in full QCD for the covariant operators given in Sec. 4, the results still exhibit some NR features. Taking   component, and should be degenerated in the NR limit. This is roughly the case in our result, the NLO MS masses of the two states are roughly equal, and from Tab. 9, one can see that they are both close to the NLO MS mass for J M-M S,4 (6.36 +0.06 −0.10 GeV), and are not consistent with that for J M-M S,3 (7.91 +0.16 −0.19 GeV). Similarly, in the NR limit, the state for J Dia T,2 is dominated by its J M-M T,1 component (see Eq. (4.31)), which leads to dimension 6 operator and can survive in the limit. Correspondingly, from Tab. 17 one can see that the mass of J Dia T,2 is close to that of J M-M T,1 , and is not consistent with that of J M-M T,2 , which leads to dimension 8 operator in the NR limit.
To see the NR behaviors of the amplitudes more explicitly, we define v = 1 − exhibit similar behaviors, especially in the near threshold region, where v is small. For the above four operators, more explicitly analysis indicates that the near threshold behaviors of ρ LO , respectively, which can be roughly seen in Fig. 4. This is to say that ρ NLO 1 is enhanced by a factor of v −2 with respect to ρ LO 1 in the near threshold region. Because of the exponential suppression of the lager v region and the threshold parameter s 0 , the dominant domain of the integration in Eq. (2.9) corresponding to v = 0.4 ∼ 0.7 for thecccc system, where the NLO contributions are comparable with the LO ones. This indicates the importance of the NLO corrections to the QCD sum rules for thecccc system.
Finally, let's compare our predictions with those given in other works [23,27,29] within the framework of QCD sum rules. The predicted masses of thecccc tetraquark states are listed in table 8. The authors of Ref. [23] adopt momentum sum rules rather than Laplace sum rules (i.e. Borel transformation) applied here. Thus, it is difficult to compare the results between theirs and ours. The Laplace sum rules were applied in Ref. [27,29]. However, the parameters (such as the renormalization scale µ) and the scheme to determine the Borel platform in Ref. [27,29] are different from ours. Moreover, only partial NLO contributions are considered in Ref. [27].    Fig. 35 (a), there is no Borel platform at LO level in MS scheme, but there is a clear and distinct platform at the NLO level. Similar phenomenon was also found in the bbb system [63]. We have checked that thebb system also has this phenomenon. This indicates that for the pure bottom system the NLO contribution is crucial to the formation of a stable Borel platform in the QCD sum rules.
Similar to the case ofcccc system, from Tabs. 19-28, one can see that NLO contributes non-negligible corrections, with mass corrections |M NLO 6 GeV in both MS and OS schemes. In addition, with the NLO corrections, the quark mass scheme dependence is improved significantly. The mass difference between the two schemes GeV at LO level, while the difference is usually smaller than 0.1 GeV at NLO level except for the J P C = 1 ++ channel.
We choose µ = k m B with k ∈ (0.8, 1.2) to explore the renormalization scale dependence of our results, because Borel platforms can not be achieved for k > 1.2 even with the NLO contributions. We find that the µ dependence is improved for the NLO results comparing with the LO ones, but the µ dependence of the NLO results for thebbbb system is more sensitive than that of thecccc system. Typical µ dependence at the LO and the NLO is shown in Fig. 54  1. × 10 -16 has been enlarged by factor 4 and 241 225 , respectively, to be compared with that of J M-M S,4 . Similar to thecccc system, the near threshold behaviors of ρ LO 1 , ρ NLO 1 and ρ GG for the above fourbbbb operators in MS scheme are of O(v 7 ), O(v 5 ) and O(v), respectively, which can be roughly seen in Fig. 6. However, the dominant domain of the integration in Eq. (2.9) corresponding to v = 0.2 ∼ 0.4 for thebbbb system, which is smaller than that for thecccc system. This is consistent with the general expectation that thebbbb system is more like a NR one than thecccc system since, in roughly speaking, m b m c Λ QCD . Due to the relative enhancement of ρ NLO 1 with respect to ρ LO 1 in the near threshold region, the NLO correction to the QCD sum rules is crucial for thebbbb system, which have been indicated by its effect on improvement of the quality of the Borel platform. On the other hand, the enhancement of the QCD perturbative correction in the near threshold region may make the perturbation convergence bad for thebbbb system within QCD sum ruls, which have been indicated by the µ-dependence of m H up to the NLO corrections showed in Fig. 54 and 55. In other words, the NNLO corrections may be important and the enhancement in the near threshold region may need to be resumed for thebbbb system within QCD sum ruls. But this is far beyond the scope of this work.
At last, we have to emphasize that there are large NLO corrections to the operator J M-M A,2 . We find that the near threshold behaviors of ρ LO 1 and ρ NLO 1 for this operator are of O(v 11 ) and O(v 7 ), respectively. This means for this operator, the ρ NLO 1 is enhanced by a factor of v −4 with respect to ρ LO 1 in the near threshold region. This enhancement is more serious forbbbb system than that forcccc system since the typical value of v is smaller for the former. Thus, the perturbation convergence is very bad for this operator, which may be indicated by the large NLO corrections and lager errors of the MS mass m H of this operator forbbbb system (see Table 23). On the other hand, since the dominant component of J Dia

Summary
In this paper, we study the NLO corrections to masses ofQQQQ states within QCD sum rules. As operators with the same J P C can mix with each other under renormalization, we diagonalize the original operators, either in meson-meson type or diquark-antidiquark type, and use the diagonalized operators in the phenomenological study. with J P C = 2 ++ , which also implies that the perturbation convergence of these operators is better than that of the others. While for thebbbb system, the difference is usually smaller than 0.1 GeV at NLO except for the J P C = 1 ++ operator. We also find that NLO corrections can significantly reduce the µ dependence. We use currents that have good perturbative convergence in our phenomenological analysis. For thecccc system, we get three J P C = 0 ++ states, with masses 6.35 +0.20 −0.17 GeV, 6.56 +0. 18 −0.20 GeV and 6.95 +0.21 −0.31 GeV, respectively. The first two may explain the broad structure around 6.2 ∼ 6.8 GeV measured by the LHCb collaboration [5], and the third one may be assigned to the observed narrow resonance X(6900). For the 2 ++ states, we find one with mass 7.03 +0. 22 −0.26 GeV, which may also be a candidate for the X(6900), and another one around 7.25 +0. 21 −0.35 GeV, which has good µ dependence but slightly large scheme dependence (with ∆M NLO H 0.7 GeV). As for thebbbb system, we find that the NLO contribution improves the quality of the Borel platform evidently in MS scheme, which is similar to the case of the bbb baryon [63]. The quark mass scheme dependence of the results are also improved significantly with the NLO contribution. However, the NLO results are still sensitive to the choice of renormalization scale µ, and we find that the Borel platforms can not be achieved for µ > 1.2m B .
Finally, we would like to emphasize the importance of NLO contributions, especially in the operator mixing or color configuration mixing for multiquark systems. (i) As a key NLO contribution, the one-gluon exchange is crucial even for charmonium and bottomonium states, because it provides the color Coulomb interaction between Q andQ, which is the most important short-range attractive force to form a heavy quarkonium. (ii) In the fully heavy tetraquark system discussed in this paper, if one starts from a colorsinglet current-current operator, the one-gluon exchange will change it to the color-octet current-current operator, therefore leads to the operator mixing. As already shown by our result, the operator mixing induced by renormalization at NLO is inevitable and has very important consequences in the QCD sum rule calculations. (iii) In the literature some works use the color-singlet current-current local operators to describe physical hadronic molecules. However, duo to the operator mixing, the color structure of the local operators must be mixed with both color-singlet and color-octet current-current configurations. It is impossible to keep the color-singlet structure unchanged if a complete NLO QCD contribution is seriously considered. In fact, a physical molecule state means that it contains two well separated color-singlet mesons at long-distances mediated by one-meson or two-meson exchanges. And a physical molecule may not be necessarily ascribed to the color-singlet current-current local operators, which only describe the very short-distance behavior of the tetraquark and are subjected to the color configuration mixing. The description for hadronic molecules needs to understand the long-distance dynamics beyond color confinement.

Acknowledgments
We thank Xiao Liu and Xin Guan for many useful and helpful discussions. We also thank Shi-Lin Zhu for helpful comments. K.T.C. thanks Cong-Feng Qiao for useful communications. The figures in this paper are drawn by using the Origin and Mathematica software.

A.1 Calculation Of Operator Renormalization Matrices
We present the calculation of operator renormalization matrices of meson-meson type operators. The operator renormalization matrices of diquark-antidiquark type operators then follow from a Fierz transformation.
A general meson-meson type operators where four quarks are different flavors, are defined as which has two independent color configurations, is called a color-singlet operator. We can also use following relation, to obtain color-octet operators, For convenience, we choose O Γ 1 ,Γ 2 , [1] and O Γ 1 ,Γ 2 , [2] as bases in our calculation. Let us first suppress the dependence of color configuration for operators. According to our operator and definition of operator renormalization matrix, denotes the bare operator replaced by renormalized fields, and O R denotes the renormalized operator. In our NLO calculation, we directly calculate O B , and thus quark self-energy diagrams are cancelled by counter term diagrams. The remaining diagrams can be divided into three parts, where A 1 denotes the contribution of gluon exchange between q 1 and q 2 , A 2 denotes the contribution of gluon exchange between q 3 and q 4 , and A 3 denotes others contributions e.g. contributions of gluon exchange between q 1 and q 3 , q 1 and q 4 and so on. Because all infrared divergences will be cancelled, we just need to consider ultraviolet (UV) divergences therein. Therefore, the mass terms in quark propagators can be discarded. Explicitly, we have After a simple manipulation, we get (A.13) According to Eq. (A.12), we get the UV divergences term (A.14) For operators with definite color configuration O Γ 1 ,Γ 2 ,[c] , we need to multiply the corresponding color configuration (δ ij δ kl or δ il δ kj ) in Eq. (A.13).
According to Eq. (A.7), to use the renormalized operator we should multiply our result where M LO = 1 is the LO amplitude and Z 2 = 1 + δZ 2 and Z O = 1 + δZ O . Demanding that final results are free of UV divergences, we get 16) in MS scheme.
A.2 J P = 0 + Operator bases are defined as (A.17) The corresponding operator renormalization matrix is given by, After renormalization, since there are identical particles in the operator of full heavy tetraquark system Q Γ 1 QQΓ 2 Q , 4-dimensional Fierz transformation can be used to related operators in different color configurations, which results in only 5 independent operators in J P C = 0 ++ channel. We can choose any 5 independent operators to perform our phenomenological study. For example, we choose J M-M S,i, [1] in this work. According to Fierz transformation we can transform J M-M S,i, [2] to J M-M S,i, [1] to get the anomalous dimension Eq. (4.4), (A.20) The operator renormalization matrix is Similay to J P C = 0 ++ , we have the Fierz transformation The operator renormalization matrix is (A.24) Similar to J P C = 0 ++ , the Fierz transformation [2] to J M-M A,i, [1] and we thus get the anomalous dimension Eq. (4.18).

(A.26)
And the operator renormalization matrix is According to Fierz transformation we can transform J M-M V,i, [2] to J M-M V,i, [1] to get the anomalous dimension Eq. (4.18).

(A.29)
And the operator renormalization matrix is According to Fierz transformation we can transform J M-M T,i, [2] to J M-M T,i, [1] to get the anomalous dimension Eq. (4.30).

B Details forcccc system
B.1 Numerical Results for J P = 0 + states Table 9. The LO and NLO Results for J P = 0 + withcccc system in the MS scheme                         (b) OS Figure 41. The Borel platform curves for J Dia P,2 with J P C = 0 −+ in the MS and On-Shell schemes 8   (b) OS Figure 43. The Borel platform curves for J Dia A,1 with J P C = 1 ++ system in the MS and On-Shell schemes 10  (a) MS       (a) MS    (a) MS  (a) MS   (a) MS