Three loop QCD corrections to heavy quark form factors

TL;DR

This paper introduces an algorithm for solving coupled differential equations in perturbative QCD, enabling the calculation of three-loop heavy quark form factors with improved efficiency and precision.

Contribution

The paper presents a novel algorithm for solving coupled linear differential systems in quantum field theory, applied here to compute three-loop heavy quark form factors.

Findings

Successfully computed master integrals for three-loop form factors

Extended the method to include light quark contributions

Enhanced the accuracy of perturbative QCD calculations

Abstract

Higher order calculations in perturbative Quantum Field Theories often produce coupled linear systems of differential equations which factorize to first order. Here we present an algorithm to solve such systems in terms of iterated integrals over an alphabet the structure of which is implied by the coefficient matrix of the given system. We apply this method to calculate the master integrals in the color-planar and complete light quark contributions to the three-loop massive form factors.

Figures

Equations

$$s=\frac{q^2}{m^2}=-\frac{(1-x{)}^2}{x}.$$

$$\frac{d}{dx}\mathcal{I}=\mathcal{M}\mathcal{I}+\mathcal{R}.$$

$${\mathcal{R}}_i=\sum _{j=-k}^{\infty }{\epsilon }^j{\mathcal{R}}_i^{(j)}.$$

$${\mathcal{I}}_i=\sum _{j=-k}^{\infty }{\epsilon }^j{\mathcal{I}}_i^{(j)}.$$

$$\frac{d}{dx}\tilde{\mathcal{I}}=\tilde{\mathcal{M}}\tilde{\mathcal{I}}+\tilde{\mathcal{R}},$$

$$\frac{d}{dx}{\cal\tilde{I}}^{(k)}={\tilde{\cal}M^{(0)}~{}\tilde{\cal}I^{(k)}+\Big{(}\tilde{\cal}M^{(1)}~{}\tilde{\cal}I^{(k-1)}+\tilde{\cal}M^{(2)}~{}\tilde{\cal}I^{(k-2)}+\cdots+\tilde{\cal}M^{(k+l)}~{}\tilde{\cal}I^{(-l)}\Big{)}+\tilde{\cal R}^{(k)}}\,.$$

$$\frac{d}{dx}{\tilde{\mathcal{I}}}^{(-l)}=\stackrel{~}{}{M}^{(0)}\stackrel{~}{}{I}^{(-l)}+{\tilde{\mathcal{R}}}^{(-l)}.$$

$$\sum _{k=0}^m{p}_k(x)\frac{d^k}{d{x}^k}{\tilde{\mathcal{I}}}_1^{(-l)}(x)=r(x).$$

$${\tilde{\mathcal{I}}}_k^{(-l)}(x)=\sum _{i=0}^{m-1}{a}_{k,i}(x)\frac{d^i}{d{x}^i}{\tilde{\mathcal{I}}}_1^{(-l)}(x)+{\rho }_k(x),$$

$$\left(\frac{d}{dx}-{\hat{p}}_1(x)\right)\left(\frac{d}{dx}-{\hat{p}}_2(x)\right)\dots \left(\frac{d}{dx}-{\hat{p}}_m(x)\right)y_1(x)=0,$$

$$\mathfrak{A}=\left\{\frac{1}{x},\frac{1}{1-x},\frac{1}{1+x},\frac{1}{1+x^2},\frac{x}{1+x^2},\frac{1}{1+x+x^2},\frac{x}{1+x+x^2},\frac{1}{1-x+x^2},\frac{x}{1-x+x^2}\right\}.$$

$$g(x)=\sum _{i=1}^m{y}_i(x){\int }_l^{x}d\tilde{x}\frac{r(\tilde{x})W_i(\tilde{x})}{W(\tilde{x})}$$

$$W_{i}(x)=(-1)^{i+m}\left|\begin{array}[]{cccccc}y_{1}&\ldots&y_{i-1}&y_{i+1}&\ldots&y_{m}\\ \frac{d}{dx}y_{1}&\ldots&\frac{d}{dx}y_{i-1}&\frac{d}{dx}y_{i+1}&\ldots&\frac{d}{dx}y_{m}\\ \vdots&&\vdots&\vdots&&\vdots\\ \frac{d^{m-2}}{dx^{m-2}}y_{1}&\ldots&\frac{d^{m-2}}{dx^{m-2}}y_{i-1}&\frac{d^{m-2}}{dx^{m-2}}y_{i+1}&\ldots&\frac{d^{m-2}}{dx^{m-2}}y_{m}\\ \end{array}\right|\,.$$

$${\bar{u}}_c(q_1){\Gamma }_{V,cd}^{\mu }{v}_d(q_2)$$

$${\bar{u}}_c(q_1){\Gamma }_{A,cd}^{\mu }{v}_d(q_2)$$

$${\bar{u}}_c(q_1){\Gamma }_{S,cd}{v}_d(q_2)$$

$${\bar{u}}_c(q_1){\Gamma }_{P,cd}{v}_d(q_2)$$

$$F_I=\sum _{n=0}^{\infty }{\left(\frac{{\alpha }_s}{4\pi }\right)}^n{F}_I^{(n)}.$$

$$q^{\mu }{\Gamma }_{A,cd}^{\mu ,\mathsf{ns}}=2m{\Gamma }_{P,cd}^{\mathsf{ns}},$$

$$F_I=Z(\mu )F_I^{\mathrm{fin}}(\mu ),$$

Full text

1 Introduction

Higher order calculations in perturbative Quantum Field Theories often produce coupled linear systems of differential equations which factorize to first order. Through the integration-by-parts identities (IBP) [1, 2], the Feynman integrals appearing in a scattering amplitude can be mapped to a smaller set of integrals, the master integrals (MIs). One method to solve these integrals is the method of differential equations [3, 4, 5, 6], where one differentiates the MIs with respect to a kinematic invariant and obtains such a coupled linear system of differential equations in which we are interested in. Such an algorithm was constructed earlier in [6, 7] in the case of a system of difference equations. Here we summarize the algorithm, presented in [8], operating on uni-variate systems of differential equations, which factorize to first order. In the cases where the systems do not factorize to first order, elliptic and even more involved structures appear, cf. e.g. [9, 10, 11, 12, 13, 14, 15, 16, 17]. Our algorithm works for any basis of master integrals, while in the one in [5] a special basis is required.

As an example, we employ this method to compute the set of MIs which appear in the color–planar and complete light quark non–singlet three-loop contributions to the heavy-quark form factors for different currents, namely the vector, axial-vector, scalar and pseudo–scalar currents. These form factors are of phenomenological importance. Being the heaviest particle of the Standard Model (SM), the top quark plays an important role in comprehending the electro-weak symmetry breaking (EWSB). Consequently, many observables like the inclusive production cross section of a pair of top quarks, the forward–backward asymmetry, etc., draw attention in the post era of the Higgs boson discovery. These form factors are basic building blocks of many such observables. The two-loop QCD contributions to these massive form factors were first computed in [18, 19, 20, 21]. Later in [22], an independent computation was performed for the vector form factors, including the ${\cal O}(\varepsilon)$ terms, where $\varepsilon$ is the dimensional regularization parameter in $D=4-2\varepsilon$ space-time dimensions.

In [23], the two-loop QCD contributions up to ${\cal O}(\varepsilon^{2})$ for all massive form factors were calculated. At the three-loop level, the color-–planar contributions to the vector form factors have been computed in [25, 24] and the complete light quark contributions in [26]. Using the method described here, we have calculated both the color–-planar and complete light quark contributions to the three–loop form factors for the axial–vector, scalar and pseudo-scalar currents in [27] and for the vector current in [8], where a detailed description of the present method along with an example is presented. In a parallel calculation, the same results have been obtained in [28]. The asymptotic behavior of these three–loop form factors has been studied in [29, 30].

2 Description of the method

We consider a set of $n$ MIs ${\cal I}=(I_{1},\ldots,I_{n})$ within the same topology. The MIs are functions of the space–time dimension $D$ and the Landau variable $x$ defined by

[TABLE]

Here $q^{2}$ is the virtuality of the boson and $m$ denotes the heavy quark mass. One obtains an $n\times n$ system of coupled linear differential equations by performing the derivative with respect to $x$ of each of the MIs followed by the IBP reduction,

[TABLE]

Here the $n\times n$ matrix ${\cal M}$ consists out of entries from the rational function field ${\mathbb{K}}(D,x)$ (or equivalently from ${\mathbb{K}}(\varepsilon,x)$) where ${\mathbb{K}}$ is a field of characteristic [math]. The inhomogeneous part ${\cal R}=({\cal R}_{1},\dots,{\cal R}_{n})$ contains MIs which are already known. In simple cases, ${\cal R}$ turns out to be just the null vector. For more involved cases, we assume that each entry ${\cal R}_{i}$ is expanded111In the following $f^{(k)}$ does not denote the $k$th derivative of $f$. into a Laurent series in $\varepsilon$ in terms of iterated integrals

[TABLE]

To proceed, we assume that the unknown integrals ${\cal I}_{i}$ can be expanded into a Laurent series in $\varepsilon$

[TABLE]

We first apply Zürcher’s algorithm [31, 33, 32, 34] (or a variant of it), implemented in the package OreSys [35], leading to a single differential equation which is then analyzed with HarmonicSums [37, 38, 39, 41, 42, 36, 40, 43, 44] for first order factorization.

The MIs can be distinguished sector–wise, where the maximal set of non–vanishing propagators in a single Feynman graph defines a sector and corresponding subsets define the sub-sectors. On the other hand, a derivative can only introduce the inverse of a new propagator. Hence, the differential equation of a MI can contain integrals from the same sector or its sub–sectors. Thus, organizing the MIs in such a way that MIs with less number of propagators are kept at the bottom of the list, provides an upper–block–triangular form of ${\cal M}$, i.e. the diagonal elements of ${\cal M}$ are square matrices of rank one or more. Each such square matrix represents a completely coupled set of MIs and we call them sub–systems of ${\cal M}$. Due to such an arrangement of the system, we can now solve it in a bottom–up approach, i.e. we first solve the last coupled sub–system, having no dependence to the other MIs and work upwards in the next step. Below we present how we solve the sub–systems.

I. Let us consider $m$ integrals ${\cal\tilde{I}}=(\tilde{I}_{1},\ldots,\tilde{I}_{m})$ which constitute a coupled sub–system,

[TABLE]

where ${\cal\tilde{M}}$ and ${\cal\tilde{R}}$ are respective parts of ${\cal M}$ and ${\cal R}$ and hence have similar expansions in $\varepsilon$. Now we exploit the fact that for a certain topology and kinematics an integral has a definitive pole structure with a fixed order ($k$) of the highest pole, as e.g. the integrals arising in three–loop heavy quark form factors can have at most a $1/\varepsilon^{3}$ pole. Keeping that in mind, we plug in Eq. (4) into Eq. (5), perform the series expansion in $\varepsilon$ and obtain the coefficient of $\varepsilon^{k}$ as follows

[TABLE]

II. Now we try to solve the sub–system order by order in $\varepsilon$. To accomplish that we start with the coefficient of the leading pole $\varepsilon^{-l}$. The corresponding sub–system is

[TABLE]

To solve Eq. (7), a natural first step is to reduce this $m\times m$ system to a higher order linear differential equation. We refer to this procedure as ‘uncoupling’. By using OreSys, we obtain

[TABLE]

Here $p_{l}(x)$ are rational functions in ${\mathbb{K}}(x)$ and $r(x)$ consists of contributions from ${\cal R}^{(-l)}_{j}(x)$ and its derivatives. OreSys also provides the additional $m-1$ relations for the remaining integrals ${\cal\tilde{I}}^{(-l)}_{k}(x)$, $k=2,\ldots,m$ in terms of linear combinations of ${\cal\tilde{I}}^{(-l)}_{1}(x)$ and its derivatives:

[TABLE]

with $a_{k,i}(x)\in{\mathbb{K}}(x)$ and $\rho_{k}(x)\in{\mathbb{K}}(x)$. We remark that the uncoupling may find several linear differential equations for several unknowns. However, the sum of the order of all the linear differential equations is $m$. For simplicity we assume in the following that only one linear differential equation for ${\cal\tilde{I}}^{(-l)}_{1}(x)$ is produced.

III. Now we consider the homogeneous solutions of Eq. (8). The differential equation can be factorized at first order as

[TABLE]

with $\hat{p}_{k}$ being rational functions in ${\mathbb{K}}(x)$ by using algorithms from [45, 46, 47]. We find the d’Alembertian solutions $y_{i}(x)$, $i=1,\ldots,m$ (see [8] for details) for the homogeneous part of Eq. (8) in terms of iterative integrals. In our application to the three–loop massive form factors, the alphabet is

[TABLE]

The arising iterative integrals in our application can be simplified to the harmonic polylogarithms (HPLs) [48] and cyclotomic HPLs [36]. We also note here that to fix the boundary condition of the differential equations, one needs to evaluate these iterative integrals at a particular value of $x$, say $x=1$ which resulted in the example on hand into multiple zeta values (MZVs) [49] and cyclotomic constants [50, 51, 36, 52]. In [53], PSLQ [54] was used to conjecture relations beyond those known from [50, 51, 36, 52].

IV. The general solution ($g(x)$) including the inhomogeneous part can now be easily given by the method of variation of constants in terms of iterative integrals as

[TABLE]

where $W(\tilde{x})$ is the Wronskian of the linear differential equation Eq. (8) and $W_{i}(\tilde{x})$ is given by

[TABLE]

V. Next we use the general solution $(g(x))$ for ${\cal\tilde{I}}^{(-l)}_{1}(x)$, plug it into Eq. (9) for $k=2,\dots,m$ and obtain solutions for the rest of the integrals and obtain ${\cal\tilde{I}}^{(-l)}(x)$ at ${\cal O}(\varepsilon^{-l})$.

VI. Now we consider the next order ($-l+1$) in the $\varepsilon$-expansion. We use the solution of ${\cal\tilde{I}}^{(-l)}(x)$ in terms of iterative integrals and plug it into Eq. (6) for $k=-l+1$. Thus we obtain a new system of the form Eq. (7) for the $\varepsilon^{-l+1}$-coefficient ${\cal\tilde{I}}^{(-l+1)}=({\cal\tilde{I}}_{1}^{(-l+1)}(x),\dots,{\cal\tilde{I}}_{m}^{(-l+1)}(x))$. Then the above procedure is repeated. From Eq. (6) it is evident that the homogeneous solution at each order remains same. The only changes happen to the inhomogeneous functions $r(x)$. As a consequence, we can reuse the already available homogeneous solutions $y_{1}(x),\dots,y_{m}(x)$ and just need to compute $g(x)$ in Eq. (12) with the updated function $r(x)$. In Ref. [8] the method has been illustrated by a detailed example.

3 Application to the massive form factors

We consider the decay of a virtual massive boson, which can be a vector ($V$), an axial-vector ($A$), a scalar ($S$) or a pseudo-scalar ($P$) of momentum $q$, into a pair of heavy quarks of mass $m$, momenta $q_{1}$ and $q_{2}$ and color $c$ and $d$, through a vertex $X_{cd}$, where $X_{cd}=\Gamma^{\mu}_{V,cd},\Gamma^{\mu}_{A,cd},\Gamma_{S,cd}$ and $\Gamma_{P,cd}$. The general form of the amplitudes can be written as

[TABLE]

Here $\bar{u}_{c}(q_{1})$ and $v_{d}(q_{2})$ are the bi–spinors of the quark and the anti–quark, respectively, with $\sigma^{\mu\nu}=\frac{i}{2}[\gamma^{\mu},\gamma^{\nu}]$. $v_{Q},a_{Q},s_{Q}$ and $p_{Q}$ are the Standard Model (SM) coupling constants for the vector, axial-vector, scalar and pseudo-scalar, respectively; $v=(\sqrt{2}G_{F})^{-1/2}$ denotes the SM vacuum expectation value of the Higgs field, with the Fermi constant $G_{F}$. More details can be found in [23]. The ultraviolet (UV) renormalized form factors ($F_{I}$) are expanded in the strong coupling constant $\alpha_{s}=g_{s}^{2}/(4\pi)$ as follows

[TABLE]

The form factors can be obtained from the amplitudes by multiplying appropriate projectors as given in [23] and performing the trace over the color and spinor indices. $n_{l}$ and $n_{h}$ are the numbers of light and heavy quarks, respectively.

The computational procedure of the three–loop massive form factors are outlined in Refs. [23, 27]. As usual the packages QGRAF [55], Color [56], Q2e/Exp [57, 58] and FORM [59, 60] have been used to generate the Feynman diagrams, calculate their color structure and to perform traces over the Dirac matrices. Since we use dimensional regularization, an appropriate description for the treatment of $\gamma_{5}$ is needed in the case of axial-vector and pseudo-scalar form factors. However, both the color–planar and complete $n_{l}$ contributions belong to the so-called non-singlet case, where the vertex is connected to open heavy quark lines. Hence, both $\gamma_{5}$-matrices appear in the same chain of Dirac matrices, which allows us to use an anti-commuting $\gamma_{5}$ in $D$ space-time dimensions, with $\gamma_{5}^{2}=1$. This is implied by the well-known Ward identity,

[TABLE]

The IBP reduction to MIs has been performed using Crusher [61]. Finally, we have obtained 109 MIs, out of which 96 appear in the color–planar case. To obtain the MIs, we have implemented the method described in the previous section. Apart from OreSys, we have intensively used HarmonicSums [37, 38, 39, 41, 42, 36, 40, 43, 44] and Sigma [62, 63]. Finally, we have performed numerical checks of all the MIs using FIESTA [64, 65, 66].

Once, we compute the MIs including the required orders in $\varepsilon$, we finally obtain the color–planar and complete light quark contributions to the unrenormalized three–loop massive form factors. In these cases, along with the strong coupling constant, the heavy quark mass and the heavy quark wave function need UV renormalization. We consider the $\overline{\rm MS}$ scheme to renormalize the strong coupling constant, where we set the universal factor $S_{\varepsilon}=\exp(-\varepsilon(\gamma_{E}-\ln(4\pi))$ for each loop order to one at the end of the calculation. However, we renormalize the heavy quark mass and wave function in the on–shell (OS) scheme. The required renormalization constants are denoted by $Z_{m,{\rm OS}}$ [67, 68, 69, 70, 71], $Z_{2,{\rm OS}}$ [67, 68, 69, 72] and $Z_{a_{s}}$ [73, 74, 75, 76, 77, 78, 79] for the heavy quark mass, wave function and strong coupling constant, respectively.

The QCD amplitudes contain infrared (IR) divergences arising from soft gluons and collinear partons. However, such behavior is universal and acts as a check of exact computations. In the case of massive form factors, the structure of IR singularities are given by [80, 81]

[TABLE]

where $F_{I}^{\mathrm{fin}}$ is finite as $\varepsilon\rightarrow 0$. $Z(\mu)$ is related to the massive cusp anomalous dimension [82, 83] through the renormalization group equation.

4 Results

The system of differential equations of the MIs appearing in the color–planar and complete light quark contributions to the massive form factors is a single scale and first order factorizable system. We apply our proposed method to this system to obtain the solutions of the MIs up to the required order in $\varepsilon$. Using these solutions for the MIs, we then obtain the form factors for vector, axial–vector, scalar and pseudo–scalar currents in the corresponding scenario. The results are lengthy and hence are provided as supplementary material along with Ref. [27, 8].

Here we illustrate in Figures 1 the behavior of the $O(\varepsilon^{0})$ parts of the vector and axial-vector form factors as a function of $x\in[0,1]$. In Figures 2, the $O(\varepsilon^{0})$ part of the scalar and pseudo-scalar form factors are presented. We also show their small– and large–$x$ expansions. The latter representations are obtained using HarmonicSums. To evaluate the HPLs and the cyclotomic HPLs numerically, we use the GiNaC package [84, 85] and the FORTRAN-codes HPOLY.f [86] and CPOLY.f [8].

4.1 Checks

We have performed a series of checks starting with maintaining the gauge parameter $\xi$ to first order and thus obtaining a partial check on gauge invariance. Fulfillment of the chiral Ward identity, Eq. (16), gives another strong check on our calculation. Considering $\alpha_{s}$–decoupling appropriately, we obtain the universal IR structure for all the UV renormalized results, confirming again the universality of IR poles. Also, in the low energy limit, the magnetic vector form factor produces the anomalous magnetic moment of a heavy quark which we cross check with [87] in this limit.

We have compared our results with those of Ref. [25, 26, 28], which have been computed using partly different methods. Both results agree.

5 Conclusion

We presented a method to solve uni-variate systems of differential equations which are first order factorizable. The system is solved in terms of iterative integrals over a finite alphabet of letters.222The present method has something in common with the method of hyperlogarithms [88], which has been used for massive systems in [89] beyond Kummer–Poincaré iterated integrals [90, 91], but applies to fully general alphabets. This method can be applied to a wide range of systems in Quantum Field Theory. We employ this method to solve such a system of differential equations appearing while solving the master integrals of the color–planar and complete light quark contributions to the three loop heavy quark form factors using the method of differential equations. Finally, we also obtain all the corresponding form factors for vector, axial–vector, scalar and pseudo–scalar currents at three loops, which play an important role in phenomenological study of top quark.