The two-mass contribution to the three-loop pure singlet operator matrix element

Available online at www.sciencedirect.com

ScienceDirect Nuclear Physics B 927 (2018) 339–367 www.elsevier.com/locate/nuclphysb

The two-mass contribution to the three-loop pure singlet operator matrix element J. Ablinger a , J. Blümlein b,∗ , A. De Freitas b , C. Schneider a , K. Schönwald b a Research Institute for Symbolic Computation (RISC), Johannes Kepler University, Altenbergerstraße 69, A-4040, Linz,

Austria b Deutsches Elektronen–Synchrotron, DESY, Platanenallee 6, D-15738 Zeuthen, Germany

Received 17 November 2017; accepted 22 December 2017 Available online 28 December 2017 Editor: Tommy Ohlsson

Abstract We present the two-mass QCD contributions to the pure singlet operator matrix element at three loop order in x-space. These terms are relevant for calculating the structure function F2 (x, Q2 ) at O(αs3 ) as well as for the matching relations in the variable flavor number scheme and the heavy quark distribution functions at the same order. The result for the operator matrix element is given in terms of generalized iterated integrals that include square root letters in the alphabet, depending also on the mass ratio through the main argument. Numerical results are presented. © 2017 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Funded by SCOAP3 .

1. Introduction Starting at 2-loop order, massive operator matrix elements (OME) Aij , which are the transition matrix elements in the variable flavor number scheme (VFNS), receive two-mass contributions [1,2]. This also applies to the Wilson coefficients in deeply inelastic scattering. The single mass PS,(3) contributions to the pure singlet (PS) OME, AQq , have been calculated in Ref. [3]. In the * Corresponding author.

E-mail address: [email protected] (J. Blümlein). https://doi.org/10.1016/j.nuclphysb.2017.12.018 0550-3213/© 2017 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Funded by SCOAP3 .

340

J. Ablinger et al. / Nuclear Physics B 927 (2018) 339–367

present paper we present the corresponding 2-mass contributions. They occur first at 3-loop order. Previously, we have calculated already the fixed moments of this OME for the Mellin variable N = 2, 4, 6, to O(η3 ln3 (η)), with η < 1 the mass ratio of the heavy quarks squared in [2,4] using the package Q2E [5,6], as well as the two-mass contributions to the 3-loop OMEs NS,(3) NS−TR,(3) (3) in the flavor non-singlet cases Aqq,Q and Aqq,Q and for Agq,Q , cf. [2], both in N - and in (3)

x-space. The 2-mass contributions to Agg,Q are in preparation [7]. Various other 3-loop single mass contributions have also been completed, cf. [8–17] and all the logarithmic contributions are (3) known to this order [18]. Furthermore, for the OME AQg , all diagrams consisting of contributions that can be obtained in terms of first order factorizing differential or difference equations have been calculated [19–21]. For all OMEs a series of moments has been calculated in the single mass case in Ref. [22] using the code MATAD [23]. From the single pole terms of all the OMEs, we have derived the contributing 3-loop anomalous dimensions, cf. e.g. [3,24]. PS,(3) We perform the calculation of the 2-mass part of the OME AQq mainly in x-space,1 using only some elements of the N -space formalism. In the present case one obtains first order factorizable expressions in x-space, but not in N -space. This implies that the N -space solution cannot be given by sum and product expressions only. The paper is organized as follows. In Section 2 we present the renormalized pure-singlet OME in the 2-mass case. Details of the calculation are given in Section 3. In Section 4 the result of the calculation is presented and numerical results are given in Section 5. Section 6 contains the conclusions. Details on new iterated integrals emerging in the representation, the G-functions, and a series of fixed moments as a function of the mass ratio of the two heavy quarks for the unrenormalized OME are given in the Appendix. 2. The renormalized 2-mass pure singlet OME The generic pole structure for the PS three-loop two-mass contribution is given by [2]   ˆˆ (3),PS 8 (0) (0) 1 1 (0) (1) 4 (0) (0) PS,(1) ˜ AQq = 3 γgq γˆqg β0,Q + 2 2γgq γˆqg β0,Q (L1 + L2 ) + γˆqg γˆgq − β0,Q γˆqq 6 3 3ε ε      1 (0) (0) 1 (0) (1) 2 2 PS,(1) + γgq γˆqg β0,Q L1 + L1 L2 + L2 + γˆ γˆ − β0,Q γˆqq (L2 + L1 ) ε 8 qg gq    1 (2),PS (2),PS (3),PS (0) (2) + γˆ˜qq − 8aQq β0,Q + γˆqg agq + a˜ Qq m21 , m22 , μ2 , (2.1) 3 where we used the short hand notation2 γˆij = γij (NF + 2) − γij (NF ) (2.2) γij (NF + 2) γij (NF ) . (2.3) − γ˜ˆij = NF + 2 NF ˆ (3),PS The tilde in A˜ˆ Qq indicates that we are considering only the genuine two-mass contributions, (l)

and the double hat is used to denote a completely unrenormalized OME. Here the γij ’s are anomalous dimensions at l + 1 loops, β0,Q = − 43 TF , and 1 This strategy has also been used in Ref. [25]. 2 In Eqs. (3.110), (3.111) of [2] unfortunately only the shift N + 1 → N has been used, which we correct here. F F

J. Ablinger et al. / Nuclear Physics B 927 (2018) 339–367



m21 L1 = ln μ2



,

m22 L2 = ln μ2

341

(2.4)

,

where m1 and m2 are the masses of the heavy quarks, and μ is the renormalization scale. Our  (3),PS 2 m1 , m22 , μ2 . goal is to compute the O(ε 0 ) term a˜ Qq In the MS-scheme, renormalizing the heavy masses on-shell, the renormalized expression is given by

1 1 1 1 2 (3),MS,PS (0) (0) = −γgq γˆqg β0,Q L2 L1 + L21 L2 + L31 + L32 A˜ Qq 4 4 3 3    1 (0) (1) 1 PS,(1) + − γˆqg γˆgq + β0,Q γˆqq L22 + L21 16 2   1 (0) (2) 1 (2),PS (0) (0) + 4aQq β0,Q − γˆqg agq − β0,Q ζ2 γgq γˆqg (L1 + L2 ) 2 4   (2),PS (3),PS (0) (2) + 8a Qq β0,Q − γˆqg a gq + a˜ Qq m21 , m22 , μ2 . (2.5) The transition relations for the renormalization of the heavy quarks in the MS-scheme is given in [3], Eq. (5.100), but it only applies to the equal mass case since for the unequal mass case (2),PS (2) and agq represent the first contributions emerge at 3-loop order. In Eqs. (2.1) and (2.5), aQq (2),PS (2) (2),PS (2) and Aˆˆ , respectively, while a and a the O(ε 0 ) terms of the two-loop OMEs Aˆˆ Qq

gq

Qq

gq

represent the corresponding O(ε) terms, cf. Refs. [26–30]. Here and in what follows, ζk = ζ (k), k ∈ N, k ≥ 2 denotes the Riemann ζ -function. 3. Details of the calculation 3.1. The basic formalism There are sixteen irreducible diagrams contributing to A˜ Qq , which are shown in Fig. 1. The unrenormalized operator matrix element is obtained by adding all the diagrams and applying the ˆ ij , quarkonic projector Pq to the corresponding Green function G Q    ij ˆˆ (3) m21 m22 ˆ ij , ˆ ij ≡ δ ( .p)−N Tr p , 2 , ε, N = Pq G (3.1) A˜ Qq /G Q Q 2 μ μ 4Nc (3),PS

where p is the momentum of the on-shell external massless quark (p 2 = 0), is a light-like D-vector, with D = 4 + ε, i and j are the color indices of each external leg, and Nc is the number of colors. The diagrams, D1 , . . . , D16 , are calculated directly within dimensional regularization in D dimensions. These diagrams have the following structure, 2 32 ε m D˜ i (η, ε, N) Di (m1 , m2 , N) = (3.2) μ2 2 32 ε 1 m = dx x N−1 Dˆ i (η, ε, x), (3.3) μ2 0

where N is the Mellin variable appearing in the Feynman rules for the operator insertions, cf. [2,22], m is the mass of the heavy quark where the operator insertion is not sitting, and

342

J. Ablinger et al. / Nuclear Physics B 927 (2018) 339–367

(3),PS Fig. 1. Diagrams for the two-mass contributions to A˜ Qq . The dashed arrow line represents the external massless quarks, while the thick solid arrow line represents a quark of mass m1 , and the thin arrow line a quark of mass m2 . We assume m1 > m2 .

η=

m22 m21

,

(3.4)

with m2 < m1 , i.e. η < 1. In previous publications where single-mass OMEs have been computed [3,8–10,12,15,16,18, 20], and in our recently published two-mass calculations [2], we have always given the results both in N - and x-space. In the case of the two-mass pure singlet OME, finding a general N -space result at three loops, turns out to be rather cumbersome. We will, therefore, present our result only in x-space, see Eq. (3.3), which is anyway all we need in order to obtain the corresponding contribution to the structure function F2 (x, Q2 ) for large values of Q2 , as well as the contribution to the variable flavor number scheme. In most of the applications one finally works in x-space. The diagrams on the first row of Fig. 1 all give the same result, and the diagrams on the second row are related to the diagrams on the first row by the exchange m1 ↔ m2 , i.e.,

J. Ablinger et al. / Nuclear Physics B 927 (2018) 339–367

343

Fig. 2. Massive bubbles appearing in the Feynman diagrams shown in Fig. 1.

Di (m1 , m2 , N) = D1 (m1 , m2 , N)

for

i = 2, 3, 4.

(3.5)

Di (m1 , m2 , N) = D1 (m2 , m1 , N)

for

i = 5, 6, 7, 8.

(3.6)

All of these diagrams vanish for odd values of N . A relation similar to (3.6) holds between the diagrams on the third row of Fig. 1 and those on the fourth row. Furthermore, the diagrams on these two last rows which differ only in the direction of a fermion arrow are related by a factor of (−1)N relative to each other. Specifically one has D10 (m1 , m2 , N) = D9 (m1 , m2 , N).

(3.7)

Di (m1 , m2 , N) = (−1) D9 (m1 , m2 , N) N

Di (m1 , m2 , N) = D9 (m2 , m1 , N)

for

Di (m1 , m2 , N) = (−1) D9 (m2 , m1 , N) N

i = 11, 12.

for

i = 13, 14. i = 15, 16.

for

(3.8) (3.9) (3.10)

We can therefore write the whole unrenormalized pure singlet operator matrix element solely in terms of diagrams 1 and 9:   N A(3),PS D9 (m1 , m2 , N) Qq (N ) = 4D1 (m1 , m2 , N) + 2 1 + (−1)   (3.11) + 4D1 (m2 , m1 , N) + 2 1 + (−1)N D9 (m2 , m1 , N). In Eq. (3.11) the use of the projection (3.1) is implicit. All of the diagrams contain a massive fermion loop with an operator insertion (Figs. 2b1 and 2b2 ) and a massive bubble without the operator (Fig. 2a1 ). The latter can be rendered effectively massless by using a Mellin–Barnes integral [31–35]   4 2 −ε/2 2 Iaμν,ab (k) = g T k − k g k (4π) F μ ν μν s 1 π +i  ∞ 2 σ  ε/2−σ (σ − ε/2) 2 (2 − σ + ε/2) (−σ ) m 2 × dσ , (3.12) −k μ2

(4 − 2σ + ε) −i ∞

where μ and ν are the respective Lorentz indices of the external legs, a and b the color indices, k√is the external momentum, m is the mass of the fermion, which can be either m1 or m2 , gs = 4παs is the strong coupling constant, and TF = 1/2 in SU (Nc ), with Nc the number of colors. Diagram 1 can then be calculated using the following expression for the massive fermion bubble containing the vertex operator insertion (Fig. 2b1 ):  ( .k)N−2 μν,ab Ib1 (k) = 16δab TF gs2

(2 − D/2) D/2 (4π)

0

1

dx x N (1 − x)

( .k) μ kν − k 2 μ ν . (m2 − x(1 − x)k 2 )2−D/2 (3.13)

344

J. Ablinger et al. / Nuclear Physics B 927 (2018) 339–367

Likewise, diagram 9 can be calculated using the following expression for the bubble containing the operator insertion on the fermion line (Fig. 2b2 ): μν,ab

Ib2

(k) = 4δab TF gs2

( .k)N−2 (4π)D/2

1

 dx x N−2 (1 − x)

0

  x 2 (3 − D/2)( .k)2 − 2 x(1 − x)(gμν k 2 − 2kμ kν ) + m2 gμν (m2 − x(1 − x)k 2 )3−D/2 x(kμ ν + kν μ )( .k) + (2 − D/2)(2N x + 1 − N ) 2 (m − x(1 − x)k 2 )2−D/2 xgμν ( .k)2 + (2 − D/2)((N − 1)(1 − 2x) − Dx) 2 (m − x(1 − x)k 2 )2−D/2  μ ν N −1 − (1 − D/2) (N (1 − x) − 1) 2 . 1−x (m − x(1 − x)k 2 )1−D/2

(3.14)

After inserting these expressions, we introduced the corresponding projector for a quarkonic OME given in Eq. (3.1), and performed the Dirac matrix algebra and the trace arising in the numerator of the diagrams using the program FORM [36]. This leads to a linear combination of integrals, the denominators of which were then combined using Feynman parameters. The momentum integrals were then performed, and we obtained expressions where one of the Feynman parameter integrals is already in the form of a Mellin transform. Therefore, we left this Feynman parameter unintegrated, and integrated the remaining ones, after which we obtained the following expression for diagram 1,  

 (3.15) D1 (m1 , m2 , N) = −128CF TF2 1 + (−1)N J1 − J2 , with CF = (Nc2 − 1)/(2Nc ) and CA = Nc , and  J1 =  J2 =

m21 μ2 m21 μ2

3ε 2

3ε 2

(N − 1) 

N + 2ε

1 dx x

N+ 2ε

1+ 2ε

(1 − x)

B1

0

(N ) 

N + 1 + 2ε

1

ε

ε

η , x(1 − x)

dx x N+ 2 (1 − x)1+ 2 B1 0



η . x(1 − x)

(3.16)

(3.17)

Here B1 (ξ ) is the following contour integral, 1 B1 (ξ ) = 2πi

i∞



 3ε ε  2 (σ + 2 − ε) dσ ξ (−σ ) (−σ + ε) σ −

σ− . 2 2 (2σ + 4 − 2ε) σ

−i∞

(3.18) For diagram 9, we get

  1 D9 (m1 , m2 , N) = 64CF TF2 − (ε + 2) − 2J3 + 2ηJ4 + (2N + 2 + ε)J5 4  



 − (N − 1)J6 + N (N − 1) J7 − J8 − (N − 1) J9 − J10 ,

(3.19)

J. Ablinger et al. / Nuclear Physics B 927 (2018) 339–367

345

where  J3 =  J4 =  J5 =  J6 =  J7 =  J8 =  J9 =  J10 =

m21 μ2 m21 μ2 m21 μ2 m21 μ2 m21 μ2 m21 μ2 m21 μ2 m21 μ2

3ε 2

3ε 2

3ε 2

3ε 2

3ε 2

3ε 2

3ε 2

3ε 2

(N + 1) 

N + 2 + 2ε

(N + 1) 

N + 2 + 2ε

(N + 1) 

N + 2 + 2ε

(N + 1) 

N + 2 + 2ε

(N − 1) 

N + 2ε

1

ε

ε

dx x N + 2 (1 − x)1+ 2 B2

dx x

N −1+ 2ε

dx x

N + 2ε



ε 2

(1 − x) B3

0

1

1+ ε2

(1 − x)

B1

0

1

ε

ε

dx x N −1+ 2 (1 − x)1+ 2 B1

η , x(1 − x)

(3.20)

η , x(1 − x)

(3.21)

η , x(1 − x)

(3.22)



0 ε

ε



0

1 dx x

N −1+ 2ε

2+ 2ε

(1 − x)

dx x

N −1+ 2ε

1+ 2ε

B4

(1 − x)

B4

0

1

ε

ε

η , x(1 − x)

η , x(1 − x)

dx x N −1+ 2 (1 − x)1+ 2 B4

0

η , x(1 − x)

η , x(1 − x)

0

1

(N) 

N + 1 + 2ε



0

1

dx x N −1+ 2 (1 − x)2+ 2 B4

(N) 

N + 1 + 2ε

(N − 1) 

N + 2ε

1



η , x(1 − x)

(3.23)

(3.24)

(3.25)

(3.26)

(3.27)

with 1 B2 (ξ ) = 2πi

i∞

 ε  2 (σ + 2 − ε) 3ε

σ +1− , dσ ξ (−σ ) (−σ + ε) σ − 2 2 (2σ + 4 − 2ε)

σ

−i∞

(3.28) B3 (ξ ) =

1 2πi

i∞ −i∞

 ε 3ε

σ +1− dσ ξ σ (−σ ) (−σ − 1 + ε) σ + 1 − 2 2 ×

1 B4 (ξ ) = 2πi

i∞

2 (σ + 3 − ε) ,

(2σ + 6 − 2ε)

(3.29)

 ε  2 (σ + 2 − ε) 3ε

σ −1− . dσ ξ (−σ ) (−σ + ε) σ − 2 2 (2σ + 4 − 2ε)

σ

−i∞

(3.30)

346

J. Ablinger et al. / Nuclear Physics B 927 (2018) 339–367

We want to compute the diagrams, and therefore the integrals (3.18), (3.28)–(3.30), as an expansion in ε up to O(ε 0 ). Notice that there is always a factor consisting on a ratio of -functions depending on N in these integrals, and there are also additional factors of N in Eq. (3.19). After the aforementioned ε expansion is performed, some of the integrals will be left with a factor of the form 1 , with l ∈ {−1, 0, 1}. (3.31) N +l In order to get our results really in terms of a Mellin transform, these factors need to be absorbed into the integral in x, which can be accomplished using integration by parts, 1 N +l

b dx x a

N−1

bN+l f (x) = N +l =

a N+l N +l

b

f (y) dy l+1 − y

b

x dx x

N+l−1

a

a

a

b

b

b

dy

f (y) + y l+1

a

dx x N+l−1 a

dy

f (y) y l+1

(3.32)

dy

f (y) . y l+1

(3.33)

x

We will see later how this is manifested in the final result. 3.2. Computation of the contour integrals We have seen that the diagrams where the operator insertion lies on the heaviest quark are related to the ones where the insertion lies on the lightest internal quark by the exchange m1 ↔ m2 , which means that in order to go from the former diagrams to the latter, we also need to do the change η → 1/η. Therefore, while in the case of diagrams D1 , . . . , D4 and D9 , . . . , D12 , we need to evaluate the contour integrals (3.18), (3.28)–(3.30) with η ξ= , (3.34) x(1 − x) in the case of diagrams D5 , . . . , D8 and D13 , . . . , D16 , we have instead 1 . (3.35) ηx(1 − x) As we will see, the contour integrals arising in the case of Eq. (3.34) will require a different treatment to those in the case of Eq. (3.35). The integrals (3.18), (3.28)–(3.30) can be computed with the help of the Mathematica package MB [37], together with the add-on package MBresolve [38]. These packages allow us to resolve the singularity structure in ε of these integrals by taking residues in σ , after which we can perform the expansion in ε. This leads to expressions of the form ξ=

(ε)

(0)

Bi (ξ ) = Bi (ξ ) + Bi (ξ ), (ε)

i = 1, 2, 3, 4,

(3.36)

where the functions Bi (ξ ) are the sum of residues. For example, in the case of B1 (ξ ), we must take residues at σ = ε/2 and σ = 3ε/2, which leads to  ε   ε  2 2 − ε  (ε) ε/2 2 B1 (ξ ) = ξ (−ε) − 2 2 (4 − ε)





2 2 + 2ε 3ε ε 3ε/2 +ξ

−

− (ε) . (3.37) 2 2

(4 + ε)

J. Ablinger et al. / Nuclear Physics B 927 (2018) 339–367

347

Once the residues are taken, it is possible to expand the original integrals in ε. This leads (0) to the contour integrals that we denote by Bi (ξ ), for which the MBresolve package finds a contour consisting of a line parallel to the imaginary axis with a real value for σ ∈ [−1, 0].3 The functions are given by (0) B1 (ξ ) =

(0) B2 (ξ ) =

(0) B3 (ξ ) =

(0) B4 (ξ ) =

1 2πi 1 2πi 1 2πi 1 2πi

−1/2+i∞ 

dσ ξ σ 2 (−σ ) 2 (σ ) −1/2−i∞

2 (σ + 2) + O(ε),

(2σ + 4)

−1/2+i∞ 

dσ ξ σ 2 (−σ ) (σ ) (σ + 1) −1/2−i∞

2 (σ + 2) + O(ε),

(2σ + 4)

−1/2+i∞ 

dσ ξ σ (−σ ) (−σ − 1) 2 (σ + 1) −1/2−i∞ −1/2+i∞ 

dσ ξ σ 2 (−σ ) (σ ) (σ − 1) −1/2−i∞

2 (σ + 3) + O(ε),

(2σ + 6)

2 (σ + 2) + O(ε).

(2σ + 4)

(3.38)

(3.39)

(3.40)

(3.41)

In order to calculate the contour integrals (3.38)–(3.41), we need to close the contour either to the right or to the left, depending on the convergence of the resulting sum of residues. The integrand always contains a -function in the denominator that can be rewritten using Legendre’s duplication formula,

1 4σ +l

(2σ + 2l) = √ (σ + l) σ + l + , l = 2, 3, (3.42) 2 2 π which means that if we close the contour to the right, the sum of residues will be convergent if and only if ξ < 4, while if we close the contour to the left, the sum of residues will be convergent if and only if ξ > 4. Therefore, in the case where ξ = (ηx(1 − x))−1 , we just need to close the contour to the left, since η < 1, and therefore ξ ≥ 4/η > 4. The calculation in the case where ξ = η/(x(1 − x)), on the other hand, is a bit more complicated, since in this case we need to split the integration range in x into the different regions where ξ can be bigger or smaller than 4, η < 4, for x ∈ (η− , η+ ) , (3.43) x(1 − x) η (3.44) > 4, for x ∈ (0, η− ) or x ∈ (η+ , 1) , x(1 − x) where

  1 (3.45) 1± 1−η . 2 Therefore, in this case, we will need to close the contour to the right when η− < x < η+ , and to the left when 0 < x < η− or η+ < x < 1. η± =

3 Here we have chosen this value to be −1/2, but any value will do, since the poles of the integrands are at integer values of σ .

348

J. Ablinger et al. / Nuclear Physics B 927 (2018) 339–367

Closing the contour to the right and summing residues, we obtain

1 3 ζ2 14 5 2 5 ζ3 82 (0) + ln(ξ ) − ln (ξ ) + ln (ξ ) + ζ2 − + B1 (ξ ) = − 6 27 36 36 18 3 81

∞ 2  2 k (k + 2) 2S1 (2k + 3) − 2S1 (k) + − ln(ξ ) , (3.46) + ξ 2 k (2k + 4) k(k + 1) k=1 ∞ 1 5 ζ2 14  k 2 (k + 2) (0) B2 (ξ ) = − ln2 (ξ ) + ξ ln(ξ ) − − + 2S1 (2k + 3) 12 18 6 27 k (2k + 4) k=1

k−1 − ln(ξ ) , (3.47) − 2S1 (k) − k(k + 1)

∞ 

2 (k + 3) 1 (0) 2S1 (k + 2) − 2S1 (2k + 5) − + ln(ξ ) , (3.48) B3 (ξ ) = ξk (k + 1) (2k + 6) k+1 k=0

1 ζ2 11 4067 1 3 ξ 107ξ (0) + ζ2 + + + ln(ξ ) − ξ+ ln (ξ ) B4 (ξ ) = − 30 9 900 6 27 13500 36

∞ 1

2 (k + 2) ζ3 49  k ξ + ln2 (ξ ) + − + 2S1 (2k + 3) ξ − 60 18 3 81 (k − 1)k 2 (2k + 4) k=2

k 2 + 3k − 2 − ln(ξ ) , (3.49) − 2S1 (k) + (k − 1)k(k + 1) while closing to the left leads to 

∞  1 4 −k (2k − 3) ln(ξ ) + ξ ln(ξ ) 4S1 (k) − 4S1 (2k − 4) − ξ k 2 2 (k − 1) k−1 k=2

1 2 2 2 + ln (ξ ) + 4S1 (k) + 4S1 (2k − 4) − 8S1 (k) S1 (2k − 4) + k−1  2 8 (3.50) S1 (2k − 4) + 2S2 (k) − 4S2 (2k − 4) + + 2ζ2 , + k−1 (k − 1)2 

∞ 1 ln(ξ )  −k (2k − 3) 2(3k − 1) (0) + ln(ξ ) 4S1 (2k − 4) − 4S1 (k) + B2 (ξ ) = − ξ ξ ξ k 2 (k − 1) (k − 1)k

(0) B1 (ξ ) =

k=2

2(3k − 2) − ln (ξ ) − 2S2 (k) − 4S12 (k) − 4S12 (2k − 4) + 4S2 (2k − 4) − (k − 1)2 k 

4(3k − 1) 3k − 1 S1 (2k − 4) + 4S1 (k) 2S1 (2k − 4) + − 2ζ2 , (3.51) − (k − 1)k (k − 1)k ζ2 5 1 14 1 1 (0) − + ln(ξ ) − ln2 (ξ ) − 2 + 2 ln(ξ ) B3 (ξ ) = − 27ξ 6ξ 18ξ 12ξ ξ ξ 

∞  2 −k (k − 2) (2k − 5) ln(ξ ) 4S1 (k − 3) − 4S1 (2k − 6) + ξ +

(k − 2) (k) k−1 k=3

1 − 2S1 (2k − 6) + 2ζ2 + 4S12 (k − 3) + 4S12 (2k − 6) + 4S1 (k − 3) k−1 2

J. Ablinger et al. / Nuclear Physics B 927 (2018) 339–367

349

 4 2 2 + ln (ξ ) , (3.52) S1 (2k − 6) + 2S2 (k − 3) − 4S2 (2k − 6) + k−1 (k − 1)2  ∞ 1 ln(ξ )  −k (k − 1) (2k − 3) (0) ξ + 4S2 (2k − 4) − S2 (k + 1) B4 (ξ ) = − − 4ξ 2ξ k (k − 1) (k + 2) k=2

4 2 − ln (ξ ) + ln(ξ ) 4S1 (2k − 4) − 2S1 (k) − 2S1 (k + 1) + k−1

8 1 2 − 4S1 (2k − 4) − S1 (2k − 4) + 4S1 (k + 1) S1 (2k − 4) + k−1 k−1

4 + S1 (k) 4S1 (2k − 4) − 2S1 (k + 1) + − S12 (k) − S12 (k + 1) k−1  2 − 2ζ2 . (3.53) − S2 (k) − (k − 1)2 −

Here Sa ≡ Sa (N ) denotes the (nested) harmonic sums [39] Sb, a (N ) =

N  (sign(b))k k=1

k |b|

Sa (k), S∅ = 1, b, ai ∈ Z\{0} .

(3.54)

In the above expressions ratios of -functions are related to special binomial coefficients, like

2 (k + 1) 1 1 . =

(2k + 2) 2k 2k k

(3.55)

All of the above sums can be performed using the Mathematica packages Sigma [40,41], HarmonicSums [42–44], EvaluateMultiSums and SumProduction [45]. The results are expressed in terms of generalized iterated integrals. z G ({f1 (τ ), f2 (τ ), · · · , fn (τ )} , z) =

dτ1 f1 (τ1 )G ({f2 (τ ), · · · , fn (τ )} , τ1 ) ,

(3.56)

0

with

 G

 1 1 1 1 , ,··· , , z ≡ lnn (z) . n! τ τ  τ

(3.57)

n times

In principle, the letters in the alphabet of these iterated integrals (i.e., the functions fk (τ )) can be any function (or distribution), for which the iterated integral exists. In the particular case 1 1 and 1+τ , these integrals correspond to the harmonic where the letters are restricted to τ1 , 1−τ polylogarithms [46], which are defined by x Hb, a (x) =

dyfb (y)Ha (y), H∅ = 1, ai , b ∈ {0, 1, −1} ,

(3.58)

0

with f0 (x) =

1 1 1 , f1 (x) = , f−1 (x) = , x 1−x 1+x

(3.59)

350

J. Ablinger et al. / Nuclear Physics B 927 (2018) 339–367

and

 H0, . . . , 0 (x) = G    n times

 1 1 1 1 , ,··· , , x = lnn (x) . τ τ τ n!   

(3.60)

n times

We will see that not only harmonic polylogarithms appear in our final result, but also iterated integrals with square roots in the letters.4 4. The massive operator matrix element We obtain the following expression for the O(ε 0 ) term of the unrenormalized 3-loop two-mass pure singlet operator matrix element 

 (3),PS 2 a˜ Qq (x) = CF TF R0 (m1 , m2 , x) + θ (η− − x) + θ (x − η+ ) x g0 (η, x)  + θ (η+ − x)θ (x − η− ) x f0 (η, x) x −



 y x dy f1 (η, y) + f2 (η, y) + f3 (η, y) x y

η−

η− y x + θ (η− − x) dy g1 (η, y) + g2 (η, y) + g3 (η, y) x y x

x − θ (x − η+ )



y x dy g1 (η, y) + g2 (η, y) + g3 (η, y) x y

η+

1 + x h0 (η, x) +



y x dy h1 (η, y) + h2 (η, y) + h3 (η, y) x y

x

η+ y x + θ (η+ − x) dy f1 (η, y) + f2 (η, y) + f3 (η, y) x y η−

1 +



 y x dy g1 (η, y) + g2 (η, y) + g3 (η, y) . x y

(4.1)

η+

Here θ (z) denotes the Heaviside function  1 z≥0 θ (z) = 0 z < 0.

(4.2)

The function R0 (m1 , m2 , x) arises from the residues taken in order to resolve the singularities in ε of the contour integrals, see Eq. (3.36). The functions fi (η, x), gi (η, x) and hi (η, x), with 4 Root-valued iterated integrals have been discussed before in Ref. [47].

J. Ablinger et al. / Nuclear Physics B 927 (2018) 339–367

351

i = 0, 1, 2, 3, arise from the sum of residues of the contour integrals that remain after the ε expansion, as described in the previous section. The functions with i = 0 are those where no additional factor depending on N needed to be absorbed. The functions with i = 1, i = 2 and i = 3 are those where a factor of 1/N , 1/(N − 1) and 1/(N + 1) was absorbed, respectively, see Eqs. (3.32), (3.33). The different Heaviside θ functions restrict the corresponding values of x to the appropriate regions. Since no contour integral needs to be performed in the case R0 (m1 , m2 , x), the easiest way to compute this function is to integrate in x and then perform the Mellin inversion using HarmonicSums. We obtain,     P 0 3 2 2 3 R0 (m1 , m2 , x) = 32 L1 + L1 L2 + L1 L2 + L2 − 2(x + 1)H0 3x 



2 1 ζ2 P0 H1 2 − + 32 L1 + L2 2(x + 1) H0,0 + H0,1 − 3 3 9x   1 x − 1

2 − (4x + 5)(7x + 5)H0 + 170x + 53x + 80 9 27x    2

x − 1

+ 128L1 L2 56x 2 + 47x + 20 + (x + 1) H0,1 − ζ2 27x 3    4

P0 H1 128ζ3

− x 2 + 7x + 4 H0 − + 64x 3 + 35x 2 − 25x + 8 9 9x 27x  29 4 2 + 64(L1 + L2 ) (x + 1) H0,1 − H0,0,0 + H0,0,1 + 2H0,1,0 9 3 3

 4 8 14 x − 1

− H0,1,1 − ζ2 H0 + ζ3 + 260x 2 + 231x + 116 3 3 3 27x

 P0 2 2ζ2 3 + H1,1 − H1,0 + 6x − 10x 2 − 19x − 6 3x 3 9x   4(x − 1) 2 1

168x 2 + 265x + 229 H0 − 5x + 23x + 5 H1 − 27 27x    2 2 64ζ2

+ 6x + 4x − 5 H0,0 + 282x 3 − 229x 2 − 85x − 120 9 81x

64P0 2 + 4H1,0,0 − 2H1,0,1 − 2H1,1,0 − H1,1,1 + ζ2 H1 9x 3 2 8 4 4 + 128(x + 1) (6x − 5)H0,1,0 + H0,0,0,0 − H0,0,0,1 − H0,0,1,0 9 9 9 3 2 4 2 2 2 + H0,0,1,1 − H0,1,0,0 + H0,1,0,1 − 2ζ3 H0 + H0,1,1,0 + H0,1,1,1 9 3 3 3 9  7 1 8 2 128

+ ζ2 H0,0 − ζ2 H0,1 + ζ2 − 12x 2 + 19x + 19 H0,1,1 9 3 15 27

 64

 256

− 813x 2 + 29x + 263 + ζ2 60x 2 + 91x + 37 H0 243 27

352

J. Ablinger et al. / Nuclear Physics B 927 (2018) 339–367

  128(x − 1)

256

22x 2 − 25x + 4 H1 + 84x 2 + 109x + 100 H0,0 81x 81   256 2 128(x − 1)

+ 6x − 5x − 5 H0,0,1 − 56x 2 − 43x + 20 H1,0 27 27x   128(x − 1)

256

2 + 40x + 49x + 40 H1,1 − 12x 2 − x − 10 H0,0,0 81x 27  128 256(x − 1)

+ (47x + 29)H0,1 + 2602x 2 − 203x + 1360 , (4.3) 81 729x +

where P0 = (x − 1)(4 + 7x + 4x 2 ),

(4.4)

L1 and L2 are the logarithms defined in Eq. (2.4), and we used the shorthand notation Ha (x) ≡ Ha . In principle (4.3) could still be reduced to a shorter basis using shuffle-algebra [48]. The fi (η, x) functions, which are defined in the range η− < x < η+ , are given by 

3/2

8P3 4x(1 − x) − η 16(x − 1) η η f0 (η, x) = K1 − K2 x(1 − x) 3x x(1 − x) 45η3/2 (x − 1)x 3 

 

 2 6η + 30x 2 − 5x P1 η − 2ζ2 + ln2 + 15x x(1 − x) 90(x − 1)3 x 5

2P2 η + ln , (4.5) x(1 − x) 45(x − 1)3 x 5

3/2

16P6 4x(1 − x) − η η 4P5 η f1 (η, x) = − K ln + 1 x(1 − x) x(1 − x) 45η3/2 (x − 1)x 3 45(x − 1)3 x 5

  64P7 η η 32(x − 1)(4x + 1) 2 + ln K2 − 2ζ + 2 3x x(1 − x) 45x 2 x(1 − x) P4 − , (4.6) 45(x − 1)3 x 5 

3/2

64P10 4x(1 − x) − η 128 η η K1 − (x − 1) K2 f2 (η, x) = 9η3/2 (x − 1)x 2 x(1 − x) 3 x(1 − x) 

  η 20 η 4 10 2 ln ln − ζ2 − − P 4P − 9 8 , 3 x(1 − x) 3 x(1 − x) 9(x − 1)3 x 4 (4.7)  √

16P13 4x(1 − x) − η 32(x − 1) η η f3 (η, x) = − K1 − K2 3/2 3 x(1 − x) 3x x(1 − x) 9η (x − 1)x 



2 η 4P12 η 2 − (3x + 5) 2ζ2 + ln ln + 3 x(1 − x) x(1 − x) 9(x − 1)3 x 5 −

P11 , 27(x − 1)3 x 5

(4.8)

where the functions K1 and K2 , which appear repeatedly in the expressions above, are given in terms of the iterated integrals defined in Eq. (3.56),

J. Ablinger et al. / Nuclear Physics B 927 (2018) 339–367

353



  √ √ √   1 √ 1

, 4 − τ τ , u + 1 − 2 ln(u) G 4 − τ τ ,u , (4.9) τ 2   √ √   √ 1 √ 2 K2 (u) = −G 4 − τ τ ,u G , 4 − τ τ , u + ln3 (u) τ 3   √ √ √ √ 1 , 4 − τ τ , 4 − τ τ , u + 4ζ2 ln(u) + 8ζ3 +G τ  √ √   1

− 1 − 2 ln(u) G2 4 − τ τ ,u . (4.10) 4 The expressions of the G-functions are presented in the Appendix in terms of harmonic polylogarithms containing square-root valued arguments, and the Pi ’s, with i = 1, . . . , 13, are polynomials in η and x given by K1 (u) = G

P1 = 1536(x − 1)4 (3x + 2)x 4 + 576(x − 1)3 (12x − 7)ηx 3 + 8(x − 1)2 (264x − 329)η2 x 2 + 16(x − 1)(12x − 37)η3 x − 45η4 , P2 = 128(x − 1) (3x − 8)x − 32(x − 1) (33x − 8)ηx 4

4

3

(4.11)

3

− 4(x − 1)2 (108x − 133)η2 x 2 − 24(x − 1)(2x − 7)η3 x + 15η4 ,

(4.12)

P3 = 4(x − 1) (6x − 1)x − 6(x − 1)(4x + 1)ηx + 15η , 2

2

2

P4 = 768(x − 1) (40x + 7)x + 576(x − 1) (20x − 1)ηx 4

4

3

(4.13)

3

− 8(x − 1)2 (260x + 197)η2 x 2 − 16(x − 1)(100x + 31)η3 x − 45(4x + 1)η4 , (4.14) P5 = 64(x − 1) (40x + 13)x + 16(x − 1) (200x + 17)ηx 4

4

3

3

− 4(x − 1)2 (100x + 79)η2 x 2 − 48(x − 1)(10x + 3)η3 x − 15(4x + 1)η4 ,

(4.15)

P6 = 8(x − 1)2 (10x + 1)x 2 − 6(x − 1)(20x + 3)ηx + 15(4x + 1)η2 ,

(4.16)

P7 = 10(x − 1)x(10x + 1) − 3η,

(4.17)

P8 = 1536(x − 1) x + 576(x − 1) ηx − 104(x − 1) η x 4 4

3

3

2 2 2

− 80(x − 1)η3 x − 9η4 ,

(4.18)

P9 = 128(x − 1) x + 160(x − 1) ηx − 20(x − 1) η x − 24(x − 1)η x − 3η , 4 4

3

3

2 2 2

3

4

(4.19) P10 = 4(x − 1)2 x 2 − 6(x − 1)ηx + 3η2 ,

(4.20)

P11 = 512(x − 1) (7x − 9)x − 1728(x − 1) (2x + 1)ηx 4

4

3

3

− 24(x − 1)2 (24x − 13)η2 x 2 + 240(x − 1)η3 x + 27η4 , P12 = 32(x − 1) (11x − 4)x − 32(x − 1) (6x + 5)ηx 4

4

3

(4.21)

3

− 4(x − 1)2 (12x − 5)η2 x 2 + 24(x − 1)η3 x + 3η4 ,

(4.22)

P13 = 16(x − 1) (3x + 1)x − 4(x − 1) (6x + 5)ηx + 6(x − 1)η x + 3η . 3

3

2

2

2

3

(4.23)

The gi (η, x) functions, defined in the ranges 0 < x < η− and η+ < x < 1, are given by 

3/2 x(1 − x) x(1 − x) 64 x − 1 64P15

+ K η − 4x(1 − x) g0 (η, x) = K 3 4 x η 3 η 45η3/2 x

354

J. Ablinger et al. / Nuclear Physics B 927 (2018) 339–367

 x(1 − x) 64ζ2 6η + 30x 2 − 35x  32  2 6η + 30x − 5x ln2 + + 45x η 45x 

  128P14 x(1 − x) 256(x − 1) 2 − ln + 3η + 24x − 34x , (4.24) 45η2 x η 45η

3/2 x(1 − x) 64P7 2 x(1 − x) 128P17

− η − 4x(1 − x) K3 ln g1 (η, x) = − 45η3/2 x 2 η 45x 2 η

x(1 − x) 256P16 128(x − 1)(4x + 1) x(1 − x) + − ln K4 3x η 45η2 x 2 η

   3η + 20x 3 − 20x 128ζ 256(x − 1) 2 3 2 3η + 80x − 36x − 44x + + , 45ηx 45x 2 (4.25) 

3/2

2 η − 4x(1 − x) x(1 − x) x(1 − x) 256 + 2K4 K3 g2 (η, x) = (x − 1) 3/2 3 3η η η

  2 16(x − 1)x 2ζ2 4 x(1 − x) − − 2 −η + 4x 2 − 4x ln − η 3η 3 3η 

x(1 − x) 5 + ln2 , (4.26) 3 η 

x(1 − x) x(1 − x) 2 64(x − 1) 2P19  + K4 η − 4x(1 − x)K3 g3 (η, x) = x η 3 η 9η3/2

 2 2P18 8x  x(1 − x) − + −7η + 36x 2 − 42x + 6 + (3x − 1)ζ2 ln 2 9η η 27η 9

 x(1 − x) 1 + (3x + 5) ln2 . (4.27) 9 η Here the functions K3 and K4 are  √   √ 

 1 1 − 4τ 1 − 4τ , K3 (u) = G , u − ln(u) + 2 G , u + ζ2 , τ τ τ  √  √   √ 1 − 4τ 1 − 4τ 1 1 − 4τ , , K4 (u) = −G , u + ζ2 G ,u τ τ τ τ  √  √ 1 − 4τ 1 − 4τ 1 , −2G , u + ζ2 ln(u) + ln3 (u), τ τ 6

(4.28)

(4.29)

and P14 = 16(x − 1)2 (6x − 1)x 3 − 8(x − 1)(9x − 4)ηx 2 − (36x − 41)η2 x − 6η3 ,

(4.30)

P15 = x(6x − 1) − 6η,

(4.31)

P16 = 32(x − 1) (10x + 1)x − 4(x − 1) (40x + 1)ηx 3

3

+ (x − 1)(20x + 23)η2 x + 3η3 , P17 = 2(x − 1)x(10x + 1) + 3η,

2

2

(4.32) (4.33)

J. Ablinger et al. / Nuclear Physics B 927 (2018) 339–367

PS,(3)

355

PS,(3)

of O(TF2 ) 2 2 2 2 2 as a function of x and μ . Dotted line (red): μ = 30 GeV . Dashed line (black): μ = 50 GeV . Dash-dotted line (blue): μ2 = 100 GeV2 . Full line (green): μ2 = 1000 GeV2 . Here the on-shell heavy quark masses mc = 1.59 GeV and Fig. 3. The ratio of the 2-mass (tm) contributions to the massive OME AQq

to all contributions to AQq

mb = 4.78 GeV [49,50] have been used.

P18 = 32(x − 1)2 (3x + 1)x 2 − 16(x − 1)(3x + 1)ηx + (2 − 7x)η2 ,

(4.34)

P19 = 4(x − 1)x(3x + 1) + (1 − 6x)η.

(4.35)

Finally, the hi (η, x) functions, defined in the full range 0 < x < 1, are just given by the gi (η, x) functions with η → 1/η, i.e.,

1 hi (η, x) = gi , x , i = 0, 1, 2, 3. (4.36) η We see that iterated integrals of up to weight three appear in our result. The alphabet of these integrals is given in terms of just three letters: √ √ √ 1 1 − 4τ 4 − τ τ, , . (4.37) τ τ In principle, we could try to calculate all of the integrals in y appearing in Eq. (4.1) and express them in terms of iterated integrals of higher weight. However, this is not really necessary or even convenient, since the expressions (4.5)–(4.8), (4.24)–(4.25) are very compact, and integrating them into higher weight iterated integrals leads to a result of considerably larger size. Furthermore, all of the iterated integrals appearing above can be written in terms of simple polylogarithms (albeit of complicated arguments), see the Appendix, for which various fastly converging numerical representations exist. Therefore, the integrals in y appearing in Eq. (4.1) can be performed numerically without problems. The convolution with parton distribution functions (which will be presented in a future phenomenological publication), in order to compute the corresponding contribution to F2 (x, Q2 ) or for the transition rate in the VNFS, is straightforward. 5. Numerical results We compare the pure singlet 2-mass contributions to the complete O(TF2 CA,F ) term as a function of x and μ2 in Fig. 3. Typical virtualities are μ2 ∈ [30, 1000] GeV2 . The ratio of the

356

J. Ablinger et al. / Nuclear Physics B 927 (2018) 339–367

2-mass contributions to the complete term of O(TF2 CA,F ) grows in this region from slightly negative contributions to ∼ 0.36 for very large virtualities in most of the x-range. The behavior of the ratio is widely flat in x, rising at very large x. 6. Conclusions PS,(3)

We have calculated the two-mass 3-loop contributions to the massive OME AQq in analytic form in x-space for a general mass-ratio η. It contributes to the 3-loop VFNS and to the massive 3-loop corrections of the deep-inelastic structure function F2 (x, Q2 ) in the region m2 Q2 . As a function of x its relative contribution to the O(TF2 CA,F ) terms of the whole matrix element PS,(3) AQq behaves widely flat and grows with the scale μ2 up to about ∼ 0.36. We used Mellin–Barnes techniques to obtain the x-space result by factoring out the N -dependence in terms of the kernel x N , and used integration by parts to absorb the N -dependent polynomial pre-factors. The result can be written as single limited integrals within the range x ∈ [0, 1] over iterated integrals containing also square-root valued letters. These integrals can be turned into polylogarithms of involved root-valued arguments depending on η. The even Mellin moments of the OME exhibit a growing number of polynomial terms in η with growing values of N . Due to this structural property and the arbitrariness of η, which enters the ground field, the method of arbitrarily large moments [51] cannot be used to find the result in the present case. Acknowledgements This work was supported in part by the Austrian Science Fund (FWF) grant SFB F50 (F5009-N15), the European Commission through contract PITN-GA-2012-316704 (HIGGSTOOLS). Appendix A. Evaluating the G-functions A.1. G-functions with support 0 < x < 1 Here, we have the argument ξ1 = x(1 − x)η ∈ (0, η/4) √ √ and therefore 4ξ1 − 1 = i 1 − 4ξ1 . We also define  ω1 = 1 − 4ξ1 . We obtain:

 √  

1 − 4τ G , ξ1 = 2ω1 + 2 ln 1 − ω1 − ln(4ξ1 ) − 2 τ

4ξ1 2 1 G ({τ |1 − 4τ |} , ξ1 ) = ξ1 − 2 3

   1 G τ 2 |1 − 4τ | , ξ1 = ξ13 − ξ1 3

(A.1)

(A.2)

(A.3) (A.4) (A.5)

J. Ablinger et al. / Nuclear Physics B 927 (2018) 339–367



 √

 1

 1 − 4τ G , ξ1 = 4ω1 − ln2 1 − ω1 − ln2 4ξ1 τ 2 



+ 4 ln 1 − ω1 − 2 ln(2) ln 1 − ω1



  − 4 ln 4ξ1 + 2 ln 1 − ω1 ln 4ξ1

1 − ω1 − 4 + ln2 (2) + 4 ln(2) + 2Li2 2  √ 

 1 − 4τ 1 , G , ξ1 = ω1 −4 − 4 ln(2) + 2 ln 4ξ1 τ τ

 √ G

357

1 , τ

(A.6)



1 + ln2 (1 − ω1 ) − ln2 (4ξ1 ) − 4 ln 1 − ω1 2 





− 2 ln(2) ln 1 − ω1 + 2 ln 4ξ1 + 2 ln(2) ln 4ξ1

1 − ω1 (A.7) − 2Li2 + 4 − ln2 (2) 2

1 − 4τ , τ

 √





 1 − 4τ , ξ1 = ω1 4 ln 1 − ω1 − 4 − 2 ln 4ξ1 − 8ξ1 τ

1 2 ln (4ξ1 ) − 4 ln (1 − ω1 ) 2 + 2 ln (4ξ1 ) − 2 ln (1 − ω1 ) ln (4ξ1 ) + 4 (A.8)   √

ω1 + 1 1 − 4τ 1 1 , , − 4 ln (1 − ω1 ) ζ2 + 2 ln (4ξ1 ) ζ2 G , ξ1 = 2Li3 τ τ τ ω1 − 1 + 2 ln2 (1 − ω1 ) +

1 3 ln (4ξ1 ) + 2 ln2 (2) 6 × ln (1 − ω1 ) − 2 ln2 (1 − ω1 ) + 2 ln(2)

+ 8ω1 − 2 ln3 (1 − ω1 ) −

× ln2 (1 − ω1 ) + 2 ln2 (1 − ω1 ) ln (4ξ1 ) − 2 ln2 (4ξ1 ) + ln(2) ln2 (4ξ1 ) + 8 ln (1 − ω1 ) − 4 ln(2) × ln (1 − ω1 ) − 8 ln (4ξ1 ) + 4 ln(2) ln (4ξ1 ) + 4 ln (1 − ω1 ) ln (4ξ1 ) − 4 ln(2) ln (1 − ω1 ) ln (4ξ1 )

1 − ω1 1 − ω1 2 + 4Li3 − 8 − ln3 (2) + 4Li2 2 2 3  G

1 , τ



√ 1 − 4τ 1 , , ξ1 τ τ



− 2 ln2 (2) + 8 ln(2) (A.9)   = −16 − 4 ln (4ξ1 ) + 2 ln2 (4ξ1 ) + 8ζ2 ln (1 − ω1 )  + ln2 (1 − ω1 ) (4 − 5 ln (4ξ1 )) + 12

 1 − ω1 + 4ω1 − 4ζ2 ln (4ξ1 ) +2Li2 2

 1 − ω1 1 − ω1 −4 −4 + 2Li2 + 2Li3 2 2

358

J. Ablinger et al. / Nuclear Physics B 927 (2018) 339–367

+Li3



  ω1 + 1 1 − ω1 + 4ω1 + −8 − 4Li2 ω1 − 1 2

−8ω1 − 2 ln2 (1 − ω1 ) + 4 ln (4ξ1 ) + 2 ln (1 − ω1 )  × ln (4ξ1 ) − ln2 (4ξ1 ) ln(2) + 4 ln3 (1 − ω1 ) 1 2 − ln3 (4ξ1 ) + ln2 (2) (ln (4ξ1 ) − 4) − ln3 (2) 6 3 (A.10)  √ √    1 1 − 4τ 1 − 4τ G , , , ξ1 = − ln (1 − ω1 ) (2 ln (4ξ1 ) − 4) − ln2 (4ξ1 ) ln(2) τ τ τ − (2 − 2 ln (1 − ω1 ) + ln (4ξ1 )) ln2 (2) − ln2 (1 − ω1 ) (−5 ln (4ξ1 ) − 6)

 1 − ω1 − 8ω1 + 8 ln (4ξ1 ) − 8 − 4Li2 2  + 2 ln2 (4ξ1 ) + 4ζ2 ln (1 − ω1 )

1 − ω1 + 4ω1 − 2ζ2 ln (4ξ1 ) − −4 + 2Li2 2

1 + ω1 1 − ω1 + 2Li2 − 2 − 4 − Li3 ω1 − 1 2

10 3 ln (1 − ω1 ) + 4ω1 + 4ξ1 − 3 1 + ln3 (4ξ1 ) + 2 ln2 (4ξ1 ) (A.11) 6  √ 

 1 − ω1 1 − ω1 1 − 4τ 1 1 , , + 2Li3 G , ξ1 = 2 −4 + 2Li2 τ τ τ 2 2

  1 + ω1 1 − ω1 + 4ω1 + 4Li2 + Li3 ω1 − 1 2 + 8ω1 + ln (1 − ω1 ) (4 − 2 ln (4ξ1 ))  − 4 (1 + ω1 ) ln (4ξ1 ) + 2 ln2 (4ξ1 ) ln(2) + (2 (2ω1 + 1) + 2 ln (1 − ω1 ) − 3 ln (4ξ1 )) ln2 (2)   + 8 − ln2 (4ξ1 ) − 4ζ2 ln (1 − ω1 ) + ln2 (1 − ω1 )

   1 − ω1 × 3 ln (4ξ1 ) − 2 + −4 − 2Li2 2  4 3 − 4ω1 + 2ζ2 ln (4ξ1 ) + ln (2) − 2 ln3 (1 − ω1 ) 3 1 + (1 + ω1 ) ln2 (4ξ1 ) − ln3 (4ξ1 ) (A.12) 6  √  √

 ω1 + 1 1 − ω1 1 − 4τ 1 1 − 4τ , , + Li2 G , ξ1 = 4 −Li3 τ τ τ ω1 − 1 2

J. Ablinger et al. / Nuclear Physics B 927 (2018) 339–367

359



 

1 − ω1 ω1 − 4ξ1 − 4 (ω1 − 1) 2  − 6 ln (4ξ1 ) ln (1 − ω1 ) − 8 (ω1 − 1)  + 4 ln2 (1 − ω1 ) + 4 ln (4ξ1 ) + 2 ln2 (4ξ1 ) ln(2)

+ Li2

− (−2 (1 + ω1 ) + 2 ln (1 − ω1 ) − ln (4ξ1 )) ln2 (2)

 1 − ω1 − − ln2 (4ξ1 ) + 4Li2 − 8ζ2 2   − 4 ln (4ξ1 ) ω1 ln (1 − ω1 ) − 2 (1 + ω1 )

  1 − ω1 + 5 ln (4ξ1 ) ln2 (1 − ω1 ) − −4 − 2Li2 2  14 3 + 4ω1 + 4ζ2 ln (4ξ1 ) + ln (1 − ω1 ) 3 1 − (ω1 − 1) ln2 (4ξ1 ) + ln3 (4ξ1 ) (A.13) 6   √ √    1 − 4τ 1 − 4τ 1 , , , ξ1 = 8 − 4 ln2 (2) − 2ζ2 − 4 ln(2)ω1 ln 4ξ1 G τ τ τ 



 + −4 1 + 2ξ1 + ω1 + 6 ln2 (2) − 4 ln 4ξ1



 

+ ln2 4ξ1 + 4 ln2 ω1 + 1 + 4ζ2 + ln 1 + ω1   



× − 4ω1 − 8 ln(2) + 4 ln(2) 1 + ω1  



× ln 1 − ω1 + 4 − 3 − 2ξ1 + ω1 + 10 ln2 (2)

  

+ 4 ln(2) 2ω1 + 1 ln 1 + ω1 + −4 ln(2)  1 − ω 

1 + 4ζ2 + 4ζ3 + 4 ln ω1 + 1 − 4ω1 Li2 2 

 

 + 8 3ξ1 + ω1 − 1 + 2 − 2 ln 4ξ1 − 2 ln(2) 



 × ln2 1 − ω1 + ln2 4ξ1 ω1 + 2 ln(2) + 1

 1 − ω1 1 − ln3 4ξ1 + 16 ln(2)ξ1 − 4Li2 2 2

ω1 + 1 1 − ω1 − 4Li3 − 2Li3 ω1 − 1 2

 

 1 + ω1 + 4 ln 1 + ω1 − 4 ln(2) − 4 − 4Li3 2

 8 1 + ω1 2 × Li2 + ln3 1 + ω1 − ln3 (2) 2 3 3

360

J. Ablinger et al. / Nuclear Physics B 927 (2018) 339–367





  − 2 1 + ω1 + 6 ln(2) ln2 ω1 + 1 − 2 ln2 (2)



  × ω1 + 2 + 2 ln3 1 − ω1 . (A.14) A.2. G-functions with support x ∈ [0, 12 −

1√ 2 1 − η]

∪ [ 12 +

1√ 2 1 − η, 1]

The integrals in this class coincide with integrals with full support (0 < x < 1), but the replacements x(1 − x) , η  ω1 → ω3 = 1 − 4ξ3 ξ1 → ξ 3 =

(A.15) (A.16)

have to be performed. A.3. G-functions with support

1 2

−

1√ 2 1−η
<

1 2

+

1√ 2 1−η

Here, we have the argument η ξ2 = ∈ (4η, 4) (A.17) x(1 − x) √ and therefore 4 − ξ2 is real. We introduce the abbreviations   ω2 = ξ2 4 − ξ2 , (A.18) √ ξ2 φ = arcsin . (A.19) 2

√ √   ξ2 4 − τ τ , ξ2 = −ω2 1 − G + 4φ (A.20) 2 

   √ 1 ξ2 1 √ , 4 − τ τ , ξ2 = − ω2 1 + + 2 − 4 ln(4 − ξ2 ) φ G 4−τ 2 2 + 4Cl2 (2φ) − 2Cl2 (4φ) (A.21) 

   √ √ ξ2 3 1 , 4 − τ τ , ξ2 = ω 2 − − 2 − 4 ln(ξ2 ) φ + 4Cl2 (2φ) G τ 4 2 (A.22)   √ √ √ 1 √ 3 2 1 3 1 4 G , 4 − τ τ , 4 − τ τ , ξ2 = 5ξ2 − ξ2 + ξ2 − ξ2 − 8ζ3 τ 2 3 32   + φ (ξ2 − 6) ω2 + 16Cl2 (2φ) + 8Cl3 (2φ) + φ 2 (8 ln(ξ2 ) − 4) . We used the Clausen function [52,53]  i Cl2 (x) = Li2 (e−ix ) − Li2 (eix ) , 2  1 Cl3 (x) = Li3 (e−ix ) + Li3 (eix ) 2 with the sum representation

(A.23)

(A.24)

J. Ablinger et al. / Nuclear Physics B 927 (2018) 339–367

Cl2 (φ) = Cl3 (φ) =

∞  sin (nφ) n=1 ∞  n=1

n2

361

, (A.25)

cos (nφ) n3

for φ ∈ (0, 2π). Appendix B. Fixed moments For fixed values of N = 2k, k ∈ N\{0}, we find the following moments: PS,(3)

(N = 2) =  1  75136 512   8192 1  11776 2048 2 2 + − L − − + + L + L − − ) (L 2 1 2 1 243 27 ε 729 9 81ε 3 ε2  125600 2944 1024 512 1024 2 H0 (η) (L2 − L1 ) − H (η) − ζ2 − − (L2 + L1 ) − 81 9 27 0 27 2187  736   128   18784 256  256  3 L2 + L31 − L22 + L21 − H0 (η) L22 − L21 − + H02 (η) − 9 27 3 243 9

256  640 496 256 ζ2 (L2 + L1 ) − H0 (η) (L2 − L1 ) + − + H1 (η) H02 (η) + 9 27 81 27 128 3 512 512 1472 1024 H0 (η) − H0 (η)H0,1 (η) + H0,0,1 (η) − ζ2 + ζ3 + 27 27 81 81

81 320 40 2 1 1 160 + + H0 (η) η + + H0 (η) η − 27 27 η 27 η 

√ 

√ 

√  

+ 32H0,−1 η + 8H0,1 (η) H0 (η) − 4 H−1 η + H1 η H02 (η)

 √

√  1 − 64H0,0,−1 η − 8H0,0,1 (η) η+ √ η 

√  40

√ 

√  160 20

H0,−1 η + H0,1 (η) H0 (η) − H−1 η + H1 η H02 (η) + 27 27 27 

√  40 1 320 3/2 (B.1) H0,0,−1 η − H0,0,1 (η) η + 3/2 − 27 27 η

AQq

PS,(3)

AQq

(N = 4) =

 1  6600284 1936   30976 1  17888 7744 − − L22 + L21 + 2 − (L2 + L1 ) + − 3 2025ε ε 6075 675 ε 455625 225 3872  24497203 4472 1936 3872 2 H0 (η) (L2 − L1 ) − H (η) − ζ2 − − (L2 + L1 ) − 2025 225 675 0 675 5467500       1650071 1118 2 484 968 3 L2 + L31 − L2 + L21 − − H0 (η) L22 − L21 − 225 675 75 151875  968 968ζ2 4294 + H0 (η)2 + H0 (η) (L2 − L1 ) (L2 + L1 ) − 225 225 3375 −

362

J. Ablinger et al. / Nuclear Physics B 927 (2018) 339–367



186109 968 484 3 1936 + H1 (η) H02 (η) + H0 (η) − H0 (η)H0,1 (η) 324000 675 2025 675

1936 2236 3872 1322 1273 2 + H0,0,1 (η) − ζ2 + ζ3 + + H0 (η) 675 2025 2025 675 5400



5239 49 1 1 49 + + − × η+ H0 (η) η − − H02 (η) η 5400 η 200 1600





√  39 39 49 1 1 H0 (η) η2 − 2 + H0,−1 η + H0,1 (η) H0 (η) × η2 + 2 − η 400 η 8 32 

√ 

√  2

√  39 39

39 1 √ − H−1 η + H1 η H0 (η) − H0,0,−1 η − H0,0,1 (η) η+ √ 64 4 32 η 

√ 

√ 

√  425 425 425

H0,−1 η + H0,1 (η) H0 (η) − H−1 η + H1 η H02 (η) + 432 1728 3456 



√ 

√  425 425 1 49 − H0,0,−1 η − H0,0,1 (η) η3/2 + 3/2 + − H0,−1 η 216 1728 400 η 

√ 

√ 

√  49 49

49 − H0,1 (η) H0 (η) + H−1 η + H1 η H02 (η) + H0,0,−1 η 1600 3200 200

 1 49 5/2 (B.2) + H0,0,1 (η) η + 5/2 1600 η

+

PS,(3)

AQq

(N = 6) =

 1  257649488 1  121472 30976 123904 − + 2 − (L2 + L1 ) + − 3 19845ε ε 297675 6615 ε 43758225  30368 7744  2 7744 15488 2 L2 + L21 − H0 (η) (L2 − L1 ) − H (η) − (L2 + L1 ) − 2205 99225 2205 6615 0   15488  18655921961 3872  3 7592  2 − ζ2 − − L2 + L31 − L2 + L21 6615 17503290000 2205 33075   64412372 3872  3872  1936 H0 (η) L22 − L21 + − − H02 (η) − ζ2 (L2 + L1 ) − 735 14586075 2205 2205

2312 3872 1936 3 78873 − H0 (η) (L2 − L1 ) + + H1 (η) H02 (η) + H (η) 9261 219520 6615 19845 0 7744 15184 15488 7744 H0 (η)H0,1 (η) + H0,0,1 (η) − ζ2 + ζ3 − 6615 6615 99225 19845



342121 2 603709 5441 27687011 1 1 + H0 (η) η + + H0 (η) η − + − + 31752000 3386880 η 1411200 η 23520



 1349 81 1 1 5261 81 − H02 (η) η2 + 2 − H0 (η) η2 − 2 + + H02 (η) 188160 η 11760 η 3136 25088





√  26939 26939 81 1 1 3 3 H0 (η) η − 3 + H0,−1 η + H0,1 (η) × η + 3 + 6272 13440 53760 η η −

× H0 (η) −

√ 

√ 

√  26939

26939 H−1 η + H1 η H02 (η) − H0,0,−1 η 107520 6720

J. Ablinger et al. / Nuclear Physics B 927 (2018) 339–367

363





√  10649 10649 26939 1 √ H0,0,1 (η) + H0,−1 η + H0,1 (η) H0 (η) − η+ √ 53760 η 24192 96768 

√ 

√  2

√  10649 10649

10649 H−1 η + H1 η H0 (η) − − H0,0,−1 η − H0,0,1 (η) 193536 12096 96768

 

√  223 1 223 3/2 × η + 3/2 + − H0,−1 η − H0,1 (η) H0 (η) η 1920 7680

√ 

√ 

√  223

223 H−1 η + H1 η H02 (η) + H0,0,−1 η 15360 960 



√  81 1 81 223 5/2 + H0,0,1 (η) η + 5/2 + H0,−1 η + H0,1 (η) H0 (η) 7680 6272 25088 η 

√ 

√  2

√  81

81 81 H−1 η + H1 η H0 (η) − H0,0,−1 η − H0,0,1 (η) − 50176 3136 25088

1 7/2 (B.3) × η + 7/2 η +

PS,(3)

AQq

(N = 8) =

 1  33262473901 1  4916 87616 21904 + 2 − (L2 + L1 ) + − 3 25515ε ε 107163 8505 ε 10126903500  1229 5476  2 5476 10952 2 L + L21 + − H0 (η) (L2 − L1 ) − H (η) (L2 + L1 ) − 2835 2 35721 2835 8505 0      2738 3 8273033473567 1229 10952 ζ2 − − L2 + L31 + L2 + L21 − 8505 27221116608000 2835 47628 2   33262473901 2738  1369 2738  H0 (η) L22 − L21 + − − H0 (η)2 − ζ2 (L2 + L1 ) − 945 13502538000 2835 2835   328686091 2738 40333 1369 3 H0 (η) (L2 − L1 ) + + H1 (η) H02 (η) + H (η) − 510300 1567641600 8505 25515 0 5476 5476 1229 10952 H0 (η)H0,1 (η) + H0,0,1 (η) + ζ2 + ζ3 − 8505 8505 71442 25515

 171113081  4243147 2 4720627 1 1 + H0 (η) η + + H0 (η) η − + 304819200 69672960 η 17418240 η

 30598577 6036587 1158389 2  2 1 1 − H0 (η) η + 2 − H0 (η) η2 − 2 + − 108864000 34836480 η 43545600 η

 2487251  825131 271091 1 1 + + H 2 (η) η3 + 3 + H0 (η) η3 − 3 47029248 41803776 0 31352832 η η

  847 847 1 1 847 + − − H 2 (η) η4 + 4 − H0 (η) η4 − 4 248832 1990656 0 η 497664 η 

√ 

√  48113 48113 48113

+ H0,−1 η + H0,1 (η) H0 (η) − H−1 η 43008 172032 344064 

√  2

√  48113 48113 1 √ + H1 η H0 (η) − H0,0,−1 η − H0,0,1 (η) η+ √ 21504 172032 η −

364

J. Ablinger et al. / Nuclear Physics B 927 (2018) 339–367

 331775 

√ 

√  331775 331775

H−1 η H0,−1 η + H0,1 (η) H0 (η) − 1161216 4644864 9289728



√  2

√  331775 331775 1 3/2 + H1 η H0 (η) − H0,0,−1 η − H0,0,1 (η) η + 3/2 580608 4644864 η 

√ 

√  1449 1449 1449

+ − H0,−1 η − H0,1 (η) H0 (η) + H−1 η 10240 40960 81920 

√  2

√  1449 1449 1 H0,0,−1 η + H0,0,1 (η) η5/2 + 5/2 + H1 η H0 (η) + 5120 40960 η 

√ 

√ 

√  95 95 95

H0,−1 η + H0,1 (η) H0 (η) − H−1 η + H1 η H02 (η) + 3584 14336 28672 



√ 

√  847 95 95 1 7/2 H0,0,−1 η − H0,0,1 (η) η + 7/2 + − H0,−1 η − 1792 14336 497664 η 

√ 

√  847 847

H0,1 (η) H0 (η) + H−1 η + H1 η H02 (η) − 1990656 3981312



√  847 1 847 (B.4) H0,0,−1 η + H0,0,1 (η) η9/2 + 9/2 + 248832 1990656 η

+

PS,(3)

AQq

(N = 10) =

 1  2689775322848 401408 1  5105152 1605632 − + + L + − ) (L 2 1 ε 1260782263125 735075ε 3 ε 2 33078375 245025  100352  2 1276288 100352 2 L 2 + L1 + H0 (η) (L2 − L1 ) − (L2 + L1 ) − 81675 11026125 81675  200704 2 50176  3 19054928458130951 200704  − H0 (η) − ζ2 − − L2 + L31 245025 245025 406677926793600000 81675  25088   672443830712 50176  319072  2 2 L 2 + L1 − H0 (η) L22 − L21 + − − H 2 (η) + 3675375 27225 420260754375 81675 0  43556878529 50176  436544 ζ2 (L2 + L1 ) − H0 (η) (L2 − L1 ) + − 81675 13476375 331195392000  25088 3 100352 100352 50176 2 H1 (η) H0 (η) + H (η) − H0 (η)H0,1 (η) + H0,0,1 (η) + 245025 735075 0 245025 245025

 226878798767 347257523 2  638144 200704 1 ζ2 + ζ3 + + H0 (η) η + + 11026125 735075 526901760000 8028979200 η

  55227289 2047449637 1 368396509553 H0 (η) η − − H 2 (η) + + − 10036224000 η 1075757760000 1405071360 0

 41219216111 27528100609 1 1 2 2 H0 (η) η − 2 + × η + 2 − 245887488000 270978048000 η η

1819513853 118201777 2  3 1 1 3 H (η) η + 3 + H0 (η) η − 3 + 9634775040 0 η 36130406400 η

 1197239  1182029 1 1 2386873 − H02 (η) η4 + 4 − H0 (η) η4 − 4 + − 100362240 802897920 401448960 η η

−

J. Ablinger et al. / Nuclear Physics B 927 (2018) 339–367

365



 507  1 1 507 507 + H02 (η) η5 + 5 + H0 (η) η5 − 5 991232 7929856 1982464 η η 

√ 

√  980747 980747 980747

+ H0,−1 η + H0,1 (η) H0 (η) − H−1 η 1351680 5406720 10813440 

√  2

√  980747 980747 1 √ + H1 η H0 (η) − η+ √ H0,0,−1 η − H0,0,1 (η) 675840 5406720 η 

√ 

√  734267 734267 734267

+ H0,−1 η + H0,1 (η) H0 (η) − H−1 η 3317760 13271040 26542080 

√  2

√  734267 734267 1 3/2 + H1 η H0 (η) − H0,0,−1 η − H0,0,1 (η) η + 3/2 1658880 13271040 η 

√ 

√  70889 70889 70889

+ − H0,−1 η − H0,1 (η) H0 (η) + H−1 η 409600 1638400 3276800 

√  2

√  70889 70889 1 5/2 + H1 η H0 (η) + H0,0,−1 η + H0,0,1 (η) η + 5/2 204800 1638400 η 

√ 

√  4179 4179

4179 H0,−1 η + H0,1 (η) H0 (η) − H−1 η + 81920 327680 655360 

√  2

√  4179 4179 1 7/2 + H1 η H0 (η) − H0,0,−1 η − H0,0,1 (η) η + 7/2 40960 327680 η 

√ 

√  39641 39641 39641

+ − H0,−1 η − H0,1 (η) H0 (η) + H−1 η 6635520 26542080 53084160 

√  2

√  39641 39641 1 + H1 η H0 (η) + H0,0,−1 η + H0,0,1 (η) η9/2 + 9/2 3317760 26542080 η 

√ 

√  507 507

507 H0,−1 η + H0,1 (η) H0 (η) − H−1 η + 1982464 7929856 15859712 

√  2

√  507 507 1 + H1 η H0 (η) − H0,0,−1 η − H0,0,1 (η) η11/2 + 11/2 . 991232 7929856 η (B.5) √ The above expressions depend on η only, but not on η, and they are symmetric for η ↔ η−1 . The expansions of the O(ε 0 ) terms for N = 2, 4, 6 up to O(η3 ln3 (η)) for η < 1 agree with the results given in Ref. [2]. With growing values of N , these expressions exhibit a growing degree of the polynomials in η. We have found that the general N formula cannot be expressed as a sum – product solution by means of difference field theory. This means that the corresponding solution will be given by a higher transcendental function depending on N and η. We will study this behavior elsewhere. +

366

J. Ablinger et al. / Nuclear Physics B 927 (2018) 339–367

References [1] J. Blümlein, A. De Freitas, C. Schneider, K. Schönwald, DESY 17-187. [2] J. Ablinger, J. Blümlein, A. De Freitas, A. Hasselhuhn, C. Schneider, F. Wißbrock, Nucl. Phys. B 921 (2017) 585, arXiv:1705.07030 [hep-ph]. [3] J. Ablinger, A. Behring, J. Blümlein, A. De Freitas, A. von Manteuffel, C. Schneider, Nucl. Phys. B 890 (2014) 48, arXiv:1409.1135 [hep-ph]. [4] J. Ablinger, J. Blümlein, S. Klein, C. Schneider, F. Wißbrock, arXiv:1106.5937 [hep-ph]; J. Ablinger, J. Blümlein, A. Hasselhuhn, S. Klein, C. Schneider, F. Wißbrock, PoS RADCOR 2011 (2011) 031, arXiv:1202.2700 [hep-ph]. [5] R. Harlander, T. Seidensticker, M. Steinhauser, Phys. Lett. B 426 (1998) 125, arXiv:hep-ph/9712228. [6] T. Seidensticker, arXiv:hep-ph/9905298. [7] J. Ablinger, et al., in preparation. [8] J. Ablinger, J. Blümlein, A. De Freitas, A. Hasselhuhn, A. von Manteuffel, M. Round, C. Schneider, F. Wißbrock, Nucl. Phys. B 882 (2014) 263, arXiv:1402.0359 [hep-ph]. [9] J. Ablinger, J. Blümlein, A. De Freitas, A. Hasselhuhn, A. von Manteuffel, M. Round, C. Schneider, Nucl. Phys. B 885 (2014) 280, arXiv:1405.4259 [hep-ph]. [10] J. Ablinger, A. Behring, J. Blümlein, A. De Freitas, A. Hasselhuhn, A. von Manteuffel, M. Round, C. Schneider, F. Wißbrock, Nucl. Phys. B 886 (2014) 733, arXiv:1406.4654 [hep-ph]. [11] J. Blümlein, G. Falcioni, A. De Freitas, Nucl. Phys. B 910 (2016) 568, arXiv:1605.05541 [hep-ph]. [12] A. Behring, J. Blümlein, G. Falcioni, A. De Freitas, A. von Manteuffel, C. Schneider, Phys. Rev. D 94 (11) (2016) 114006, arXiv:1609.06255 [hep-ph]. [13] J. Ablinger, J. Blümlein, S. Klein, C. Schneider, F. Wißbrock, Nucl. Phys. B 844 (2011) 26, arXiv:1008.3347 [hep-ph]. [14] J. Blümlein, A. Hasselhuhn, S. Klein, C. Schneider, Nucl. Phys. B 866 (2013) 196, arXiv:1205.4184 [hep-ph]. [15] A. Behring, J. Blümlein, A. De Freitas, A. von Manteuffel, C. Schneider, Nucl. Phys. B 897 (2015) 612, arXiv:1504.08217 [hep-ph]. [16] A. Behring, J. Blümlein, A. De Freitas, A. Hasselhuhn, A. von Manteuffel, C. Schneider, Phys. Rev. D 92 (11) (2015) 114005, arXiv:1508.01449 [hep-ph]. [17] J. Ablinger, et al., DESY 15-112. [18] A. Behring, I. Bierenbaum, J. Blümlein, A. De Freitas, S. Klein, F. Wißbrock, Eur. Phys. J. C 74 (9) (2014) 3033, arXiv:1403.6356 [hep-ph]. [19] J. Ablinger, A. Behring, J. Blümlein, G. Falcioni, A. De Freitas, A. Hasselhuhn, A. von Manteuffel, M. Round, C. Schneider, F. Wißbrock, PoS LL 2016 (2016) 065, arXiv:1609.03397 [hep-ph]. [20] J. Ablinger, A. Behring, J. Blümlein, A. De Freitas, A. von Manteuffel, C. Schneider, Comput. Phys. Commun. 202 (2016) 33, arXiv:1509.08324 [hep-ph]. [21] J. Ablinger, A. Behring, J. Blümlein, A. De Freitas, A. von Manteuffel, C. Schneider, DESY 17-199. [22] I. Bierenbaum, J. Blümlein, S. Klein, Nucl. Phys. B 820 (2009) 417, arXiv:0904.3563 [hep-ph]. [23] M. Steinhauser, Comput. Phys. Commun. 134 (2001) 335, arXiv:hep-ph/0009029. [24] J. Ablinger, A. Behring, J. Blümlein, A. De Freitas, A. von Manteuffel, C. Schneider, Nucl. Phys. B 922 (2017) 1, arXiv:1705.01508 [hep-ph]. [25] J. Blümlein, A. De Freitas, W. van Neerven, Nucl. Phys. B 855 (2012) 508, arXiv:1107.4638 [hep-ph]. [26] M. Buza, Y. Matiounine, J. Smith, R. Migneron, W.L. van Neerven, Nucl. Phys. B 472 (1996) 611, arXiv:hepph/9601302. [27] M. Buza, Y. Matiounine, J. Smith, W.L. van Neerven, Eur. Phys. J. C 1 (1998) 301, arXiv:hep-ph/9612398. [28] I. Bierenbaum, J. Blümlein, S. Klein, Nucl. Phys. B 780 (2007) 40, arXiv:hep-ph/0703285. [29] I. Bierenbaum, J. Blümlein, S. Klein, C. Schneider, Nucl. Phys. B 803 (2008) 1, arXiv:0803.0273 [hep-ph]. [30] I. Bierenbaum, J. Blümlein, S. Klein, Phys. Lett. B 672 (2009) 401, arXiv:0901.0669 [hep-ph]. [31] E.W. Barnes, Proc. Lond. Math. Soc. (2) 6 (1908) 141. [32] E.W. Barnes, Q. J. Math. 41 (1910) 136. [33] H. Mellin, Math. Ann. 68 (3) (1910) 305. [34] E.T. Whittaker, G.N. Watson, A Course of Modern Analysis, Cambridge University Press, Cambridge, 1927, reprinted 1996. [35] E.C. Titchmarsh, Introduction to the Theory of Fourier Integrals, Calendron Press, Oxford, 1937, 2nd edition, 1948. [36] J.A.M. Vermaseren, arXiv:math-ph/0010025. [37] M. Czakon, Comput. Phys. Commun. 175 (2006) 559, arXiv:hep-ph/0511200.

J. Ablinger et al. / Nuclear Physics B 927 (2018) 339–367

367

[38] A.V. Smirnov, V.A. Smirnov, Eur. Phys. J. C 62 (2009) 445, arXiv:0901.0386 [hep-ph]. [39] J.A.M. Vermaseren, Int. J. Mod. Phys. A 14 (1999) 2037, arXiv:hep-ph/9806280; J. Blümlein, S. Kurth, Phys. Rev. D 60 (1999) 014018, arXiv:hep-ph/9810241. [40] C. Schneider, Sémin. Lothar. Comb. 56 (2007) B56b. [41] C. Schneider, Computer algebra in quantum field theory: Integration, summation and special functions, in: C. Schneider, J. Blümlein (Eds.), Texts and Monographs in Symbolic Computation, Springer, Wien, 2013, p. 325, arXiv:1304.4134 [cs.SC]. [42] J. Ablinger, PoS LL (2014) 019; Computer Algebra Algorithms for Special Functions in Particle Physics, Ph.D. Thesis, J. Kepler University, Linz, 2012, arXiv:1305.0687 [math-ph]; A Computer Algebra Toolbox for Harmonic Sums Related to Particle Physics, Diploma Thesis, J. Kepler University, Linz, 2009, arXiv:1011.1176 [math-ph], 2009. [43] J. Ablinger, J. Blümlein, C. Schneider, J. Math. Phys. 52 (2011) 102301, arXiv:1105.6063 [math-ph]. [44] J. Ablinger, J. Blümlein, C. Schneider, J. Math. Phys. 54 (2013) 082301, arXiv:1302.0378 [math-ph]. [45] J. Ablinger, J. Blümlein, S. Klein, C. Schneider, Nucl. Phys. Proc. Suppl. 205–206 (2010) 110, arXiv:1006.4797 [math-ph]; J. Blümlein, A. Hasselhuhn, C. Schneider, PoS RADCOR (2011) 032, arXiv:1202.4303 [math-ph]; C. Schneider, J. Phys. Conf. Ser. 523 (2014) 012037, arXiv:1310.0160 [cs.SC]. [46] E. Remiddi, J.A.M. Vermaseren, Int. J. Mod. Phys. A 15 (2000) 725, arXiv:hep-ph/9905237. [47] J. Ablinger, J. Blümlein, C.G. Raab, C. Schneider, J. Math. Phys. 55 (2014) 112301, arXiv:1407.1822 [hep-th]. [48] J. Blümlein, Comput. Phys. Commun. 159 (2004) 19, arXiv:hep-ph/0311046. [49] S. Alekhin, J. Blümlein, K. Daum, K. Lipka, S. Moch, Phys. Lett. B 720 (2013) 172, arXiv:1212.2355 [hep-ph]. [50] K.A. Olive, et al., Particle Data Group, Chin. Phys. C 38 (2014) 090001. [51] J. Blümlein, C. Schneider, Phys. Lett. B 771 (2017) 31, arXiv:1701.04614 [hep-ph]. [52] L. Lewin, Dilogarithms and Associated Functions, Macdonald, London, 1958. [53] L. Lewin, Polylogarithms and Associated Functions, North-Holland, New York, 1981.