NLO Effects for Doubly Heavy Baryon in QCD Sum Rules

With the QCD sum rules approach, we study the newly discovered doubly heavy baryon $\Xi_{cc}^{++}$. We analytically calculate the next-to-leading order (NLO) contribution to the perturbative part of $J^{P} = \frac{1}{2}^{+}$ baryon current with two identical heavy quarks, and then reanalyze the mass of $\Xi_{cc}^{++}$ at the NLO level. We find that the NLO correction significantly improves both scheme dependence and scale dependence, whereas it is hard to control these theoretical uncertainties at leading order. With the NLO contribution, the baryon mass is estimated to be $m_{\Xi_{cc}^{++}} = 3.66_{-0.10}^{+0.08} \text{~GeV}$, which is consistent with the LHCb measurement.

Introduction -Quark model predicts rich structures of hadronic states according to symmetries of quarks. Numerous predicted states have been observed experimentally, which indicates the validity of quark model. Yet, for a class of states, which contain more than one heavy quark, the discovery has not been confirmed for decades. Recently, LHCb collaboration observed a highly significant structure in the Λ + c K − π + π + mass spectrum, which is interpreted as the doubly charmed baryon Ξ ++ cc [1] with mass 3621 ± 0.72 ± 0.27 ± 0.14 MeV. Early experimental studies of Ξ + cc were performed by SELEX [2], Babar [3] and Belle [4] collaborations.
The discovery of Ξ ++ cc demands more rigorous theoretical studies. A number of methods have been used in the literature [5][6][7][8][9][10][11][12][13]. Among them, the QCD sum rules, which are based on the first principle of QCD, are powerful tools to study various properties of hadronic states [14,15]. Many works have devoted to the study of doubly heavy baryons within QCD sum rules [16][17][18][19][20][21][22], and got very interesting results. But in all these works, only leading-order (LO) in α s expansion of perturbative contribution and Wilson coefficients of vacuum condensates are considered. Without higher order contributions, it is hard to control theoretical uncertainties in QCD sum rules, which limits the predictive power. It was in fact known long time ago that next-to-leading order (NLO) correction has sizable contributions to meson and nucleon sum rules [23][24][25]. Therefore, the study of NLO effect for doubly heavy baryons in QCD sum rules is badly needed.
Higher order calculations in QCD sum rules become harder and harder when more particles or more massive particles are involved. For mesons, the state-of-art calculation has developed to O(α 4 s ) in terms of mass expansion [26][27][28][29][30][31]. While for baryons, only the O(α s ) correction (or NLO) is available in the literature for nucleon and singly heavy baryon [24,25,32].
In this Letter, we calculate the NLO correction to perturbative contribution of doubly heavy J P = 1 2 + baryon, and show its important effects in QCD sum rules. With the help of integration-by-parts [33,34] and differential equations [35,36] methods, we get a fully analytical expression. We reproduce the massless result in the literature when we set all quark masses to zero. Based on this calculation, we reanalyze the newly discovered Ξ ++ cc in QCD sum rules.
QCD Sum Rules -The central object in QCD sum rules is the following two-point correlation function [14,37] where Ω denotes the QCD vacuum, and η is the baryon current to be defined later.
On the one hand, one can calculate Π(q 2 ) using operator product expansion, which gives where C 1 is the perturbative contribution and C i is the Wilson coefficient of a gauge invariant Lorentz scalar operator O i . Both C 1 and C i are perturbatively calculable. O i is a shorthand for the vacuum condensates Ω|O i |Ω , which is a nonperturbative but universal quantity. It means that the value of O i determined from other processes should be the same as its value in the process considered in this Letter. On the other hand, Π(q 2 ) satisfies the dispersion relation where ρ is the spectrum density. Based on the optical theorem, one assume spectrum density ρ(q 2 ) to be [37] ρ(q 2 ) = λ 2 where s th is the threshold of continuum spectrum, λ H is defined by λ H u(p, s) = 0|η(0)|H(p, s) , where u(p, s) is the Dirac spinor of hadron. Define and employ quark-hadron duality and Borel transformation, we obtain a sum rule corresponding to Π 1 (q 2 ) [37] where s 0 is the threshold parameter, and m B is the Borel parameter, which are introduced in quark-hadron duality and Borel transformation respectively. One can also obtain a similar sum rule corresponds to Π 2 (q 2 ), but we will not discuss it in this Letter.
To obtain the baryon mass, we differentiate both sides of Eq. (7) with respect to −m −2 B and solve for m 2 H , which results in In this Letter, we keep vacuum condensates up to dimension 4, and evaluate ρ 1,qq up to O(m q ). Contributions from even higher dimension operators are highly power suppressed and thus can be neglected up to a desired precision. Baryon Currents -The most general current of baryon containing two identical heavy quarks is where Q is the heavy quark with mass m Q , while q is the light quark with mass m q . ǫ abc is the antisymmetric matrix in color space, C is the charge conjugation matrix, and Γ 1 and Γ 2 are Dirac matrices with possible Lorentz indices suppressed. Spinor indices are contracted within the bracket, and therefore transposing the bracket part should keep the current intact. Note that C T = −C, one can see that Γ 1 can only be γ µ or σ µν [37]. For a J P = 1 2 + baryon, there are only two possible currents where η 1 corresponds to Ioffe current [37] if we take Q as u quark and q as d quark. It is well-known that η 1 and η 2 are renormcovariant [38], Thus it is advantageous to work with these currents when calculating the NLO correction. There exist other choices of current [16,39,40], which can be expressed by η 1 and η 2 with the help of Fierz identity, where b is a complex mixing parameter. NLO Correction to C 1 -It is known that C 1 and C i can be calculated perturbatively, and results at LO are available in Refs. [16,41]. Among them, the most important one is C 1 , because all other coefficients will be multiplied by higher dimensional operators which are power suppressed. Thus the main theoretical uncertainty is due to NLO correction to C 1 .
In order to perform NLO calculation for C 1 , we use FeynArts [42,43] to generate all Feynman diagrams (see Fig. (1)), and FeynCalc [44,45] to manipulate resulting amplitude. After these steps, we are left with some threeloop-like scalar integrals. These integrals can be further simplified by integration-by-parts (IBP) method [33,34]. FIRE [46] and LiteRed [47] are used to reduce the full amplitude to linear combination of a complete set of 29 master integrals (see Fig. (2)), where ε is defined by dimension D = 4 − 2ε, v = 1 − 4m 2 Q q 2 , and all coefficients c k are purely imaginary. Note that here I k is defined to be dimensionless. corresponding cut diagrams of I k . But evaluating fourbody phase space in the presence of two massive particles is still a formidable task.
To proceed, we employ differential equation method [35,36] by first differentiating I k with respect to v, then reducing the resulting integrals using IBP, and obtaining a system of differential equations, where I represents the vector of master integrals I k , and A is a 29 × 29 matrix. To solve this differential equation, we implement algorithm proposed in [48] to transform the equation into so-called ε-form [35], where v i ∈ 0, ±1, ± √ 3i , B i are constant matrices, and I ′ is related to I with an invertible linear transformation. The virtue of this ε-form is that the right hand side of Eq. (17) is proportional to ε, which can be easily solved iteratively in terms of Goncharov polylogarithms [49]. The boundary value of I(ε, v) at v = 1, i.e. m Q = 0, is nothing but massless four-body phase space integral, which is very easy to work out. By evaluating the boundary value I(ε, 1), and solving the equation iteratively, we finally obtain I k and finish our calculation.
We find that the Coulomb divergence, which appears as v → 0, does not present at this order. Then by combining all terms together, infrared divergences are canceled out, so we only need to deal with ultraviolet divergences. After performing wavefunction and mass renormalization of quarks (m Q is renormalized either in MS scheme or onshell scheme), the remained ultraviolet divergences can be removed by operator renomalization of η 1 and η 2 . We renormalize them in MS scheme, anomalous dimensions of which are which confirm the results in Refs. [25,50]. We then get the finite result at NLO. Our NLO result confirms the massless result [24,25] in the limit of m Q → 0. Our analytical result is provided as ancillary file of the arXiv preprint.
Phenomenology -In our analysis, we use with θ a complex mixing parameter. We choose following parameters [16,[51][52][53]: and α s (m Z = 91.1876 GeV) = 0.1181. According to Eq. (8), the evolution of the current η is irrelevant to the estimation of hadron mass, thus we do not include it. We use two-loop running of coupling constant α s and heavy quark mass m Q . We evolve vacuum condensates according to their one-loop anomalous dimensions: γ qq = −γ mq and γ g 2 s GG = 0 [54]. By default choice in the following, we choose central values for all parameters, set renormalization scale µ = m B [14,55], and choose MS scheme for heavy quark mass renormalization.
In Eq. (8), the baryon mass m H depends on two parameters: m B and s 0 . In order to obtain a reliable result, we should keep m B inside the so-called Borel window to ensure the validity of OPE, and the choice of s 0 should ensure the ground-state pole contribution domination. Since m H is a property of hadron, it does not depend on m B and s 0 , thus within the valid parameter space (we shall call it "window" hereafter), we should find a region in which m H depends weakly on m B and s 0 . m H in this region is considered to be the estimated hadron mass in QCD sum rules.
We define relative contributions of condensates and continuum spectrum as and impose the following constraints on our sum rule We find that with mixing parameter θ = 0.018i, we can obtain a very stable plateau of m B and s 0 , as shown in Fig. (3). Note, however, that QCD sum rules alone cannot tell which mixing current is the physical one. For example, there is a family of mixing parameters that can yield similar good plateau of m B and s 0 . We thus also provide another set of results by choosing θ = − i 3 , which corresponds to the mixing used in [16]. The relative importance of each term in OPE is shown in Fig. (4) values shown in Tab. (I). We find that the NLO correction has important contribution. In m MS Q scheme, the ratio of NLO correction to LO is about 29% (19%) if θ = 0.018i (θ = − i 3 ). While in m on-shell Q scheme, these ratios reaches 233% (146%), signaling the bad convergence of perturbative expansion, which is the reason why we choose MS scheme by default. Nevertheless, with NLO correction, the difference of predicted m Ξ ++ cc between MS scheme and on-shell scheme for m Q is substantially reduced. As shown in Tab. (I), the mass differences obtained from LO and LO + NLO results are 0.27 GeV and 0.01 GeV, respectively. Thus NLO correction largely reduces the scheme dependence.
To study the renormalization scale µ dependence, we fix all other parameters by their default choices (or central values) and freely vary µ. The variation of m Ξ ++ cc with respect to µ is shown in Fig. (5). We find the scale dependence is much weaker when NLO correction is included. More precisely, the error of m Ξ ++ cc . We find that our NLO result is consistent with the LHCb measurement. As a comparison, we also list the results with m on-shell Q renormalization scheme or with θ = − i 3 . We find that all plots above are almost unchanged when changing m q from m u to m d , thus our prediction of mass of Ξ + cc (ccd) is almost the same as that of Ξ ++ cc (ccu).
Summary -The NLO calculation for hadrons with massive quarks in QCD sum rules is important but hard to do. With the help of recent development of multi-loop calculation technique, we are able to analytically calculate NLO perturbative correction to the imaginary part of the two-point correlation function of J P = 1 2 + baryon current with two identical heavy quarks. We apply our result to the QCD sum rules analysis of newly discovered baryon Ξ ++ cc by LHCb [1]. The QCD sum rules estimation of m Ξ ++ cc is 3.67 +0.09 −0.10 GeV, which is consistent with the LHCb measurement within uncertainties. By comparing LO with LO + NLO results, we find the NLO perturbative correction substantially reduces m Q renormalization scheme dependence and renormalization scale µ dependence, thus makes the theoretical uncertainties under better control.
We thank H. X. Chen and S. L. Zhu for many useful communications and discussions. The work is supported in part by the National Natural Science Foundation of