Two-loop triangle integrals with 4 scales for the $HZV$ vertex

We calculate analytically the two-loop triangle integrals entering the $\mathcal{O}(\alpha\alpha_s)$ corrections to the $HZV$ vertex with $V=Z^*,\gamma^*$ using the method of differential equations. Our result provides a prototype to study the analytic properties of multi-loop multi-scale Feynman integrals, and also allows fast numeric evaluation for phenomenological studies. We apply our results to the leptonic decay of the Higgs boson and to $ZH$ production at electron-positron colliders. Besides the top quark loop, we include also the bottom quark loop contributions, whose evaluation takes a lot of time using purely numeric methods, but is very efficient with our analytic results.


INTRODUCTION
The HZV vertex, where V = Z * , γ * , is relevant for several important processes in Higgs physics, such as the leptonic decay H → Z * Z * → 4l which is crucial in the discovery of the Higgs boson and is also important for measuring its mass and width. The vertex also enters the production of the Higgs boson through vector boson fusion (ZZ → H) at the Large Hadron Collider (LHC), as well as the associate production of the Higgs boson with a Z boson (pp → ZH) at the LHC and at future Higgs factories (e + e − → ZH). Theoretically, the HZV vertex probes the electroweak symmetry breaking mechanism of the standard model, and is sensitive to alternative models of the Higgs sector. For example, if the Higgs is a composite Goldstone boson, the tree-level HZZ coupling is modified by O(v 2 /f 2 ) which can be as large as a few percent if the symmetry breaking scale f is in the TeV range. Therefore, precise knowledge of the HZV vertex is essential for the quest to better understand the properties of the Higgs boson and the mechanism of electroweak phase transition.
Experimentally, the HZZ vertex can be measured with an uncertainty of 1.5% at the High Luminosity LHC (HL-LHC) with an integrated luminosity of 3000 fb −1 [1] (in the κ-interpretation). At future Higgs factories such as the Circular Electron Positron Collider (CEPC), the Future Circular Collider (FCC-ee) and the International Linear Collider (ILC), the experimental precision for the HZZ vertex can reach the sub-percent level to about 0.2% [2][3][4].
In order to match the accuracy of future experimental measurements, it is necessary to provide high precision theoretical calculations for the HZV vertex within the standard model perturbation theory. The next-toleading order O(α) corrections were given in [5,6]. At the two-loop level, we are concerned with the O(αα s ) mixed QCD-electroweak corrections. A typical diagram is depicted in Fig. 1, where the heavy quark Q running in the loop can be the top quark t or the bottom quark b. To fix the notation, we consider the momenta assign- ments as V (q) → Z(p Z ) + H(p H ). The integrals we need to calculate are of the form where d = 4 − 2 in dimensional regularization, k 1 and k 2 are loop momenta, {a i } ≡ {a 1 , a 2 , a 3 , a 4 , a 5 , a 6 , a 7 } are the powers of the propagators D i in the integrand, and the propagators are given by The results of the integrals are functions of the 4 Lorentz invariants m 2 Q , q 2 , p 2 Z and p 2 H . We employ the method of integration-by-parts to reduce all these integrals into 41 master integrals, using the program packages FIRE5 [7], LiteRed [8] and Reduze2 [9]. These master integrals can be computed numerically using the method of sector decomposition [10], as performed in [11,12]. Note that the numeric integration is highly resource-demanding if any of q 2 , p 2 Z and p 2 H is larger than 4m 2 Q , which is the case if one considers collider energies above the tt threshold, or if one wants to take into account the contribution from bottom quark loops.
In the case where q 2 , p 2 Z , p 2 H < 4m 2 Q , one may perform a series expansion of the integrals in 1/m 2 Q . This approach has been taken in [11]. The benefit of this method arXiv:1905.11463v1 [hep-ph] 27 May 2019 is that after the expansion (which can be done at the amplitude level), the remaining integrals are single-scale ones which can be easily evaluated to analytic expressions. Therefore, one can implement the result straightforwardly into any event generators for phenomenological analyses.
If q 2 > 4m 2 Q but p 2 Z , p 2 H < 4m 2 Q , the 1/m 2 Q expansion fails but one may instead employ an expansion in powers of p 2 Z and p 2 H . This is in spirit very similar to the method of Ref. [13], where a small m H expansion is performed for Higgs boson pair production at the LHC. This kind of expansion is generically valid for top quark loops, but will not work for bottom quark loops.
Given the disadvantage of the purely numeric method and the limited applicability of various approximations, the goal of this paper is to provide an exact analytic solution to the master integrals appearing in the O(αα s ) corrections to the HZV vertex. The analytic expressions are valid for arbitrary values of the internal mass and the external momenta. This will allow fast numeric evaluations at any phase-space point, and will also serve as a prototype for analyzing the structure of loop integrals with many scales.

ANALYTIC SOLUTION
To obtain the analytic solution for the master integrals, we employ the method of differential equations [14,15]. We define the dimensionless variables We are able to find a basis f (x, y, z, ) of the master integrals which satisfies the canonical-form differential equation [16] d f (x, y, z, ) = dA(x, y, z) f (x, y, z; ) where A i are constant matrices independent of kinematic variables and the dimensional regulator. The "letters" α i ≡ α i (x, y, z) are rational functions of x, y, z and the following 4 kinds of square roots: with the Kallen function given by The solution to the canonical-form differential equation (3) can be formally written as a path-ordered integral from the boundary point r 0 = (x 0 , y 0 , z 0 ) to the point r = (x, y, z) We choose the boundary point to be r 0 = (0, 0, 0), where the values of the master integrals are simply given by In practice, we are interested in the master integrals as power series in where the n-th order coefficient function with n > 0 can be represented as a linear combination of Chen's iterated integrals [17] with transcendental weight n of the form (9) These integrals can be mapped to the symbol representation which can be manipulated using its algebraic properties. In general, these iterated integrals can be converted to linear combinations of multiple polylogarithms (MPLs) [18]. If all the letters are rational functions, this conversion can be handled straightforwardly, at least for low weights. However, in our case the presence of the 4 kinds of square roots in Eq. (4) makes the conversion much more difficult. Therefore, we need to describe the procedure in more detail. Only 3 letters appear at weight 1, which are β i ≡ β(r i ) with r 1 = x, r 2 = y, r 3 = z and The inversion of the above formula gives Noting that β(t) = 1 for t = 0, we conclude that the weight-1 part of the master integrals can be entirely written in terms of log(β(x)), log(β(y)) and log(β(z)) which satisfy the right boundary condition (7). At weight 2, another 13 letters come into play, and the independent integrable symbol tensors can be chosen as where i, j and k take distinct values from 1, 2 and 3. Except for the last one, the above symbol tensors can be converted to MPLs with the help of the known symbol maps of logarithm and dilogarithm functions. The last symbol tensor in Eq. (13) involves all 4 kinds of square roots, which cannot be rationalized simultaneously. To derive its functional form, we exploit its iterated integral representation. Being an integrable symbol tensor, the corresponding combination of integrals is path-independent. We can therefore choose a specific path to evaluate the 3 terms separately, while keeping in mind that only the sum of them is meaningful (i.e., we should choose the same path for all 3 terms). We parametrize the path as t r, with t going from 0 to 1. As an example, the first term can be written as where we have exploited the fact that λ(x, y, z) is a homogeneous quadratic polynomial of x, y and z, so that R 2 (tx, ty, tz) = tR 2 (x, y, z) = tR 2 .
According to Eq. (12), the square root R 1 (tx) can be rationalized via a change of variable to u ≡ 1 − β(tx), and the integral I 1 is then given by (16) Note that in the above formula, x, y and z should be treated as constants when taking the differential (i.e., only u is the active variable). The integration over u can be performed using the definition of MPLs, and the first term then becomes and the functional form of the whole symbol tensor can be obtained via permutation. We now turn to the weight-3 part of the solution. With the presence of non-rational functions, it is a highly complicated task to determine integrable symbol tensors at high weights. Nevertheless, since the master integrals as a whole correspond to a path-independent combinations of iterated integrals, we can still choose a specific path to evaluate them. At weight 3, there can be two R 1 functions appearing together with R 2 in one symbol. Without loss of generality, we take them to be R 1 (x) and R 1 (y). Similar to the weight-2 case, we choose the path parameterized by t r. To obtain the solution in terms of MPLs, we now need to rationalize R 1 (tx) and R 1 (ty) simultaneously, which can be done by the change of variable (18) The integration over v can then be handled using the definition of MPLs. Note that this variable change may lead to a proliferation of terms, and should only be used when necessary.
Using the above method, we are able to express all the master integrals in analytic form in terms of MPLs up to weight-3. Note that it is very rare that such two-loop integrals with many scales can be solved in closed form, and our result serves as a prototype to study the analytic structure of multi-loop multi-scale Feynman integrals. For example, it is interesting to investigate the behaviors of the scattering amplitude in various asymptotic limits, such as the massless limit and the threshold limit. These will be presented in a forthcoming article [19].
Finally, in the weight-4 part of the solution, all the 4 square roots in Eq. (4) can appear in a single symbol. We are not able to simultaneously rationalize R 1 (tx), R 1 (ty) and R 1 (tz) via a variable change (with respect to t). It is therefore not possible to convert the symbols to MPLs using the above method. It remains possible that one can construct the solution in terms of polylogarithms, using the function arguments at weight-2 and weight-3 as hints. We leave this for future investigation. In the current work, we evaluate the weight-4 part as a one-fold numeric integral. This is sufficiently fast for practical purposes.

NUMERIC RESULTS
We now use our analytic result to calculate the O(αα s ) corrections to the e + e − → ZH production cross section and the H → ZZ * → 4l partial width. We consider both the top quark and the bottom quark in the loop. In our numeric calculations, we take m t = 173.3 GeV, Γ t = 1.35 GeV, m b = 4.78 GeV, m Z = 91.1876 GeV, m W = 80.385 GeV, m H = 125.1 GeV, α s (m Z ) = 0.118 [20].
In Fig. 2, we show the exact top quark loop contributions to the e + e − → ZH cross sections for a centerof-mass energy √ s ranging from 220 GeV to 500 GeV. Thanks to our analytic results, it is much faster to perform the numeric computation compared to purely numeric methods such as sector decomposition. Especially for the tt threshold region √ s ∼ 2m t , the purely numeric integration is badly convergent, while it poses no difficulty for our analytic formulas. In the plot we also show the result of large m t expansion derived in [11] up to  O(1/m 4 t ). As expected, the expansion behaves well for low energies, but ceases to be valid near or above the tt threshold.
In Fig. 3, we show the impact of adding the contributions from bottom quark loops to the ZH cross section. Again, computing this contribution is rather timeconsuming with sector decomposition, but is much faster with the analytic results at hand. Phenomenologically this contribution only amounts to a few percent of the O(αα s ) corrections (which is below 1 per mille of the total cross section), and is therefore not important.
We now turn to the leptonic decay H → 4l. The leading m 2 t enhanced contributions at O(αα s ) have been considered in [21]. Here we give the result for the exact O(αα s ) corrections including bottom quark loops. For simplicity, we consider the process H → Zl + l − and treat the leptons as massless. A more dedicated study, including the decay of both Z bosons and the lepton mass effects, will be presented in [19]. In Fig. 4 we show the O(αα s ) corrections to the differential decay rate dΓ/dM , where M is the invariant mass of the lepton pair. We incorporate both top quark and bottom quark loop contributions. Note in particular the kink at M ≈ 2m b , which is due to the Coulomb singularity at the bb threshold. A proper treatment of this region would require resumming the Coulomb exchanges as well as dealing with non-perturbative effects, which is beyond the scope of this work.

SUMMARY AND OUTLOOK
In this paper, we calculate analytically the two-loop triangle integrals entering the O(αα s ) corrections to the HZV vertex. We derive the canonical-form differential equations for the 41 master integrals appearing in the calculation. For integrals with 4 mass scales, these differential equations are not easy to solve due to the presence of many non-rational functions. We are able to find fully analytic solutions up to weight 3 in terms of multiple polylogarithms. We apply our results to the e + e − → ZH production cross section and the H → ZZ * → 4l decay width, including both top quark loops and bottom quark loops. For the bottom quark loop contributions, and for cases when the collider energy is near the tt threshold, the integrals are rather time-consuming using purely numeric methods such as sector decomposition. This poses no difficulty for our analytic results, whose numeric evaluation is efficient for all phase-space points.
Loop integrals with many mass scales are very common in electroweak physics, Higgs physics and top quark physics. However, it is not easy to evaluate them in closed form, especially at high orders in . Our result serves as a prototype to study the analytic structure of multi-loop multi-scale Feynman integrals. For the HZV vertex, it is interesting to study its behaviors in various asymptotic limits. For example, near the tt threshold √ s ∼ 2m t , it is expected that the amplitude can be factorized in the framework of non-relativistic effective field theory [22][23][24]. It is interesting to see how such a factorization is achieved at O(αα s ) using our analytic results. Another interesting region is the small internal mass limit, which is relevant for bottom quark loops. In this limit, logarithms of the internal mass may develop at each order in perturbation theory, which have been studied in [25,26] for the Hgg vertex. It is interesting to investigate the same limit for the HZV vertex, where more scales come into play. These studies will be presented in a forthcoming article [19].
Finally, although we have obtained the analytic results at weight-3 in terms of MPLs, their expressions are rather lengthy, mainly due to the variable change Eq. (18). It is possible that we can find a suitable basis of trilogarithm functions to shorten the expressions. This may also help us to construct an analytic form for the weight-4 part of the solution. It is interesting to investigate these in the future.