Two-Loop integrals for CP-even heavy quarkonium production and decays: Elliptic Sectors

By employing the differential equations, we compute analytically the elliptic sectors of two-loop master integrals appearing in the NNLO QCD corrections to CP-even heavy quarkonium exclusive production and decays, which turns out to be the last and toughest part in the relevant calculation. The integrals are found can be expressed as Goncharov polylogarithms and iterative integrals over elliptic functions. The master integrals may be applied to some other NNLO QCD calculations about heavy quarkonium exclusive production, like $\gamma^*\gamma\rightarrow Q\bar{Q}$, $e^+e^-\rightarrow \gamma+ Q\bar{Q}$,~and~$H/Z^0\rightarrow \gamma+ Q\bar{Q}$, heavy quarkonium exclusive decays, and also the CP-even heavy quarkonium inclusive production and decays.


Introduction
Precision physics in colliders requires more higher-order corrections in perturbation theory. Unravelling the mathematical structure of Feynman integrals in multiloop calculation is somehow critical to handle the complexity of higher order calculations and may help us to obtain a better control of the perturbative expansion. In recent years, the corresponding research achieved some breakthroughs and becomes now one of the hot topics in physics and mathematics.
The heavy quarkonium production and decay are one of the hot topics in particle physics ever since the first discovery in 1974, especially with the advent of Nonrelativistic Quantum Chromodynamics (NRQCD) factorization formalism [41]. Up to date there still γ * γ γ * γ Figure 1. Typical two-loop Feynman diagrams for CP-even heavy quarkonium production. exist some discrepancies between experimental data and theoretical expectations [42][43][44][45], which appeal for precision calculations. In one of our previous works [46] we gave out a set of 86 two-loop master integrals about heavy quarkonium production and decay, which can be cast into the canonical form and expressed in terms of multiple polylogarithms. However, for those Feynman integrals with functions beyond the realm of multiple polylogarithms the calculation is not done yet. In fact, to date, only a limited number of similar calculations have been performed in the literature.
In this work, we calculate analytically all remaining integrals with different mathematical structures from multiple polylogarithms in CP-even heavy quarkonium production and decays. The master integrals will be classified into two sectors, one with integrals containing sub-topologies related to the two-loop massive sunrise integrals and the other involving non-planar two-loop three-point integrals. Following the strategy suggested in Ref. [39] and with properly chosen basis, we cast the differential equations of those integrals in the first sector into a proper form that can be solved recursively. Of the second sector, the key point is to find the homogeneous solutions for the second-order differential equations of the two-loop non-planar three-point massive integrals, with that the full solutions can then be obtained by constant variation.
The paper is organized as follows. In section 2, the kinematics is discussed and the derivatives with respect to kinematic variables will be given. In section 3, the iterative integrals and complete elliptic integrals are introduced. In section 4, the elliptic type integrals will be separated into two sectors, and the calculation procedure for them will be elucidated respectively. For illustration, specific examples will be given. Section 5 is remained for conclusions and outlooks. The definition of master integrals is given in appendix A, and several simple but typical analytical results are presented in appendix B.

Notation and kinematics
The heavy quarkonium exclusive production in electron-positron collision has a relatively low background, and has played an important role in the study of quarkonium production mechanism. Here we calculate the CP-even quarkonium production in two correlated processes, that is in γ * γ collision and in electron-position annihilation associated with a photon, where k 2 1 = 2ss, k 2 2 = 0 and k 2 q = k 2 q = m 2 q . The typical Feynman diagrams are showed in Fig. 1. The process (2.1) is in Euclidean region with ss < 0, and the momenta satisfy the following relations Whereas, the process (2.2) is in Minkowski region with 2ss > 4m 2 q , and Note, in the threshold expansion approach, quark and anti-quark momenta are taken to be equal, i.e. k q = kq.
In order to express the results compactly, here we introduce three dimensionless variables x, y and z as follows: The NNLO QCD corrections to processes (2.1) and (2.2) are calculated in light of Feynman diagrams. As a routine, with some algebraic manipulations, the amplitudes can be reduced to a set of scalar integrals. We use the Mathematica package FIRE [47][48][49] to reduce the scalar integrals to a minimum set of independent master integrals. The calculation of these master integrals is the central issue, and normally turns out to be a nontrivial work. In our calculation, we apply the method of differential equations to calculate the master integrals.
The first step of deriving differential equations is taking derivatives of the Lorentz invariant kinematic variables, and expressing them as linear combinations of master integrals. The FIRE is also employed in the derivation of differential equations. The derivatives of the external momenta can be expressed as the derivatives of ss and m 2 q , like with i(j) = 1 or 2. And in reverse, the derivative ∂ ∂ss can be expressed as a linear combination of derivatives k i · ∂ ∂k j , i.e., The derivative transform can be readily obtained according to equation (2.5). With the variables chosen in above, analytical results of the integrals can then be formulated in a compact form, in terms of iterative integrals and elliptic integrals.

Iterated integrals and complete elliptic integrals
The Goncharov polylogarithms (GPLs) [50] are defined as which in fact are special cases of a more general type of integrals, named Chen-iterated integrals [51]. If all indices a i belong to set {0, ±1}, the Goncharov polylogarithms can then be transformed into the well-known Harmonic polylogarithms (HPLs) [52]

3)
H a 1 ,a 2 ,...,an (x) = (−1) k G a 1 ,a 2 ,...,an (x) , (3.4) where k equals to the number of times the element (+1) appearing in (a 1 , a 2 , . . . , a n ) . The GPLs satisfy the following shuffle rules: In above equation, aXb is composed of the shuffle products of a i (i = 1, 2 . . . m) and b i (i = 1, 2 . . . n), which is defined as the set of lists containing all elements of a i and b i , with the order of elements a i and b i preserved. The GPLs and HPLs can be numerically evaluated by implementing the GINAC [53,54], and the Mathematica package HPL [55,56] is applicable to the HPLs reduction and evaluation. Both GPLs and HPLs can be transformed into functions ln, Li n and Li 22 up to weight four in light of the method described in Ref. [57]. In our calculation, the complete elliptic integrals are necessary to express the integrals encountered. The first and second kinds of complete elliptic integrals are defined as and They satisfy the following derivative relations: The Legendre relation is useful in simplifying the complete elliptic integrals, i.e., (3.9) -4 -

Elliptic integral sectors
The symbols and canonical basis in the calculation of elliptic integrals keep the same as in the preceding work [46], where the linear differential equations can be expressed, via a suitable basis choice of master integrals, as canonical form [6] d with F being the vector of canonical master integrals F i (i = 1 . . . 86) [46]. Whereas, the two-loop massive Feynman integrals concerned in this work may involve elliptic functions, and hence the calculation of the integrals should be further explored. We separate them into two elliptic sectors: one with integrals containing sub-topologies related to the two-loop massive sunrise integrals, the other with two-loop non-planar three-point integrals. In the following we elucidate the calculation procedures of these integrals.

Sector I : integrals with massive sunrise integrals as subtopology
The 39 Feynman integrals E i (i = 1 . . . 39) belonging to this subsection are shown in Fig.  2, which contain sub-topologies related to the two-loop massive sunrise integrals. The expressions of master integrals without numerators can be readily read off from the figure, and those with numerators are given in appendix A. Note, the massive sunrise integrals are composed of the complete elliptic integrals and cannot be expressed as pure Goncharov polylogarithms. The two-loop massive sunrise integrals (E 1 , E 2 ) have been widely studied.
Here, the bases (A 1 , A 2 ), which contain (E 1 , E 2 ), are of the same as their first appearance in Ref. [39]: , (4.2) In the following we sketch the calculation of this sector. With a suitable choice of the basis in the high topologies (E 3 . . . E 39 ), the homogeneous part of the differential equations for integrals (E 3 . . . E 39 ) can be cast into the canonical form, whereas depending on the inhomogeneous terms of massive sunrise integrals (E 1 , E 2 ), or (A 1 , A 2 ). To be more specific, after a proper selection of bases Here, A ′ is a 37-dimensional basis vector containing integrals E i (i = 3 . . . 39) and F i (i = 1 . . . 86); F is a 86-dimensional basis vector that was given in Ref. [46]; W and Y are 37×37 Notice that in equation (4.4) the inhomogeneous term that contain A 2 is free of ǫ, and the differential equation for A 1 given in Ref. [39] can be reexpressed as Since in above equation the inhomogeneous term containing A 2 is also ǫ free, we therefore are legitimate to perform a basis shift as With the basis shift, Q 3 A 2 will be removed from the differential equation (4.4). Here b i (ss) are algebraic functions to be determined. Moreover, the basis shift may also simplify the inhomogeneous term containing A 1 , considerably. For illustration, we take the differential equations for (E 4 , E 5 , E 6 ) as an example, which have the same topology. By properly choosing the basis, the differential equations for Here, e(ss, ǫ), a 3-dimensional basis vector containing integrals (E 4 , E 5 , E 6 ) and F 12 , may be expressed as and A 1 and A 2 are scalar functions defined as (4.3). To remove the A 2 dependence from the inhomogeneous part of the differential equations, we perform the basis shift where b i (ss) are algebraic functions to be determined. By virtue of the differential equation for A 1 , one can figure out the shift functions b i (ss) in (4.10), which may be formulated in a 3-dimensional vector form (4.11) The differential equation for e 1 (ss, ǫ) is in canonical form, and hence no need to make the shift. After the basis shift, Λ 0 (ss)A 2 and Ω 0 (ss)A 1 terms in differential equation for e 3 (ss, ǫ) vanish, and the differential equation for e 3 (ss, ǫ) turns to be canonical. Of the differential equation for e 2 (ss, ǫ), though Λ 0 (ss)A 2 term does not exist, Ω 0 (ss)A 1 term remains. Note, with the basis shift the inhomogeneous part of the differential equations for e 2 (ss, ǫ) will be greatly simplified, and the differential equations turn to be solvable recursively.
The method described above is also applicable to high sectors with more propagators. Except for integrals (E 1 , E 2 , E 5 , E 9 ), differential equations for the remaining 35 integrals can be transformed into the canonical form (4.1), with the method employed in this work. The basis vector A is built up with 39 functions A i (ss, m q , ǫ), the linear combinations of master integrals E i and F i with the latter given in Ref. [46]. Explicitly, the 39 bases that contain planar and non-planar two-loop integrals can be formulated as With the basis chosen above, the differential equations for (A 3 . . . A 39 ) then turn to the canonical form, except for A 5 and A 9 . The differential equations for A 5 and A 9 with respect to x write as: Notice that the above two equations are not in canonical form, and they both have the ǫ free A 1 terms, by a factor of 2 difference. Those terms without A 1 can be expressed in d-log form. By using the method described in above, different from casting all terms into canonical form via (non-algebraic) basis change in Ref. [40], the obtained differential equations are greatly simplified and are suitable for solving recursively. Taking the known result on A 1 [39] as an input, the differential equations for (A 3 . . . A 39 ) can be integrated straightforwardly order by order in ǫ. The corresponding lengthy expressions is given as an auxiliary file in arXiv version of this paper. After determining the bases, to fix the boundary conditions is necessary for solving the differential equations. Here, we apply the regularity conditions as in Ref. [4] to assist the determination of boundary conditions. Noticing that the integrals (E 3 , E 4 , E 5 , E 7 , E 8 , E 9 , E 11 , E 13 , E 15 , E 16 , E 17 , E 19 . . . E 23 , E 25 , E 26 , E 28 . . . E 39 ) are regular at ss = 2m 2 q and multiplying the normalization factor (ss − 2m 2 q ) to A i , one may find that the corresponding bases A i turn to be zero at ss = 2m 2 q . The boundary condition for A 6 at ss = 2m 2 q can be fixed in a similar way, that is (4.14) Here, the integral F 12 is known, and the boundary condition for A 1 may be determined from its definition in (4.12), i.e. A 1 | ss=2m 2 q = 3 4 F 12 . The integral E 10 is regular at ss = m 2 q with the normalization factor (ss − m 2 q ). Multiplied by this normalization factor, we then find A 10 = 0 at ss = m 2 q . Since the integrals (E 12 , E 14 , E 18 , E 27 ) are also regular at ss = 2m 2 q , the boundaries of corresponding bases A i can be determined by differential equations. For instance, the differential equation for A 12 reads where ellipses stand for less singular terms at y = 0, i.e. ss = 2m 2 q . Since all integrals in (4.15) have finite limits at y → 0, the following relation between different integrals exists: Because (F 24 , A 7 , A 11 ) are zero at y = 0 (ss = 2m 2 q ), we then have Similarly, from those boundaries for integrals E 14 , E 18 , E 24 and E 27 , one can fix all boundary conditions for bases (A 1 . . . A 39 ), of which the none-zero ones up to weight-4 write as: ) − 8 ln 2 (2)(π 2 + ln 2 (2))) + O(ǫ 5 ),

Sector II : non-planar two-loop three-point integrals
In this subsection we consider the non-planar two-loop three-points integrals that appear in the massive light-by-light Feynman diagrams. There are eight master integrals, as shown in Fig. 3, with the corresponding bases B i as Note, here the integrals (C 1 . . . C 6 ) were first calculated in Ref. [58], and the left two non-planar two-loop integrals (B 7 , B 8 ) cannot be cast into the canonical form via algebraic change of basis. A similar topology of Feynman diagram as that of (C 7 , C 8 ), but with different kinematics and outgoing momentum squared, was handled in Ref. [59].
In order to get expressions for B 7 and B 8 we first derive two coupled first-order differential equations with the evolution of variable ss, and then transform them to a second-order differential equation for B 7 . That is: 20) with N (ǫ, ss, m 2 q ) denoting the non-homogeneous term. Here, the tough issue is how to determinate the homogeneous solution. To this aim, we make a variable transformation of ss to v = −i(ss−2m 2 q ) 4m 2 q , then the homogenous part of the differential equation turns to The solutions of equation (4.21) can be readily obtained. The two homogeneous solutions (y 1 (v), y 2 (v)) read with K(x) being the first kind complete elliptic integral. Note that the recently development on maximal-cut [60][61][62] is also applicable to the determination of the homogeneous solution.
The Wronskian of the homogeneous solution reads .
With the homogeneous solutions and Wronskian, a particular solution can be obtained by means of the constant variation. The general solution is then where i refers to the order of ǫ in B 7 .
Since the integral C 7 has no singularity at ss = 2m 2 q , and the normalization for C 7 in B 7 is (ss − 2m 2 q ) 2 , we know Hence, the constants c 1 and c 2 can be fixed to Once B 7 is obtained, we can then determine the B 8 from the first order differential equation with respect to B 7 straightforwardly. Before calculating the differential equations for integrals in this sector, still the corresponding boundary conditions should be fixed. Since the integrals (B 1 , B 2 , B 3 , B 4 , B 6 ) are regular at ss = 2m 2 q , by multiplying their normalization factor (ss − 2m 2 q ) to B i , the corresponding bases B i then turn out to be zero at ss = 2m 2 q . Considering that the master integrals in basis B 5 are regular as ss = 0 and have a common normalization factor ss, we readily know B 5 = 0 when ss = 0. With these discussions, all necessary boundary conditions to fix the solutions of differential equations are ready.

Analytic continuation and discussions
With the analytical results obtained in above, the next necessary step is to determinate the analytic continuation of the master integrals, which is similar to the procedure in our previous work [46]. The correct analytic continuation can be achieved by the replacement of ss → ss + i0 at fixed m 2 q , which corresponds to x → x + i0, y → y + i0 and z → z + i0. The canonical bases in (4.12) contain 4 independent square roots ( √ ss, ss − 2m 2 q , ss + 2m 2 q , ss + 6m 2 q ) , (4.27) which cannot be simultaneously rationalized via one variable change. This means it is not possible to integrate the differential equations directly in terms of Gongcharov polylogarithms. It is worth mentioning that Refs. [58,63] proposed some novel ways to express the results of canonical bases for non-elliptic sectors in terms of multiple polylogarithms, without considering the existence of rational parametrization of the alphabet. However the results tend to be rather lengthy when expressed in multiple polylogarithms. In order to calculate the integrals numerically in a faster and convenient way, we construct a one-fold integral representation for the integrals that can be cast into the canonical form by means of what proposed in Ref. [58]. For integrals in elliptic sectors we need the two-fold integral representation to express the results up to weight four. The one fold and two fold integral representations we adopted are suitable for fast and precise numerical evaluation with Mathematica program on a single core computer. The analytic calculation in this work is performed by our own developed Mathematica code, and in order to guarantee the correctness of our results, we ask all analytical expressions for master integrals experiencing at least one independent examination. We check all results in contrast to those obtained via numerical programs Fiesta [64,65] and SecDec [66,67]. We have achieved an excellent agreement in analytical and numerical approaches with kinematics in both Euclidean and Minkowski regions.

Conclusions and outlooks
The integrals involving elliptic functions in the NNLO QCD corrections to heavy quarkonium exclusive production and decays are calculated, which turns out to be a tough issue. Those integrals are classified into two sectors, one with integrals containing sub-topologies related to the two-loop massive sunrise integrals and another with two massive two-loop non-planar three-points integrals. We find the simple example studied in Ref. [39] is in fact applicable to more general cases, that is, the expressions for two master integrals composed of two-loop massive sunrise integrals are still suitable for our case. In order to compute the first sector Feynman integrals under consideration we exploit the result for the two-loop massive sunrise integrals in Ref. [39]. We find a suitable linear combination of Feynman integrals such that only one of the master integrals about the solutions of two-loop massive sunrise integrals is required. By properly choosing canonical basis, we transform the differential equations into a simple and compact form that can be solved recursively. For another elliptic sector, the key point is to solve the homogeneous equation, with that inhomogeneous solutions can be obtained by means of constant variation.
Together with those 86 integrals calculated in our previous work [46], all master integrals appearing in the calculation of NNLO QCD correction to CP-even heavy quarkonium exclusive production and decays, such as γ * γ → QQ and e + e − → γ + QQ [68], are ready. The master integrals take the form of mutilple polylogarithms, iterative integrals over complete elliptic integrals and multiple polylogarithms. It is noteworthy that the integrals calculated in this work may also appear in the calculation of NNLO corrections in other processes, such as the exclusive decay of Higgs or Z 0 boson to quarkonium plus a photon and the inclusive hadronic production or decay of η c /η b , which are also phenomenologically meaningful. Moreover, we tend to believe that the calculation procedure and results in this work might be helpful to the mater integrals calculation of processes beyond the scope of heavy quarkonium physics, for instance the NNLO corrections to top quark pairs hadronic production, and NNLO corrections to heavy quark pair production plus a jet in electron-positron collision.
Note, only simple results are given in the appendix, however the full but lengthy results will be provided upon request.

A The definition for integrals
The integral A 1 is defined as where the measure of the integration is For master integrals without numerators, their definition can be read off from Fig. 2 and Fig. 3, with the normalization defined in above. For master integrals with numerators, we can define a series of propagators as P 9 = −(q 2 + k 2 − k q ) 2 , P 10 = −(q 1 + q 2 + k q ) 2 + m 2 q , P 11 = −(q 1 + k q ) 2 , P 12 = −(q 2 + k q ) 2 , P 13 = −(q 1 + k 1 + k q ) 2 , P 14 = −(q 1 − k q ) 2 , Then, the master integrals with numerators can be expressed as M 24 = D D q 1 D D q 2 P 11 P 1 P 2 P 2 4 P 7 P 10 , M 34 = D D q 1 D D q 2 P 7 P 1 P 2 P 3 P 4 P 9 P 10 , M 36 = D D q 1 D D q 2 P 12 P 1 P 2 P 3 P 6 P 8 P 10 , M 38 = D D q 1 D D q 2 P 12 P 1 P 2 P 3 P 4 P 5 P 10 , C 2 = D D q 1 D D q 2 P 11 P 2 P 9 P 10 P 13 , C 5 = D D q 1 D D q 2 P 12 P 2 P 10 P 11 P 14 P 15 . (A.4)