The Three-Loop Splitting Functions $P_{qg}^{(2)}$ and $P_{gg}^{(2, N_F)}$

TL;DR

This paper presents the calculation of three-loop unpolarized splitting functions in quantum chromodynamics, providing precise results for parton evolution with advanced computational techniques.

Contribution

It introduces a novel calculation of the three-loop splitting functions $P_{qg}^{(2)}$ and $P_{gg}^{(2, m N_F)}$ using massive operator matrix elements and recurrence relations.

Findings

Confirmed previous results in the literature.

Calculated $P_{gg}^{(2, m N_F)}(x)$ directly.

Derived $P_{qg}^{(2)}(x)$ from 1200 moments.

Abstract

We calculate the unpolarized twist-2 three-loop splitting functions $P_{qg}^{(2)}(x)$ and $P_{gg}^{(2,\rm N_F)}(x)$ and the associated anomalous dimensions using massive three-loop operator matrix elements. While we calculate $P_{gg}^{(2,\rm N_F)}(x)$ directly, $P_{qg}^{(2)}(x)$ is computed from 1200 even moments, without any structural prejudice, using a hierarchy of recurrences obtained for the corresponding operator matrix element. The largest recurrence to be solved is of order 12 and degree 191. We confirm results in the foregoing literature.

Tables

Table 1. Table 1: Recurrence and run-time parameters to determine the anomalous dimension γ q ​ g ( k ) ​ ( N ) superscript subscript 𝛾 𝑞 𝑔 𝑘 𝑁 \gamma_{qg}^{(k)}(N) .

	color/${\zeta }_l$	order	degree	${\tau }_{\mathrm{guess}}$ [sec]	${\tau }_{\mathrm{solve}}$ [sec]
${\gamma }_{qg}^{\left(0\right)}$	$T_F{N}_F$	1	4	0.16	0.12
${\gamma }_{qg}^{\left(1\right)}$	$C_A{T}_F{N}_F$	3	29	0.40	3.35
	$C_F{T}_F{N}_F$	3	26	0.44	2.73
${\gamma }_{qg}^{\left(2\right)}$	$C_A^2{T}_F{N}_F$	12	191	36.41	1038.31
	$C_A^2{T}_F{N}_F{\zeta }_3$	2	14	0.19	1.18
	$C_F^2{T}_F{N}_F$	10	140	12.14	487.35
	$C_F^2{T}_F{N}_F{\zeta }_3$	1	7	0.16	0.19
	$C_F{C}_A{T}_F{N}_F$	12	185	34.12	998.08
	$C_F{C}_A{T}_F{N}_F{\zeta }_3$	2	14	0.28	1.19
	$C_A{T}_F{N}_F^2$	5	59	1.62	28.95
	$C_F{T}_F{N}_F^2$	5	68	3.46	22.21

Equations

$${\hat{\hat{A}}}_{Qg}^{(3)}$$

$${\hat{\hat{A}}}_{gg,Q}^{(3)}$$

$$\hat{f}(N_F)$$

$$\tilde{f}(N_F)$$

$$(\Delta .p{)}^N\to \sum _{k=0}^{\infty }{x}^N(\Delta .p{)}^N=\frac{1}{1-x\Delta .p}\text{leavevmode}.$$

$$S_{b,\vec{a}}(N)=\sum _{k=1}^N\frac{(\mathrm{sign}(b){)}^k}{k^{|b|}}{S}_{\vec{a}}(k),\text{leavevmode}\text{leavevmode}b,a_i\in \mathbb{Z}\backslash \{0\},\text{leavevmode}\text{leavevmode}N\in \mathbb{N}\backslash \{0\},S_{\varnothing }=1.$$

$${\gamma }_{qg}^{(0)}$$

$${\gamma }_{qg}^{(1)}$$

$$P_1$$

$$P_2$$

$${\gamma }_{qg}^{(2)}$$

$$P_3$$

$$P_4$$

$$P_5$$

$$P_6$$

$$P_7$$

$$P_8$$

$$P_9$$

$$P_{10}$$

$$P_{11}$$

$$P_{12}$$

$$P_{13}$$

$$P_{14}$$

$$P_{15}$$

$$P_{16}$$

$$P_{17}$$

$$P_{18}$$

$$P_{19}$$

$$P_{20}$$

$$P_{21}$$

$$P_{22}$$

$$P_{23}$$

$$P_{24}$$

$$P_{25}$$

$$P_{26}$$

$$P_{27}$$

$$P_{28}$$

$$P_{29}$$

$$P_{30}$$

$$P_{31}$$

$$P_{32}$$

$$P_{33}$$

$$P_{34}$$

$$P_{35}$$

$$P_{36}$$

$$P_{37}$$

$$P_{38}$$

$$P_{39}$$

$${\gamma }_{gg}^{(0)}$$

$${\gamma }_{gg}^{(1)}$$

$$P_{40}$$

Full text

1 Introduction

The three-loop anomalous dimensions and splitting functions $P_{ij}^{(2)}(x)$ govern the scale-evolution of the parton distribution functions in Quantum Chromodynamics (QCD) at next-to-next-to-leading order (NNLO) and are of importance for precision predictions at $ep$- and hadron colliders, such as HERA, the Tevatron and the LHC. They are instrumental for the precise prediction of the Higgs-, top-quark, and jet-production cross sections, concerning both, the parton distribution functions [1, 2] and the strong coupling constant $\alpha_{s}(M_{Z})$ [3], thus allowing the exploration of the hitherto heaviest systems in the Standard Model in greatest possible detail. High precision predictions are also of importance for planned $ep$-facilities as EIC and LHeC [4], operating at high luminosity. A first computation of the three-loop splitting functions was performed in Refs. [5, 6], after the calculation of fixed moments for these quantities in Refs. [7, 8, 9]. In these references the anomalous dimensions were calculated using massless processes.

In the present paper we compute the splitting function $P_{qg}^{(2)}$ and the $N_{F}$-dependent part of $P_{gg}^{(2)}$ independently using massive three-loop operator matrix elements (OMEs). They are necessary for the computation of the heavy flavor contributions to deep-inelastic scattering in the region of virtualities $Q^{2}$ much larger than the heavy quark mass squared $m^{2}$. While we obtain the complete expression for $P_{qg}^{(2)}$, only the $N_{F}$-part of $P_{gg}^{(2)}$ emerges, since the respective process is proportional to the color factor $T_{F}$.111In $SU(N_{c})$ the color factors contributing to the three-loop anomalous dimensions are given by $T_{F}=1/2,C_{A}=N_{c},C_{F}=(N_{c}^{2}-1)/(2N_{c})$, with $N_{c}$ the number of colors and $N_{c}=3$ in case of QCD.

We calculate $P_{gg}^{(2,\rm N_{F})}$ directly by applying the techniques described in Ref. [10]. In the case of $P_{qg}^{(2)}$, we apply the recently proposed method of large Mellin moments [11]. In previous papers we have obtained the corresponding contributions $\propto T_{F}$ for the splitting functions $P_{gq}^{(2,\rm N_{F})},P_{qq}^{(2,\rm NS,N_{F})},P_{qq}^{(2,\rm NS,TR,N_{F})}$ [12, 13] and the complete splitting function $P_{qq}^{(2),\rm PS}$ [14]. A series of Mellin moments for these quantities has been obtained before in Ref. [9]. Our results agree with those of the previous calculation in Refs. [5, 6]. Recently, parts of the three-loop anomalous dimensions $\propto T_{F}$ have also been obtained for a series of anomalous dimensions in two-mass calculations [15].

The paper is organized as follows. In Section 2 we describe the calculation of the anomalous dimensions $\gamma_{qg}^{(2)}$ and $\gamma_{gg}^{(2),\rm N_{F}}$ in Mellin-$N$ space. Section 3 contains the conclusions. In the Appendices we present the splitting functions in $x$-space and give some technical details in calculating master integrals. In all representations we project onto the algebraic basis [16], both for the harmonic sums [17, 18] and the harmonic polylogarithms (HPLs) [19].

2 The Anomalous Dimensions

In the following, we will use the massive OMEs $A_{Qg}^{(3)}$ and $A_{gg,Q}^{(3)}$, cf. Ref. [9], to calculate the anomalous dimensions. The unrenormalized OMEs $\hat{\hat{A}}_{Qg}^{(3)}$ and $\hat{\hat{A}}_{gg,Q}^{(3)}$ were given in [9]

[TABLE]

Here, $m$ denotes the heavy quark mass, $\varepsilon=D-4$ the dimensional parameter, $\mu$ the factorization scale, $\zeta_{l},l\in\mathbb{N},l\geq 2$ denotes the values of the Riemann $\zeta$ function at integer argument, $N_{F}$ is the number of massless quark flavors, $\beta_{i}$ the expansion coefficients of the QCD $\beta$-function, $\beta_{i,Q}$ related expansion coefficients associated to heavy quark effects, $\gamma_{ij}^{(k)}$ the expansion coefficients of the anomalous dimensions, and $\delta m_{k}^{(l)}$ the expansion coefficients of the unrenormalized quark mass, cf. [9]. The coefficients $a_{ij}^{(k)}$ denote the constant terms of the OMEs at $k$-loop order and $\bar{a}_{ij}^{(k)}$ the corresponding terms at $O(\varepsilon)$, [20, 21, 22, 23, 24, 9].

Furthermore, we use the convention

[TABLE]

The pole terms $O(1/\varepsilon^{3})$ and $O(1/\varepsilon^{2})$ contain the complete $O(a_{s})$ and $O(a_{s}^{2})$ anomalous dimensions, and the single pole term allows to obtain the terms $\propto T_{F}N_{F}$ of the $O(a_{s}^{3})$ anomalous dimensions $\gamma_{qg}^{(2)}$ and $\gamma_{gg}^{(2)}$, which in the case of $\gamma_{qg}^{(2)}$ is the complete anomalous dimension.

The Feynman diagrams contributing to the OMEs $A_{Qg}^{(3)}$ (2.1) and $A_{gg,Q}^{(3)}$ (2.2) were generated by the code QGRAF [25]222See Ref. [9] for the implementation of the local operators.. The color configurations were calculated using the package Color [26] and the Feynman integrals were reduced to master integrals using the package Reduze 2 [27, 28].333The package Reduze 2 uses the packages FERMAT [29] and Ginac [30].. There are different techniques available to calculate the master integrals, which have been summarized in Refs. [11] and [10].

2.1 The Anomalous Dimension $\gamma_{qg}^{(2)}$

In calculating $\gamma_{qg}^{(2)}$ we have used the method of large moments [11]. First the local operator insertions $(\Delta.p)^{N}$ are formally resummed into propagator like terms by

[TABLE]

The system of differential equations in the variable $x$ obeyed by the master integrals can be obtained by integration by parts identities through Reduze 2. The master integrals can then be expanded into a formal Laurent-Taylor series, which can be inserted in the differential equations, leading to associated difference equations, with known initial conditions. One can then exploit these difference equations to generate a large set of moments for each color and $\zeta$-value coefficient, expressed as a sequence of rational numbers. In the present approach it is required that the master integrals have to be known to a higher order in the dimensional parameter $\varepsilon$. Usually standard techniques allow to derive the corresponding expansion. In one particular case we had to spend some effort to obtain the initial conditions for some master integrals, which is described in Appendix B.

With the corresponding series of rational numbers at hand, we use the algorithm of guessing [31] to derive their associated minimal recurrence. Using the package Sigma [32, 33] this recurrence is then solved in terms of nested sums and the solutions are afterwards transformed to special function formats using the package HarmonicSums [34, 35, 36, 37, 38].

The recurrences derived for the anomalous dimension are first order factorizable and have therefore sum- and product solutions according to difference field theory [39, 40, 41, 42, 43, 44, 45, 46]. Moreover, the solutions turn out to be given even in terms of harmonic sums only. To obtain a solution conform with a later analytic continuation one should consider even moments only, which would require a twice as large set of moments as input. In many problems this number is so high that it might be difficult to reach. Alternatively, one can start by taking all moments and derive a first expression at general values of $N$, from which only the even moments are then projected. This is algebraically possible by considering $(-1)^{N}$ as a variable $y$ and performing division with remainder w.r.t. the polynomial $1-y$ (or $1+y$). In this way we obtain the $N$th moment of the single pole term. The use of also odd moments is possible in other approaches and recommended if a lower number of moments allows to obtain the difference equation to be determined. It yet corresponds to a different problem in general. However, given a recurrence for odd and even values, one can now generate easily a very high number of moments using this recurrence, can extract afterwards the needed even moments and can repeat the guessing process with the now sufficiently large number of even moments.

In the calculation we computed only the irreducible contributions and added the other terms known already as a function of $N$, cf. Refs. [20, 21, 22, 23, 24, 9]. Here the master integrals are not solved in terms of sums, but only their moments are determined and those are inserted into the expressions of $A_{Qg}^{(3)}$, for which, again, only moments are calculated at first. We have then projected the irreducible contributions, using the known terms, on to $\gamma_{qg}^{(2)}$ using Eq. (2.1).

The order and degree of the recurrences associated to the different color and $\zeta$-factors, as well the times to guess the corresponding recurrences using Guess [31] and to solve them using the packages Sigma [32, 33] and HarmonicSums [34, 35, 36, 37, 38] are summarized in Table 1. We have generated 2000 even moments for the unrenormalized expression $\hat{\hat{A}}_{qg}^{(3)}$ up to its $1/\varepsilon$ term using the package SolveCoupledSystem in which the method of large moments is implemented [11]. The computation of the moments for the master integrals took 9.30 days per kernel, where we parallelized to 15 kernels. The creation of all pole terms in $\hat{\hat{A}}_{qg}^{(3)}$ took 10.67 days/kernel where we parallelized to 25 kernels. The moments of the OME formed the input of the subsequent calculation. We used the even moments to guess the corresponding recurrences using Guess [31]. The implementation of the guesser in Sage [47, 48] needs a much shorter time to find a recurrence than an earlier one written in Mathematica [31]. The total guessing time of $\sim 89$ sec was smaller than the corresponding one using the Mathematica-version by a factor 75. The Sage-version of the guesser does currently not report the number of moments used. For the most demanding example, the $C_{A}^{2}T_{F}N_{F}$-term, 1000 moments were not sufficient, but 1200 moments were. In our earlier study [49], slightly different recurrences were obtained, since there we determined them from even and odd moments. The present recurrences are somewhat shorter.

The computation time for the last two steps to obtain $\gamma_{qg}^{(2)}$ amounted to $44.6$ min, which included the time to create the different recurrences through guessing and their solution to obtain the closed expressions as a function of the Mellin variable $N$.

This computational time and memory requirement are much shorter than that one using the traditional method, where first all master integrals are calculated in explicit form and are then inserted to obtain the final result. The costs to determine the different master integrals as functions of $N$ would be much higher, if possible at all, because elliptic structures may arise at higher orders in $\varepsilon$.

In any case, the time needed to reduce the OME to master integrals is the longest part of the calculation. Using Reduze 2 [27, 28] the reduction time amounted to several months. Similar or larger run times are expected also for other reduction programs.

The $O(a_{s})$ and $O(a_{s}^{2})$ anomalous dimensions are read from the $O(1/\varepsilon^{3})$ and $O(1/\varepsilon^{2})$ terms and $\gamma_{qg}^{(2)}$ is determined from the $O(1/\varepsilon)$ term. The expressions are reduced to the polynomial basis in the harmonic sums [17, 18] $S_{\vec{a}}(N)\equiv S_{\vec{a}}$,

[TABLE]

One obtains

[TABLE]

with the polynomials

[TABLE]

The anomalous dimensions $\gamma_{qg}^{(0,1)}$ agree with results given in Refs. [50, 51, 52, 53, 54, 55, 56, 57], or with their Mellin transform.

For the three-loop anomalous dimension $\gamma_{qg}^{(2)}$, we obtain

[TABLE]

where

[TABLE]

The anomalous dimension $\gamma_{qg}^{(2)}$, Eq. (2.11), agrees with the corresponding expression given in Ref. [6]. The terms $\propto N_{F}^{2}$ were confirmed before in [58].

2.2 The Anomalous Dimension $\gamma_{gg}^{(2),\rm N_{F}}$

The anomalous dimensions $\gamma_{gg}^{(0,1)}(N)$ can be obtained from the pole terms $O(1/\varepsilon^{3})$ and $O(1/\varepsilon^{2})$ of massive OME $A_{gg,Q}^{(3)}$, while the single pole term $O(1/\varepsilon)$ contains $\gamma_{gg}^{(2),\rm N_{F}}(N)$, by virtue of being $\propto T_{F}$. We compute the contributions to $A_{gg,Q}^{(3)}$ up to the $1/\varepsilon$ terms in the traditional way, i.e. the contributing master integrals are first calculated in explicit form. A subset of master integrals can be obtained using representations through (generalized) hypergeometric functions [59, 60] or using Mellin-Barnes integrals [61, 62]. The resulting definite sums have been simplified afterwards with symbolic summation using the packages Sigma [32, 33], EvaluateMultiSums [63], and HarmonicSums [34, 35, 36, 37, 38]. A further large set is calculated using the method of differential equations [64, 10]. Here we use the algorithm described in Ref. [10], which is available in SolveCoupledSystem using the packages Sigma, HarmonicSums, EvaluateMultiSums, SumProduction [63] and OreSys [65]. This algorithm allows to solve the differential equations analytically in the one-parameter case without the need to use any specific basis of master integrals. In the present calculation of anomalous dimensions, the corresponding difference equations are factorizing in first order and can be solved analytically using the package Sigma. In a single case the computation time turned out to be large. Therefore, we used the Almkvist-Zeilberger theorem [66] and our packages MultiIntegrate [35] to derive a linear recurrence. Solving the recurrence with Sigma and HarmonicSums provided then an immediate solution. The master integrals are obtained as expressions in the variable $x$, see Eq. (2.5). The solution in Mellin $N$-space is then given by the $N$th expansion coefficient of these expressions, which is found by a routine of HarmonicSums.

We obtain the following expressions for $\gamma^{(0,1),N_{F}}_{gg}$:

[TABLE]

with

[TABLE]

The anomalous dimensions $\gamma_{gg}^{(0,1)}$ agree with results given in [50, 51, 52, 53, 54, 55, 56, 57], or with their Mellin transform.

For the contribution $\propto N_{F}$ to $\gamma^{(2)}_{gg}$ at three-loop order, $\gamma^{(2),N_{F}}_{gg}$, we obtain:

[TABLE]

where the polynomials are

[TABLE]

The anomalous dimension $\gamma_{gg}^{(2),\rm N_{F}}$, Eq. (2.55), agrees with the corresponding expression given in Ref. [6]. We have confirmed the leading $N_{F}$-terms $\propto N_{F}^{2}$ in [68, 67], which were also given in [69].

3 Conclusions

We have computed the 3-loop anomalous dimension $\gamma_{qg}^{(2)}(N)$, the contributions $\propto N_{F}$ for $\gamma_{gg}^{(2)}(N)$ and the associated splitting functions in a massive calculation, which is fully independent of the earlier computation in Ref. [6]. We agree with the previous results. In the case of $\gamma_{qg}^{(2)}(N)$, we have performed a fully automatic calculation generating the Feynman diagrams, reducing them to master integrals and calculating them using differential and difference equations. We used formal Taylor series representations in terms of the subsidiary parameter $x$ resumming the local operator insertions. The new method of high Mellin moments [11] has been used. Here the moments are calculated recursively using the difference equation systems associated to the differential equations. We did not compute the master integrals as explicit functions of $N$, but kept them as vectors of moments. The method of guessing has then been used to obtain one difference equation for each color/$\zeta$-value factor. All difference equations obtained factorize at first order and one finally obtains a representation in harmonic sums.

The contributions $\propto N_{F}$ to $\gamma_{gg}^{(2)}(N)$ have been calculated using more traditional methods. After the reduction to the master integrals, we have calculated these as functions of the Mellin variable $N$. The corresponding techniques are discussed in detail in Ref. [10] and include (generalized) hypergeometric and associated higher transcendental function representations [59, 60], the method of hyperlogarithms [70, 71, 72], differential equation techniques [64, 10], and integration using the Almkvist-Zeilberger theorem [66, 35]. The final solution has then been obtained by inserting the master integrals and extracting $\gamma_{gg}^{(2),\rm N_{F}}(N)$ from the single pole term of $\hat{\hat{A}}_{gg}^{(3)}$. In this process, quite a series of sums cancel, including all non-harmonic sums.

The following Mathematica, Fortran, and FORM [73] packages, see also Ref. [9]444Here the renormalization of the massive operator matrix elements has been carried out in a maple code., were instrumental in performing the present calculation: QGRAF [25, 9], Color [26], MB [62], Reduze 2 [27, 28], FERMAT [29], Ginac [30], Sigma [32, 33], EvaluateMultiSums and SumProduction [63], SolveCoupledSystems [64, 10, 11], OreSys [65], HarmonicSums [34, 35, 36, 37, 38] and Guess [31, 48].

The universality of the QCD anomalous dimension allows to compute them within various setups. In the present calculation, they were obtained from the pole structure of massive OMEs. The calculation of these OMEs is part of an ongoing project with the final goal to compute the massive Wilson coefficients for deep-inelastic scattering in the region $Q^{2}\geq m^{2}$.

The anomalous dimensions and splitting functions presented in this paper, together with our earlier results in Refs. [12, 13, 58, 68, 67, 14], are given in Mathematica.m format in the attachment, together with the expressions for $I^{\rm B5a}_{\rm 101}$ for $D=4+\varepsilon$ to $D=30+\varepsilon$ up to $O(\varepsilon^{6})$ presented in appendix B. Their conversion to maple or Form-inputs [73] is straightforward.

Appendix A The Splitting Functions

In the following we present the splitting functions $P_{qg}^{(i)}$ and $P_{qg}^{(i)},i=0,1,2$ in $x$-space. They are obtained as a Mellin-inversion from the even moments of $\gamma_{qg}^{(i)}$ and $\gamma_{gg}^{(i)}$ via

[TABLE]

The splitting functions to 3-loop order can be expressed in terms of harmonic polylogarithms $H_{\vec{a}}(x)\equiv H_{\vec{a}}$ over the alphabet $f_{b}(x)\in\{1/x,1/(1-x),1/(1+x)\}$ [19]. They are defined by

[TABLE]

The splitting functions are :

[TABLE]

[TABLE]

and

[TABLE]

[TABLE]

Appendix B Calculation of the initial values

A sufficient number of initial values has to be provided to solve the recurrence relations which determine the master integrals. These are the moments of the master integrals for low values of the Mellin variable $N$. Up to $\mathcal{O}(\varepsilon^{0})$ they can be calculated using the MATAD [74] setup from [9, 75]. Beyond this order, other methods have to be used. If a multi-sum representation has already been derived to the necessary order in $\varepsilon$, the required initial values can be derived from it. However, since this involves some degree of manual intervention and since an appropriate representation cannot always be found, a fully algorithmic method is desirable. In the following, we describe an algorithmic method which works for all master integrals arising in the ongoing project to calculate the massive OMEs to three-loop order. The method has been briefly described in [67] and applied in several calculations of this project [67, 14, 10]. It is based on the $\alpha$-parameterization for Feynman integrals, see, e.g., [76, 77, 78], and the idea [79, 80] of rewriting integrals with scalar products in the numerator in terms of scalar integrals in shifted dimensions with raised propagator powers. This allows us to map fixed moments of the master integrals with operator insertions to operator-less tadpole integrals. These can be reduced using integration-by-parts (IBP) identities [81, 27, 28] and lead to linear combinations of only three tadpole master integrals. Once these are calculated in shifted dimensions to sufficient order in $\varepsilon$, we can assemble the moments of the master integrals with operator insertions. The steps described here can be implemented efficiently in Mathematica.

B.1 Structure of the integrals

The scalar integrals with operator insertions for which we need initial values are two-point integrals with an on-shell, massless external momentum $p^{2}=0$. The Feynman rules of the operators involve a light-like vector $\Delta$ and they introduce symbolic powers of scalar products $\Delta.k$, where $k$ is a linear combination of loop momenta and $p$. Due to the structure of the Feynman rules, we can deal with these scalar products elegantly by introducing an ordinary generating function for the Mellin moments of the integrals. In the simplest case of a scalar product $(\Delta.k)^{N}$, we introduce a tracing variable $x$, multiply by $x^{N}$ and sum over $N$, cf. Eq. (2.5). Similar terms appear also for more complicated operator insertions. The operator insertions always lead to products of terms like $(1-x\,\Delta.k)^{-1}$, which resemble propagators that are linear in the momentum $k$.

All scalar integrals of this project can be expressed in terms of 24 integral families with nine quadratic propagators (which depend quadratically on a linear combination of loop and external momenta) and three linear propagators that encode the structure of the operator insertion. We use the notation

[TABLE]

where $P_{i}$ denotes the propagator denominators of integral family $f$. The first nine indices refer to the quadratic propagators and the last three to the linear propagators. The 24 twelve-propagator integral families are based on six nine-propagator integral families (denoted B1, B3, B5, C1, C2 and C3), which are supplemented by different sets of linear propagators (indicated by an additional letter in the family identifier). For example, the integral family B5a has the propagators

[TABLE]

We first discuss the properties of the quadratic propagators. The corresponding nine-propagator families differ by the placement of the mass $m$ and by the routing of the external momentum $p$. If no external linear propagator is present, they can only depend on the mass $m$ and the external momentum $p$. Thus, they reduce to massive tadpole integrals, since $p^{2}=0$, and we can remove $p$ from the propagators altogether. Then three quadratic propagators become linearly dependent of the remaining six ($P_{2}$, $P_{4}$ and $P_{9}$ in the example above) and the six nine-propagator families reduce to only two six-propagator tadpole integral families.

Let us briefly recapitulate the basic facts about the $\alpha$-parameterization and its use for expressing integrals with irreducible numerators in terms of scalar integrals in shifted dimensions. For more details, see for example [76, 77, 78, 79, 80]. The $\alpha$-parameterization for the above integrals without linear propagators reads after applying a Wick rotation

[TABLE]

where we have explicitly pulled out a factor $(-1)$ from each propagator $P_{i,E}$ due to the Wick rotation. Below we suppress the subscript $E$ and work with Euclidean momenta unless explicitly stated otherwise. Each propagator is at most quadratic in the loop momenta. Hence, we can write any linear combination of propagators in the form

[TABLE]

The matrix $A$ is real and symmetric, and the vectors $b_{i}$ will be linear combinations of external momenta. Performing the loop integrals leads to

[TABLE]

The elements of the matrix $A$ are the coefficients of the scalar products of the loop momenta, while the elements of the vector $b$ are the coefficients of the terms linear in the loop momenta – i.e. scalar products with the external momenta (or other vector quantities). The scalar term $c$ finally contains the masses of the propagators. For our purposes, we have to consider the following form of the propagators:

[TABLE]

The $\delta$ coefficients are either $1$, [math] or $-1$ depending on how these terms enter the propagator under consideration. Each propagator enters the exponential multiplied by its corresponding $\alpha$ parameter if it is present in the integral. We can model this by replacing $\alpha_{i}$ with $\beta_{i}=\alpha_{i}\Theta(\nu_{i}-\frac{1}{2})$, which vanishes if the propagator is not present. Thus, we obtain

[TABLE]

Here, we have already used that $p^{2}=0$.

For the integral families of this project, the explicit form of the matrix $A$ and the vector $b$ depend only on whether the integral belongs to one of the families B1, B3 or B5 or to C1, C2 or C3. Using the notation introduced above, the quantities read

[TABLE]

The scalar term $c$ is different for each of the six integral families. However, we skip it here since its specific form does not play a role in the following.

B.2 Treatment of the linear propagators

If the integrals contain scalar products involving the loop momenta in the numerator, we can include them in the $\alpha$-parameterization by introducing an exponential generating function for them, i.e.

[TABLE]

We introduce one auxiliary parameter $r$ for each scalar product. Thus, auxiliary parameters associated to scalar products of two loop momenta, $k_{i}.k_{j}$, end up as additional terms in the matrix elements $A_{ij}$ (e.g., $r_{1}\,k_{1}.k_{2}$ in the exponential leads to the replacement $A_{12}\to A_{12}+r_{1}$). The auxiliary parameters which accompany scalar products with one loop momentum, $p.k_{i}$, are added to the parameter $b_{i}$ (e.g., the scalar product $r_{2}\,p.k_{1}$ yields $b_{1}\to b_{1}+r_{2}p$).

To derive a parameter representation for fixed moments of the master integrals with operator insertions, we first derive the $\alpha$-parameterization for the corresponding operator-less integral. Afterwards, we add in the linear propagators. Each of them has the structure $(1-x\,\Delta.k)^{-1}$, where $k$ is a linear combination of loop and external momenta. As discussed above, it is introduced as an ordinary generating function for the moments of the master integrals given in Eq. (2.5). If there are multiple linear propagators, all of them have the same tracing parameter $x$, such that the expansion involves finite sums arising from the Cauchy products. For example, two linear propagators yield

[TABLE]

which provides the generating function of the Feynman rule for an operator insertion on a vertex. Analogously, for three linear propagators, we obtain

[TABLE]

In case a linear propagator is raised to a higher power, we use

[TABLE]

Thus, given a number of linear propagators, each of them raised to a given power $n_{i}$, i.e. $(1-x\,y_{i})^{-n_{i}}$, we can extract the coefficient of $x^{N}$ of the series expansion in $x$ as follows: We introduce a summation index $j_{i}$ and a factor

[TABLE]

for each linear propagator, multiply them together and sum over $j_{i}$ from [math] to $j_{i+1}$. Of course, there is no sum over $j_{0}$ and the outermost index (which is never summed over) must be renamed to $N$.

To include the operators into our description, we think of the sums we just introduced as a fully expanded polynomial in the scalar products of the linear propagators. Each term has a number of scalar products each raised to a definite power and for each term, we use the exponential generating function representation, described in Eq. (B.10). To summarize: For each linear propagator we get a binomial coefficient which encodes the combinatorics of multiple, identical linear propagators, and we replace each power of a scalar product by a corresponding power of the derivative operator w.r.t. the auxiliary parameter $r_{i}$ (which belongs to the $i$th linear propagator).

Given all these ingredients, we can assemble a general formula for the master integrals with operator insertions for the case of three linear propagators

[TABLE]

Here, $A$ is completely determined by the operator-less part of the integral ($A^{\text{B}}$ and $A^{\text{C}}$ from Eqs. (B.7), (B.8)). $b$ has one term, proportional to $p$, which is contributed by the operator-less part ($b^{\text{B}}$ and $b^{\text{C}}$ from Eq. (B.9)) and one part, proportional to $r\Delta$, which stems from the linear propagators. We describe this part in the next paragraph below. The scalar term $c$ contains the mass term, but it remains unchanged by the operators, unless the linear propagator involves the external momentum $p$.

The coefficients $b$ of the linear term are modified by the introduction of the exponential generating functions for the operators: To each $b_{i}$ we have to add a term $r_{j}\Delta$ where the linear propagators contain the loop momentum $k_{i}$. Let us clarify this with an example: For family B1a we have the linear propagators

[TABLE]

We introduce three parameters, $r_{1}$, $r_{2}$ and $r_{3}$, for the three propagators and arrive at the exponential generating function

[TABLE]

and thus we obtain for $b$

[TABLE]

Similar terms arise for the other integral families.

If the linear propagators contain the external momentum $p$, we have to also modify the scalar term $c$. In that case we have to add a term proportional to $\Delta.p$. As an example, we use family B3b, which has the linear propagators

[TABLE]

and leads to the exponential generating function

[TABLE]

and hence

[TABLE]

The central idea of the dimensional shifts is now that derivatives with respect to $r_{i}$ will act on the exponential function and introduce factors of $(A^{-1})_{ij}=(\det(A))^{-1}\tilde{A}_{ij}$, where $\tilde{A}_{ij}$ is the $(i,j)$ cofactor of the matrix $A$. Since the dimension $D$ enters the parametric representation Eq. (B.15) only in the exponent of the determinant (except for constant prefactors), we can reinterpret the integrals arising from taking the derivatives as integrals in shifted dimensions, according to the exponent of the determinant. The cofactors $\tilde{A}_{ij}$ are polynomials in the $\alpha$ parameters and can be absorbed into the powers of $(\alpha_{k})^{\nu_{k}-1}$ which encode the powers to which the quadratic propagators are raised. More concretely, the powers of the propagators enter the $\alpha$ representation in the form

[TABLE]

which means that an additional factor of $\alpha_{j}^{k}$ in the integrand corresponds to raising the power of propagator $P_{j}$ by $k$ (i.e., $\nu^{\prime}_{j}=\nu_{j}+k$) and multiplying the expression by, cf. e.g. [60],

[TABLE]

to compensate for the changed $(-1)^{\nu^{\prime}_{i}}$ and $\Gamma(\nu^{\prime}_{i})$. Here $(\nu_{i})_{k}$ denotes the Pochhammer symbol. Thus, each derivative shifts the dimension and raises the powers of some propagators. After applying all derivatives, we obtain a linear combination of scalar integrals in shifted dimensions with dotted propagators and no operator insertions. These can then be mapped to tadpole integrals.

Let us now have a closer look at how the derivatives act on the parameter representation. For the master integrals we are concerned with, the auxiliary parameters $r_{i}$ only enter $b$ and $c$. The vector $b$ has two terms—a term proportional to $p$ and one proportional to $\Delta$. Since $b$ only appears in the bilinear form $b^{\mathsf{T}}A^{-1}b$ and since $p^{2}=0$ and $\Delta^{2}=0$, only the mixed terms proportional to $\Delta.p$ survive. Therefore, all terms are linear in the $r_{i}$. Thus, each derivative acting on this part of the exponential brings down exactly one matrix element of $A^{-1}$ and one power of $\Delta.p$. This means that each derivative shifts the dimension of the integral from $D$ to $D+2$. For the scalar term $c$ we also have a linear dependence on the $r_{i}$. Here, however, we do not get any shift of dimension from a derivative acting on this part of the exponential, since it only brings down a power of $2\Delta.p$. This leads to linear combinations of integrals with different dimensional shifts. In all cases, each derivative introduces a factor $\Delta.p$, so that the $N$th moment of any master integral will always be proportional to $(\Delta.p)^{N}$.

Finally, we have to undo the Wick rotation. The only scalar product left in the final expression is $(\Delta.p)^{N}$ which introduces an additional, global factor of $(-1)^{N}$ into (B.15).

B.3 Reduction of the tadpole integrals

The procedure described above maps a scalar integral with an operator insertion in $D=4+\varepsilon$ dimensions to a (possibly large) linear combination of operator-less tadpole integrals in shifted dimensions. To cope with these, we use IBP relations derived with Reduze 2 to map them onto linear combinations of master integrals. Their coefficients are rational functions in the dimension $D$ and the mass squared $m^{2}$. The reductions are done for general values of $D$ so that they can be used for integrals in any dimension. Particular care must be taken if the initial values are mapped to a linear combination of tadpole integrals in different dimensions.

As already observed in [82], all tadpole integrals can be reduced to just three master integrals. They correspond to

[TABLE]

and their diagrams are shown in Figure 1.

B.4 Required tadpole integrals

The first two tadpole integrals are straightforwardly obtained and can be expressed completely in terms of $\Gamma$-functions, for example using their Feynman parameterization. The first one is just the one-loop tadpole cubed, which reads

[TABLE]

where we use the normalization factor

[TABLE]

and the second one is the three-loop banana graph with two massless lines, yielding

[TABLE]

These expressions can be expanded in $\varepsilon$ in any dimension $D=2n+\varepsilon$. When expanding in $\varepsilon$, we keep the normalization $\mathcal{N}$, the mass scale $m^{3D-2\nu}$ and an additional factor of $\exp(\frac{3}{2}\varepsilon\gamma_{\text{E}})$ unexpanded to simplify the expressions. The latter is the same regardless of the integer dimension around which we expand. The third integral is the three-loop banana graph with fully massive lines and reads

[TABLE]

In [82] this integral was found to correspond to an expression involving $\Gamma$-functions and a generalized hypergeometric function ${}_{3}F_{2}$ for general values of $D$. Here, we are interested in the expansion around $\varepsilon=D-2n$ for different values of $n\in\mathbb{N}$ up to $2n=30$. To calculate it, we follow the method used in [83, 84, 85], which we briefly describe below.

We calculate this integral as a special case of the same integral with two different masses, by deriving a differential equation in the ratio of the masses and then specializing to the case of equal masses. The two-mass diagram is depicted in Figure 2a.

[TABLE]

where we define $x=\frac{m_{2}}{m_{1}}$. Since the overall factor of $m_{1}^{\alpha}$, $\alpha\in\mathbb{R}$ can be reconstructed from dimensional arguments, we set $m_{1}=1$ in the remaining discussion. For $x=0$ the integral with $a=1$ corresponds to $J_{1}(0)=I^{\text{B1a}}_{101}$ while for $x=1$ it corresponds to $J_{1}(1)=I^{\text{B5a}}_{101}$. Moreover, we can derive the property

[TABLE]

from the unevaluated integral.

The system of differential equations with respect to $x$ can be derived using for example Reduze 2 [27, 28]. In order to close the system, four integrals are necessary. Besides $J_{1}(x)$, also $J_{2}(x)$, i.e. the same integral with the first propagator raised to the second power (see Figure 2b), contributes. Moreover, the integrals $T_{1}(x)$ and $T_{2}(x)$ appear (see Figure 2c and 2d, respectively). They correspond to products of the massive one-loop tadpole integral and can be expressed in terms of $\Gamma$-functions for any value of $D$. Expressed in terms of the single scale tadpole they evaluate to (again setting $m_{1}=1$)

[TABLE]

Therefore, we treat them as known inhomogeneous parts of the differential equation. The resulting system of first order differential equations reads

[TABLE]

These equations can be uncoupled into a single, second order differential equation – either by using the package OreSys [86, 65] or here even by hand. We obtain

[TABLE]

which is a second order differential equation with regular singular points at $x\in\{-1,0,1,\infty\}$. Thus, we can only use the initial condition $J_{1}(0)=I^{\text{B1a}}_{101}$, while the derivative of $J_{1}$ at $x=0$ becomes linearly dependent. We could use the value of $J_{1}$ at $x=1$ to fix the second constant, however, this is just the quantity we want to calculate here. Therefore, we use the relation between $J_{1}(x)$ and $J_{1}(x^{-1})$ given in Eq. (B.32) instead.

The differential equation can be expanded in $\varepsilon=D-2n$, $n\in\mathbb{N}$ and solved order by order. The homogeneous equation is the same at every order and therefore also the homogeneous solutions are the same. They can be found for example using the Mathematica routine DSolve. The solutions contain polynomials and logarithms in $x$. For $D=4$ they are given by

[TABLE]

The inhomogeneous part $r(x)$ of the differential equation arises from the expansion of $T_{1}(x)$ and $T_{2}(x)$ in $\varepsilon$ as well as lower orders of the solution to $J_{1}(x)$. We then use the variation of constants to calculate a particular solution to the inhomogeneous differential equation, cf. [87],

[TABLE]

where $W(x)=h_{1}(x)h_{2}^{\prime}(x)-h_{2}(x)h_{1}^{\prime}(x)$ is the Wronskian of the homogeneous solutions. For $D=4$ it is $W(x)=x^{2}(1-x)(1+x)$. The structure of these integrals tells us that they only lead to iterated integrals over the letters $x$, $1-x$ and $1+x$, i.e. the harmonic polylogarithms [19]. Performing these integrals can be done efficiently using the HIntegrate command from the package HarmonicSums. Finally, we fix the constants of integration using the expansion of $J_{1}(0)=I^{\text{B1a}}_{112}$ along with the relation Eq. (B.32). Since only HPLs occur, the relation maps them to linear combinations of HPLs again. If both sides are algebraically reduced, we can solve for the constants by comparing coefficients of the HPLs. Of course, this requires the algebraic relations for HPLs of argument $x$ to be available. In intermediate steps of the calculation at hand, we needed these relations up to weight ${\sf w=9}$, which can be calculated, e.g., using HarmonicSums, cf. also [17, 88]. The expansion of the solution for $J_{1}(x)$ is given completely in terms of HPLs with polynomial coefficients in $x$. For $D=4+\varepsilon$ we obtain

[TABLE]

where we suppressed the higher order terms in $\varepsilon$ due to their length.

To obtain the solution for $I^{\text{B5a}}_{101}=J_{1}(1)$, we express the HPLs at $x=1$ in terms of multiple zeta values and use their relations [88] to reduce the occurring constants to an algebraically independent basis. The full result for $D=4+\varepsilon$ reads

[TABLE]

Besides the usual zeta values $\zeta_{k}$ and the polylogarithms evaluated at $\frac{1}{2}$, $\mathrm{Li}_{n}(\frac{1}{2})$, also the following constants appear [17, 88]

[TABLE]

where the $\sigma_{\vec{a}}$ are multiple zeta values, defined by $\sigma_{\vec{a}}=\lim_{N\to\infty}S_{\vec{a}}(N)$. Except for the zeta values $\zeta_{k}$ and $B_{4}$ [82] the remaining multiple zeta values cancel in the final results for the operator matrix elements.

Similar results are obtained for this integral in all even dimensions $D=4+\varepsilon$ to $D=30+\varepsilon$ up to and including the $\mathcal{O}(\varepsilon^{6})$ terms. This allows the calculation of moments up to $N=13$ for the master integrals with operator insertions. Since this exceeds the number of required initial values for the recurrences, we can use the remaining moments as a cross check for their solutions.

Acknowledgment. The diagrams were drawn using Axodraw [89]. We would like to thank C. Anastasiou, M. Kauers, P. Marquard, M. Round, and Y. Schröder for discussions. This work was supported in part by the Austrian Science Fund (FWF) grants P20347-N18 and SFB F50 (F5009-N15) and the European Commission through contract PITN-GA-2012-316704 (HIGGSTOOLS). We gratefully acknowledge support by the MITP and the Institute for Physics at Johannes Gutenberg University, Mainz, where part of the research by A.v.M. presented in this article was performed.