Z c and Z cs systems with operator mixing at NLO in QCD sum rules

: We study the mass spectra of hidden-charm tetraquark systems with quantum numbers ( I G ) J P = (1 + )1 + using QCD sum rules. The analysis incorporates the complete next-to-leading order (NLO) contribution to the perturbative QCD part of the operator product expansions, with particular attention to operator mixing effects due to renormalization group evolution. For the ¯ dc ¯ cu system, the masses of two mixed operators, J Mixed1 , 5 and J Mixed2 , 6 , are determined to be 3 . 89 +0 . 18 − 0 . 12 GeV and 4 . 03 +0 . 06 − 0 . 07 GeV, respectively, closely matching those of Z c (3900) and Z c (4020). Similarly, for the ¯ sc ¯ cu states, the masses of J Mixed1 , 5 and J Mixed2 , 6 are found to be 4 . 02 +0 . 17 − 0 . 09 GeV and 4 . 21 +0 . 08 − 0 . 07 GeV, respectively, in close proximity to Z cs (3983)/ Z cs (4000) and Z cs (4220), consistent with the expectation that they are the partners of Z c (3900) and Z c (4020). Our results highlight the crucial role of operator mixing, an inevitable effect in a complete

In recent years, a large number of new hadronic states containing heavy quarks (the charm quark c or bottom quark b) have been observed at hadron colliders and e + e − colliders [1].
They are expected to be candidates of tetraquark states, pentaquark states, and other exotic states which contain two heavy quarks [2][3][4][5][6][7][8].These findings have opened up a new stage for the study of hadron physics and QCD.Among these exotic states, the Z c (3900) is a good candidate of tetraquark state, and this state can help us to better and deeper understand tetraquark structure and properties.
The QCD Sum Rules [68,69] approach is a powerful tool to study hadronic properties [70][71][72][73].Currently, there have been many leading order (LO) in α s calculations for the Z c(s) system [20, 21, 27, 33, 36-41, 74, 75].But it is well known that the LO results have large theoretical uncertainties, e.g. the ad hoc choice of quark masses and QCD strong coupling constant, and that the next-to-leading order (NLO) QCD corrections to the perturbative part C 1 may significantly reduce the theoretical errors and bring about sizable corrections to the results.This has been emphasized in many works, e.g. for the proton [76,77], the singly heavy baryons [78], the doubly heavy baryon Ξ ++ cc [79], the fully heavy baryons Ω ++ ccc and Ω − bbb [80], the fully heavy tetraquarks QQ QQ (Q = c, b) [81], and other systems (based on partial NLO contributions) [39,40,82,83].
In general, the NLO contribution of C 1 may lead to large corrections to the results and substantially reduce the dependence of renormalization scale and scheme as well as other parameters.This has been shown in our studies at NLO in α s for the doubly heavy baryon Ξ ++ cc [79], and fully heavy baryons Ω ++ ccc and Ω − bbb [80].Moreover, according to the NLO study for the fully heavy tetraquark ( QQ QQ) [81], we find that the NLO contribution can lead to the operator-mixing (color-configuration mixing) and its effect can be vitally important in understanding the color-structure of the states and phenomenological results.In contrast, at LO in α s , there exist infinite operator-mixing schemes, so in the previous LO work [20, 21, 27, 33, 36-38, 41, 74, 75] it is hard to deal with this problem.Whereas at NLO, we can naturally determine the operator-mixing scheme by diagonalizing the anomalous dimension matrix of those operators.Based on our previous studies, we expect that the NLO contribution and operator-mixing effect are also important for the Z c system.So far in the literature there exist no complete calculations of the NLO contribution of C 1 , and in Ref. [39,40] only partial NLO contribution originating from the so-called factorized diagrams, was considered.In contrast, we will perform a complete NLO calculation for C 1 , and study the induced operator-mixing or color configuration mixing effects, and do phenomenological analysis.This will be the aim of our present paper.
The rest of the paper is organized as the following.In Sec. 2, QCD sum rules for calculation of the mass of uc dc are given.In Sec. 3, we present our methods to calculate perturbative coefficients.In Sec. 4, we introduce the current operators and the operatormixing scheme.Phenomenological results and discussions are given in Sec. 5. Some details of our calculations and results are given in Appendices A, B.

QCD Sum Rules
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. [70] for more details.We start with a two-point correlation function where D denotes the space-time dimension, Ω denotes the QCD vacuum and J is a tetraquark current operator 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 [70], 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 the 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 five, which gives the approximated expression of OPE as where ⟨g s qqG⟩ denotes ⟨g s qσ • Gq⟩.
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 where s th = 4m 2 c is the QCD threshold for the dccu system (light quarks are considered massless in Wilson coefficients), and the integral in the second line has been assumed to be convergent.Then by employing the quark-hadron duality and Borel transformation [70], we obtain a sum rule for Π(q 2 ), , where ρ 1 = 1 π ImC 1 and ρ i = 1 π ImC i .Similar to Eq. (2.1), for the axial vector tetraquark currents J µ (to be defined later), one can introduce two-point correlation functions as For J P = 1 + axial vector particle, the correlation function Π A µν can be decomposed as .11)In this paper we use Π A 1 to construct sum rules, as it contains the genuine spin-1 degrees-of-freedom we are interested in.The calculation of the corresponding ground state masses is similar to that in Eq. (2.9).

Calculation of Wilson coefficients C 1 and C i
In QCD Sum Rules, there are two kinds of expansions: the OPE and the perturbative expansion in α s .For the OPE, we only consider vacuum condensates up to dimension five (see Eq. 2.5) because other higher dimensional operators are further power suppressed in the OPE, and those condensates values are not well determined.According to Eq. (2.9), we need to calculate the imaginary parts of C 1 and C i perturbatively.If the OPE is valid for the correlation function under consideration, then the perturbative contribution C 1 should dominate the whole expansion.And the next important contribution can be the NLO corrections for C 1 or the LO contribution of condensates C i .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 condensates C i as the LO contribution and the NLO corrections to C 1 as the NLO contribution in the following.
We use FeynArts [84,85] to generate Feynman diagrams and Feynman amplitudes of C i .Some representative Feynman diagrams at the LO and the NLO are shown in Fig. 1 and Fig. 2, respectively.Here we emphasize that Fig. 2 includes all NLO diagrams, in particular those one gluon exchange diagrams between cd and cu, without which the gauge invariance can no longer be maintained.This is a crucial point that distinguishes our work from Ref. [39,40], in which only the so-called factorized diagrams are considered to intentionally keep the (c d) 1 -(cu) 1 color singlet-singlet configuration for the c dcu system.
Our calculation procedure for C 1 and C i are summarized below: • 1.We use FeynCalc [86,87] to simplify spinor structures of Feynman amplitudes with the Larin scheme [88].
• 2. We use Reduze [89] to reduce all loop integrals to linear combinations of a set of simpler integrals, which are called master integrals (MIs).
• 4. Renormalization.There is no infrared divergence in the NLO amplitude of C 1 .
After performing the wave-function and mass renormalization of quarks (m Q is renormalized in either the MS scheme or the on-shell scheme), the remaining ultraviolet divergences can be removed by the renomalization of the current operators.When there are more than one current operators sharing the same quantum numbers J P C , they are usually mixed with each others under the renormalization.We get operator renormalization matrices in the MS scheme, which are shown explicitly in Appendix A.
Because the expressions of C 1 and C i are too complicated to be shown in the paper, we attach the imaginary part of C 1 and C i , which are the only needed information in phenomenological studies, as ancillary files.

Mixed operators
A physical state can couple to any current that shares the same quantum numbers, therefore, it can also couple to a mixture of these currents.And the operator mixing has significant effects on the QCD sum rules calculations, say, for the heavy baryon spectrum [79] and fully heavy tetraquark spectrum [81].However, the LO calculation of the Wilson coefficients alone provides no guidance to pin down the mixing scheme of the operators.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, the dependence on the renormalization scale µ tends to be cancelled out on the righthand side of Eq. (2.9), which is desirable since the left-hand side M 2 H is a physical quantity.
In this mixing scheme, no matter with which type (Eq.(4.1) or (4.2)) operators basis we start from, we always get the same mixed operators.In this paper, we choose the meson-meson type currents in Eq. (4.1) as the operator basis, the operator anomalous dimension matrix γ M-M can be obtained according to the operator renormalization matrix Eq. (A.17), where δ = − αs 16π .To diagonalize the matrix in Eq. (4.4), one needs the following transformation matrix .5) So, we get the new basis And the anomalous dimension matrix of ⃗ J Dia is diagonal, which is given by From Eq. (4.7), one can see that the eigenvalues of anomalous dimension matrix are degenerate.Since the operators with the same eigenvalue can be further mixed with each other, we defined the following operators as where the parameters θ i (i = 1, 2, 3, 4) and ϕ i (i = 1, 2, 3, 4) are real.
In the following, we call ⃗ J Mixed as mixed operators.

Phenomenology
In our numerical analysis, we choose the following parameters [79,[100][101][102][103], where q denotes light quark u, d.Moreover, it is worth emphasizing that consistent with our calculation, α s (µ) and the heavy quark mass m MS c (µ) are obtained through two-loop running, where µ is the renormalization scale.⟨qq⟩(µ) and ⟨g s qqG⟩(µ) are obtained through one-loop running, and their anomalous dimensions are given by [104] Since the anomalous dimension vanishes at one-loop level, we don't need to consider the running of the GG condensate ⟨g2 s GG⟩.As a typical choice, we set µ = M B in our phenomenological analysis [68,105], but the renormalization scale dependence will also be discussed.
According to Eq. (2.9), numerical result M H also depends on other two parameters: s 0 and M 2 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 2 B is weak.On the other hand, the choice of M 2 B and s 0 should ensure the validity of the OPE and ground-state contribution dominance, which constrain the two parameters to be within 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 2 B weakly.
To search for the Borel window, we define the relative contributions of the condensate and continuum as where ⟨O i ⟩ ∈ {⟨qq⟩, ⟨g 2 s GG⟩, ⟨g s qqG⟩} 2 .And we impose the following constraints: To guarantee the validity of OPE and the ground-state contribution dominance, there are constraints for t r and t c , that t r ≤ 1 and t c ≤ 0.5.In addition to the conditions given in Eq. ( 5.4), we also impose the following constraint on s 0 : since, roughly speaking, s 0 denotes the energy scale where the continuum spectrum begins to contribute and that the binding energy for a hadron is usually smaller than 1 GeV.
To find the Borel platform, we search for the point where the parametric dependence of M H is weakest within the Borel window.More explicitly, we choose the variables as x = s 0 and y = M 2 B and define the function (5.6) By minimizing the function ∆(x, y) within the Borel window and with the constraint Eq. (5.5), we get a point (x 0 , y 0 ), which will be used to calculate the central value of M H .To estimate errors of M H , we vary the values of s 0 and M 2 B 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 for meson-meson type and diquark-antidiquark type operators of Z c system
The comprehensive results for both meson-meson type (see Eq. (4.1)) and diquark-antidiquark type (see Eq. (4.2)) operators are shown in Tab.9-12 in Appendix B, where we include both LO and NLO, both MS scheme and on-shell scheme for the results of Z + c mass M H .In these tables, we set µ = M B and thus the errors of M H are due to choices of s 0 and M 2 B .Further information of s 0 and M 2 B dependence is shown in Fig. 9-24 in Appendix B. In these plots, a black dot denotes the central point (x 0 , y 0 ), and shadows denote the Borel window determined by Eq. (5.7).
We also summarize the results of Tab.9-10 (for meson-meson type operators) and Tab.11-12 (for diquark-antidiquark type operators) in Fig. 3 and Fig. 4, respectively.Through these figures, one can see that the NLO corrections are important, especially for the on-shell scheme.However, we should emphasize that to obtain the Borel windows for both meson-meson type and diquark-antidiquark type operators, one needs very loose constraints for the Borel conditions given in Eq. ( 5.3), which are corresponding to (5.7) Thus, the convergence of the OPE may be bad, and the results shown in Fig. 3 and Fig. 4 may suffer severe theoretical uncertainties.Furthermore, the qualities of these Borel platforms are bad (see Fig. 9-24), which cause large errors of M H 's from their dependence on s 0 and M 2 B (see Tab. 9-12).In addition, although the NLO corrections tend to reduce the mass gap between MS scheme and on-shell scheme, the NLO mass difference M NLO-OS H − M NLO-MS H > 0.5 GeV for most of these operators (see Fig. 3 and Fig. 4), which indicates that the perturbative convergence of the results for these operators may be bad.
Due to large uncertainties of these results, we will not use the masses shown in Fig. 3 and Fig. 4 to do further phenomenological analysis, although we have shown the mass of Z c (3900) as a horizontal line in these figures.
Z c ( 3 9 0 0 ) Z c D i q u a r k -A n t i d i q u a r k T y p e O p e r a t o r s

Numerical results for mixed operators of Z c system
As have been mentioned in the subsection 4.2, the mixed operators given in Eq. (4.8) are more adequate to be used in the QCD sum rules.On the one hand, due to the µdependence tends to be cancelled out in the righthand side of Eq. (2.9) for the mixed operator, the corresponding results are expected to have better perturbative convergence compared to those given in last subsection.On the other hand, one can reduce the s 0 -and M 2 B -dependence of the results by choosing ideally mixing parameters θ i and ϕ i defined in Eq. (4.8) to get Borel platforms with high quality.
We scan in the parameter spaces of θ i and ϕ i , and determine the following mixing parameters through minimizing s 0 -and M 2 B -dependence of the results: Thus we determine four ideally mixed operators by inserting Eq. (5.8) into Eq.(4.8).Comparing with the un-mixed operators given in Eq. (4.1) and Eq.(4.2), for these ideally mixed operators we can find better Borel windows even if we choose relative strict constraints (5.9) Thus, both the convergence of the OPE and the validity of the ground state dominance are guaranteed better after considering mixed effect.
The results for these ideally mixed operators are listed in Tab.13-14 in Appendix B, where we include both LO and NLO, both MS scheme and on-shell scheme for the results of Z + c mass M H .In these tables, again, we set µ = M B and thus errors of M H are only due to choices of s 0 and M 2 B .Further information of s 0 and M 2 B dependence is shown in Fig. 25-28 in Appendix B. In these plots, a black dot denotes the central point (x 0 , y 0 ), and shadows denote the Borel window determined by Eq. (5.7).Form Fig. 25-28, one can see the quality of the Borel platforms are evidently better than those for the un-mixed operators, thus, the errors of the results listed in Tab.13-14 are vastly smaller than those in Tab.9-12.
We also explore the m c -and the µ-dependence of the above results, and the corresponding errors coming from these dependence are listed in Tab.1-4 for the four ideally mixed operators, respectively.To explore the µ-dependent of the MS mass M MS H , we set µ = kM B and vary the value of k from 0.8 to 2.0.The µ-dependent of both the LO and NLO MS masses are shown in Fig. 6.Furthermore, we also summarize the M H results of Tab.1-4 with errors in Fig. 5.

Current
Order Error from s 0 and M

Current
Order Error from s 0 and M    , one can find that the NLO corrections are important for both MS and on-shell schemes, and the NLO corrections largely reduce the mass gap between the two schemes.To be specific, for these two operators (see Tab. 1-2), the LO mass central value difference of the two schemes 05 GeV.This indicates that the perturbative convergence for the results of these two operators are very good, which is also reflected in the weaker µ-dependence of the NLO MS masses than that of the LO ones of these two operators (see Fig. 6 (a) and (b)).
As for the last two ideally mixed operators J Mixed 3,7 and J Mixed 4,8 , with the NLO corrections the mass gap between the two schemes are not reduced, and is even enlarged for the J Mixed 4,8 operator.Correspondingly, the µ-dependence is not improved with the NLO corrections for these two operators (see Fig. 6 (c) and (d)), and is even worse for the J Mixed 4,8 operator.All the above indicate that the perturbative convergence for the results of these two operators are bad, and the higher order corrections may be important for these two operators, which are beyond the scope of this paper.
Phenomenologically, the NLO mass of J Mixed 1,5 in the MS scheme is given by 3.89 +0.18 −0.12 GeV (see Tab. 1), which is very close to that of Z c (3900), and the NLO mass of J Mixed 2,6 in the MS scheme is given by 4.03 +0.06 −0.07 GeV (see Tab. 2), which is very close to that of Z c (4020).
Thus our results supports that there are two Z c states with the same quantum numbers (I G )J P = (1 + )1 + .But it is not sure if there are Z c states other than Z c (3900) and Z c (4020) through the NLO calculations in QCD sum rules, since the results other than those for J Mixed 1,5 and J Mixed 2,6 suffer large uncertainties and thus are not reliable to make any predictions.Even for the operators J Mixed , the uncertainties of the NLO masses due to renormalization scheme dependences still remain significant, at approximately 0.05 GeV and 0.03 GeV (see Tab. 1 and Tab. 2), respectively.This may introduce potential systematic errors into our predictions of mass at the NLO level.Nevertheless, it is conceivable that the discrepancy between the MS scheme and the on-shell scheme could be further mitigated by higher-order QCD corrections.
One may think that the threshold parameter s 0 is a little bit large for the NLO result for J Mixed 1,5 ) in the MS scheme.That is, for the central values of M NLO-MS H and the corresponding s 0 (see Tab. 1), √ s 0 − m H = 0.80 GeV, which is larger than the typical energy gap ∆(∼ 0.5 GeV) between the ground state and the first radial exited state for a hadronic system.However, one can see from Tab. 1 that our prediction on the NLO-MS mass is not very sensitive to the threshold parameter, which can also be seen from Fig. 25(a).That means, if we choose s 0 = 19.8GeV 2 , then the central value of the mass is about M H = 3.87 GeV, and √ s 0 − m H = 0.58 GeV, which approaches to the value of ∆.

Numerical results for mixed operators of Z cs system
In this paper, we regard Z + cs (3985)/Z + cs (4000) and Z + cs (4220) as the partners of Z + c (3900) and Z + c (4020), respectively.Thus, all the calculation details about Z cs system are almost the same as Z c system.One needs just replace down quark d with the strange quark s in both the operators and the quark condensate.For numerical analysis associating to s quark, we choose the following parameters [103]  (5.10) where ⟨qq⟩(2GeV) has been listed in Eq. (5.1).Similar to the case of Z c system, the QCD sum rules for the un-mixed operators of the Z cs system are not reliable even with NLO corrections.Firstly, to obtain the Borel windows, one needs very loose constraints as those given in Eq. (5.7), and therefore the convergence of OPE is bad and the theoretical uncertainties may be very large for these operators.Secondly, the quality of the Borel platforms are bad and the results are sensitive to the choice of s 0 and M 2 B .Thirdly, the perturbative convergence for the results of these operator are bad.Thus, we will not list the results for these unmixed operators here.
As for the mixed operators, we choose the same mixing parameters as those given in Eq. (5.8) and the same Borel conditions as those given in Eq. (5.9).The results are listed in Tab.5-8 for the four ideally mixed operators, respectively.For µ = M B , the curves of the Borel platforms are shown in Fig. 29-32 in Appendix B. In these plots, again, a black dot denotes the central point (x 0 , y 0 ), and shadows denote the Borel window determined by Eq. (5.9).The quality of the platforms are quite good, thus, the errors of M H form s 0 and M 2 B are very small (see Tab. 5-8).In Tab.5-8, we also list the errors coming from m c and µ.Again, to explore the µ-dependence of M H , we set µ = kM B and vary the value of k from 0.8 to 2.0.The µ-dependent of both the LO and NLO MS masses are shown in Fig. 8. Furthermore, we also summarize the M H results with errors in Fig. 7.

Current
Order M H (GeV) s 0 (GeV  of Z cs system in MS and On-Shell schemes.Here the errors for M H are from s 0 , M B , the charm quark mass, and the renormalization scale µ with µ = kM B and k ∈ (0.8, 1.2) (the central values correspond to µ = M B ).

Current
Order M H (GeV) s 0 (GeV   , perturbative convergence of the results are bad, and the µ-dependence is not improved evidently with the NLO corrections for these two operators (see Fig. 6 (c) and  (d)), and is even worse for the J Mixed are reliable for the Z cs system.Phenomenologically, the NLO mass of J Mixed 1,5 in the MS scheme is given by 4.02 +0.17 −0.09GeV (see Tab. 5), which is very close to that of Z cs (3985)/Z cs (4000), and the NLO mass of J Mixed   2,6   in the MS scheme is given by 4.21 +0.08 −0.07 GeV (see Tab. 6), which is very close to that of Z c (4220) 3 .Thus, our results support to assign Z cs (3985)/Z cs (4000) and Z cs (4220) as the partners of Z c (3900) and Z c (4020), respectively.
It should be mentioned here that, in our calculation, we set m s = 0 in evaluation of the perturbation part C 1 , the non-zero m s given in Eq. (5.10) is only used to evaluate the condensate parts C ss ⟨ss⟩ and C ssG ⟨g s ssG⟩.And we choose r sq = ⟨ss⟩/⟨qq⟩ = 0.6 ± 0.1 to approximately fit the mass of Z cs (3985)/Z cs (4000) or that of Z cs (4220).In Ref. [108], the ratio r sq is fitted to be 0.8 ± 0.1 by the mass difference between Σ * and ∆ baryons with a fixed strange quark mass m s = 0.15 GeV, which is different from our choice in Eq. (5.10).If we choose r sq = 0.8 ± 0.1, the central values of the M H for the mixed operators will be increased by around 0.15 GeV.
In the above calculation, we have taken Z cs as the partner of Z c .If we deal with Z cs independently, the number of operators in the basis will be doubled since there are no constraint from G-parity as that for the Z c system (see Eq. (4.1) or Eq.(4.2)).Thus, for J P = 1 + cusc system, there are 16 operators in the basis, and the diagonalized anomalous dimension matrix is fourfold degenerated.As a result, the mixing of the operators will be complex for this system.For example, there will be two extra operators which can be mixed with J Mixed 1,5 further.However, the operator J Mixed 1,5 is still corresponding to a special choice of the mixing parameters in the four operator mixing space, which can give result with good perturbative convergence as one had seen.In this aspect, our calculation still supports to assign Z cs as the partner of Z c .

Summary
In this paper, we study the NLO corrections to the spectrum of dccu tetraquark system with (I G )J P = (1 + )1 + , i.e., the Z c system, in QCD sum rules.As operators with the same quantum numbers can mix with each other under renormalization, we recombine the original operators, either in meson-meson type or diquark-antidiquark type, by diagonalizing the anomalous dimension matrix.However, different from the case of fully heavy tetraquark system [81], the α s -order anomalous dimension matrix in Z c system is degenerated and the mixing of operators can not be determined uniquely by diagonalization of the matrix.We then determine the mixing parameters by minimizing the s 0 -and M 2 B -dependence to give the four so-called ideally mixed operators.And we call the operator in the original basis (either in meson-meson type or diquark-antidiquark type) as the unmixed one.
Numerically, we find the results for almost all the unmixed operators may suffer from large theoretical uncertainties caused by both higher order power corrections in OPE and higher order QCD corrections.Thus, these results are not reliable to use in the phenomenological analysis.As for the four ideally mixed operators, we find that the convergence of OPE and the quality of Borel platform are good.But for the QCD perturbative convergence, it is good only for the first two ideally mixed operators J Mixed .In particular, the first two ideally mixed operators lead to results with much weaker µ-dependence at NLO than LO.Phenomenologically, the NLO MS masses of the first two ideally mixed operators are 3.89 +0.18  −0.12 GeV and 4.03 +0.06 −0.07 GeV, which are very close to those of Z c (3900) and Z c (4020), respectively.As for the Z cs system with quark content dccu and J P = 1 + , we study it as a partner of the Z c system in QCD sum rules.Thus, one needs to just replace the down quark d with the strange quark s in both the operators and the quark condensate.The calculations and conclusions are almost the same as those for the Z c system.We find that both the convergence of OPE and the convergence of QCD perturbation corrections are good for the ideally mixed operators J Mixed 1,5 and J Mixed   2,6   in the Z cs system, and they give NLO MS masses 4.02 +0.17 −0.09GeV and 4.21 +0.08 −0.07 GeV in turn, which are very close to Z cs (3985)/Z cs (4000) and Z cs (4220), respectively.Thus, our results support to assign Z cs (3985)/Z cs (4000) and Z cs (4220) as the partners of Z c (3900) and Z c (4020), respectively.
Moreover, we also study the Z b system, however, both the perturbative convergence and the quality of the Borel platforms are bad even for the mixed operators.This case is very similar to those of fully bottom baryon Ω bbb [80] and fully bottom tetraquark bb bb [81].Thus, we can not give credible results for Z b system at present.
Finally, as having been emphasized in fully heavy tetraquark system [81], we would like to emphasize again 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 and Q, which is the most important short-range attractive force to form a heavy quarkonium.(ii) In the dccu tetraquark system discussed in this paper, if one starts from a color-singlet 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, due 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.

A.2 The operator renormalization matrix
For meson-meson type operator bases Eq.(4.1), the corresponding operator renormalization matrix is given by, , (A.17) where δ MS = 1 ϵ + ln(4π) − γ E .For diquark-antidiquark type operator bases Eq.(4.2), the corresponding operator renormalization matrix is given by,                   (a) MS (a) MS               with MS and OS renormalization schemes of Z cs system.with MS and OS renormalization schemes of Z cs system.with MS and OS renormalization schemes of Z cs system.with MS and OS renormalization schemes of Z cs system.

3 3 6 4 1 26 B
Calculation of Wilson coefficients C 1 and C i Current operators of Z c system ( dccu) with (I G )J P = (1 + )1 + 7 4.1 Meson-meson type and diquark-antidiquark type operators 7 Numerical results for meson-meson type and diquark-antidiquark type operators of Z c system 12 5.2 Numerical results for mixed operators of Z c system 14 5.3 Numerical results for mixed operators of Z cs system 18 Details for Various Operators 27 B.1 Numerical results for meson-meson type operators of Z c system 27 B.2 Numerical results for diquark-antidiquark type operators of Z c system 32 B.3 Numerical results for mixed operators of Z c system 37 B.4 Numerical results for mixed operators of Z cs system 40 1 Introduction

Figure 2 .
Figure 2. NLO Feynman diagrams (including counter terms) of C 1 .Z c denotes the interpolating current.

Figure 3 .
Figure 3.The LO and NLO mass spectra of meson-meson type operators of Z c system in the MS and OS schemes.The errors of masses shown in this figure just come from the parametric dependence on s 0 and M 2 B , listed in Tab.9-10 in Appendix B.

Figure 4 .
Figure 4.The LO and NLO mass spectra of diquark-antidiquark type operators of Z c system in the MS and OS schemes.The errors of masses shown in this figure just come from the parametric dependence on s 0 and M 2 B , listed in Tab.11-12 in Appendix B.

Figure 5 .
Figure 5.The LO and NLO mass spectra of mixed operators of Z c system in the MS and OS schemes.The errors of masses shown in this figure just come from the parametric dependence on s 0 , M 2 B , heavy quark mass m Q and renormalization scale µ, which are shown in the Tab.1-4.

3 , 7 of
Z c system in MS and On-Shell schemes.Here the errors for M H are from s 0 , M B , the charm quark mass, and the renormalization scale µ with µ = kM B and k ∈ (0.8, 1.2) (the central values correspond to µ = M B ).

4 , 8 of
Z c system in MS and On-Shell schemes.Here the errors for M H are from s 0 , M B , the charm quark mass, and the renormalization scale µ with µ = kM B and k ∈ (0.8, 1.2) (the central values correspond to µ = M B ).

Figure 7 .
Figure 7.The LO and NLO mass spectra of mixed operators of Z cs system in the MS and OS schemes.The errors of masses shown in this figure just come from the parametric dependence on s 0 , M 2 B , heavy quark mass m Q and renormalization scale µ, which are show in the Tab.5-8.

Figure 8 . 5 , 8 .in
Figure 8.The µ-dependence of the LO and NLO results for mixed operators J Mixed i of Z cs system in the MS scheme.(a) for J Mixed 1,5 , (b) for J Mixed 2,6 , (c) for J Mixed 3,7

4 , 8
operator.Thus, based on our NLO calculations, only the results of operators J Mixed 1

18 )B
Details for Various Operators B.1 Numerical results for meson-meson type operators of Z c system .

Figure 9 .
Figure 9. LO and NLO Result of J M-M 1

Figure 10 .
Figure 10.LO and NLO Result of J M-M 2

Figure 15 .
Figure 15.LO and NLO Result of J M-M 7

Figure 16 .
Figure 16.LO and NLO Result of J M-M .

Figure 17 .
Figure 17.LO and NLO Result of J Di-Di 1

Figure 18 .
Figure 18.LO and NLO Result of J Di-Di 2

Figure 19 .
Figure 19.LO and NLO Result of J Di-Di 3

Figure 20 .Figure 21 .
Figure 20.LO and NLO Result of J Di-Di

Figure 22 .
Figure 22.LO and NLO Result of J Di-Di

Figure 23 .
Figure 23.LO and NLO Result of J Di-Di 7

Figure 24 .
Figure 24.LO and NLO Result of J Di-Di

Table 1 .
The LO and NLO results for the mass of J Mixed 1,5 of Z c system in MS and On-Shell schemes.Here the errors for M H are from s 0 , M B , the charm quark mass, and the renormalization scale µ with µ = kM B and k ∈ (0.8, 1.2) (the central values correspond to µ = M B ).

Table 2 .
The LO and NLO results for the mass of J Mixed 2,6 of Z c system in MS and On-Shell schemes.Here the errors for M H are from s 0 , M B , the charm quark mass, and the renormalization scale µ with µ = kM B and k ∈ (0.8, 1.2) (the central values correspond to µ = M B ).

Table 3 .
The LO and NLO results for the mass of J Mixed

Table 4 .
The LO and NLO results for the mass of J Mixed

Table 5 .
The LO and NLO results for the mass of J Mixed 1,5 of Z cs system in MS and On-Shell schemes.Here the errors for M H are from s 0 , M B , the charm quark mass, and the renormalization scale µ with µ = kM B and k ∈ (0.8, 1.2) (the central values correspond to µ = M B ).

Table 6 .
The LO and NLO results for the mass of J Mixed 2,6 of Z cs system in MS and On-Shell schemes.Here the errors for M H are from s 0 , M B , the charm quark mass, and the renormalization scale µ with µ = kM B and k ∈ (0.8, 1.2) (the central values correspond to µ = M B ).

Table 7 .
The LO and NLO results for the mass of J Mixed

Table 8 .
The LO and NLO results for the mass of J Mixed 4,8 of Z cs system in MS and On-Shell schemes.Here the errors for M H are from s 0 , M B , the charm quark mass, and the renormalization scale µ with µ = kM B and k ∈ (0.8, 1.2) (the central values correspond to µ = M B ).

Table 9 .
LO and NLO results of meson-meson type operators of Z c system with MS renormalization scheme.The errors of masses shown in this table just come from the parametric dependence on s 0 and M 2 B .

Table 10 .
LO and NLO results of meson-meson type operators of Z c system with OS renormalization scheme.The errors of masses shown in this table just come from the parametric dependence on s 0 and M 2 B

Table 11 .
LO and NLO results of diquark-antidiquark type operators of Z c system with MS renormalization scheme.The errors of masses shown in this table just come from the parametric dependence on s 0 and M 2 B .

Table 12 .
LO and NLO results of diquark-antidiquark type operators of Z c system with OS renormalization scheme.The errors of masses shown in this table just come from the parametric dependence on s 0 and M 2 B

Table 13 .
LO and NLO results of mixed operators of Z c system with MS renormalization scheme.The errors of masses shown in this table just come from the parametric dependence on s 0 and M 2 B .

Table 14 .
LO and NLO results of mixed operators of Z c system with OS renormalization scheme.The errors of masses shown in this table just come from the parametric dependence on s 0 and M 2 B .

Table 15 .
LO and NLO results of Z cs system with MS renormalization scheme.The errors of masses shown in this table just come from the parametric dependence on s 0 and M 2 B .

Table 16 .
LO and NLO results of Z cs system with OS renormalization scheme.The errors of masses shown in this table just come from the parametric dependence on s 0 and M 2 B .