Color-octet nonrelativistic QCD matrix elements for heavy quarkonium decays in the refined Gribov-Zwanziger theory

We determine color-octet nonrelativistic QCD matrix elements for quarkonium decays from moments of the two-point correlation function of the QCD field-strength tensor computed in the refined Gribov-Zwanziger theory. We find that a tree-level calculation in the refined Gribov-Zwanziger theory can give a suitable description of the QCD field-strength correlation function at both short and long distances, which leads to moments that are infrared finite and can be properly renormalized. By using the color-octet matrix elements we obtain, we quantitatively improve the nonrelativistic effective field theory description of quarkonium decay rates, especially for the χ QJ and η Q states, where Q = c or b .


I. INTRODUCTION
The two-point correlation function of the QCD fieldstrength tensor has been considered an important quantity in phenomenological studies of the strong interaction [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15].In particular, it has been known that they can be used to compute nonperturbative matrix elements that arise in decays of heavy quarkonium states [16][17][18].The contributions to the decay rates of heavy quarkonia from probabilities for a heavy quark Q and antiquark Q pair to be in a color-octet state are encoded in the nonrelativistic QCD (NRQCD) color-octet matrix elements, which can be expressed as products of quarkonium wave functions at the origin and moments of field-strength correlators as have been shown in the potential NRQCD (pNRQCD) effective field theory formalism [17][18][19][20][21].The color-octet contributions can have significant effects on decay rates of heavy quarkonia.Most notably, in inclusive decays of P -wave heavy quarkonia the color-octet contribution appears at leading order in the nonrelativistic expansion in powers of v, the typical velocity of the heavy quark inside the quarkonium [22].Color-octet contributions also appear in two-photon decays of P -wave heavy quarkonia as corrections of order v 2 [23].Even in the case of S-wave quarkonium decays, inclusion of coloroctet contributions are necessary for improving the theory description of decay rates and two-photon branching fractions of η c [24,25].It is also known that color-octet contributions are enhanced by an inverse power of α s in the case of J/ψ and Υ decays [26,27].As the quarkonium wave functions at the origin can usually be determined from potential models and electromagnetic decay rates of heavy quarkonia, and even be computed accurately from first principles [28,29], quarkonium decays can provide useful probes of the moments of QCD field-strength correlators.
The moments of field-strength correlators are sensitive to both the short-distance and long-distance behaviors of the correlators.The correlators as functions of the separation of the field strengths have been studied in perturbative QCD [30], in operator-product expansion [31,32], and on the lattice [33][34][35].The calculations in both perturbative QCD and operator-product expansion approaches only reproduce the short-distance behaviors of the correlators, so they cannot be used to compute the moments; these calculations give results for the correlators that are given by inverse powers of the separation of the two field strengths, whose moments are scaleless divergent and vanish in dimensional regularization.On the other hand, lattice QCD calculations have been shown to reproduce both the nonperturbative long-distance behavior and the power-law behavior at short distances.However, the short-distance behaviors of the lattice results are not well understood in terms of perturbative QCD calculations, and this makes it difficult to renormalize the ultraviolet (UV) divergences in the moments.Renormalization of the moments is important because it is directly related to renormalization of the UV divergences in the color-octet matrix elements.Hence, currently available results for the field-strength correlators do not immediately lead to quantitative results for their moments.It would be desirable to have calculations for the correlators that are compatible with perturbative QCD at short distances, and also describe the nonperturbative longdistance behaviors at least approximately.
It has been known that non-Abelian gauge theories quantized in the standard Fadeev-Popov method contain Gribov copies [36]; that is, a gauge-fixing condition such as ∂ • A = 0, where A is a gauge field, does not uniquely fix the gauge and there are distinct gauge-field configurations that satisfy the gauge-fixing condition that are related by large gauge transformations.In order to remove this ambiguity, the functional integral can be restricted to a region free of Gribov copies called the fundamental modular region.A semiclassical calculation of the gluon propagator in the Landau gauge leads to a result that is modified in the infrared region by a dimensionful quantity called the Gribov parameter in a way that the poles of the propagator are shifted to nonphysical locations, while the large-momentum behavior of the propagator coincides with perturbative QCD.The Gribov parameter appears in the modified propagator in the form of a complex-valued gluon mass, so that infrared divergences are regulated while single-gluon states do not appear in the physical spectrum.Later it was shown that by introducing auxiliary fields the restriction to the fundamental modular region can be incorporated into a local and renormalizable action, hereafter referred to as the Gribov-Zwanziger (GZ) theory [37], whose tree-level gluon propagator reproduces the semiclassical result in Ref. [36].We refer readers to Refs.[36][37][38][39][40] and the review in Ref. [41] for details of the Gribov ambiguity and the GZ theory.Because the tree-level propagator follows from a semiclassical calculation, one would hope that it would at least approximately describe the nonperturbative behavior of the gluon propagator.Unfortunately, lattice QCD calculations have shown that this is not the case; while the tree-level gluon propagator at zero momentum vanishes in the GZ theory, lattice data shows that it is nonzero [42,43].This discrepancy suggests that more nonperturbative effects will need to be included in order to describe the nonperturbative behavior of non-Abelian gauge theories.
The refined Gribov-Zwanziger (RGZ) theory has been obtained by adding effects of dimension-two condensates associated with the gauge field and the auxiliary fields to the GZ action [44][45][46][47][48][49].The dimension-two condensate of the gauge field has been considered an important quantity in the study of non-Abelian gauge theories [50][51][52]; a gauge-invariant definition of the dimension-two condensate is given by its minimum, which occurs in the Landau gauge.Analyses in the effective potential formalism suggest that a nonvanishing value of the dimension-two condensate is favored [49,52,53]; this is supported by lattice QCD calculations of the gluon propagator [54], the quark propagator [55], and also by a study based on resummation of Feynman diagrams [56].As such a condensate leads to a dynamically generated gluon mass, it would modify the infrared behavior of the gluon propagator.Remarkably, tree-level calculation in the RGZ theory leads to a good description of the lattice measurement of the Landau-gauge gluon propagator [48,57].This suggests that a perturbative calculation in the RGZ theory may be able to account for a good part of the nonperturbative dynamics of gluons in the SU (3) gauge theory.It would therefore be interesting to compute other quantities involving gauge fields in the RGZ theory in perturbation theory, such as the QCD field-strength correlators and their moments.
In this work, we compute the two-point correlation functions of the QCD field-strength tensor and their moments in the RGZ theory at tree level.At tree level, perturbative calculations in the RGZ theory simply amounts to using the modified gluon propagator in usual perturbative QCD calculations.Due to the dimensionful parameters associated with the Gribov parameter and the condensates, perturbative calculations in the RGZ theory are infrared finite.When the field-strength correlators are computed perturbatively in the RGZ theory, their long-distance behaviors are modified from the power-law behaviors in the usual perturbative QCD calculations so that their moments are infrared finite, while at short distances they reproduce the leading UV behavior in perturbative QCD.We thus obtain finite values of the moments by subtracting the UV divergences through renormalization, and obtain renormalized color-octet matrix elements.From this we quantitatively determine the coloroctet contributions to quarkonium decay rates, which can be compared directly with experiment.
This paper is organized as follows.In Sec.II, we compute the two-point correlation function of the QCD fieldstrength tensor in the RGZ theory, and compare the results with lattice QCD.In Sec.III, we compute the moments of the field-strength correlation functions that appear in decay rates of heavy quarkonia.We use the results for the moments to obtain the color-octet matrix elements for quarkonium decays and compute the coloroctet contributions to decay rates of heavy quarkonia in Sec.IV.We conclude in Sec.V.

II. QCD FIELD-STRENGTH CORRELATORS
We define the two-point correlation function of the QCD field-strength tensor in Euclidean space as where ν is the QCD field-strength tensor, A µ is the gauge field, g is the strong coupling, T F = 1/2, and Φ ab (x, 0) is a straight Wilson line defined by Φ(x, 0) = P exp ig where P is the path ordering and A adj is the gauge field in the adjoint representation.The definition in Eq. ( 1) is consistent with Refs.[33][34][35], while it contains an additional factor of g 2 T F compared to what was used in the perturbative calculation in Ref. [30].The general form of the correlator can be written as where D(x 2 ) and D 1 (x 2 ) are invariant functions of x 2 .The form (3) follows from the fact that D(x 2 ) vanishes in the Abelian gauge theory.Because of this, D(x 2 ) vanishes at tree level in the non-Abelian gauge theory.At tree level we can compute D 1 (x 2 ) at order α s while D(x 2 ) = 0 + O(α 2 s ), where α s = g 2 /(4π).Useful representations of the invariant functions can be obtained by defining the two linear combinations where d = 4 − 2ϵ is the number of spacetime dimensions, and we work in a frame where x µ is aligned to the time (µ = 0) direction.Note that in this frame, D E (x 2 ) involves only the chromoelectric fields and D B (x 2 ) involves only the chromomagnetic fields.We work in arbitrary spacetime dimensions to make possible the use of dimensional regularization in the computation of the moments in the following section.The invariant functions can be conveniently computed in perturbation theory by using the first equalities.At tree level, the correlation function can be computed by differentiating the tree-level gluon propagator.The Landau-gauge gluon propagator in the RGZ theory is given by [48] where m and M are dimension-one parameters associated with the condensates of the gauge and auxiliary fields in the RGZ action, respectively, and λ 4 = 2g 2 N c γ 4 +m 2 M 2 with γ the Gribov parameter.Note that the tree-level propagator in the GZ theory is obtained by setting m and M to zero, while keeping γ nonzero [37].The values of m, M , and λ that give a satisfactory description of the gluon propagator have been found in Ref. [48] based on lattice QCD calculations of the gluon propagator with coupling β = 6/g 2 set to 5.7, 6.0, and 6.2.They read M 2 = 2.15 ± 0.13 GeV 2 , m 2 = −1.81± 0.14 GeV 2 , and λ 4 = 0.26 GeV 4 .By using Eq. ( 6) we obtain where Note that this expression is valid in arbitrary dimensions.We can integrate over k 0 by using contour integration.
The poles in k 0 are located at ±i Note that Q 2 is purely imaginary for the values of the parameters determined in Ref. [48].Because √ x 2 > 0, we must close the contour on the upper half plane.We obtain where Note that the quantity in the square brackets can be written as twice the real part of the first term.To obtain D B (x 2 ) at a given value of x 2 , the remaining integral in d − 1 = 3 dimensions can be computed numerically.Note that the small-x 2 behavior can be computed by setting the dimensionful parameters m, M , and λ to zero, which leads to the perturbative QCD result ), in agreement with Ref. [30].We show the numerical result for D B (x 2 )| RGZ in Fig. 1 compared to the lattice measurement from Ref. [35].For the numerical calculation we set the strong coupling to be β = 6.0, which is close to the average value used in the determination of the parameters of the RGZ theory in Ref. [48], and is also same as the value used for the main results of the lattice study in Ref. [35].Note that our definition of D B (x 2 ) corresponds to D ⊥ (x 2 ) in Ref. [35].We display the lattice measurement by using the functional form given by D lat ⊥ (x 2 ) = Ae −|x|/λa with A = 0.94 +0.32  −0.16 GeV 4 and λ a = 0.120 +0.009 −0.012 fm determined in Ref. [35].Because this result was obtained from lattice data for √ x 2 > 0.3 fm, and lattice data are only shown for √ x 2 < 0.6 fm in Ref. [35], we only display the result from D lat ⊥ (x 2 ) for √ x 2 between 0.3 and 0.6 fm.We can see that the RGZ result agrees with the lattice QCD result within uncertainties.Note that similarly to the lattice study, D B (x 2 )| RGZ exhibits an exponentially decaying behavior at large x 2 with a finite correlation length, as can be seen from Fig. 1 for √ x 2 larger than about 2 GeV −1 .A fit of the same functional form Ce −|x|/λc for |x| larger than 2 GeV −1 (≈ 0.4 fm) to the RGZ result leads to C = 0.58 ± 0.03 GeV 4 and λ c = 0.77 ± 0.01 GeV −1 .This correlation length reads in distance units λ c ≈ 0.15 fm.The lattice QCD study in Ref. [35] also found that at small √ x 2 the function D B (x 2 ) shows a power-law behavior given by 0.4/x 4 , which agrees very well with the tree-level perturbative QCD result if we set the coupling using β = 6.0 as was done in the lattice study in Ref. [35].For comparison, in Fig. 1 we also show the result from the GZ theory that we obtain by setting m = M = 0 while keeping the parameter γ the same as the one obtained in the RGZ theory.We find that the GZ theory results in the longdistance behavior of D B (x 2 ) that falls off too quickly compared to the RGZ theory and is incompatible with lattice QCD data.In fact, we find that the GZ theory result for D B (x 2 ) always falls off faster than the perturbative QCD result for any nonzero value of the Gribov parameter, so that it would not be possible to obtain a result that is compatible with lattice QCD measurements.
Now we look into D E (x 2 ).Because D(x 2 ) vanishes at tree level, we could obtain x 2 ∂D 1 (x 2 )/∂x 2 at tree level by differentiating and multiplying by x 2 our result for D B (x 2 ).However, for computing moments of the correlators in the following sections it is necessary to obtain a dimensionally regulated momentum-integral representation for D E (x 2 ).By using the definition for D E (x 2 ) and the tree-level gluon propagator in the RGZ theory we obtain We can again integrate over k 0 by using contour integration, closing the contour on the upper half plane.We obtain Similarly to the case of D B (x 2 ), we can obtain D E (x 2 ) at a given value of x 2 by evaluating the remaining integral in d − 1 = 3 dimensions numerically.
In order to compare with lattice data, we consider the combination D B (x 2 ) − D E (x 2 ) = −x 2 ∂D 1 (x 2 )/∂x 2 which corresponds to −D * (x 2 ) in the lattice QCD study in Ref. [35].We show the result in Fig. 1 against the lattice measurement from Ref. [35].We display the lattice measurement by using the functional form given by −D lat * (x 2 ) = Be −|x|/λ b with B = 0.47 +0.20 −0.06 GeV 4 and λ b = 0.189 +0.013 −0.029 fm in Ref. [35].Because the fit in Ref. [35] for D lat * (x 2 ) was restricted to √ x 2 > 0.35 fm, and lattice data are only shown for √ x 2 < 0.6 fm, we only display the lattice result for √ x 2 between 0.35 and 0.6 fm.Similarly to the case of D B (x 2 ), the result for D B (x 2 ) − D E (x 2 ) in the RGZ theory agrees with the lattice result in Ref. [35] within uncertainties.The authors of Ref. [35] did not present a functional form of the short-distance behavior of D B (x 2 ) − D E (x 2 ).Nevertheless, the lattice QCD data available in Ref. [35] at distances shorter than about 0.35 fm seem to agree very well with what is expected from perturbative QCD at tree level given by D B (x 2 ) − D E (x 2 )| pQCD ≈ 0.8/x 4 for β = 6.0.As we have done for D B (x 2 ), we also show the results from the GZ theory in Fig. 1, which we obtain by setting m = M = 0 and keeping the Gribov parameter γ the same as the one obtained in the RGZ theory.We see that also in the case of D B (x 2 ) − D E (x 2 ) the result from the GZ theory has a long-distance behavior that falls off too quickly and is incompatible with lattice QCD results.
As we have discussed earlier, in perturbation theory D(x 2 ) appears from next-to-leading order accuracy, and so, at tree level only the D 1 (x 2 ) and its derivative contribute to D B (x 2 ) and D E (x 2 ).Based on the fact that the tree-level results for the invariant functions in the RGZ theory are compatible with the lattice QCD study in Ref. [35], we may assume that D(x 2 ) is indeed small and is negligible compared to the uncertainties in the lattice results.The agreement in the short-distance behaviors of the invariant functions between the lattice QCD study and the RGZ results supports the suppression of D(x 2 ) even at short distances.We may estimate the effects from a nonzero D(x 2 ) by using the lattice QCD results in Ref. [35] in our numerical analysis.
We note that there are other lattice QCD studies of the two-point field-strength correlation function based on the cooling method in Refs.[33,34].These lattice studies result in a much larger D(x 2 ) compared to D 1 (x 2 ) at both short and long distances, which contradicts the naïve expectation from perturbation theory that D(x 2 ) would be suppressed compared to D 1 (x 2 ), and is also in conflict with the smallness of D(x 2 ) that can be inferred from the lattice study in Ref. [35].As was discussed in Ref. [30], the two-point field-strength correlation function involves UV divergences at loop level, and after renormalization the functions D(x 2 ) and D 1 (x 2 ) mix under change of the renormalization scheme and scale; since the cooling method used in Refs.[33,34] removes short-range fluctuations, it would not be possible to make direct comparisons until the results are converted to a common scheme.In our case, because we work at tree level and the lattice results in Ref. [35] are renormalized nonperturbatively, it is not possible to quantify the effect of scheme dependence; nonetheless, the good agreement between the tree-level RGZ results and the lattice results in Ref. [35] at both long and short distances may indicate that the effect of the scheme change is small and may even be negligible between the two results.
So far the calculation in this section has been done in the Landau gauge.Because we work with the gaugeinvariant definition of the field-strength correlator (1), we expect that the results for D E and D B we obtain to be valid in any gauge.We may check this by, for example, computing the correlators from the gluon propagator in a covariant gauge.We first need to replace the dimension-two condensate term A 2 in the action with a gauge-invariant expression [58,59] where the corrections from order g involve products of three or more gauge fields and can be computed systematically as formal power series in A. At the current level of accuracy, we only need to keep the lowest-order nonvanishing contribution.We can then obtain the gluon propagator in the RGZ theory in a generic covariant gauge; the result in momentum space contains an extra term to the Landau-gauge expression in Eq. ( 6), with α the gauge parameter.It is easy to check explicitly that this additional term makes vanishing contributions in D E and D B .This may be understood as a consequence of the Ward identity, arising from the gauge invariance of the field-strength correlator.
In this section we have established the two-point fieldstrength correlation functions at tree level in the RGZ theory as dimensionally regulated momentum integrals in Eqs. ( 8) and (10).We will use these results in the following sections to compute moments of the correlation functions in dimensional regularization.

III. MOMENTS OF FIELD-STRENGTH CORRELATORS
In this section we compute the moments of the twopoint field-strength correlator.Rather than considering all possible moments of the correlator, we focus on the ones that appear in heavy quarkonium decay rates.Roughly speaking, the moments of the correlator represent the squared amplitude for a heavy quark-antiquark pair in a color-octet state to evolve into a color-singlet state through insertions of chromoelectric or chromomagnetic fields on the heavy quark or antiquark lines.We refer the readers to Refs.[17,18,60] for details of the pNRQCD formalism that were used to obtain expressions for color-octet NRQCD matrix elements in terms of moments of field-strength correlators.
The following dimensionless moment of D E (x 2 ) defined by appears in inclusive decay rates of P -wave heavy quarkonia at leading order in v.The factor 1/N c appears from the projection of the Q Q color, and the sign arises from the fact that E 3 was defined in Minkowski space in Ref. [18] while the right-hand side is computed in Euclidean space.Because the D E (τ ) scales like 1/τ 4 at small τ , the integral over τ contains a logarithmic UV divergence.After renormalization, E 3 acquires a scale dependence given by where Λ is the renormalization scale for E 3 .Throughout this paper, we use the superscript (Λ) to denote that the quantity is renormalized in the MS scheme and depends on the MS scale Λ.It has been known that this scale dependence directly reproduces the renormalization-scale dependence of the color-octet matrix element that appears in the P -wave quarkonium decay rates [16,17].
Note that E 3 also appears in the order-v 2 correction to the leading-order color-singlet matrix elements for S-wave quarkonium decays [18].
Similarly, the following dimension-two moments of D E (x 2 ) and D B (x 2 ) defined by appear in inclusive decay rates of S-wave heavy quarkonia at relative order v 4 and v 3 , respectively [18].Although the contributions from these moments to the decay rates are suppressed by powers of v, their effects can be enhanced by large short-distance coefficients.In the case of spin-1 S-wave heavy quarkonium decays, the short-distance coefficients associated with color-octet contributions are enhanced by 1/α s compared to the one at leading order in v, which can make the coloroctet contributions numerically significant [26,27].In the decays of spin-0 S-wave states, such as the η c and η b , the lack of knowledge of the color-octet matrix element arising from B 1 is a significant source of uncertainties [24,25].The quantity E 1 also appears in the orderv 2 correction to two-photon decay rates of P -wave heavy quarkonia [18,60].Hence, accurate determinations of E 1 and B 1 are phenomenologically important.Note that both E 1 and B 1 contain power UV divergences at small τ , which must be subtracted in dimensional regularization consistently with the calculation of short-distance coefficients in perturbation theory.
We also compute the dimension-one moment iE 2 defined by where the phase is chosen to recover the Minkowski space definition in Ref. [18].Note that E 2 is purely imaginary, so that iE 2 is real.Although this moment does not appear directly in color-octet matrix elements, it appears in the corrections to the quarkonium wave functions at the origin associated with the velocity-dependent potential at zero separation [28,29], and can be useful in heavy quarkonium decay phenomenology.We note that the overall phases of the right-hand sides of the definitions of the moments can be checked against the Minkowski space definitions in Ref. [18] by computing them in perturbative QCD and comparing the integrands of the momentum integral.
We now proceed to compute the moments in the RGZ theory.

A. Dimensionless moment E3
We first compute the dimensionless moment E 3 .By using Eq. ( 10) and the identity ∞ 0 dτ τ 3 exp(−Aτ ) = 6/A 4 , which is valid when ReA > 0, we obtain the following ex-pression which is valid in arbitrary dimensions.The d − 1dimensional integral over k is straightforwardly evaluated by using where µ is a scale associated with dimensional regularization.We obtain where in the last equality we expanded in powers of ϵ = (4 − d)/2 and set µ 2 = Λ 2 e γE /(4π), so that Λ is a MS scale.Here, γ E is the Euler-Mascheroni constant.Note that the result in the last equality is invariant under simultaneous rescalings of the denominator factors M 2 in the arguments of the logarithms.The pole in ϵ is exactly what we expect from the order-α s scale dependence in Eq. ( 13).After renormalization we then have in the MS scheme with Λ the MS renormalization scale.This is our result for E 3 in the RGZ theory at tree level.Note that the explicit dependence on log Λ satisfies the evolution equation in Eq. ( 13).
B. Dimension-two moments E1 and B1 Now we consider the dimension-two moments E 1 and B 1 .We first compute E 1 .From Eq. ( 10) and the identity ∞ 0 dτ τ exp(−Aτ ) = 1/A 2 we obtain which is valid for arbitrary d.The (d − 1)-dimensional integral gives where the terms in the first square brackets come from the part of the integrand in Eq. ( 20) that is proportional to k 2 , while the terms in the second square brackets come from the part that is proportional to d − 1.Note the appearance of the pole at d = 2 from Γ(ϵ − 1), indicating that the individual integrals are quadratically power divergent.However, the contributions from the two parts of the integrand cancel each other, so that E 1 | RGZ vanishes through order g 2 : which is valid for arbitrary d.
We can understand this result for E 1 in terms of the functions D(x 2 ) and D 1 (x 2 ) by rewriting the definition as where we changed the integration variable to τ 2 , and used integration by parts to obtain the last equality.
Note that this expression is valid to all orders.While τ 2 D 1 (τ 2 )| τ 2 =∞ = 0, the vanishing of the last term in the last equality at τ 2 = 0 follows from the fact that the the UV divergences of E 1 must be regularized dimensionally.For example, in dimensional regularization the short-distance behavior of D 1 (τ 2 ) is proportional to 1/τ d , which must be computed with d < 2 in order to regulate power UV divergences; hence we obtain τ 2 D 1 (τ 2 )| τ 2 =0 = 0. Therefore, in dimensional regularization E 1 comes solely from D(τ 2 ), and because this appears only from order α 2 s , E 1 vanishes at tree level.Note that this result is valid not just in the RGZ theory, but generally holds in a non-Abelian gauge theory.Now we compute B 1 .We use Eq. ( 8) and the identity which is valid in arbitrary dimensions.A straightforward evaluation of the (d−1)-dimensional integral over k yields Note the appearance of a power divergence at d = 2.
The result also contains a logarithmic UV divergence; if we expand in powers of ϵ, we obtain Note that the UV pole is proportional to m 2 .Since m 2 is the parameter associated with the dimension-two condensate of the gauge field, this UV pole is of nonperturbative origin and does not have a counterpart in the NRQCD factorization formalism.Hence, it is not possible to properly renormalize this UV divergence within the NRQCD formalism.However, this logarithmic UV divergence can be related to the dimension-two condensate at tree level, which contains the same divergence.In the Landau gauge the dimension-two condensate ⟨(gA) 2 ⟩ ≡ ⟨Ω|g 2 A a µ A a µ |Ω⟩ can be computed at tree level as The first equality follows from the Landau-gauge gluon propagator at zero separation, and the second equality is obtained by integrating over k 0 using contour integration.Note that, in a general gauge, A a µ A a µ must be replaced by a gauge-invariant expression given in Eq. (11).In this case, any additional term in the gluon propagator involving k µ or k ν makes a vanishing contribution to Eq. ( 27), so that the result is gauge invariant.We compute the remaining integral to obtain By comparing this result with Eq. ( 25) we find that Note that this relation holds for arbitrary d.The relation in Eq. ( 29) arises in a similar fashion as the vanishing of E 1 at tree level; note that the integrand of B 1 in Eq. ( 24) is identical to the term proportional to k 2 in the integrand of E 1 in Eq. ( 20), and the integrand of ⟨(gA) 2 ⟩ in Eq. ( 27) is just the term proportional to d − 1 in the integrand of E 1 in Eq. ( 20).As we have seen from the evaluation of E 1 , the two contributions are equal, leading to the relation in Eq. ( 29).Similarly to the case of E 1 , this relation may be modified at higher orders in α s , notably from a nonzero D(τ 2 ).
In fact, the relation in Eq. ( 29) can be understood in terms of gauge fields by rewriting the condensate in a manifestly gauge-invariant form in terms of the fieldstrength tensor, given by [58] where D is the covariant derivative.The corrections from order g involve products of three or more fieldstrength tensors.Then, by using , which follows from rotational invariance and taking a d-dimensional angular average, and by exponentiating (D 2 0 ) −1 in momentum space at order g 0 , we obtain This, together with the vanishing of E 1 at order g 2 , leads to the the relation in Eq. ( 29).
As can be obtained from Eq. ( 28), or from the relation in Eq. ( 29) and the result for B 1 in Eq. ( 26), the UV pole in ⟨(gA) 2 ⟩ that we obtain reproduces the lowestorder result in Eq. ( 22) of Ref. [52] for the counterterm of the generating functional proportional to m 2 .While perturbative calculations of B 1 and ⟨(gA) 2 ⟩ both yield logarithmic UV divergences, a finite result for the ⟨(gA) 2 ⟩ has been obtained in the effective potential approach [52]: This result is obtained by minimizing the effective potential in the presence of the dimension-two condensate.Since this finite result has been obtained after subtracting the 1/ϵ pole, the ϵ-dependent coefficient in Eq. ( 29) produces an additional finite piece.That is, From this and Eq. ( 32) we obtain Note that since m 2 is negative, the first term is positive while the second term is negative.This is our finite result for B 1 in the RGZ theory.

C. Dimension-one moment iE2
From Eq. ( 10) and the identity The integral can be evaluated straightforwardly: We can see from the first equality that the integral contains a linear power divergence from the pole at ϵ = 1/2.In the second equality, we expanded in powers of ϵ, which subtracts this power divergence according to dimensional regularization and we obtain a finite result.

D. Numerical results
Now we are ready to present numerical results for the moments.We first compute E 3 in the MS scheme.Because we work at tree level, there is some ambiguity in the choice of the strong coupling.Since the parameters of the RGZ theory are obtained from the lattice data with the lattice coupling around β = 6.0, it seems appropriate to compute the coupling from the relation β = 6/g 2 when computing low-energy operator matrix elements.We take the numerical values of the parameters M , m, and λ from Ref. [48] as we have done in the previous section.By using Eq. ( 19) we obtain at scale Λ = 1 GeV the value for E 3 given by where the uncertainties come from the parameters in the RGZ theory.Note that our result for E 3 at tree level misses the contribution from D(τ 2 ), which only occurs from order α 2 s .We can estimate the effect of the inclusion of D(τ 2 ) at long distances by using the lattice measurement of D B (τ 2 ), which includes D(τ 2 ), and comparing it with the result from the RGZ theory containing only D 1 (τ 2 ).We neglect the short-distance contribution from D(τ 2 ) because it is of higher orders in the strong coupling, and the short-distance behavior approximately cancels in the difference between the lattice and the RGZ results for D B (τ 2 ).Hence, we use the long-distance functional forms we obtained in Sec.II to estimate the longdistance contribution from D(τ 2 ).The contribution to where Ae −τ /λa comes from the lattice QCD result for D B (τ 2 ) in Ref. [35] and Ce −τ /λc is the long-distance behavior of the RGZ result for D B (τ 2 ) determined in Sec.II.Alternatively, we could compute the same quantity by obtaining D 1 (τ 2 ) from lattice measurements of ), whose long-distance behavior is given in Ref. [35] by −Be −τ /λ b .From this we obtain an expression for D 1 (τ 2 ) in terms of the exponential integral z dt e −t /t.Note that this functional form behaves at large τ like D lat 1 (τ 2 ) ≈ 2Bλ b e −τ /λ b /τ , which is different from the functional form used for D B (τ 2 ).From these we obtain the contribution to E 3 from D(τ 2 ) given by Remarkably, this result has a central value that is very close to the result in Eq. ( 38), but with much larger uncertainty.We take the result in Eq. ( 38) as our estimate for the contribution from D(τ 2 ) to E 3 .We combine the results in Eqs.(37) and (38) to obtain where the uncertainties are added in quadrature.Note that the uncertainties are dominated by the unknown D(τ 2 ).This is our final numerical result for E 3 at the MS scale Λ = 1 GeV.
In the case of E 1 , we have established that the contribution from D 1 (τ 2 ) vanishes not just in the RGZ theory, but in general, as long as power divergences are properly subtracted in dimensional regularization.Hence, the only contribution to E 1 comes from D(τ 2 ).We can again estimate this contribution by comparing our result for D B (τ 2 ) in the RGZ theory with the lattice measurement.
We have If instead we use the lattice measurement for D * (τ 2 ) in Ref. [35] to obtain D 1 (τ 2 ), then we obtain which is compatible with Eq. ( 41), but with a larger uncertainty.Note that neither results are precise enough to determine the sign.We take Eq. ( 41) as our estimate for E 1 because the uncertainty is smaller.
Next, B 1 in the RGZ theory at tree level can be computed from Eq. (34).We obtain where the uncertainty comes from the uncertainty in m 2 .
Similarly to E 3 and E 1 , this result is missing the contribution from D(τ 2 ).By adding this contribution, which is equal to our estimate for E 1 in Eq. ( 42), we obtain which is our final numerical result for B 1 .
Finally, we compute iE 2 in the RGZ theory: where the uncertainties come from the parameters in the RGZ theory.Similarly to E 3 , E 1 , and B 1 , this result is missing the contribution from D(τ 2 ).We estimate the contribution from D(τ 2 ) to iE 2 in the same way as we have done for E 3 .We have If we instead use only the lattice QCD results in Ref. [35] to estimate this contribution, then we obtain −0.13 ± 0.30, which is consistent with above but with larger uncertainties.We combine Eqs. ( 45) and ( 46) to obtain our final numerical result for iE 2 : We may compare our numerical results based on the RGZ theory with estimates based on the lattice QCD analysis in Ref. [35].In the case of E 1 , we would obtain the same result as in Eq. ( 42), since the contribution from D 1 (τ 2 ) cancels in E 1 and we obtain the contribution from D(τ 2 ) by using the lattice results in Ref. [35].In the case of B 1 , we may use the long-distance functional form of D B (τ 2 ) given in Ref. [35] by D ⊥ (τ 2 ) = Ae −τ /λa to obtain This result is much smaller than what we obtained in Eq. ( 44) based on the RGZ theory.Note that because we used the long-distance functional form used in the lattice QCD analysis, the quadratic power divergence expected from the perturbative QCD calculation is completely missing in Eq. ( 48).In the case of the dimensionless moment E 3 , the perturbative power-law contribution is both UV and IR divergent, so that it is not possible to obtain a dimensionally regulated result simply from the long-distance functional forms used in the lattice QCD analysis.That is, if we compute the moment E 3 from the long-distance functional forms obtained in the lattice QCD analysis, the result will be missing the perturbative order-α s scheme-dependent finite pieces as well as the logarithmic dependence on the renormalization scale.Simply computing the moment from the longdistance functional form of D E (τ 2 ) = D ⊥ (τ 2 ) + D * (τ 2 ) from Ref. [35] gives Although the central value is similar in size to the MSrenormalized result in the RGZ theory at Λ = 1 GeV in Eq. ( 40), it is completely ambiguous which to renormalization scheme and scale this result corresponds.Also the uncertainty in this lattice-based estimate is too large to be phenomenologically useful.Similarly, a lattice estimate of iE 2 gives which is compatible with our result in Eq. ( 47), but with much larger uncertainties.Now that we have obtained our numerical estimates for the dimensionless moment E 3 in the MS scheme and the dimension-two moments E 1 and B 1 in dimensional regularization, we can proceed to discuss phenomenological applications in the following section.

IV. QUARKONIUM DECAYS
Based on the results for the moments of two-point field-strength correlators in the previous section, we now compute the color-octet matrix elements that appear in heavy quarkonium decay rates and discuss phenomenological applications.
We begin with the inclusive decays of χ QJ states (Q = c or b), which involve the dimensionless moment E 3 .The relation between E 3 and the color-octet matrix element is given by [17,18] where m Q is the heavy quark pole mass, and we take the definitions of the NRQCD operators O 8 ( 3 S 1 ) and O 1 ( 3 P J ) in Refs.[16,18].The subscripts 1 and 8 stand for the color state of the Q Q that is created and annihilated by the NRQCD operator, and the spectroscopic notation 2S+1 L J is used to indicate the spin, orbital, and total angular momentum state of the Q Q.The factor m 2 Q is included on the left-hand side in order to make the ratio dimensionless.This relation has been obtained in Refs.[17,18] in three spatial dimensions1 .Since both ⟨χ QJ |O 8 ( 3 S 1 )|χ QJ ⟩ and E 3 contain logarithmic UV divergences, the relation must be obtained in d = 4−2ϵ spacetime dimensions in order to facilitate correct renormalization of both quantities in the MS scheme; the above result valid for any d can be obtained by rederiving the relation in arbitrary spacetime dimensions.It is easy to see that the denominator factor d − 1 comes from the tensor δ ij /(d − 1) which arises from taking an average over spatial directions of a rank-2 Cartesian tensor in Eq. ( 73) of Ref. [18].Due to the UV poles in the unrenormalized color-octet matrix element and the unrenormalized moment E 3 , the following relation holds between the MS-renormalized color-octet matrix element and the moment: where the extra finite piece comes from the product of the order-ϵ term in the d-dependent coefficient 3/(d − 1) = 1 + 2 3 ϵ + O(ϵ 2 ) and the 1/ϵ pole in the unrenormalized E 3 .Following Refs.[61,62] we refer to this ratio as ρ (Λ) 8 .In phenomenological determinations of the quantity E 3 , such extra finite terms were not taken into account, and so, when comparing with the MS-renormalized moment E (Λ) 3 , the phenomenologically determined E 3 must be compared with the quantity in the parenthesis of Eq. ( 52) that contains the extra finite term.By using our result for E (Λ=1 GeV) 3 in Eq. ( 40), we obtain This agrees within uncertainties with a recent phenomenological determination 2.05 +0.94 −0.65 in Ref. [60], but with much smaller uncertainties.It is straightforward to compute χ QJ decay rates from this result.The NRQCD factorization formula for χ QJ decay rates to light hadrons at leading order in v is given by [16] Γ where 2 Imf (Λ) 1 ( 3 P J ) and 2 Imf 8 ( 3 S 1 ) are short-distance coefficients which begin at order α 2 s and are known to order-α 3 s accuracy [26]; exceptionally 2Imf (Λ) 1 ( 3 P 1 ) vanishes at order α 2 s because a spin-1 state cannot decay into two gluons at tree level.The short-distance coefficients for the color-singlet channels contain logarithmic dependences on Λ that cancel the scale dependence of the color-octet matrix element.While the color-octet matrix element ⟨χ QJ |O 8 ( 3 S 1 )|χ QJ ⟩ (Λ) can be reexpressed in terms of E (Λ) 3 and the color-singlet matrix element ⟨χ QJ |O 1 ( 3 P J )|χ QJ ⟩ by using Eq. ( 52), the color-singlet matrix element can be computed in terms of quarkonium wave functions: where R χ QJ (r) is the radial wave function of the χ QJ state, which can be obtained by solving a Schrödinger equation.Rather than using potential models to compute |R ′ χ QJ (0)| 2 , which can vary greatly depending on the choice of the model (see, for example, Ref. [60]), we take the first-principles calculations of the P -wave wave functions in Ref. [29].Because we include corrections of relative order α s in the short-distance coefficients, we include the one-loop correction to the wave function which comes from the radiative correction to the static potential, as has been done in Ref. [29].We neglect higherorder corrections that were included in Ref. [29] that are associated with two-loop corrections to the wave function, because that would require including corrections of relative order α 2 s to the short-distance coefficients that are currently unknown.We also include a part of the relativistic correction that comes from the phase space by replacing an overall factor of 1/(2m Q ) by 1/M χ QJ as have been done in Ref. [29].Consistently with Ref. [29], we set m c = 1.316GeV and m b = 4.743 GeV, and compute α s in the MS scheme by using RunDec at four-loop accuracy at scales µ = 2.5 GeV for Q = c and µ = 5 GeV for Q = b.We use the quarkonium masses M χ cJ = 3.47 GeV, M χ bJ (1P ) = 9.94 GeV, M χ bJ (2P ) = 10.26GeV, and M χ bJ (3P ) = 10.53GeV, which were computed in Ref. [29] and agree within 0.1 GeV with the Particle Data Group (PDG) values in Ref. [63].In the short-distance coefficients, we set the number of light quark flavors n f to be n f = 3 for Q = c and n f = 4 for Q = b, and fix Λ = 1 GeV.We obtain the following decay rates for χ cJ : Γ(χ c0 → LH) = 13.5 where the first uncertainties come from the uncertainty in E 3 , the second uncertainties come from varying µ between 1.5 and 4 GeV, and the third uncertainties come from uncalculated corrections of order v 2 , which we estimate by 0.3 times the central values, reflecting that typically v 2 ≈ 0.3 for charmonia.In the last equalities we add the uncertainties in quadrature.These results may be compared directly with experiment; however, the total widths of χ cJ include sizable contribution from radiative decays into J/ψ + γ, especially for J = 1 and 2, which will not be captured by the calculation of Γ(χ cJ → LH).After subtracting these radiative decay rates by using the measured branching fractions, we obtain from Ref.When compared with the calculation in Ref. [60] based on the phenomenological determination of E 3 , the precision has improved for the decay rate of χ c1 , while the uncertainty in the decay rate of χ c0 is increased.While the improvement in precision has mostly to do with the improved determination of E 3 and the calculation of the wave functions from first principles, the increased uncertainty mainly comes from the fact that we included an uncertainty from variation of the QCD renormalization scale, while in Ref. [60] the uncertainties from uncalculated corrections of higher orders in α s were estimated to be α 2 s times the central value.The size of the uncertainty in the decay rate of χ c2 is comparable to the result in Ref. [60].
Similarly, we obtain the following decay rates for χ bJ (1P ): where the first uncertainties come from the uncertainty in E 3 , the second uncertainties come from varying µ between 2 GeV and 8 GeV, and the third uncertainties come from uncalculated corrections of order v 2 , which we estimate by 0.1 times the central values, reflecting that typically v 2 ≈ 0.1 for bottomonia.In the last equalities we add the uncertainties in quadrature.Unfortunately, experimental results for the total widths of χ bJ states have not been made available yet.There are, however, determinations based on the theoretical calculations of the radiative decay rates into Υ + γ and the measured branching fractions of the same process.In Ref. [64], the radiative decay rates were computed in potential NRQCD at weak coupling, from which the total decay rates were computed by using the measured branching fractions.After subtracting the radiative decay rates we obtain from Ref. [64] the values Γ(χ b0 (1P ) → LH)| Ref. [ 57).In Ref. [65], the authors took the results for the radiative decay rates computed in a nonrelativistic model based on the Cornell potential available in Table 4.16 in Ref. [66] and used the measured results for the branching fractions in Ref. [63] to compute the total decay widths.We again subtract the radiative decay rates from the results listed in Ref. [65] to obtain Γ(χ b0 (1P ) → LH)| Ref. [65] = 1.121 +0.185 −0.140 MeV, Γ(χ b1 (1P ) → LH)| Ref. [65] = 0.051 +0.005 −0.004 MeV, and Γ(χ b2 (1P ) → LH)| Ref. [65] = 0.144 +0.010 −0.009 MeV.These are in good agreement with our results in Eq. ( 57).Note that, however, the uncertainties in the results in Ref. [65] may be underestimated because they reflect the ones in the measured branching fractions only, as no uncertainty is given in the model calculation of the radiative decay rates in Ref. [66].
We also compute decay rates for χ bJ (2P ) and χ bJ (3P ) states: Compared to the results in Ref. [60] based on the phenomenological determination of E 3 , the decay rates for defined in Eq. ( 52) at Λ = 1 GeV in this work (RGZ), compared to determinations based on χcJ decay rates in Ref. [17] (BEPSV01) and Ref. [60] (BCMV20), and the determination based on measurements of Br(χ bJ (1P ) → D0 +X) in Ref. [62] (CLEO).J = 0 and 2 are in agreement within uncertainties.Although the results for J = 1 also agree within uncertainties with Ref. [60], we obtain smaller central values for the decay rates.The uncertainties in the χ bJ decay rates are comparable in size with the results in Ref. [60]; however, the uncertainties would have been much smaller if we estimated the uncertainties from uncalculated corrections of higher orders in α s in the same way as Ref. [60], instead of varying the QCD renormalization scale.
In Refs.[61,62], the dimensionless ratio ρ (Λ) 8 defined in Eq. ( 52) for the bottomonium state at the scale Λ = 4.6 GeV was investigated by computing the decays χ bJ → c + X and comparing them with measurements of branching fractions Br(χ bJ → D 0 + X).If we take into account the evolution of the color-octet matrix element, then our results lead to ρ (Λ=4.6 GeV) 8 = 0.17 ± 0.01.This is compatible with the result 0.160 +0.071 −0.047 obtained in Ref. [62] from the measured branching fractions for the 1P states.In contrast, Ref. [62] obtained a smaller result 0.074 +0.010 −0.008 from measurements involving 2P states, while we expect Eq. ( 52) to be approximately independent of the radial excitation; however, we note that the quality of the fit for the 2P states in Ref. [62] is much worse compared to the 1P states, resulting in a value of χ 2 that is more than an order of magnitude larger than that of the 1P case.
In Fig. 2 we compare our result for the χ QJ matrix element ratio ρ (Λ) 8 defined in Eq. ( 52) at Λ = 1 GeV with the phenomenological determinations in Ref. [17] (BEPSV01) and Ref. [60] (BCMV20) based on χ cJ decay rates, and the determination in Ref. [62] (CLEO) from measurements of Br(χ bJ (1P ) → D 0 + X).Note that because the result in Ref. [62] has been obtained at the scale of the bottom quark mass, we solved the evolution equation for E 3 in Eq. ( 13) to obtain the result at Λ = 1 GeV.We can see that from Fig. 2, our result based on the RGZ theory has much smaller uncertainties compared to the phenomenological determinations based on measured decay rates, although the results are consistent within uncertainties.We now list our numerical results for the color-octet matrix elements for χ cJ and χ bJ states obtained from the RGZ result for the moment E 3 and the calculation of the P -wave quarkonium wave functions in Ref. [29] at one-loop level: Here we use the shorthand ) .These results are much more precise compared to previous phenomenological determinations in Ref. [60], with uncertainties reduced by a factor of about 1/3.Note that all of the color-octet matrix elements in Eq. ( 60) are computed in the MS scheme at scale Λ = 1 GeV, and any direct comparison with other results must be done at the same scheme and scale.
Next we look into the decays of S-wave states into light hadrons.In the case of η Q , where Q = c or b, the following relations hold for the color-octet matrix elements and the dimension-two moments [18]: where c F = 1 + O(α s ) is the short-distance coefficient in the NRQCD Lagrangian associated with the spin-flip interaction [67][68][69]; even though c F is known to three loops, we only need the result at tree level consistently with our tree-level evaluation of B 1 .We take the definitions of the operators O 8 ( 3 S 1 ), O 1 ( 1 S 0 ), and O 8 ( 1 P 1 ) in Refs.[16,18,27].We have included a factor of 1/m 2 Q on the left-hand side of Eq. (61b) in order to make the ratio dimensionless.Similarly to the χ QJ case, these ratios do not depend on the radial excitation of the η Q state, and the only dependence on the heavy quark flavor comes from the heavy quark mass m Q .Note that Eq. (61b) is taken from Eq. (155) in Ref. [18], where one needs to set J = 1 for the relation to hold for the spin-singlet case.There is an additional color-octet matrix element ⟨η Q |O 8 ( 1 S 0 )|η Q ⟩ that enters the decay rate of spin-singlet S-wave quarkonium into light hadrons; we do not consider this matrix element because it is proportional to a moment of the four-point field-strength correlation function [see Eqs. ( 80) and (153) of Ref. [18]], which is outside the scope of this work.Before we consider the decay rate of η Q let us compare Eqs.(61) with the estimates of the sizes of the color-octet matrix elements given in Ref. [27].In Ref. [27], the ratio in Eq. (61a) has been estimated to be the order of v 3 /(2N c ), which corresponds to about 0.03 for charmonium, and 0.005 for bottomonium.By using our result for B 1 in Eq. ( 44), we obtain −0.1 and −0.008 for the right-hand side of Eq. (61a) for charmonium and bottomonium, respectively.These results are larger in size than the estimates given in Ref. [27], and the signs are negative due to the positivity of B 1 .In the case of Eq. (61b), this ratio was estimated in Ref. [27] to be about v 4 /(2N c ), which gives 0.015 and 0.0017 for charmonium and bottomonium, respectively.If we use our result for E 1 in Eq. ( 41), then the ratio ranges between −0.014 and 0.008 for charmonium, and between −0.0010 and 0.0006 for bottomonium, which are smaller than the estimates in Ref. [27].
The sizes of the contributions from the color-octet matrix elements in Eqs. ( 61) to the decay rates Γ(η Q → LH) can be computed by multiplying the corresponding shortdistance coefficients.Numerical sizes of the tree-level short-distance coefficients compared to the one at leading order in v are listed in Ref. [27].Compared to the decay rate at leading order v, which only involves the color-singlet matrix element ⟨η Q |O 1 ( 1 S 0 )|η Q ⟩, the relative contribution from the color-octet matrix element ⟨η Q |O 8 ( 3 S 1 )|η Q ⟩ is given by Eq. (61a) times 0.75n f .If we put n f = 3 for charmonium and n f = 4 for bottomonium, then the correction to the decay rate from B 1 is about −22 ± 3% for η c , and about −2.3 ± 0.3% for η b .In the case of the matrix element ⟨η Q |O 8 ( 1 P 1 )|η Q ⟩, the contribution from E 1 relative to the leading-order decay rate is given by 1.13 times Eq.(61b).Because of the small size of E 1 , the contribution from E 1 to the decay rate is at most of the order of 1% for η c , and is at most of the order of 0.1% for η b .The order-v 4 correction from the matrix element ⟨η Q |O 8 ( 1 S 0 )|η Q ⟩ involving the fourpoint correlation function would be of the order of 3% for charmonium and 0.3% for bottomonium according to the estimate given in Ref. [27].Hence, the only significant color-octet contribution to the η Q decay rate comes from B 1 .
In order to properly quantify the color-octet contribution to the η Q decay rate, we need to include radiative corrections to the short-distance coefficients, which are known to converge slowly [24].In Ref. [25], the authors computed the ratio R η Q = Γ(η Q → LH)/Γ(η Q → γγ) by resumming a class of radiative corrections in the largen f limit, based on an earlier resummation calculation in Ref. [24] and a fixed-order calculation at two-loop accuracy in Ref. [70].It was found in this work that the lack of knowledge of the color-octet matrix element ⟨η Q |O 8 ( 3 S 1 )|η Q ⟩ is a significant source of uncertainty.
Since we can determine this matrix element from our result for B 1 , we can remove this uncertainty.Since the calculation in Ref. [25] was done by cutting off the quadratic power divergence in the color-octet matrix element with a hard cutoff, first we need to convert our result for B 1 from dimensional regularization to cutoff regularization.The power divergence in B 1 is simply given by a scaleless power-divergent integral which can be read off from Eq. ( 24) by setting all dimensionful parameters to zero: where Λ is the hard cutoff on the spatial loop momentum, consistently with what was done in Ref. [25].At Λ = 1 GeV, which was used in Ref. [25] for charmonium, Eq. ( 62) amounts to 0.15 GeV 2 if we set α s (m c ) = 0.35.At Λ = 2 GeV, which was used in Ref. [25] for bottomonium, the cutoff-regulated B 1 is larger than the dimensionally regulated one by 0.37 GeV 2 if we use α s (m b ) = 0.22.We can now improve the result in Ref. [25] by including the color-octet contribution we compute with our result for B 1 .We first compute R ηc for the charmonium states.We obtain, from Eq. ( 66) of Ref. [25], where we shifted the central value by the correction from the color-octet matrix element.The first uncertainty comes from the variation of the QCD renormalization scale between 1 and 2 GeV, and the second uncertainty comes from B 1 .Because we have now included the effect of the color-octet matrix element in the correct scheme, we have removed the uncertainties arising from varying the cutoff Λ and from the unknown color-octet matrix element.Here, NNA refers to the naïve non-Abelianization scheme used for the resummation of radiative corrections in the large-n f limit [71].The size of the correction from the color-octet matrix element can be easily obtained by replacing the estimate ⟨η 3 used in Ref. [25] with the pNRQCD result Eq. (61a) and by using our result for B 1 .Here, we take m c = 1.5 GeV consistently with Ref. [25].The uncertainties are slightly reduced compared to the result in Ref. [25], because the uncertainty in the NNA result is dominated by the variation of the QCD renormalization scale.Similarly, from Eq. (67) of Ref. [25] we obtain which again we obtain by shifting the central value by the correction from the color-octet matrix element.Here, BFG refers to the background-field gauge method that was used for the resummation of radiative corrections in the large-n f limit [72].The sources of uncertainties are same as in Eq. ( 63).The uncertainty in the BFG result [73] (HPQCD23), large-n f resummation calculations in Ref. [25] (BCK18) and Ref. [24] (BC01), a fixed-order calculation at two-loop accuracy in Ref. [70] (FJS17).Results for the η b case are not available in HPQCD23 [73] and BC01 [24].
The PDG value for the ηc case from Ref. [63] is shown as a blue band.
|⟨O 8 ( 1 P 1 )⟩ ηc(1S) | ≲ 6.5 × 10 −3 GeV 5 , (70a) Here we again use the shorthand In obtaining these results, we included the order-α s corrections to the wave function arising from loop corrections to the static potential, but neglected corrections of order α 2 s because short-distance coefficients for the color-octet channel are generally known up to next-to-leading order in α s .In the case of ⟨η Q |O 8 ( 1 P 1 )|η Q ⟩, we only show the upper estimates of the absolute values because our result for E 1 is not precise enough to determine its sign; nevertheless the upper estimates in Eq. ( 70) are smaller than the power-counting estimates in Ref. [27].
Next, we consider the decays of spin-triplet S-wave quarkonium V , which includes J/ψ and Υ.The following relations hold for the color-octet matrix elements for the V state and the dimension-two moments [18]: where J = 0, 1, and 2. We again take the definitions of the NRQCD operators O 8 ( 1 S 0 ), O 1 ( 3 S 1 ), and O 8 ( 3 P J ) in Refs.[16,18,27].Note that these relations can also be obtained from Eqs. (61) by using heavy-quark spin symmetry (see Ref. [27]).Similarly to the η Q case, we do not consider the matrix element ⟨V |O 8 ( 3 S 1 )|V ⟩, because it is proportional to a moment of the four-point fieldstrength correlator, which is outside of the scope of this work.The values of these color-octet matrix elements can be computed in the same way as the η Q matrix elements, by using ⟨V |O 1 ( 3 S 1 )|V ⟩ = 2N c |Ψ V (0)| 2 and the firstprinciples calculations of the quarkonium wave functions Ψ V (r) in Ref. [28].We obtain Again we use the shorthand These results are same as what could be obtained from the η Q matrix elements in Eqs. ( 69) and ( 70) and heavy-quark spin symmetry.Similarly to the η Q case, we only show the upper estimates of the absolute values of ⟨V |O 8 (3 P J )|V ⟩.
The results for the matrix elements we obtain can be used to estimate sizes of the color-octet contributions to the decay rates Γ(V → LH).As have been known from Refs.[26,27], many of the short-distance coefficients associated with the color-octet matrix elements are enhanced by a factor of π/α s , because a heavyquark antiquark pair in color-octet states can decay into two gluons or a light-quark pair at tree level, while a color-singlet S-wave spin-triplet state can only decay into three or more gluons at tree level.For example, the size of the color-octet contribution from the matrix element ⟨V |O 8 ( 1 S 0 )|V ⟩ to the decay rate Γ(V → LH) relative to the one at leading order in v is given by 11.64 × (π/α s ) times Eq.(71a) [27].Similarly, the contributions from the matrix elements ⟨V |O 8 ( 3 P 0 )|V ⟩ and ⟨V |O 8 ( 3 P 2 )|V ⟩ relative to the leading-order decay rate are given by Eq. (71b) times 34.93 × (π/α s ) and 9.3 × (π/α s ), respectively2 .Due to the large coefficients, the color-octet contributions are significant for Γ(V → LH), and can even exceed the color-singlet contribution at leading order in v in the charmonium case.For example, for J/ψ the correction from B 1 is about −3.1 times the leading-order decay rate, and the correction from E 1 ranges from −1.0 to +0.6 times the leading-order decay rate.The size of the contribution from the matrix element ⟨J/ψ|O 8 ( 3 S 1 )|J/ψ⟩ that involves the four-point correlation function would also be order one if we use the estimate for the matrix element given in Ref. [27].Even in the case of Υ, the color-octet contributions from ⟨Υ|O 8 ( 1 S 0 )|Υ⟩ and ⟨Υ|O 8 ( 3 S 1 )|Υ⟩ can be as large as 50% of the leading-order decay rate, although the correction from ⟨Υ|O 8 ( 3 P J )|Υ⟩ is expected to be at most ±0.15 times the leading-order decay rate.Since the contributions from B 1 to the decay rate is negative, it is possible that large cancellations may occur between color-octet contributions.Nonetheless, even in the case of Υ decays, knowledge of the matrix element ⟨V |O 8 ( 3 S 1 )|V ⟩ would be necessary to compute the decay rate to some precision.
Aside from the exceptional case of the decays of spintriplet S-wave quarkonium into light hadrons, the moment E 1 may be neglected in most cases due to its tiny size.For example, E 1 appears in order-v 2 corrections to the electromagnetic production and decays of χ QJ (see Refs. [18,60]).The correction from E 1 to the two-photon decay amplitude of χ QJ is given by 4  9 E 1 /m 2 Q for J = 0 and 2  3 E 1 /m 2 Q for J = 2 relative to the leading-order result [60].These are at most of the order of 5% compared to the tree-level leading-order result for the charmonium case, and less than 1% relative to the leading-order result for the bottomonium case, so they may be neglected in the phenomenology of P -wave quarkonium decays at the current level of accuracy.
Finally, we discuss phenomenological applications of iE 2 .This quantity appears in the correction to the relation between the color-singlet matrix element and the wave function at the origin for P -wave states [Eq.( 55)], which was computed in Ref. [60].This correction was later found to be canceled by the correction to the P -wave wave function at the origin coming from the velocitydependent potential at zero separation (see Ref. [29]).Hence, iE 2 does not have a phenomenological significance in P -wave quarkonium decays or production.However, iE 2 does appear in the correction to the wave function at the origin for S-wave quarkonia that arises from the velocity-dependent potential at zero separation [28,29], which unlike the P -wave case is not canceled by corrections to the relations The correction to the S-wave wave function at the origin Ψ(0) relative to the leading-order result is given by −iE 2 /(3m Q ), which can be obtained from the calculation in Ref. [28] and the identification of the velocity-dependent potential at zero separation in terms of iE 2 in Ref. [29].In Ref. [28] this correction was considered only in the uncertainties, which were estimated by assuming that |2iE 2 /3| ≲ 0.5 GeV.This estimate is more than twice as large compared to the central value of our result for iE 2 in Eq. ( 50), and hence, our result for iE 2 can be used to significantly reduce this uncertainty.By including the correction from the velocity-dependent potential at zero separation, the S-wave charmonium wave functions at the origin are enhanced by 7 ± 4%, and the S-wave bottomonium wave functions at the origin are enhanced by 2 ± 1%.In Ref. [28], the uncertainty from the unknown iE 2 was estimated to be about 19% and 5% for S-wave charmonium and bottomonium wave functions at the origin, respectively; by using our result for iE 2 , these uncertainties would be greatly reduced to only about 4% and 1%, respectively.For example, using our result for iE 2 would improve the first-principles calculation of the J/ψ leptonic decay rate in Ref. [28] from 4.5 +2.5 −1.9 to 5.1 +1.6 −1.4 keV, which not only reduces uncertainties but also brings the central value closer to the PDG value 5.53±0.10keV [63].In the case of the two-photon decay rate of η b (1S), we obtain only a small improvement in uncertainty, from 0.433 +0.165  −0.065 to 0.450 +0.149 −0.047 keV.When combined with our results for R η b in Eq. ( 66), we obtain the total decay rate of η b (1S) given by 11.3 +3.8 −1.2 MeV, which is slightly improved compared to Eq. (68a).

V. SUMMARY AND DISCUSSION
In this work we computed color-octet nonrelativistic QCD matrix elements that appear in heavy quarkonium decay rates based on the refined Gribov-Zwanziger theory [44][45][46][47][48][49].The color-octet matrix elements correspond to the probabilities for a heavy quark and antiquark pair in a color-octet state to evolve into a color-singlet state, and can be related to the moments of correlation functions of QCD field-strength tensors by using the potential nonrelativistic QCD formalism [17][18][19][20][21].We found that tree-level calculations in the refined Gribov-Zwanziger theory, which can reproduce the nonperturbative gluon propagator very well, can also give a satisfactory description of the two-point correlation functions of the QCD field-strength tensor that agrees with a lattice QCD study in Ref. [35].This allowed us to compute moments of the correlation functions that are infrared finite, while also reproducing the correct leading ultraviolet divergences that are expected from the nonrelativistic QCD factorization formalism [16] and can be properly renormalized.In particular, the dimensionless moment E 3 defined in Eq. ( 12) contains a logarithmic ultraviolet divergence, which we renormalize in the MS scheme consistently with calculations in the nonrelativistic QCD factorization formalism.Similarly, in the case of the dimension-two moments E 1 and B 1 defined in Eq. ( 14), which involve quadratic power divergences, correct results in dimensional regularization are obtained by computing them in arbitrary spacetime dimensions and later setting the number of dimensions to 4. Such a calculation in dimensional regularization would be difficult to carry out in a numerical study, such as a calculation based on lattice QCD data, especially for the dimensionless moment E 3 where the conversion to dimensional regularization would involve not only ultraviolet divergences but also infrared divergences.
The numerical result for E 3 that we obtain in the refined Gribov-Zwanziger theory is given in Eq. ( 40).This result is much more precise than the results from phenomenological analyses in Refs.[17,60].We have used this result to compute the color-octet matrix element for P -wave quarkonium decays and to compute decay rates of χ cJ and χ bJ into light hadrons.The uncertainties that arise from color-octet matrix element have significantly improved compared to previous phenomenological stud-ies, although the results are generally compatible with previous results and experimental data.
In the case of the dimension-two moment B 1 , we found that a straightforward calculation in the refined Gribov-Zwanziger theory leads to not only a power ultraviolet divergence, as expected in perturbative QCD, but also a logarithmic ultraviolet divergence that is proportional to a nonperturbative parameter associated with the dimension-two condensate of the gauge field.While appearance of subdivergences is not completely unexpected given that the refined Gribov-Zwanziger theory contains dimensionful parameters, due to the nonperturbative origin of the subdivergences, they could not be subtracted within the nonrelativistic QCD factorization formalism.Fortunately, we found that at tree level, this result is exactly proportional to the dimension-two condensate of the gauge field which also contains a logarithmic ultraviolet divergence as was found in Refs.[49,52].By using the finite result for the condensate obtained in the effective action formalism in Ref. [52], we obtain a finite result for B 1 in dimensional regularization, which is given in Eq. ( 44).This result lead to the calculation of the color-octet matrix elements that appear in S-wave quarkonium decays associated with the spin-flip interaction.We have used this result to improve a previous calculation of the two-photon branching fraction of η c and η b in Ref. [25], where the color-octet matrix elements have been left unknown.By including the color-octet matrix elements computed from the calculation of B 1 in the refined Gribov-Zwanziger theory, the central values of the branching fractions have increased, and the uncertainties associated with the color-octet matrix elements have reduced.
Finally, we found that the dimension-two moment E 1 vanishes at tree level in the refined Gribov-Zwanziger theory.We also found that in general, the contribution from the invariant function D(x 2 ) vanishes in E 1 , and only D 1 (x 2 ) contributes to E 1 , which in perturbation theory appears from next-to-leading order in the strong coupling.We have estimated the size of the contribution from D 1 (x 2 ) to E 1 by using lattice QCD measurements of the invariant functions in Ref. [35] to obtain the result in Eq. (41).The result suggests that E 1 is at most about an order of magnitude smaller than B 1 , and so, in most phenomenological applications, contributions from E 1 may be neglected, unless it is multiplied by an unusually large short-distance coefficient.
The results for the two-point field-strength correlation function in this work are based on tree-level calculations of the gluon propagator in the refined Gribov-Zwanziger theory.Similarly to the case of the gluon propagator, tree-level results in this theory seem to be able to reproduce bulk of the long-distance nonperturbative behavior of the correlation function, as we find agreements with the lattice QCD analysis in Ref. [35].We also found that the Gribov-Zwanziger theory alone [37], without the refinements from the dimension-two condensates, cannot reproduce the long-distance behavior found in lattice studies.This suggests that inclusion of the dimensiontwo condensates is necessary in describing the nonperturbative nature of the field-strength correlation functions, similarly to the case of the gluon propagator.Although the results in this work are based on calculations at tree level, in obtaining the numerical results we included effects of contributions at higher orders in the strong coupling estimated by comparing with lattice QCD results.It would still be interesting to compute the correlation functions at next-to-leading order accuracy, which may also help reduce the uncertainties in the numerical results.We also note that in this work we have not computed the color-octet matrix element that arises from the four-point field-strength correlation function.In principle, we could compute moments of the four-point correlation function at tree level similarly to our calculation of the two-point correlation function.However, we found that such a calculation leads to logarithmic ultraviolet divergences, which, unlike the case of B 1 , involve not only m 2 but also M 2 and λ, which makes it impossible to reexpress the divergence in terms of the dimension-two condensate of the gluon field.Extending the analysis to the case of the four-point correlation function and also to next-to-leading order accuracy would be important in the phenomenology of J/ψ and Υ decays, as the shortdistance coefficients associated with color-octet matrix elements are very large for these processes.

10 FIG. 1 .
FIG.1.Results for the invariant functions DB(x 2 ) (left panel) and DB(x 2 ) − DE(x 2 ) (right panel) from the RGZ theory (black bands) shown against results in the GZ theory (black dashed lines) and lattice QCD results from Ref.[35] (orange bands).