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

The two-point correlation function of the QCD field-strength 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 $\bar{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 [19, 20, 17, 18, 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 color-octet contributions are necessary for improving the theory description of decay rates and two-photon branching fractions of $\eta_{c}$ [24, 25]. It is also known that color-octet contributions are enhanced by an inverse power of $\alpha_{s}$ in the case of $J/\psi$ and $\Upsilon$ 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 long-distance 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 $\partial\cdot 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, 38, 39, 37, 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 [52, 49, 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 color-octet 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 field-strength 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 color-octet contributions to decay rates of heavy quarkonia in Sec. IV. We conclude in Sec. V.

We define the two-point correlation function of the QCD field-strength tensor in Euclidean space as

where $|\Omega\rangle$ is the QCD vacuum, $G_{\mu\nu}^{a}=\partial_{\mu}A_{\nu}^{a}-\partial_{\nu}A_{\mu}^{a}+gf^{abc}A_{\mu}^{b}A_{\nu}^{c}$ is the QCD field-strength tensor, $A_{\mu}$ is the gauge field, $g$ is the strong coupling, $T_{F}=1/2$, and $\Phi_{ab}(x,0)$ is a straight Wilson line defined by

where $P$ is the path ordering and $A^{\rm 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 ${\cal D}(x^{2})$ and ${\cal D}_{1}(x^{2})$ are invariant functions of $x^{2}$. The form (II) follows from the fact that ${\cal D}(x^{2})$ vanishes in the Abelian gauge theory. Because of this, ${\cal D}(x^{2})$ vanishes at tree level in the non-Abelian gauge theory.

At tree level we can compute ${\cal D}_{1}(x^{2})$ at order $\alpha_{s}$ while ${\cal D}(x^{2})=0+O(\alpha_{s}^{2})$, where $\alpha_{s}=g^{2}/(4\pi)$. Useful representations of the invariant functions can be obtained by defining the two linear combinations

where $d=4-2\epsilon$ is the number of spacetime dimensions, and we work in a frame where $x^{\mu}$ is aligned to the time $(\mu=0)$ direction. Note that in this frame, ${\cal D}_{E}(x^{2})$ involves only the chromoelectric fields and ${\cal 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 $\lambda^{4}=2g^{2}N_{c}\gamma^{4}+m^{2}M^{2}$ with $\gamma$ the Gribov parameter. Note that the tree-level propagator in the GZ theory is obtained by setting $m$ and $M$ to zero, while keeping $\gamma$ nonzero [37]. The values of $m$, $M$, and $\lambda$ 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 $\beta=6/g^{2}$ set to $5.7$, $6.0$, and $6.2$. They read $M^{2}=2.15\pm 0.13$ GeV², $m^{2}=-1.81\pm 0.14$ GeV², and $\lambda^{4}=0.26$ GeV⁴. By using Eq. (II) we obtain

where $C_{F}=(N_{c}^{2}-1)/(2N_{c})$. 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 $\pm i\sqrt{\bm{k}^{2}+(M^{2}+m^{2}\pm Q^{2})/2}$ and $\pm i\sqrt{\bm{k}^{2}+(M^{2}+m^{2}\mp Q^{2})/2}$, where $Q^{2}=i\sqrt{4\lambda^{4}-(M^{2}+m^{2})^{2}}$. Note that $Q^{2}$ is purely imaginary for the values of the parameters determined in Ref. [48]. Because $\sqrt{x^{2}}>0$, we must close the contour on the upper half plane. We obtain

where $\kappa^{\pm}=\sqrt{\bm{k}^{2}+(M^{2}+m^{2}\pm Q^{2})/2}$. Note that the quantity in the square brackets can be written as twice the real part of the first term. To obtain ${\cal 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 $\lambda$ to zero, which leads to the perturbative QCD result ${\cal D}_{B}(x^{2})|_{\rm pQCD}=g^{2}N_{c}C_{F}/(\pi^{2}x^{4})+O(\alpha_{s}^{2})$, in agreement with Ref. [30].

We show the numerical result for ${\cal D}_{B}(x^{2})|_{\rm RGZ}$ in Fig. 1 compared to the lattice measurement from Ref. [35]. For the numerical calculation we set the strong coupling to be $\beta=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 ${\cal D}_{B}(x^{2})$ corresponds to $D_{\perp}(x^{2})$ in Ref. [35]. We display the lattice measurement by using the functional form given by ${D}_{\perp}^{\rm lat}(x^{2})=Ae^{-|x|/\lambda_{a}}$ with $A=0.94^{+0.32}_{-0.16}$ GeV⁴ and $\lambda_{a}=0.120^{+0.009}_{-0.012}$ fm determined in Ref. [35]. Because this result was obtained from lattice data for $\sqrt{x^{2}}>0.3$ fm, and lattice data are only shown for $\sqrt{x^{2}}<0.6$ fm in Ref. [35], we only display the result from ${D}_{\perp}^{\rm lat}(x^{2})$ for $\sqrt{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, ${\cal D}_{B}(x^{2})|_{\rm RGZ}$ exhibits an exponentially decaying behavior at large $x^{2}$ with a finite correlation length, as can be seen from Fig. 1 for $\sqrt{x^{2}}$ larger than about 2 GeV^-1. A fit of the same functional form $Ce^{-|x|/\lambda_{c}}$ for $|x|$ larger than 2 GeV^-1 ($\approx 0.4$ fm) to the RGZ result leads to $C=0.58\pm 0.03$ GeV⁴ and $\lambda_{c}=0.77\pm 0.01$ GeV^-1. This correlation length reads in distance units $\lambda_{c}\approx 0.15$ fm. The lattice QCD study in Ref. [35] also found that at small $\sqrt{x^{2}}$ the function ${\cal 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 $\beta=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 $\gamma$ the same as the one obtained in the RGZ theory. We find that the GZ theory results in the long-distance behavior of ${\cal 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 ${\cal 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 ${\cal D}_{E}(x^{2})$. Because ${\cal D}(x^{2})$ vanishes at tree level, we could obtain $x^{2}\partial{\cal D}_{1}(x^{2})/\partial x^{2}$ at tree level by differentiating and multiplying by $x^{2}$ our result for ${\cal 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 ${\cal D}_{E}(x^{2})$. By using the definition for ${\cal 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 ${\cal D}_{B}(x^{2})$, we can obtain ${\cal 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 ${\cal D}_{B}(x^{2})-{\cal D}_{E}(x^{2})=-x^{2}\partial{\cal D}_{1}(x^{2})/\partial 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}_{*}^{\rm lat}(x^{2})=Be^{-|x|/\lambda_{b}}$ with $B=0.47^{+0.20}_{-0.06}$ GeV⁴ and $\lambda_{b}=0.189^{+0.013}_{-0.029}$ fm in Ref. [35]. Because the fit in Ref. [35] for ${D}_{*}^{\rm lat}(x^{2})$ was restricted to $\sqrt{x^{2}}>0.35$ fm, and lattice data are only shown for $\sqrt{x^{2}}<0.6$ fm, we only display the lattice result for $\sqrt{x^{2}}$ between $0.35$ and $0.6$ fm. Similarly to the case of ${\cal D}_{B}(x^{2})$, the result for ${\cal D}_{B}(x^{2})-{\cal 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 ${\cal D}_{B}(x^{2})-{\cal 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 ${\cal D}_{B}(x^{2})-{\cal D}_{E}(x^{2})|_{\rm pQCD}\approx 0.8/x^{4}$ for $\beta=6.0$. As we have done for ${\cal 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 $\gamma$ the same as the one obtained in the RGZ theory. We see that also in the case of ${\cal D}_{B}(x^{2})-{\cal 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 ${\cal D}(x^{2})$ appears from next-to-leading order accuracy, and so, at tree level only the ${\cal D}_{1}(x^{2})$ and its derivative contribute to ${\cal D}_{B}(x^{2})$ and ${\cal 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 ${\cal 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 ${\cal D}(x^{2})$ even at short distances. We may estimate the effects from a nonzero ${\cal 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 ${\cal D}(x^{2})$ compared to ${\cal D}_{1}(x^{2})$ at both short and long distances, which contradicts the naïve expectation from perturbation theory that ${\cal D}(x^{2})$ would be suppressed compared to ${\cal D}_{1}(x^{2})$, and is also in conflict with the smallness of ${\cal 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 ${\cal D}(x^{2})$ and ${\cal 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 gauge-invariant definition of the field-strength correlator (1), we expect that the results for ${\cal D}_{E}$ and ${\cal 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 $\alpha k^{2}/[k^{4}+2\alpha g^{2}N_{c}\gamma^{4}k^{2}/(k^{2}+M^{2})]\times k_{\mu}k_{\nu}/k^{2}$ compared to the Landau-gauge expression in Eq. (II), with $\alpha$ the gauge parameter. It is easy to check explicitly that this additional term makes vanishing contributions in ${\cal D}_{E}$ and ${\cal 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 field-strength correlation functions at tree level in the RGZ theory as dimensionally regulated momentum integrals in Eqs. (II) and (II). We will use these results in the following sections to compute moments of the correlation functions in dimensional regularization.

In this section we compute the moments of the two-point 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 ${\cal 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\bar{Q}$ color, and the sign arises from the fact that ${\cal E}_{3}$ was defined in Minkowski space in Ref. [18] while the right-hand side is computed in Euclidean space. Because the ${\cal D}_{E}(\tau)$ scales like $1/\tau^{4}$ at small $\tau$, the integral over $\tau$ contains a logarithmic UV divergence. After renormalization, ${\cal E}_{3}$ acquires a scale dependence given by

where $\Lambda$ is the renormalization scale for ${\cal E}_{3}$. Throughout this paper, we use the superscript $(\Lambda)$ to denote that the quantity is renormalized in the $\overline{\rm MS}$ scheme and depends on the $\overline{\rm MS}$ scale $\Lambda$. 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 ${\cal 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 ${\cal D}_{E}(x^{2})$ and ${\cal 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/\alpha_{s}$ compared to the one at leading order in $v$, which can make the color-octet contributions numerically significant [26, 27]. In the decays of spin-0 $S$-wave states, such as the $\eta_{c}$ and $\eta_{b}$, the lack of knowledge of the color-octet matrix element arising from ${\cal B}_{1}$ is a significant source of uncertainties [24, 25]. The quantity ${\cal E}_{1}$ also appears in the order-$v^{2}$ correction to two-photon decay rates of $P$-wave heavy quarkonia [18, 60]. Hence, accurate determinations of ${\cal E}_{1}$ and ${\cal B}_{1}$ are phenomenologically important. Note that both ${\cal E}_{1}$ and ${\cal B}_{1}$ contain power UV divergences at small $\tau$, 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 $i{\cal E}_{2}$ defined by

where the phase is chosen to recover the Minkowski space definition in Ref. [18]. Note that ${\cal E}_{2}$ is purely imaginary, so that $i{\cal E}_{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.