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

By employing the method of differential equations, we compute the various types of two-loop master integrals involved in CP-even heavy quarkonium exclusive production and decays. All the integrals presented in this paper can be casted into canonical forms and expressed in terms of Goncharov polylogarithms and Harmonic polylogarithms. These master integrals are frequently used in the calculation of NNLO corrections of the heavy quarkonium production processes, as $\gamma^*\gamma\rightarrow Q\bar{Q}$, $e^+e^-\rightarrow \gamma+ Q\bar{Q}$, and $H/Z^0\rightarrow \gamma+ Q\bar{Q}$, and decay processes. They are also applicable in the calculation of NNLO corrections to CP-even quarkonium inclusive production and decay processes.


Introduction
The production and decay mechanism of heavy quarkonium has been a longstanding topic since the discovery of heavy quarkonium [1,2]. Large discrepancies between the experimental data and the leading order calculations stimulate great amount of theoretical researches [3,4]. The complexity of the calculation for higher order QCD corrections of these processes made the investigation of this topic very challenging. Fortunately, the advance of Nonrelativistic Quantum Chromodynamics (NRQCD) factorization formalism enables people to study this mechanism more and more reliably [5]. The progress in NRQCD calculations has deepen our understanding of strong interactions. It was found that the discrepancies can be remedied by introducing next-to-leading order (NLO) Quantum Chromodynamics (QCD) corrections [6,7]. The calculations according to NRQCD recognize significant NLO corrections comparing to the LO results. However, people noticed that even after introducing the relatively large NLO corrections, the renormalization scale dependence and uncertainties can still be quite significant. Therefore, NNLO corrections will most likely be necessary if one hopes to further classify the issue.
The calculation of multi-loop corrections for heavy quarkonium processes is always considered to be difficult primarily due to the massive two-loop integrals it involves. In 1997, the first complete analytical calculation of the NNLO corrections to the processes of J/ψ/Υ → e + e − and e + e − → J/ψ/Υ was achieved [8,9] (the NNNLO corrections to the same processes have not been obtained until 2014 [10]). Since then, NNLO corrections to the heavy quarknoium processes have been studied extensively [11][12][13]. Recently, the NNLO corrections to the process γγ * → η c were calculated [14]. We notice that most of the master integrals computed in that work are in non-physical region. Moreover, the calculations of the integrals were performed purely numerically. The accuracy of the numerical approach is limited by the applicability of the numerical packages available. Higher accuracy is expected, especially for the processes in physical region. For example, the master integrals involved in the computation of e + e − → γ * → γ + η c /η b are all in physical region. Due to the limitation of today's numerical packages, a full numerical calculation of these integrals will most certainly be unreliable. Analytical calculations of the master integrals are highly desired.
Threshold expansion [15] is widely used in the calculations of loop-corrections for heavy quarkonium production and decay. This method is not only very convenient, but also provides an implicit definition of NRQCD in dimensional regularization [8]. Within threshold expansion, only contributions coming from the hard loop momenta are needed, the contributions from soft, potential, and ultrasoft loop momenta can all be factored out. Thanks to this method, the calculations of heavy quarkonium processes become much more efficient than the calculations of corresponding heavy quark processes.
The Feynman diagrams of the NNLO corrections to CP-even heavy quarkonium exclusive production γ * γ → QQ and e + e − → γ * → γ + QQ share the same topology (both massless and massive "light-by-light" diagrams are included). By performing reducitons of the amplitudes, the scalar integrals involved are reduced to 133 master integrals. We notice that 86 out of the 133 master integrals can be expressed in terms of Harmonic polylogarithms and Goncharov polylogarithms. In this work, we will focus on the analytical solutions of these integrals.
The paper is organized as follows. In section 2, we introduce the kinematics and notations for CP-even quarkonium exclusive production processes. We also present the generic form of the differential equations with respect to the kinematics variables. In section 3, the Goncharov polylogarithms as well as Harmonic polylogarithms are introduced. In section 4, the canonical basis is explicitly presented, followed by the discussion of their solutions. In sections 5, the determination of the boundary conditions, as well as the analytical continuation are explained. Discussions and conclusions are made in section 6. All the analytical results up to weight four from our computation are collected in an ancillary file that we submit to the arXiv, the results up to weight three are listed in appendix.A. Figure 1. Sample of Feynman diagrams contributing to the two-loop QCD corrections of γ * + γ → QQ and QQ → γ * /γ + γ.

Notations and Kinematics
We consider the kinematics for the production of heavy quarkonium through γ * γ collision, and the process associated with a photon via a virtual photon.
where k 2 1 = 2ss, k 2 2 = 0 and k 2 q = k 2 q = m 2 q . Sample of Feynman diagrams contributing to the NNLO corrections for CP-even heavy quarkonium exclusive production and decay are showed in Fig.(1). For processes (2.1) in Euclidean region with ss < 0, we have the following relation For processes (2.2) in Minkowski region with 2ss > 4m 2 q , we have Within the threshold expansion approach, we take the momentum of quark and anti-quark to be equal k q = kq. In order to express the integrals more compactly, we introduce three dimensionless variables x, y and z as The QCD corrections to the processes (2.1) and (2.2) are calculated using Feynman diagram approach. After some manipulations, the amplitudes can be expressed in terms of a set of scalar integrals. The calculation of these scalar integrals always turns out to be the most difficult part in the whole work. We first use packages FIRE [36][37][38] to reduce the group of scalar integrals into a minimum set of independent master integrals. FIRE is also adopted in the following derivations of differential equations.
The first step of deriving differential equations is taking derivatives of lorentz invariant, and writing them down as linear combinations of master integrals. The derivatives of the external momenta can be transformed to the derivatives of ss and m 2 where i = 1, 2. Using these linear equations, we can write the derivatives ∂ ∂ss in terms of the derivatives The transform of the derivatives can readily be obtained according to (2.5). With the variables chosen above, we express the integrals in terms of the so-called Goncharov polylogarithms and Harmonic polylogarithms, which will be discussed in the next section.

Goncharov polylogarithms and Harmonic polylogarithms
The Goncharov polylogarithms (GPLs) [39] are defined as follow They can be viewed as a special case belonging to a more general type of integrals called Chen-iterated integrals [40]. If all the index a i belong to the set {0, ±1}, the Goncharov polylogarithms turns into the well-known Harmonic polylogarithms (HPLs) [41]

3)
H a 1 ,a 2 ,...,an (x) = (−1) k G a 1 ,a 2 ,...,an (x), (3.4) where k equals to the times of element (+1) taken in (a 1 , a 2 , . . . , a n ) . The GPLs fulfil the following shuffle rules Here, aXb is composed of the shuffle products of list a and b. It is defined as the set of the lists containing all the elements of a and b, with the ordering of the elements of a and b preserved. The GPLs and HPLs can be numerically evaluated within the GINAC implementation [42,43]. A Mathematica package HPL [44,45] is available to reduce and evaluate the HPLs. Both the GPLs and HPLs can be transformed to the function of ln, Li n and Li 22 up to weight four, with the methods and packages described in [46]. Figure 2. Set of master integrals that can be cast into canonical form and expressed in terms of GPLs. The thin lines denote massless propagators and on-shell massless external particles; the thick lines denote massive quark propagators and on-shell external quarks. The dash lines denote off-shell external particles with squared momentum equal to 2ss. A dot on a propagator indicates that the power of the propagator is raised to 2. Two dots means that the propagator is raised to power 3.

The canonical basis
It is widely known that the two-loop sunrise integrals with non-zero masses and off-shell external legs cannot be expressed only in terms of multiple polylogarithms. They can be evaluated and expressed, however, in terms of elliptic multiple polylogarithms [47,48]. All the integrals of the processes that we concern can be reduced to a minimum set of 133 master integrals. 41 of these master integrals cannot be cast into canonical form and then expressed in terms of GPLs. Among the 41 integrals, the solutions of 2 full massive sunrise integrals with off-shell external legs and 1 kite integral have been discussed [47][48][49][50][51][52]. 6 of the 133 master integrals were calculated in terms of Chen-iterated integrals (integrals 32) in [53]). Here we focus on the remaining 86 integrals. They all can be cast into canonical form and expressed in terms of HPLs and GPLs. The remaining 41 integrals will be discussed in detail elsewhere [54]. The corresponding of 86 linear differential equations can be expressed, via a suitable basis choice of the master integrals, into the canonical form [21] Here F is a vector of 86 canonical master integrals. The ǫ dependence is then completely factorized from the matrix dA. The vector F depends only on the dimensionless variables {x, y, z} , defined in Eq. (2.5). The matrix dA represents a total differential as follow where R k are constant matrices that contain only rational numbers, and the l k are linear functions of {x, y, z}. The so-called alphabet of the GPLs is preserved in order to obtain the result where x n ∈ {x, y, z} . The constant matrices R k with rational numbers are collected in an ancillary file called "Matrix.txt" that we submit to the arXiv.
The search of the canonical basis is based on an experimental approach with trial and error, with the aid of Mathematica code written by ourselves. We notice that packages for choosing a canonical basis appeared recently [55,56].
The vector of basis F is built up with 86 functions F i (ss, m q , ǫ), defined in terms of the master integrals M i drawn in Fig.(2), 11) F 10 = ǫ 2 ss M 10 , (4.13) 14) , (4.24) The integrals M 1 is defined as follow where the measure of the integration is defined as For master integrals without numerators, their definition can be read off from Fig.2, with the normalization defined above. For master integrals with numerators, we first define a series of propagators Then, the master integrals with numerators can be expressed as and . (4.94)
Since For M 12 whose results are constants, we found that its results could be obtained by manipulating F 19 . First the base F 19 equal to ǫ 2 ss (ss − 2m 2 q ) M 19 + 3ss/2 M 12 , considering the fact that M 19 does not have singularity at ss = 2m 2 q , and the normalization factor of M 19 in F 19 contains (ss − 2m 2 q ), we can take the limit at ss → 2m 2 q for F 19 and obtain the results of F 12 . With the assistance of PSLQ algorithm [61], the results of F 12 can be simplified and expressed as follow the results of weight four of F 12 is new to the best of our knowledge. The result of F 17 can be obtained by taking the limit ss → 2m 2 q for F 6 and expressed as Because all the master integrals appeared in F 40 are regular at ss = 8m 2 q 3 , the boundary of F 40 is calculated by taking the limit ss → 8m 2 q 3 for F 40 . Considering the fact that M 49 does not have singularity at ss = 0, the boundary condition of F 52 can be determined from the differential equation of F 49 . To further illustrate it, we consider the differential of F 49 with respect to variable x and find that where ellipses stand for less singular terms. All the integrals that appear in Eq.(5.3) have finite limits at x → 1(ss → 0). This consistency leads to a relation between different integrals lim x→1 We can then obtain the boundary condition of F 52 at x = 1 from the above equation.
A more general calculation of the kinematics of F 86 was performed before [53,62,63]. By taking the limit p 2 = 4m 2 for F 10101 in [63], and by adopting the PSLQ algorithm, we obtain the results of F 86 By now, all the boundary conditions are fixed. The next step is to determinate the analytic continuations of the master integrals. When we consider the processes (2.1,2.2), the variables of the master integrals lie either in Euclidean region or in Minkowski region. Their analytic continuations should be considered carefully. The proper analytic continuation can be achieved by the replacement ss → ss+i0 at fixed m 2 q . This transfer corresponds to x → x + i0, y → y + i0 and z → z + i0. For single scale integrals that depend only on m 2 q , the replacement m 2 q → m 2 q + i0 is sufficient for obtaining the correct results.
The calculations are performed with our self-written Mathematica code. All the analytical expressions of the master integrals require an independent examination. We check all the results against the results obtained from numerical programs Fiesta [64,65] and SecDec [66,67]. Good agreement has been achieved between the analytical and numerical approaches with kinematics in both Euclidean region and Minkowski region.

Discussions and Conclusions
In this work, we obtained the analytic results of 86 out of the 133 master integrals involved in the calculation of NNLO corrections to CP-even heavy quarkonium production processes such as γ * γ → QQ and e + e − → γ + QQ. By choosing a proper canonical basis, the differential equation group is cast into a canonical form. All of the 86 master integrals are then expressed in terms of Harmonic polylogarithms and Goncharov polylogarithms. The integrals obtained here may also be applied to the calculation of NNLO corrections of other process, such as the exclusive decay of higgs boson or Z 0 boson into CP-even quarkonium plus a photon and the inclusive hadron production or decay of η c /η b . The remaining integrals which cannot be expressed in Goncharov polylogarithms require a further investigation.