Matching coefficients in NRQCD to two-loop accuracy

TL;DR

This paper computes two-loop matching coefficients in NRQCD up to order 1/m^2, providing essential components for high-precision calculations of top quark pair production and heavy quarkonium decay.

Contribution

It presents the first two-loop calculations of key matching coefficients in NRQCD, advancing precision in theoretical predictions for heavy quark processes.

Findings

Analytic results for matching coefficients in a general covariant gauge.

Foundation for higher-order corrections in top quark and quarkonium physics.

Techniques for two-loop calculations in effective field theories.

Abstract

We consider the Lagrange density of non-relativistic Quantum Chromodynamics expanded up to order $1/m^2$, where $m$ is the heavy quark mass, and compute several matching coefficients up to two-loop order. Our results are building blocks for next-to-next-to-next-to-leading logarithmic and next-to-next-to-next-to-next-to-leading order corrections to the threshold production of top quark pairs and the decay of heavy quarkonia. We describe the techniques used for the calculation and provide analytic results for a general covariant gauge.

Tables

Table 1. Table 1: The coefficients R j k superscript subscript 𝑅 𝑗 𝑘 R_{j}^{k} introduced in Eq. ( 34 ) for the matching of the scattering amplitude.

$k$ 0 $j$ 1 $-1$ 2 $-1$ 3 $1$ 4 $1$ 5 $1$ 6 $1$ 7 $-1$ 8 $-1$ $k$ 0 1 9 $-d$ $2$ 10 $-d$ $2$ 11 $d$ $-2$ 12 $d$ $-2$ 13 $3d-2$ $-6$ 14 $3d-2$ $-6$ 15 $-3d+2$ $6$ 16 $-3d+2$ $6$ $k$ 0 1 2 17 $-d^2-4d+4$ $4d+8$ $-4$ 18 $-d^2-4d+4$ $4d+8$ $-4$ 19 $d^2+4d-4$ $-4d-8$ $4$ 20 $d^2+4d-4$ $-4d-8$ $4$ 21 $5d^2-4$ $-20d$ $20$ 22 $5d^2-4$ $-20d$ $20$ 23 $-5d^2+4$ $20d$ $-20$ 24 $-5d^2+4$ $20d$ $-20$

Table 2. Table 2: The coefficients R j c , k superscript subscript 𝑅 𝑗 𝑐 𝑘 R_{j}^{c,k} introduced in Eq. ( 37 ) for the matching of the annihilation amplitude.

$k$ 0 $j$ 1 2 3 4 5 $-1$ 6 $-1$ 7 $1$ 8 $1$ $k$ 0 1 9 $2$ 10 $2$ 11 $-2$ 12 $-2$ 13 $-d-2$ $2$ 14 $-d-2$ $2$ 15 $d+2$ $-2$ 16 $d+2$ $-2$ $k$ 0 1 2 17 $4d$ $-8$ 18 $4d$ $-8$ 19 $-4d$ $8$ 20 $-4d$ $8$ 21 $-d^2-8d+4$ $4d+16$ $-4$ 22 $-d^2-8d+4$ $4d+16$ $-4$ 23 $d^2+8d-4$ $-4d-16$ $4$ 24 $d^2+8d-4$ $-4d-16$ $4$

Equations

$${\mathcal{L}}_{\mathrm{NRQCD}}={\mathcal{L}}_g+{\mathcal{L}}_l+{\mathcal{L}}_{\psi }+{\mathcal{L}}_{\chi }+{\mathcal{L}}_{\psi \chi },$$

$${\mathcal{L}}_g=-\frac{1}{4}{G}^{\mu \nu a}{G}_{\mu \nu }^a+\frac{1}{4}\frac{c_1^g}{m^2}g{f}_{abc}{G}_{\mu \nu }^a{G}_{\alpha }^{\mu b}{G}^{\nu \alpha c},$$

$${\mathcal{L}}_l=\sum _{i=1}^{n_l}{\overline{q}}_{i}i/D{q}_i+\mathcal{O}\left(\frac{1}{m^2}\right),$$

$$\displaystyle{\cal L}_{\psi}=\psi^{\dagger}\Bigg{[}iD_{0}+\frac{c_{k}}{2m}\vec{D}^{2}+g_{s}\frac{c_{F}}{2m}\vec{\sigma}\cdot\vec{B}+g_{s}\frac{c_{D}}{8m^{2}}(\vec{D}\cdot\vec{E}-\vec{E}\cdot\vec{D})$$

$$\displaystyle\qquad\qquad+ig_{s}\frac{c_{S}}{8m^{2}}\vec{\sigma}\cdot(\vec{D}\times\vec{E}-\vec{E}\times\vec{D})+\mathcal{O}\left(\frac{1}{m^{3}}\right)\Bigg{]}\psi\,,$$

$${\mathcal{L}}_{\chi }=-{\mathcal{L}}_{\psi }\mathrm{with}\psi \to \chi ,i{D}^0\to -i{D}^0,E^i\to -E^i,$$

$${\mathcal{L}}_{\psi \chi }$$

$$d_{ss}$$

$$d_{sv}$$

$$d_{vs}$$

$$d_{vv}$$

$$\frac{{\alpha }_s^0}{{\alpha }_s(\mu )}{\left(\frac{{\mu }^2{e}^{{\gamma }_{\mathrm{E}}}}{4\pi }\right)}^{-ϵ}=Z_{{\alpha }_s}=1-\frac{{\beta }_0}{ϵ}\frac{{\alpha }_s(\mu )}{\pi }+\left(\frac{{\beta }_0^2}{{ϵ}^2}-\frac{{\beta }_1}{2ϵ}\right){\left(\frac{{\alpha }_s(\mu )}{\pi }\right)}^2+\mathcal{O}({\alpha }_s(\mu {)}^3),$$

$${\beta }_0=\frac{11}{12}{C}_A-\frac{1}{3}(n_l+n_h)T_F,{\beta }_1=\frac{17}{24}{C}_A^2-\left(\frac{5}{12}{C}_A+\frac{1}{4}{C}_F\right)(n_l+n_h)T_F,$$

$$T_F=\frac{1}{2},C_F=\frac{N_c^2-1}{2N_c},C_A=N_c.$$

$$\frac{{\alpha }_s^{(n_l+1)}(\mu )}{{\alpha }_s^{(n_l)}(\mu )}=1-\frac{{\alpha }_s^{(n_l)}(\mu )}{\pi }\frac{1-ϵI_0}{3ϵ}{T}_F+{\left(\frac{{\alpha }_s^{(n_l)}(\mu )}{\pi }\right)}^2{T}_F[T_F\frac{(1-ϵI_0{)}^2}{9{ϵ}^2}+C_A(-\frac{5}{24ϵ}$$

$$+\frac{ϵ\left(4{ϵ}^3+4{ϵ}^2-11ϵ-10\right)I_0^2}{8(ϵ-2)(2ϵ+1)(2ϵ+3)})+C_F\left(\frac{-ϵ\left(4{ϵ}^3-7ϵ-1\right)I_0^2}{4(ϵ-2)(2ϵ-1)(2ϵ+1)}-\frac{1}{8ϵ}\right)]+\mathcal{O}({\alpha }_s^3),$$

$$l_{\mu }$$

$$q_1(p)+{\bar{q}}_2(p)\to q_1(p)+{\bar{q}}_2(p)$$

$${\mathcal{M}}_{\mathrm{QCD}}^{scat.}=\sum _{j=1}^{24}\left(C_{\mathrm{s},j}\bar{u}{B}_j^{(1)}u\bar{v}{B}_j^{(2)}v+C_{\mathrm{o},j}\bar{u}{T}^a{B}_j^{(1)}u\bar{v}{T}^a{B}_j^{(2)}v\right),$$

$${\mathcal{M}}_{\mathrm{QCD}}^{anni.}=\sum _{j=1}^{24}\left(C_{\mathrm{s},j}^c\bar{v}{B}_j^{(1)}u\bar{u}{B}_j^{(2)}v+C_{\mathrm{o},j}^c\bar{v}{T}^a{B}_j^{(1)}u\bar{u}{T}^a{B}_j^{(2)}v\right).$$

$$\{{\gamma }^{\mu },{\gamma }^{\nu }\}=2g^{\mu \nu },g^{\mu }{}_{\mu }=d.$$

$$B_1^{(1)}\otimes B_1^{(2)}$$

$$B_2^{(1)}\otimes B_2^{(2)}$$

$$B_3^{(1)}\otimes B_3^{(2)}$$

$$B_4^{(1)}\otimes B_4^{(2)}$$

$$B_5^{(1)}\otimes B_5^{(2)}$$

$$B_6^{(1)}\otimes B_6^{(2)}$$

$$B_7^{(1)}\otimes B_7^{(2)}$$

$$B_8^{(1)}\otimes B_8^{(2)}$$

$$B_9^{(1)}\otimes B_9^{(2)}$$

$$B_{10}^{(1)}\otimes B_{10}^{(2)}$$

$$B_{11}^{(1)}\otimes B_{11}^{(2)}$$

$$B_{12}^{(1)}\otimes B_{12}^{(2)}$$

$$B_{13}^{(1)}\otimes B_{13}^{(2)}$$

$$B_{14}^{(1)}\otimes B_{14}^{(2)}$$

$$B_{15}^{(1)}\otimes B_{15}^{(2)}$$

$$B_{16}^{(1)}\otimes B_{16}^{(2)}$$

$$B_{17}^{(1)}\otimes B_{17}^{(2)}$$

$$B_{18}^{(1)}\otimes B_{18}^{(2)}$$

$$B_{19}^{(1)}\otimes B_{19}^{(2)}$$

$$B_{20}^{(1)}\otimes B_{20}^{(2)}$$

$$B_{21}^{(1)}\otimes B_{21}^{(2)}$$

$$B_{22}^{(1)}\otimes B_{22}^{(2)}$$

$$B_{23}^{(1)}\otimes B_{23}^{(2)}$$

$$B_{24}^{(1)}\otimes B_{24}^{(2)}$$

$$\displaystyle\gamma^{0}=\left(\begin{array}[]{rr}1&0\\ 0&-1\end{array}\right),\quad\vec{\gamma}=\left(\begin{array}[]{rr}0&\vec{\sigma}\\ -\vec{\sigma}&0\end{array}\right)$$

Full text

1 Introduction

Non-Relativistic Quantum Chromodynamics (NRQCD) [1] has proven to provide accurate predictions for systems of two heavy quarks, which move with a small relative velocity. Among them are decay rates and binding energies of quarkonia and the threshold production of top quark pairs in electron positron annihilation. For comprehensive compilations of results we refer to the review articles [2, 3, 4] and restrict ourselves here to recent next-to-next-to-next-to-leading order (N3LO) results. These include predictions for top quark pair production [5],111 In Ref. [6] next-to-next-to-leading logarithmic (NNLL) corrections have been obtained, see also [7]. the decay of the $\Upsilon(1S)$ meson [8], and energy levels of heavy quarkonia ground and excited states [9, 10, 11] together with phenomenological applications [12, 13].

Despite the high accuracy reached for a number of observables, it is desirable to extend the precision of the predictions. For example, the perturbative uncertainty of the N3LO top quark threshold prediction of about 3% will constitute the main uncertainty in the top quark mass value extracted from the comparison with future cross section measurements (see, e.g., Ref. [14]). Furthermore, the dominant source of uncertainty in the determination of the charm and bottom quark masses from bound state energies originates from the renormalization scale dependence, due to unknown higher order corrections [12, 11]. Currently a complete N4LO calculation is out of reach, note, however, that the completion of the ingredients necessary for the N3LO predictions took more than ten years and the combined effort of several groups (see, e.g., Ref. [4]). It is thus reasonable to proceed in a similar way at N4LO and gradually provide the individual building blocks required. In this work we compute two-loop matching coefficients which are building blocks of the NRQCD Lagrange density at N4LO.

A further and more short-term motivation of our work is the construction of logarithmically enhanced contributions which complement the N3LO predictions. The potential NRQCD (pNRQCD) Lagrange density relevant for $S$-wave states with next-to-next-to-next-to-leading logarithmic (N3LL) accuracy has been constructed in Ref. [15] up to a few missing contributions to the so-called soft running. Among them are the coefficients $d_{ss}$ and $d_{vs}$ (see the next section for a precise definition) which are computed in this work. Note that for $P$-wave states the N3LL pNRQCD Lagrange density is complete and can be found in Ref. [16].

The main purpose of this paper is the computation of the matching coefficients between QCD and NRQCD to two-loop order. We concentrate on the four-fermion operators but also compute the matching coefficients for gluon-quark interactions ($c_{D}$, $c_{F}$ and $c_{S}$) which are needed to obtain gauge invariant results. The corresponding one-loop results have been obtained in Refs. [17] and [18], respectively (see also Refs. [4]). The gauge dependence has its origin in the non-minimality of the operators entering the NRQCD Lagrange density. If fact, some of the effective operators can be absorbed into other operators by using the equation of motion or field redefinitions. The relevant equation of motion in our calculation is that which relates some of the four-fermion operators and the gluon-quark interaction [19] and thus only a particular combination is gauge invariant (see, e.g., Ref. [20]). In this paper, we perform our calculations in the general covariant gauge and present results for an arbitrary gauge parameter $\xi$. We check the cancellation of $\xi$ in the proper combination of the matching coefficients entering physical quantities. The computation of $d_{xy}$ requires a precise definition of the Pauli matrices in $d=4-2\epsilon$ dimensions, which we discuss in detail.

The calculation of the matching coefficients for four-fermion operators is naturally divided into two parts, which we call the annihilation and the scattering channel. The tree-level contribution of the former originates from the diagrams where a quark-anti-quark pair annihilates into a (virtual) gluon which subsequently “decays” into a quark-anti-quark pair (cf. Fig. 1). The corresponding one- and two-loop sample diagrams are shown in Figs. 1 and 3. In the case of the scattering channel one considers the scattering of a quark and an anti-quark, which may have different flavours and thus also different masses.

The remainder of the paper is organized as follows: In the next section we provide the relevant parts of the NRQCD Lagrange density and define the matching coefficients which we want to compute. In Section 3 we concentrate on the four-fermion matching coefficients and provide details of our two-loop calculation. Section 4 is devoted to the computation of the gluon fermion form factor and the extraction of the corresponding matching coefficients. The main results of the paper are presented in Section 5 where we provide analytic expressions for the four-fermion matching coefficients. In the appendix we provide additional material such as the matching coefficients needed for the redefinition of the gluon operators. Furthermore, analytic results for all two-loop master integrals are given in Appendix A.

2 ${\cal L}_{\rm NRQCD}$

The NRQCD Lagrange density to order $1/m^{2}$ which we use for our calculations is given by (see, e.g., Refs. [2, 4])

[TABLE]

where $i\vec{D}=i\vec{\nabla}+g_{s}\vec{A}$, $E^{j}=G^{j0}$, $B^{j}=-\varepsilon_{jkl}G^{kl}/2$, with $G^{ij}$ being the field strength tensor, and $n_{l}$ is the number of light quarks. In order to arrive at the canonical kinetic term of the gluon (2), one has to apply the field redefinition and the rescaling [21] (see also Appendix B). The main purpose of this work is the computation of the matching coefficients of ${\cal L}_{\psi\chi}$ (see below). However, in order to construct a gauge invariant combination we also need $c_{D}$, which we discuss in Section 4. Results for $c_{F}$ and $c_{S}$ are presented in Appendix C.

The interaction of four heavy quarks is given by

[TABLE]

where $\psi_{1}$ ($\psi_{2}$) are Pauli spinors annihilating a heavy quark with mass $m_{1}$ ($m_{2}$), and $\chi_{1}$ ($\chi_{2}$) are Pauli spinors creating a heavy anti-quark with mass $m_{1}$ ($m_{2}$). In this work we will identify the two masses and write $m=m_{1}=m_{2}$. We furthermore use the notation for the subscripts which is usually used in the literature: The first index in the matching coefficients $d_{xy}$ refers to the colour (“s” for singlet and “v” for octet) and the second denotes the singlet (“s”) and triplet (“v”) quark-anti-quark state.

The effective Lagrange density in Eq. (6) can be rewritten with the help of Fiertz transformations to arrive at

[TABLE]

which is better suited for the annihilation part of the matching calculation whereas we prefer version (6) for the scattering part. The relations between the coefficients in Eqs. (6) and (7) are given by [17]

[TABLE]

where $N_{c}=3$ corresponds to QCD. We compute the one- and two-loop four-quark amplitudes in Section 3 and provide results for $d_{xy}$ in Section 5.

Let us now describe the procedure which is used to obtain the NRQCD matching coefficients. We consider QCD with $n_{h}=1$ heavy quarks and $n_{l}$ light quarks, and compute the four quark scattering amplitudes (see Eqs. (15) and (16) below), the vertex corrections (see Eq. (39)), and the corrections to the matching coefficients in the gluon sector (see Eq. (60)). The ultra-violet (UV) renormalization is done in the $(n_{l}+n_{h})$-flavor theory. The relation between the bare coupling constant $\alpha_{s}^{0}$ and the $\overline{\mathrm{MS}}$ renormalized coupling constant $\alpha_{s}(\mu)$ reads

[TABLE]

where $\mu$ is the renormalization scale, and the colour factors for the $\mathrm{SU}(N_{c})$ gauge group are given by

[TABLE]

The heavy quark mass and wave function are renormalized on-shell. The renormalization constants are well known in the literature (see, e.g., Refs [22, 23]). We recompute them here in order to retain the exact $\epsilon$-dependence. Note that the wave function renormalization of the gluon is given by $1/\sqrt{Z_{\alpha_{s}}}$ because we use the background field method [24].

We first compute $F_{1}^{\prime}(0)$, $F_{2}(0)$ (see Section 4), and $d_{1}$, $d_{2}$ (see Appendix B). After UV renormalization, we convert the four-component Dirac spinors to the two-component Pauli spinors, and the Dirac matrices $\gamma^{\mu}$ to the Pauli matrices $\sigma^{j}$ assuming the non-relativistic limit. We then canonicalize the gluon sector (see Appendix B) and simultaneously decouple the heavy quark in the gluon wave function. Finally, we express $\alpha_{s}^{(n_{l}+n_{h})}(\mu)=\alpha_{s}^{(n_{l}+1)}(\mu)$ in terms of $\alpha_{s}^{(n_{l})}(\mu)$ by using the relation (for the bare version see Ref. [25])

[TABLE]

with $I_{0}=(\epsilon-1)I_{1}^{a}$, where $I_{1}^{a}$ is given in Eq. (59). Equation (12) is exact in $\epsilon$; $\epsilon$-expanded versions can be found in Refs. [26, 27]. In order to keep the expressions in this paper simple we provide the results in terms of $\alpha_{s}(m)$, which means that the renormalization scale $\mu$ is set to $m$. Using the renormalization group equations it is possible to reexpress $\alpha_{s}(m)$ by $\alpha_{s}(\mu)$. After expanding Eq. (12) in $\epsilon$ one obtains $\log\mu^{2}/m^{2}$ terms which we abbreviate by

[TABLE]

3 Four-fermion matching coefficients

In this section we describe the calculation of the full-QCD amplitudes which are needed for the matching coefficients $d_{xy}$ and $d_{xy}^{c}$ defined in Eqs. (6) and (7). They are obtained from the four-quark amplitude

[TABLE]

with the special kinematics indicated in the arguments of the quark fields $q_{1}$ and $q_{2}$. Sample Feynman diagrams, which one has to consider at one- and two-loop order, are shown in Fig. 1. In general one can sub-divide them into “annihilation” (top row) and “scattering” contributions (bottom row). Note that in the case that the two heavy quarks have different flavours (and thus also different masses) only scattering diagrams contribute whereas in the equal-mass case also the annihilation diagrams are needed. In this paper we consider only the limit that both quarks have equal masses. Nevertheless we discuss the two contributions separately.

3.1 Matching

Let us in the following briefly describe the individual steps which are necessary to perform the matching between QCD and NRQCD. The general idea is to consider the four-fermion amplitude in QCD in the limit of a heavy quark mass and compare to the corresponding expression in NRQCD, which provides results for $d_{xy}$ and $d_{xy}^{c}$.

We start with the QCD amplitudes which for the scattering and annihilation channel have the form

[TABLE]

where $u$ ($v$) is the quark (anti-quark) spinor and $2T^{a}$ are the Gell-Mann matrices. The superscript “c” in Eq. (16) denotes that the result is matched to the Lagrange density (7), whereas in the scattering channel we match our expressions to Eq. (6). The coefficients $C_{\mathrm{s/o},j}$ and $C_{\mathrm{s/o},j}^{c}$, where “s” and “o” refer to singlet and octet colour states, are determined by an explicit calculation of the amplitude in Eq. (14). In calculating the QCD amplitude, we treat the $\gamma$ matrices as $d$-dimensional objects which satisfy

[TABLE]

Unlike the case of 4-dimensional $\gamma$ matrices, products of more than four $d$-dimensional $\gamma$ matrices can not be expressed in terms of simpler products of $\gamma$ matrices, and we have to treat all such products as independent basis elements. Taking into account this fact, we consider the following basis elements222Note that $B_{22}^{(1)}\otimes B_{22}^{(2)}$, $B_{23}^{(1)}\otimes B_{23}^{(2)}$ and $B_{24}^{(1)}\otimes B_{24}^{(2)}$ do not enter our calculation since, up to two-loop order, at most five $\gamma$ matrices are present in one fermion line. Nevertheless, for symmetry reasons, we provide also these basis elements.

[TABLE]

where ${/\!\!\!v}={/\!\!\!p}/m$ and the superscript refers to the fermion line. We have explicitly introduced the external momentum $p$ since we do not use the Dirac equation in the course of the computation of the Feynman diagrams.

In matching to the NRQCD amplitude, we use the following representation of the $\gamma$ matrices

[TABLE]

in terms of $(d-1)$-dimensional Pauli matrices which satisfy

[TABLE]

In particular, we do not use the commutation relation of the Pauli matrices at this point.

The NRQCD amplitudes for the scattering and annihilation channels can be written as

[TABLE]

where $\phi$ and $\eta$ are two-component spinors which in the limit of vanishing 3-momentum are related to the $u$ and $v$ spinors in full QCD via

[TABLE]

The factor $\sqrt{2m}$ for each external quark appears due to our convention for the normalization of the non-relativistic quark fields [4]. Note that in Eqs. (25) and (26) different bases have been introduced for the scattering and annihilation channels (see also Eqs. (6) and (7)). In $d=4-2\epsilon$ dimensions the basis elements are related to the Pauli matrices as

[TABLE]

For the two-loop calculation of $d_{xy}$ and $d^{c}_{xy}$ only $\Sigma_{i}$ and $\Sigma_{i}^{c}$ with $i=0,1,2$ are needed. At three loops basis elements constructed from products of more than five Pauli matrices are necessary.

In order to obtain the matching coefficients in Eqs. (6) and (7), one has to reduce the structure of the Pauli matrices to $1\!\!1\otimes 1\!\!1$ and $\sigma^{j}\otimes\sigma^{j}$ instead of those in Eqs. (32). In other words, one has to take the limit $d\to 4$. There are different prescriptions to do this; one can use the commutation relation $[\sigma_{j},\sigma_{k}]=2i\varepsilon_{jkl}\sigma_{l}$ assuming $\varepsilon^{jkl}\varepsilon^{jkl^{\prime}}=(d-2)\delta^{ll^{\prime}}$ [17], or $\varepsilon^{jkl}\varepsilon^{jkl^{\prime}}=2\delta^{ll^{\prime}}$. Since it is unclear which prescription should be used, we provide the $d$-dimensional results in the basis of Eqs. (32). Nevertheless, it is useful to have the conventional matching coefficients $d_{xy}$. For this purpose we adopt $\varepsilon^{jkl}\varepsilon^{jkl^{\prime}}=2\delta^{ll^{\prime}}$ and obtain

[TABLE]

In the following, we refer to this prescription as “taking the limit $d\to 4$”.

At this point it is convenient to discuss the scattering and annihilation channel separately. In the former case one has to consider $\gamma^{\mu_{1}}\cdots\gamma^{\mu_{n}}$ sandwiched between $\bar{u}$ and $u$ or $\bar{v}$ and $v$, which means that only diagonal parts of $\gamma^{\mu_{1}}\cdots\gamma^{\mu_{n}}$ contribute. Then we obtain

[TABLE]

where the $R_{j}^{k}$ are given in Tab. 1. In order to obtain the table entries one can use the equation of motion for the external fermions

[TABLE]

Afterwards, we insert the explicit expressions for the spinors $u$ and $v$ in terms of $\phi$ and $\eta$ (cf. Eq. (31)). After substituting Eq. (34) into Eq. (15) and comparing with Eq. (25), we obtain the relations between NRQCD coefficients $c_{\mathrm{s/o},k}$ and QCD coefficients $C_{\mathrm{s/o},j}$:

[TABLE]

In the case of the annihilation channel $\gamma^{\mu_{1}}\cdots\gamma^{\mu_{n}}$ is sandwiched between $\bar{v}$ and $u$ or $\bar{u}$ and $v$ and thus only the off-diagonal parts contribute, which means that one needs an odd number of $\vec{\gamma}$ matrices. In analogy to Eq. (34) we can write

[TABLE]

where $R_{j}^{c,k}$ are given in Tab. 2. Substituting Eq. (37) into Eq. (16) and comparing with Eq. (26) leads to the relations between NRQCD coefficients $c^{c}_{\mathrm{s/o},k}$ and QCD coefficients $C^{c}_{\mathrm{s/o},j}$:

[TABLE]

Results up to two loops for $c_{\mathrm{s/o},k}$ and $c^{c}_{\mathrm{s/o},k}$ are presented in Section 5.

3.2 Loop integrals

In the following we briefly describe the workflow of our calculation. We first generate the full QCD amplitudes with qgraf [28] and map the output to general four-point families which have four and nine independent propagators at one and two loops, respectively. Next, we apply projectors to obtain the coefficients of the basis elements $B_{i}$ which leads us to scalar expressions. Afterwards, we specify the kinematics given in Eq. (14). At two loops this leads to five (instead of nine) linearly independent propagators. One has to apply a partial fraction decomposition in order to obtain integral families which can be reduced to master integrals using FIRE [29] and LiteRed [30].

In an alternative approach, which we use for some of the integral families, we specify only some of the kinematic relations such that the propagators are still linearly independent. Then we perform an integration-by-parts reduction, apply the full kinematic information of Eq. (14) to the resulting master integrals, perform a partial fraction decomposition to these masters, and a further (very simple) reduction in order to arrive at the same set of master integrals as in our standard approach. Note that in all cases the reduction problem is quite simple and takes at most, even for general QCD gauge parameter, a few minutes on a desktop computer.

Our final result for the QCD amplitude can be expressed in terms of two one-loop and ten two-loop master integrals (cf. Fig. 2). We retain the exact $\epsilon$-dependence up to this point and provide the corresponding results in an ancillary file [31]. Most of the master integrals are available in the literature [32, 33, 34]. However, not all of them are known analytically, and for some higher orders in $\epsilon$ are needed. Furthermore, to our knowledge the box-type integral $I_{2}^{g}$ is not available in the literature so far. For this reason we (re)compute those integrals analytically and present the results in Appendix A.

After inserting the master integrals into the four-fermion amplitudes we use Eqs. (36) and (38), expand in $\epsilon$ and thus obtain the matching coefficients $c_{s/o,k}$ and $c^{c}_{s/o,k}$. Analytic results are presented in Section 5. Let us mention that the colour and Lorentz part of the QCD amplitude factorizes such that they can be computed independently.

4 Gluon fermion matching coefficients

The purpose of this section is the computation of $c_{D}$ which has to be combined with $d_{vs}$ in order to cancel the $\xi$ dependence. Since the calculation of $c_{F}$ and $c_{S}$ proceeds among similar lines we compute all three matching coefficients simultaneously and present results up to two loops.

The matching coefficients $c_{D}$, $c_{F}$ and $c_{S}$ can be extracted from the gluon-quark vertex function which we parameterize as

[TABLE]

where $p$ ($p^{\prime}$) is the outgoing (incoming) quark (anti-quark) momentum and $q=p-p^{\prime}$. The quark momenta are on-shell, i.e. $p^{2}=(p^{\prime})^{2}=m^{2}$ and we have $\sigma^{\mu\nu}=i[\gamma^{\mu},\gamma^{\nu}]/2$. The fundamental indices in the matrix $T^{a}$ are suppressed.

The calculation is performed in the background field method [24] where the gauge parameter $\xi$ enters via the gluon propagator

[TABLE]

and the vertex of the background gluon and two quantum gluons, which contains a factor $1/(1-\xi)$. Note that the $\xi$-dependence is treated exactly throughout the calculation.

For the matching calculation it is sufficient to consider $\Gamma^{a}_{\mu}$ in the limit of small gluon momentum $q$. In fact, after considering the non-relativistic limit in Eq. (39) the comparison to the tree-level Feynman rules from ${\cal L}_{\psi}$ in Eq. (4) leads to

[TABLE]

where the prime indicates the derivative w.r.t. the argument and $d_{1},d_{2}$ can be found in Appendix B. The tilde in Eq. (41) indicates that no rescaling of the gluon field has been performed. Thus, in order to obtain the matching coefficients present in the Lagrange density (4) one has to apply Eq. (64) in Appendix B. Note that $d_{1}=1+\mathcal{O}(\alpha_{s})$ and $d_{2}=\mathcal{O}(\alpha_{s})$, and thus $d_{2}/d_{1}\to d_{2}$ at one-loop order. We can Taylor-expand the form factors $F_{1}$ and $F_{2}$ in the gluon momentum and are left with one- and two-loop on-shell integrals which are well studied in the literature (see, e.g., Refs. [35, 36]).

In the following we provide results for the form factors and their derivatives for $q^{2}=0$. We parametrize the form factors as

[TABLE]

Note that the $F_{i}$ still contain poles and also have an explicit $\mu$ dependence. Below we show the $\epsilon$-expanded expressions and provide the $\epsilon$-exact results in an ancillary file [31]. Our results for $F_{i}^{\prime}(0)$ and $F_{2}(0)$ read

[TABLE]

Our two-loop result for $F_{2}(0)$ agrees with Refs. [37, 27] and the QED part333The QED result is obtained for $C_{A}=0$, $C_{F}=1$, $T_{F}=1$, $n_{l}=0$, $n_{h}=1$ and the coupling constant renormalized in the on-shell scheme. of $F^{\prime}_{1}(0)$ can be found in [38, 39]. The two-loop QCD corrections to $F^{\prime}_{1}(0)$ are new.

We can now use Eq. (41), apply the rescaling of Eq. (64) and decouple the heavy quark in the gluon wave function and the coupling constant444Note that we apply the decoupling also to the factor $g_{s}$ in Eq. (4). in order to compute $c_{D}$, $c_{F}$ and $c_{S}$. In the following we present one- and two-loop expressions for $c_{D}$ and postpone $c_{F}$ and $c_{S}$ to Appendix C. By parameterizing the matching coefficients $c_{X}$ as

[TABLE]

we obtain for $c_{D}$

[TABLE]

Note the $\xi$ dependence in the second last line which is inherited from $F_{1}^{\prime(2)}(0)$ and $d_{2}$ according to Eq. (41).

5 Results for the four-fermion matching coefficients

In this section we present first our results in $d$ dimensions and afterwards take the limit $d\to 4$. We discuss both the scattering and the annihilation channel.

5.1 NRQCD four quark coefficients in $d$ dimensions

We parametrize the matching coefficients as follows

[TABLE]

and use an analogous expansion for $c_{\mathrm{s/o},k}^{c}$. At tree level we have

[TABLE]

and all the other coefficients are zero. We have obtained exact results in $d$ dimensions both at one and two loops and provide the corresponding results in an ancillary file [31]. Below we show the $\epsilon$-expanded expressions.

5.1.1 One-loop results

Our one-loop results for the scattering channel are given by

[TABLE]

Note that $c_{s,2}^{(1)}=0$ and $c_{o,2}^{(1)}=0$ since at one-loop order at most two $\sigma$ matrices are present in a spinor line. In the literature, the factor $1/N_{c}$ in the colour singlet matching coefficients are expressed as $(C_{A}-2C_{F})$. Here and in the following, we use $1/N_{c}$ in order to have more compact expressions.

For the annihilation channel we have

[TABLE]

where we have again $c_{s,2}^{c,(1)}=0$ and $c_{o,2}^{c,(1)}=0$.

5.1.2 Two-loop results

At two-loop order the matching coefficients obtained form the scattering process read

[TABLE]

All six coefficients are new and not yet present in the literature. This is also true for the following six matching coefficients obtained from the annihilation-type diagrams

[TABLE]

Note that for the annihilation channel, products of two one-loop diagrams also have to be taken into account. Furthermore, two-loop vertex corrections as shown in Fig. 3(a) contribute to the colour-octet vector current. After adapting the colour factors, we have cross-checked these contributions against the explicit results provided Ref. [33].

In the next subsection we use the results presented above in order to obtain the four-quark matching coefficients present in ${\cal L}_{\rm NRQCD}$.

5.2 NRQCD four quark coefficients in four dimensions

In the following we use the expressions from the previous subsection and apply $[\sigma_{i},\sigma_{j}]=2i\varepsilon_{ijk}\sigma_{k}$ and $\varepsilon^{jkl}\varepsilon^{jkl^{\prime}}=2\delta^{ll^{\prime}}$. Using Eq. (33) one obtains the following linear combinations of $c_{s/o,k}$ which provide the matching coefficients present in the NQRCD Lagrange density of Eq. (6):

[TABLE]

Note that at one-loop order we have $c_{s/o,2}^{(1)}=0$ and thus the relations are trivial. The $\epsilon$-exact one-loop expressions agree with Ref. [4]. Note that in Ref. [17] a different prescription for $\varepsilon^{ijk}$ in three dimensions has been used (cf. discussion between Eqs. (32) and (33)) which leads to different relations compared to those in Eq. (52).555At one-loop order one has $d_{xv}^{\mbox{\tiny\cite[cite]{[\@@bibref{}{Pineda:1998kj}{}{}]}}}=(1-\epsilon)d_{xv}^{\mbox{\tiny Eq.(\ref{eq::c2d})}}$ with $x\in\{s,v\}$.

By denoting the loop corrections as

[TABLE]

the two-loop scattering coefficients are given by

[TABLE]

The relations between $c^{c}_{s/o,k}$ and $d^{c}_{xy}$ are also obtained from Eq. (33) and are given by

[TABLE]

At tree-level, $c^{c}_{s/o,1}=c^{c}_{s/o,2}=0$ and the relations are trivial.

We define the coefficients $d_{xy}^{c,(j)}$ in analogy to Eq. (53) and obtain for the one-loop annihilation matching coefficients

[TABLE]

The $\epsilon$-exact expressions agree with Ref. [4] and the expanded expressions with Ref. [17]. The two-loop annihilation matching coefficients read

[TABLE]

Note that all two-loop coefficients are $\xi$ independent except $d_{vs}^{(2)}$. In fact, the gauge parameter dependence cancels in the combination $(\alpha_{s}/\pi)c_{D}^{(2)}+d_{vs}^{(2)}$ which enters physical quantities.

The imaginary parts of $d_{ss}^{c,(2)}$, $d_{vs}^{c,(2)}$, and $d_{vv}^{c,(2)}$ are calculated in the context of the heavy quarkonium inclusive decays [40], and our results agree with the literature.

All the matching coefficients from the annihilation process are finite after the UV renormalization except $d_{vv}^{c,(2)}$. The remaining divergences originate from diagrams shown in Fig. 3(a) and (b). They are well studied in the literature [41] where it is shown that the divergences from the purely hard regions, which are contained in our expressions, are canceled against contributions from the potential region. We have confirmed this cancellation for the contribution from Fig. 3(b) where explicit results for the different regions are given in Ref. [41].

6 Conclusions and outlook

In this paper we compute two-loop corrections to the matching coefficients $d_{ss}$, $d_{sv}$, $d_{vs}$, $d_{vv}$, $d^{c}_{ss}$, $d^{c}_{sv}$, $d^{c}_{vs}$ and $d^{c}_{vv}$ of the operators in the NRQCD Lagrange density involving four heavy quarks. We carefully discuss the treatment of the Pauli matrices in a non-integer number of dimensions which leads to an enlargement of the basis and six (instead of four) two-loop coefficients in intermediate steps (see Section 5.1). The results for $d_{xy}$ and $d^{c}_{xy}$, which are obtained after using the usual commutation relations between the Pauli matrices, are given in Section 5.2.

Our calculation is performed in the covariant $R_{\xi}$ gauge with a general gauge parameter $\xi$. One observes that starting from two loops the coefficient $d_{vs}$ is $\xi$ dependent which arises from our non-minimal choice of the operator basis in $\mathcal{L}_{\text{NRQCD}}$. We check the $\xi$ dependence by computing two-loop corrections to the heavy-quark-gluon vertex functions. We extract the related matching coefficients, in particular $c_{D}$, and show that the combination $(\alpha_{s}/\pi)c_{D}^{(2)}+d_{vs}^{(2)}$ is independent of $\xi$. Note that in Feynman gauge the one-loop results $c_{D}^{(1)}$ and $d_{vs}^{(1)}$ are individually $\xi$ independent. However, the gauge dependence can be observed by comparing to the results in Coulomb gauge [20].

The results obtained in this paper enter as building blocks various physical quantities involving two slowly moving heavy quarks at the N3LL and N4LO accuracy.

The annihilation channel only contributes to the case where the two heavy quarks in ${\cal L}_{\phi\chi}$ (cf. Eq. (6)) have the same flavour. On the other hand, for different quark flavours the matching coefficients $d_{xy}$ receive contributions only from the scattering channel. We use the same mass for quarks and anti-quarks and provide only results for this equal-mass case. A possible next step would thus be the extension of our calculation of the scattering contribution to the case of different quark masses. A further next step is the computation of two-loop corrections to the matching coefficient of the operator with two heavy and two light quarks usually denoted by $c_{1}^{hl}$ (see, e.g., Ref. [15]).

Acknowledgements

We would like to thank Alexander Penin, Jan Piclum and Antonio Pineda for many useful discussions and communications. We thank Florian Herren for technical help in connection to the partial fraction decomposition and Joshua Davies for carefully reading the manuscript. This research was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under grant 396021762 — TRR 257 “Particle Physics Phenomenology after the Higgs Discovery” and the Graduiertenkolleg “Elementarteilchenphysik bei höchster Energie und höchster Präzision”.

Appendix A Master integrals

In this appendix we collect analytic results for the master integrals which we need for the computation of the matching coefficient. Most of them are already needed for two-loop matching coefficients between QCD and NRQCD of the vector, axial-vector, scalar and pseudo-scalar currents [32, 33, 42] and the integrals have been studied in the literature [34] (see also Refs. [43, 44]). Note, however, that for $I_{2}^{e}$ the $\epsilon$ expansion was not sufficiently deep and the $\epsilon^{0}$ was only known numerically. Furthermore $I_{2}^{g}$ was (to our knowledge) not available in the literature.

The master integrals are defined as (cf. Fig. 2)

[TABLE]

where $\mathcal{N}=(\mu^{2}e^{\gamma_{\mathrm{E}}})^{\epsilon}$. We normalize the master integrals such that they have the mass dimension zero. Our results read

[TABLE]

For the integral $I_{2}^{e}$ we derive a Mellin-Barnes representation with non-zero parameter $\epsilon$ and use MB.m [45] to analytically continue to $\epsilon\to 0$. The resulting (at most) two-dimensional Mellin-Barnes integrals are reduced to one-dimensional Mellin-Barnes integrals with a help of the generalized Barnes lemma [46, 47]. The one-dimensional integrals can be evaluated numerically with a very high precision, and we apply the PSLQ algorithm [48] to obtain the analytic results.

Using the Mellin-Barnes method for $I_{2}^{g}$ leads to a complicated four-dimensional Mellin-Barnes integral, and we adopt a different strategy for its computation. Note that $I_{2}^{g}$ is a finite integral and we require only the $\epsilon^{0}$ term. This means we can set $\epsilon=0$ from the very beginning of our computation. We use the Lee-Pomeransky representation [49] which turns out to be useful since the integrand is now a simple rational function. We can perform most of the integrations analytically and remain only with a two-dimensional integral with good covergence properties. Thus, numerical integration leads to sufficiently high precision such that the PSLQ algorithm can be applied. We cross-check all master integrals with the help of FIESTA [50].

Appendix B Gluon field redefinition

In Ref. [3] the NRQCD Lagrange density has been defined such that the kinetic term of the gluon field has a canonical normalization which has been achieved by a redefinition of the gluon field. The procedure is presented in Ref. [21]. As a consequence the constants $d_{1}$ and $d_{2}$ appear on the right-hand side of the formula for $c_{D}$ in Eq. (41). In this section we provide analytic expressions for $d_{1}$ and $d_{2}$ up to two-loop order.

Our starting point is the following Lagrange density which describes the interaction of the heavy quarks with a gluon before the redefinition of the gluon field

[TABLE]

where $G_{\mu\nu}$ is the gluon field strength tensor and $a,b,c$ are colour indices. The matching coefficients $d_{1}$ and $d_{2}$ can be computed from the hard contribution of the gluon two-point function. For convenience we provide the results which we parametrize by

[TABLE]

and $d_{1}^{(0)}=1$, $d_{2}^{(0)}=0$. Up to two-loop order our results read

[TABLE]

Note that the external gluon fields have been renormalized in the $\overline{\rm MS}$ scheme.

It is common practice to perform a redefinition of the gluon field as

[TABLE]

which eliminates the second term in Eq. (60). A subsequent rescaling of the form

[TABLE]

leads to the canonical factor “$-1/4$” in the first term of Eq. (60).

Appendix C Results for $c_{F}$ and $c_{S}$

In this appendix we provide analytic results for $c_{F}$ and $c_{S}$ up to two loops. Our results read

[TABLE]

The one-loop results agree with Refs. [18] and [4]; the two-loop results are new.