Master integrals for ${\cal O}(\alpha \alpha_s)$ corrections to $H \to ZZ^*$

We present analytic results for all the Feynman integrals relevant for ${\mathcal O}(\alpha \alpha_s)$ virtual corrections to $H \rightarrow ZZ^*$ decay. We use the method of differential equations to solve the master integrals while keeping the full dependence on the masses of all the particles including internal propagators. Due to the presence of four mass scales we encounter multiple square roots. We argue that all the occurring square roots can not be rationalized at the same time as a simultaneous rationalization brings us to integrals over $CY_3$ manifolds. Hence we rationalize only three square roots simultaneously and construct suitable ans\"atze to obtain dlog-forms containing the square root, after obtaining an epsilon-factorised form for the differential equations. We present the alphabet and the analytic form of all the boundary constants that appear in the solutions of the differential equations. The results for master integrals are expressed in terms of Chen's iterated integrals with dlog one-forms.

semi-leptonic (H → 2 2q) and hadronic (H → 4q) final states appear at the one-loop level via correction to the Zqq vertex. In the H → 4 decay mode with charged leptons in the final state, the QCD virtual corrections, however, involve HZZ vertex corrections at the two-loop level with massive quarks in the loop. In perturbation theory, it is an O(αα s ) effect. The relevant two-loop triangle Feynman diagrams are similar to those which appear in the QCD corrections to H → γγ, γZ decays. Due to the presence of many mass scales, the analytic calculation of two-loop three-point Feynman diagrams relevant to QCD corrections in H → 4 decay is quite non-trivial. Efforts are going on to provide these corrections by evaluating the two-loop Feynman integrals numerically [11].
State-of-the-art for analytic studies of two-loop three-point functions including massive propagators are as follows. For two-loop master integrals which appear in the virtual QCD corrections of H → γγ, the analytic results in terms of harmonic polylogarithms have been known for quite some time [12]. Also, the analytic results of master integrals for massless 2-loop 3-point functions relevant for QCD corrections to decay H → V * V * (V = Z or W ) with three off-shell legs have been obtained in terms of Goncharov polylogarithms [13]. The results for the master integrals relevant to NLO QCD corrections to H → Zγ with a massive top-quark loop are also available [14], and the results are expressed in terms of Goncharov polylogarithms. In [15], two-loop master integrals for QCD correction to neutral massive boson coupling to a pair of W bosons were computed. The master integrals in these cases depend on two or three mass scales at maximum. The full analytic results for the master integrals appearing in O(αα s ) corrections to HW + W − vertex are also known in terms of multiple polylogarithms [16]. A few examples of one-loop triangle diagrams with arbitrary mass dependence are given in [17][18][19][20]. A canonical set of master integrals applicable to O(αα s ) corrections in H → ZZ * decay with four mass scale dependence was studied in [21], however the analytic results for these integrals are not public.
In this article, we provide for the first time the full analytic result for all the twoloop canonical master integrals that contribute to O(αα s ) corrections to H → ZZ * decay, keeping the full dependence on the mass of the top quark (m t ) in the loop. We also present the analytic results of all the boundary constants along with a compact set of letters called the alphabet. These are important for phenomenological applications as a minimal set of letters is needed for a fast numerical evaluation, and analytic boundary constants help in evaluating the master integrals with higher precision. The results also cover a more general decay H → Z * Z * as the number of mass scales remains the same. In addition, we identify the hypersurface obtained during a simultaneous rationalization of all the square roots, occurring because of the multiple mass scales, to be a CY 3 manifold. This establishes a connection between the relevant Feynman integrals for this process and integrals that occur in more symmetric quantum field theories. The choice of a set of canonical master integrals presented here is also useful for similar 2-loop massive triangles with three offshell legs and a massive internal loop. For instance, this set of master integrals, apart from contributing to H → ZZ * , is also useful for mixed electroweak-QCD corrections to Higgs production via the process e + e − → HZ and e + e − → Hµ + µ − , which will be relevant at future e + e − colliders [21][22][23]. The results can also be used for precision studies in processes like µ + µ − → HZ, H + − at future muon collider [24].
The standard technique of evaluating multi-loop Feynman diagrams consists of first obtaining the set of master integrals using integration-by-part (IBP) identities [25,26], then setting up a system of differential equations for the master integrals, and finally solving the differential equations using appropriate boundary conditions [27][28][29][30][31][32][33][34]. In general, after bringing the system of differential equations into an epsilon-factorized form, their solutions can be written in terms of iterated integrals. The iterated integrals with rational integrands can be expressed in terms of Goncharov polylogarithms [35][36][37]. However, the presence of multiple square roots in the system of differential equations poses technical challenges in the analytic computations of these master integrals [38,39].
The study on the effects of the presence of multiple square roots in analytic calculations of Feynman integrals has gained a lot of momentum recently [40][41][42][43][44][45][46][47][48]. On the formal side, a lot of efforts have also been going on to understand the geometrical aspects of Feynman integrals. Many recent analytic computations of Feynman integrals show that the geometry is manifested by the square roots that appear in the differential equations of the master integrals or in their evaluation via direct integration [49,50]. Characterizing Feynman integrals by their geometries helps to understand the connection among various theories like QCD, electroweak and super Yang-Mills (SYM) theories. In our case, we study all the occurring square roots and comment on their rationalizability while establishing a connection with studies on geometrical aspects of Feynman integrals. We then construct ansätze to obtain a dlog-form of the system of differential equations with a minimal set of independent one-forms involving a non-rationalizable square root. We also present the results of all the master integrals in terms of Chen's iterated integrals [51] and perform numerical checks on our results by matching them against those obtained using publicly available numerical tools.
This article is organized as follows. In section 2, we discuss the notations of the Feynman integrals and the kinematic invariants. In section 3, we present all the pre-canonical as well as canonical master integrals. In section 3.2, we present the details of the square root(s) and discuss their rationalization. In section 3.3 we discuss all the one-forms present in the differential equations and present the alphabet. In section 4, we present the analytic results of all the boundary constants for the canonical master integrals and perform checks.

Two-loop Feynman diagrams and auxiliary topology
The representative Feynman diagrams contributing to H(q) → Z(p 1 )+Z * (p 2 ) at O(αα s ) are shown in figure 1. There are two independent external momenta p 1 and p 2 corresponding to Z and Z * respectively, with p 2 1 = m 2 Z and p 2 2 = m 2 Z * . The external momenta corresponding to Z * is off-shell implying p 2 1 = p 2 2 , in general. The Higgs boson is on-shell with q 2 = (p 1 + p 2 ) 2 = m 2 H . There are seven linearly independent scalar products involving loop momenta k 1 and k 2 , and external momenta p 1 and p 2 (k 2 1 , k 2 2 , k 1 .k 2 , k 1 .p 1 , k 1 .p 2 , k 2 .p 1 and k 2 .p 2 ) which appear in our Feynman integrals. The relevant two-loop Feynman diagrams, however, have a maximum of six propagators. In order to express Feynman integrals in terms of master integrals, we introduce an auxiliary topology with seven propagators as shown in figure 2. The corresponding integral family is given by Here γ E denotes the Euler-Mascheroni constant, D is the space-time dimension, µ is an arbitrary scale introduced to render the Feynman integral dimensionless, the quantity ν is given by The propagators P j are given by We obtain the set of master integrals using the IBP program Fire [52][53][54] combined with LiteRed [55,56], which gives us 41 master integrals. The set of 41 master integrals forms a basis which we denote by I. For convenience, the master integrals can be classified into sectors where a sector is defined by a common set of non-zero propagators. Master integrals belonging to each sector may have different powers of the common propagators. In our case, some sectors have more than one master integral. The list of all 41 master integrals, classified sector-wise and following the notation in equation 2.1, is shown in the column 2 of table 1. We may set µ 2 = m 2 t , after which our master integrals depend kinematically on three dimensionless quantities defined by

An epsilon-form for the master integrals
We set up a linear system of differential equations of the form d I = A I for all the master integrals of basis I by taking their derivatives with respect to the kinematic invariants (u, v and w) and using IBP identities. The 41×41 matrix A depends on kinematic invariants and space-time dimensions (D). The system of differential equations for the master integrals is solved iteratively at each order in the dimensional regularization parameter after switching from basis I to a canonical basis J such that the system of differential equations has an epsilon-form [33] given by Here matrixÃ is independent of the dimensional regulator . We aim to find a transformation matrix that brings the pre-canonical basis I to a canonical basis J leading to a differential system given in equation 3.1.

Finding the canonical basis
Arranging the master integrals sector-wise gives a block-diagonal form for their differential equations. We also use a "bottom-up" approach where we bring the differential equations of a lower sector to a canonical form using a coordinate system suitable to that sector before moving on to a higher sector. For some of the master integrals, it is convenient to work in D = 2 − 2 dimensions. We use the dimension shifting operator D − to shift the integrals from 4 − 2 to 2 − 2 . We express such integrals as a linear combination of master integrals in D = 4 − 2 dimensions using the dimensional shift relations [55]. To find a transformation matrix that brings our pre-canonical basis to a canonical form, we use the information of the maximal cuts of each of the sector [57][58][59][60][61]. Following the same method, as described in [44,47,62], we first construct an ansatz for a canonical basis using the information from the maximal cut of the sector in consideration. This fixes the diagonal part of the differential equations for that sector. We then include the contribution from the sub-sectors number of master integrals master integrals kinematic propagators basis I basis J dependence In the above, λ(u, v, w) is the Källen function defined by Note that in terms of the coordinates of equation 2.4, the system of differential equations contains multiple square roots. In particular, the following square roots appear A non-trivial task is to find a coordinate system which rationalizes all these square roots simultaneously. We introduce the following set of transformations that takes us from (u, v, w) to (x, y, z) and rationalizes the first three square roots in equation 3.3 These transformations are ubiquitous in literature [63]. After using these set of transformation we are left with only one square root λ (u, v, w), which in the new coordinate system becomes r = (y 2 z 2 + x 4 y 2 z 2 − 2xyz(z + y(1 + z(−2 + y + z))) − 2x 3 yz(z + y(1 + z(−2 + y + z))) + x 2 (−2y(−1 + z) 2 z − 2y 3 (−1 + z) 2 z + z 2 + y 4 z 2 + y 2 (1 + z(4 + z(−6 + z(4 + z)))))). (3.5) This square root r cannot be rationalized by doing further coordinate transformations. We discuss the appearance of square roots and their rationalizability in more detail in the next section.
Sector-wise classification of the master integrals of canonical basis and their kinematic dependence is given in columns 3 and 4, respectively, of table 1. The master integrals from higher sectors have more propagators and are usually more complicated to solve than the master integrals from the lower sectors. We also look at the topology of each sector, where multiple sectors may correspond to the same topology. The form of canonical basis can be recycled within sectors belonging to the same topology since these sectors are usually obtained by a permutation of the external legs. In our case, we have 9 master topologies which are shown in figure 3. Apart from being a guiding principle in the construction of a canonical basis, the classification of integrals in terms of master topologies also helps in fixing the boundary constants while solving the differential equations. The boundary constants are discussed in section 4.

The square root
Due to the presence of many mass scales in our system, the differential equations contain 4 square roots in the coordinate system given in equation 2.4. However, our choice of variables given in equation 3.4 rationalizes 3 out of 4 square roots present in the system. The fourth square root is r = P (x, y, z) given in equation 3.5. This square root is present in the analytic expressions of the master integrals J 8 , J 11 , J 18 , J 20 − J 41 . In general, one can try to perform a change of variables using with rational functions ψ 1 , ψ 2 and ψ 3 , such that P (ψ 1 (x 1 ), ψ 2 (x 2 ), ψ 3 (x 3 )) becomes a perfect square. Several algorithms have been devised in the past to identify rational parametrizations for square roots that appear in certain Feynman integral calculations [64], or to  demonstrate that rationalization is not achievable [39,[65][66][67][68]. In particular, it is simpler to find a transformation for rationalizing the square roots with one-variable dependence. The condition to find such a parametrization is that the degree of the polynomial involved must be less than or equal to 2. However, in some cases, the presence of a single square root of a multi-variable polynomial or a polynomial with a degree greater than 2 makes the subject of rationalizability more complicated [39,64,[66][67][68]. There exists a number of examples where the analytic computation of Feynman integrals involve non-rationalizable square roots of degree ≥ 3, for which the results are expressed in terms of a more general class of functions such as elliptic multi-polylogarithms (eMPLs) (iterated integrals over moduli space of torus) [42-44, 48, 69-89] and K3 surfaces [64,65,[90][91][92]. Recent studies show that even more complicated functions are expected to appear for massive cases at higher orders in scattering amplitudes in various QFTs which can be expressed as integrals over manifolds in higher dimensions. In many cases, these manifolds are found to be Calabi-Yau [65,90,[92][93][94][95][96]. In order to study the square root r and to understand its connection with a Calabi-Yau threefold variety, discussed before in the literature [94,96], we proceed as follows. We start by noting that P (x, y, z) in the square root r we obtain after simultaneously rationalizing three out of the four square roots of equation 3.3 P (x, y, z) = r 2 = y 2 z 2 + x 4 y 2 z 2 − 2xyz(z + y(1 + z(−2 + y + z))) − 2x 3 yz(z + y(1 + z(−2 + y + z))) + x 2 (−2y(−1 + z) 2 z − 2y 3 (−1 + z) 2 z + z 2 + y 4 z 2 + y 2 (1 + z(4 + z(−6 + z(4 + z))))) (3.7) is a degree 8 in-homogeneous polynomial with no repeated roots. Then we assign weight one to all the three coordinates x, y and z to realize this surface as a projective variety. In the next step, we introduce an additional auxiliary coordinate 'l' to homogenize P (x, y, z). As a result, we obtain a homogeneous polynomial P 8 (x, y, z, l) of overall degree 8 in 4 variables. Writing this hypersurface as Q(x, y, z, l, r) = r 2 − P 2 8 (x, y, z, l) = 0, (3.8) we can identify it as a degree 8 hypersurface in 4-dimensional weighted projective space WP 1,1,1,1,4 , with weight 4 assigned to r. The sum of the weights of this weighted projective space is 8 which is exactly equal to the degree of the polynomial Q, known as a Calabi-Yau threefold CY 3 in WP 1,1,1,1,4 [96] 1 . With this identification, one can obtain the Hodge structure and Euler characteristic that characterize the Calabi-Yau manifold. After establishing this interesting connection, we proceed to bring the epsilon-factorized differential equation system into a dlog-form, and express our results for the master integrals in terms of Chen's iterated integrals with dlog one-form kernels. Even though from a mathematical point of view, it might be possible to obtain a representation of these iterated integrals in terms of elliptic multiple polylogarithms (eMPLs) [48] with such kernels, we would like to keep this discussion for the future and proceed keeping the phenomenological applications of our results in mind. More details on one-forms and the motivation behind using iterated integrals over dlog one-forms for phenomenological applications are provided in the following section.

dlog one-forms
We Taylor expand the canonical basis integrals J k around = 0 as Putting it in the -factorized differential equation of 3.1, allows us to express each J (j) k in terms of Chen's iterated integrals up to a boundary constant. Chen's iterated integrals are defined as follows. Let there be a path γ on an n-dimensional manifold M , γ : [0, 1] → M , where γ(0) is the starting point and γ(1) is the endpoint. Let the set of differential oneforms be denoted by ω i and their pullback to the interval [0, 1] be given as γ * ω j = f j (λ)dλ. The k-fold iterated integral over the one-forms is defined as (3.10) and, the 0-fold iterated integral satisfies I γ (; λ) = 1 [97]. These iterated integrals obey the shuffle algebra. Very often Feynman integrals give rise to one-forms ω's which are of the form dlog p i (x 1 , .., x m ) where x i , ..., x n are the coordinates and p i (x 1 , ..., x m ) is an algebraic function, also known as a letter. Multiple polylogarithms (MPLs) are special cases of iterated integrals when p i is rational and linear-reducible in the variables x i . For more details regarding the conditions of expressing iterated integrals as MPLs, one may refer to [98]. For non-rational p i (algebraic letters with square roots), there is no general algorithm to express the iterated integrals in terms of MPLs. Nevertheless, in physics, many examples of iterated integrals with dlog-forms with non-rational arguments appear, which can also be evaluated in terms of MPLs [48]. In order to present the analytic results in terms of iterated integrals over algebraic dlog one-forms, we need to bring the epsilon-factorized differential equation system to a dlog-form. In the dlog-form, the entries of the matrixÃ are Q-linear combination of dlog one-forms so that From the matrixÃ, a set of rational letters can be trivially obtained for the rational oneforms. For the one-forms with a dependence on r, getting a dlog one-form and therefore the list of non-rational letters, is not a trivial exercise. To construct the non-rational letters, we implement the algorithm from [38,99] explained in the following: 1. We find the list of all the rational letters and the list of square roots. In our case, we have only one square root r with a degree-8 polynomial.
2. We construct all possible monomials upto degree 8 using the rational letters including r 2 . Further, we identify those monomials which can be factorized as (q a + r)(q a − r). This factorization gives the list of q a 's.
3. Using these q a 's, we construct ansätze l a for dlog one-forms as l a = q a + r q a − r . (3.12) 4. Using these ansätze, we fit all the non-rational one-forms which gives us a minimal set of dlog one-forms p a .
In recent times, iterated integral representation with dlog one-forms has been shown to be very efficient for phenomenological applications, see for example [100,101]. Iterated integrals with logarithmic one-forms have a clear branch cut structure, and the results expressed in terms of these functions have a compact analytic form. For a numerical evaluation of these functions, we can use a local series expansion of the one-forms combined with path-decomposition property satisfied by iterated integrals, see for example [102], to integrate from one phase-space point to another. Logarithmic one-forms can be series expanded, and since they have a power-log expansion a numerical evaluation is easy to implement. A representation in terms of iterated integrals is even useful for cases where public tools cannot be used due to an absence of power-log representation of the oneforms [103].

Results and checks
Regarding the results presented in this article, we numerically evaluate the iterated integrals using an in-house implementation in Mathematica. For the reader's convenience, we explain the steps that can be taken to do the numerical evaluation of iterated integrals in the following: 1. we parametrize the one-forms on a path, 2. we series expand the one-dimensional one-forms around a point on the path, 3. we perform the iterated integration of the expanded one-forms.
We might need multiple path segments if the phase-space point is far from the boundary point. In these cases, we perform the steps mentioned above on each path segment and use the path-decomposition formula [97,103] to obtain the final numerical result. For phenomenological applications, the numerical evaluation of these functions can also be combined with several publicly-available tools for fast evaluation. For example, one can set up a differential system just for the iterated integrals and evaluate them with generalized power series expansions [104] using tools like DiffExp [105] or parametrize the one-forms on a path and use GiNaC [106].
As explained in the previous section, the result of master integrals can be written as Taylor expansion in . Since our results are applicable to calculations involving two-loop integrals, we calculate master integrals only up to O( 4 ). We provide results for all the 41 master integrals of canonical basis in terms of the iterated integrals with the dlog one-forms in the supplementary material attached to this paper. We want to emphasize once again that the choice of the coordinate system given in equation (3.4) allows us to write the iterated integrals appearing in J 2 − J 7 , J 9 , J 10 , J 12 − J 17 , J 19 in terms of MPLs.
The complete results for the master integrals depend on the boundary terms. To obtain the analytic form of these boundary constants, we first evaluate the integrals in a regular limit [62,107]. We then match these values against the functional part of the results by first evaluating the iterated integrals up to many (100) digits and then using PSLQ [108] to extract the analytic constants. The analytic expressions of all the boundary constants B i = (J i ) |(x=0,y=0,z=0) in equation 3.9 upto O 4 are as follows: To verify the correctness of our results, we numerically evaluate the master integrals at multiple phase-space points. For this, we parametrize the one-forms on a path, use a series expansion of the integrands around (x = 0, y = 0, z = 0), and integrate the iterated integrals at some phase-space point. We then match the results against the numerical values obtained from pySecdec [109] and AMFlow [107] at the same phase-space point and find good agreement (∼ 100 digits). We present the numerical results of one of the master integrals from the top sectors J 41 up to O( 4 ) evaluated at point (x, y, z) = (0.5, 0.5, 0.5) up to 90 decimal places in table 2.  Table 2. Numerical values of J 41 evaluated at (x, y, z) = (0.5, 0.5, 0.5) for first five orders of using analytic results.
-17 -We have presented the analytic results for the master integrals that contribute to the mixed electroweak-QCD corrections in H → ZZ * decay. We obtain a canonical form of the differential equations for all the master integrals. Despite the simultaneous nonrationalizability of all the square roots that appear in the differential equations, we construct ansätze to find dlog-forms having algebraic dependence on the square roots. Our results are obtained in terms of Chen's iterated integrals order-by-order in , the parameter of dimensional regularization. Numerical evaluations of the result have been performed, and checks have been carried out using publicly available tools. The results presented here are straightforwardly applicable for computing the partial decay width of H → ZZ * → 4 at O(αα s ) accuracy and, in any production or decay process involving two-loop threepoint Feynman diagrams with four mass scales. Our analytic results can be used in event generation codes for Higgs production and decay to carry out phenomenological studies and data analysis for future collider experiments such as e + e − , µ + µ − and FCC-hh colliders.