On the evaluation of two-loop electroweak box diagrams for e + e − → HZ production

U.S.A


Introduction
Since the discovery of the Higgs boson in 2012 [1,2], experiments at the Large Hadron Collider (LHC) have measured many properties of this new particle, in overall good agreement with the predictions of the Standard Model (SM) [3]. However, in most models beyond the SM, one would expect deviations of the Higgs boson couplings at the per-cent level [4], which is beyond the achievable precision at the LHC.
For this reason, several proposals have been made for so-called e + e − Higgs factories, operating at center-of-mass energies of 240-250 GeV: the International Linear Collider (ILC) [5,6], the Future Circular Collider (FCC-ee) [7], and the Circular Electron-Positron Collider (CEPC) [8]. These machines would be able to study the Higgs boson through the process e + e − → HZ in a clean environment and produce per-cent level precision measurements of the dominant Higgs couplings. The HZZ coupling can be extracted from the HZ production cross-section itself, σ HZ , which is anticipated to be measurable with a precision of about 1.2% at ILC, 0.4% at FCC-ee, and 0.5% at CEPC.
The interpretation of σ HZ in terms of the HZZ coupling requires precise theoretical predictions for the process e + e − → HZ, including radiative corrections. The next-toleading order (NLO) corrections within the SM have been known since a long time for unpolarized beams [9][10][11], and more recently for polarized beams [12]. Mixed electroweak-QCD O(αα s ) corrections have been computed in refs. [13][14][15], which required the evaluation of two-loop self-energy and vertex diagrams. Given the relatively large decay width of the Z boson, the predictions can be further refined by including the Z-boson decay at the same order, i.e. by computing corrections to the process e + e − → Hff . The NLO electroweak contributions to this process for the final states f = ν e and f = e have been studied in refs. [16][17][18], and the O(α) and O(αα s ) corrections for f = µ have also become available recently [19].
The numerical impact of the one-loop corrections is at the level of 5-10%, with a dominant contribution stemming of initial-state radiation (ISR) of soft and collinear photons.

JHEP04(2021)179
These ISR effects are enhanced by logarithmic terms of the form log(s/m 2 e ). In the soft and collinear limit, these logarithmic terms are process-independent, and higher-order ISR contributions can be included through the structure function method [20,21]. The impact of ISR on e + e − → HZ has been recently studied [22,23]. It was found that, when including third-order corrections in the structure function [24][25][26], the uncertainty from missing higher ISR orders is at the level of 10 −4 and thus negligible [23].
The O(αα s ) contributions modify the HZ cross-section by about 1.5% when parametrizing the elecroweak couplings in terms of α, and about 0.4% using the Fermi constant G µ instead. These corrections are sizeable and must be taken into account for analyses at future e + e − Higgs factories. The largest unknown higher-order contribution stems from electroweak two-loop effects, which are expected to have an impact at the level of O(1%) [27], and thus also have to be included. Given that the width-to-mass ratio of the Z boson, Γ Z /m Z ∼ 2.7%, is comparable to one order in electroweak perturbation theory, it is sufficient to compute these NNLO electroweak corrections for the process e + e − → HZ with an on-shell Z boson, whereas the full process e + e − → Hff should be treated at the NLO level. These two contributions can be consistently combined by performing an expansion about the pole of the Z boson [28][29][30][31].
Among the NNLO electroweak corrections, diagrams with closed fermion loops are typically dominant, due to the large top-quark Yukawa coupling and the large multiplicity of fermions in the SM. 1 Within this class of diagrams, the most challenging piece are planar and non-planar two-loop box graphs with top quarks inside one sub-loop. 2 Even when neglecting all fermion masses besides the top quark, these contributions depend on up to four independent mass scales (m H , m Z , m W , m t ), as well as two additional momentum scales (which can be represented by the Mandelstam variables s and t). Therefore it is difficult to find analytical solutions, since the expressions will be impractically large and may require the development of new special functions. On the other hand, generic numerical methods (such as numerical integration over Feynman parameters [33]) are highly computationally intensive.
In this paper, a more efficient numerical method for the evaluation of two-loop box integrals is proposed. It is based on a combination of a dispersion relation and Feynman parameters for one of the two sub-loops [34]. The method of ref. [34] is extended to enable the direct evaluation of tensor integrals (rather than attempting to reduce them to a set of master integrals). 3 This approach leads to three-dimensional numerical integrals for the two-loop boxes, which can be evaluated with about four-digit precision within minutes on a single CPU core.
In the following section, the derivation of the numerical integral representations for the planar and non-planar two-loop box diagrams is discussed in detail. Section 3 describes the application of this method to the evaluation of two-loop box diagrams contributing to the process e + e − → HZ, including several important aspects of the numerical implemen-JHEP04(2021)179 tation, as well as numerical results for these diagrams. The main findings of this paper are summarized in section 4.

Planar diagrams
To illustrate the method, let us initially consider a basic scalar planar integral, which contains the propagators of the diagram in figure 1 (top-left) but simply 1 in the numerator: The extension to non-trivial tensor structures in the numerator will be discussed below.
The following approach is based on the technique used in ref. [37], which is makes use of the basic dispersion relation for the one-loop self-energy function B 0 , where D is the space-time dimension and λ(a, b, c) = (a − b − c) 2 − 4bc. This dispersion relation is derived from the analytical properties of the B 0 function: for complex p 2 , B 0 (p 2 , m 2 1 , m 2 2 ) has a branch point at p 2 = (m 1 + m 2 ) 2 , with the branch cut on the real-axis interval ((m 1 + m 2 ) 2 , ∞). When using Cauchy's theorem,

JHEP04(2021)179
, one must choose a contour C that circumvents the branch cut, as illustrated in figure 2 (left). The discontinuity ∆B 0 accounts for the difference of B 0 (σ, m 2 1 , m 2 2 ) for values of σ just below and just above the branch cut. The contour is closed with a circle at infinity, which gives vanishing contribution for sufficiently small dimension D.
In order to apply this relation to the planar box diagram, it is useful to introduce Feynman parameters for the propagators that depend only on loop momemtum q 2 [34]: Then the q 2 loop can be expressed as

JHEP04(2021)179
where we have introduced the short-hand notation and used the fact that ∆B 0 (σ 0 , m 2 , m 2 q ) = 0. Unfortunately, the σ integral blows up at the lower boundary, and the term in [ ] is also divergent for σ → σ 0 , whereas only the sum of the two is finite. To circumvent this problem, one can modify the integrand according to (2.8) Here the extra term in the integrand of dσ is added back in integrated form, where the function ∂ 2 m B 0 can be expressed in terms of basic logarithms (see appendix). With the modified integrand, the boundary term in eq. (2.6) evaluates to zero.
Inserting eq. (2.8) into the remainder of the q 1 loop integral, one obtains and D 0 is the well-known scalar one-loop box function [38][39][40].
Since the double box diagrams are UV finite, all expressions in eq. (2.9) can be computed for D=4 dimensions.
The full diagram respresented by figure 1 (top-left) contains additional terms with momenta q 1,2 in the numerator stemming from the Dirac propagators and vertex structures. For terms depending on q 2 , it is convenient to perform a Passarino-Veltman decomposition [40,41] of ∂ 2 m B 0 (q 2 1 , m 2 , m 2 q ) after introduction of the Feynman parameters. As a first step, let us shift the integration momentum to q 2 ≡ q 2 + k : . (2.10) The terms with powers of q 2 in the numerator can be decomposed according to can then be represented through a dispersion relation in the same manner as above: (2.12) Similarly, the q 1 loop will in general contain terms with different powers of q 1 in the numerator, some of which in fact originate from eq. (2.11). These lead to Passarino-Veltman functions D 1 , D 2 , D 3 , D 00 , etc. [41], which can be evaluated numerically by using, for example, the techniques introduced in refs. [39,42]. In some cases, there are cancellations between terms in the numerator and denominator, resulting in C 0 , C 1 , C 2 , C 00 , . . . and B 0 , B 1 , B 00 , . . . functions.
For the double-box diagrams in figure 1 tensors up to rank 3 in at least one of the two loops are encountered.

Non-planar diagrams
The approach in the previous sub-section can be adapted also to the case of non-planar box diagrams. Let us begin with the scalar non-planar integral corresponding to the diagram in figure 1 (top-right): . (2.13) Introducting two Feynman parameters, the q 2 loop can be written as where

JHEP04(2021)179
In the last step of eq. (2.14), a threshold subtraction has again been utilized to ensure that the σ integral is convergent at the lower boundary. Together with the q 1 integral, I np then becomes where all components of the integrands are well-known analytical functions. As before, the extension to tensor integrals can be realized by using dispersion relations for ∂ m 1 ∂ m 2 B 1 , ∂ m 1 ∂ m 2 B 00 , etc. (see appendix for explicit formulas). Similarly, the usual Passario-Veltman tensor functions D 1 , D 2 , etc. can be used for tensor structures in the q 1 loop integrals.
An additional complication arises for the non-planar diagram with W bosons, in which case m q = m b ≈ 0, so that m 2 1 is negative. As a result, the branch point σ 0 of the ∂ m 1 ∂ m 2 B 0 function is in the lower complex half-plane rather than on the real axis (see figure 2), so that the dispersion relation (2.2) must be modified. One option is to choose a contour along the real axis, which is closed via a semi-circle in the upper complex half-plane, leading to Using this relation, I np can be expressed as The i in the last mass parameter of the D 0 function is important to properly define its result for all values of σ. In fact, as also discussed in the next section, it turns out to be necessary to include a small numerical value for when using LoopTools [46] for the evaluation of certain Passarino-Veltman functions.

Implementation and numerical results
In this section we describe how the approach described in the previous section has been applied to the calculation of all box diagrams of the form in figure 1. The results presented in this section are based on two independent realizations of the calculation, in order to enable cross-checks between the two. Both implementations employ Mathematica [43] as the framework for algebraic manipulations and FeynArts 3 [44] for the generation of diagrams and amplitudes in Feynman gauge. One implementation uses FeynCalc 9 [45] for carrying out the Lorentz and

JHEP04(2021)179
Dirac algebra and then divides the expressions into individual tensor integral terms, as discussed at the end of section 2.1. Each of these terms is integrated separately within C++, using the LoopTools 2.15 [46] package for the Passarino-Veltman functions and the adaptive Gaussian quadrature integration routine from the Boost library [47]. The integration results are then added up to obtain full diagram results.
The second implementation performs the Lorentz and Dirac algebra with in-house routines and then tranforms the expressions for complete diagrams into a single integrand each. The numerical integration is carried out in C++ using the adaptive Gauss-Kronrod integration routine from TVID [48], which is based on the Quadpack library [49]. It also uses LoopTools for the Passarino-Veltman functions in the integrand.
In light of the fact that the double-box integrals are UV-finite, it is advantageous to perform the Lorentz and Dirac algebra in 4 dimensions, thus avoiding any ambiguities in the treatment of γ 5 . Even though the sum of all box diagrams considered here is IR finite, individual diagrams with photons are IR divergent, and thus an IR regulator is required. A convenient choice is the use of a small photon mass, m γ , since it is trivially compatible with the 4-dimensional Lorentz and Dirac algebra.
It is advantageous to implement the three-dimensional numerical integrals in a nested structure, with the σ integral being the inner-most integral, since this makes the adaptive integration algorithms most effective. The achievable precision is limited by the double precision floating point algebra used in the default compilation of LoopTools. In fact, numerical instabilities are typically encountered near the lower and upper limits of the σ integration. These can be mitigated by introducing cut-offs at both ends, where δ 1 and Λ should be much larger than all mass and momentum scales in the matrix element. The error due to these cut-offs can be further mitigated by observing that the integrand approximately behaves like ∼ (σ − σ 0 ) −1/2 near the lower threshold and ∼ (A + B log σ)/σ 2 for large σ. Thus one can introduce additional correction terms, For the non-planar diagrams, additional instabilities occur for x ≈ y, when the Gram determinants for some Passarino-Veltman tensor functions vanish. Our two implementations use two different strategies for mitigating this problem: (a) splitting one of the two integration intervals, such that none of the Gaussian points of the x integration lies too close to the ones for the y integration; or (b) interpolating the y integration across a small interval, y ∈ [x − ∆x, x + ∆x]. A reasonable comprise between accuracy and stability is achieved for ∆x ∼ O(10 −2 ). Both methods yield consistent results, and the impact of varying ∆x by a factor of a few can be used as a contribution to the final error estimate. Finally, the evaluation of the non-planar W W box requires an explicit value for the Feynman i , see eq. (2.18), to avoid instabilities in LoopTools for negative σ. A value of = 10 −9 |σ| is chosen for the results presented below. The results are not significantly affected when increasing this value by a factor 10 or using a constant value ∼ 10 −5 .
In the following, numerical results will be presented for the different classes of box diagrams, which are distinguished by the gauge bosons V 1,2 appearing inside the loops. The numbers are obtained by contracting the matrix elements for the two-loop box diagrams, M 2 , with the tree-level matrix element M 0 , averaging over e ± helicities and summing over the final-state Z-boson polarization states.
Using the inputs in table 1 (a), we obtain the numbers shown in table 1 (b). For the diagrams with photons, the dependence on the photon mass regulator only drops out when adding planar and non-planar diagrams, as illustrated in figure 3. Also shown in table 1 (b) is an esimate of the precision, as obtained by varying the lower and upper cut-off of the σ by one order of magnitude each. For the non-planar diagrams, the impact of varying the width ∆x of the window around y = x by a factor 2 is also considered. The integration times for each line in table 1 (b) range from a few minutes up to about half an hour on one CPU core. 4 Figure 4 shows the dependence on the scattering angle θ. The differential distributions are symmetric, since each subset of box diagrams has a t ↔ u crossing symmetry, where t, u are the usual Mandelstam variables.
One can see that the diagrams with W bosons produce results that are about one order of magnitude larger than the ones with neutral bosons. This may be explained by the fact than the effective W W ZH interaction (corresponding to the fermion loop in our

Summary
Box diagrams in 2 → 2 process can be efficiently evaluated with a numerical method that combines Feynman parametrization and a dispersion relation for one sub-loop, while standard analytical expressions are used for the other sub-loop. Tensor structures in the numerator can be handled by adjusting the dispersion relation for the first loop and using Passarino-Veltman reduction for the second sub-loop. The resulting three-dimensional numerical integrals can be efficiently evaluated using nested adaptive one-dimensional integration algorithms.
The efficacy of the technique has been demonstrated by computing planar and nonplanar box diagrams with top quarks contributing to the two-loop electroweak corrections for the process e + e − → HZ. Infrared divergencies from QED can be controlled with a JHEP04(2021)179 photon mass, without loss of numerical precision. Results with a relative uncertainty of about 0.1% can be obtained in a few minutes on a single CPU core. The longest run-time (about half an hour) is required for diagrams with a physical cut in the fermion sub-loop, which occurs for the non-planar topology in the top-right of figure 1 with V 1,2 = W and q = b. In this case a modified version of the dispersion relation is used, with an integration contour along the entire real axis instead of just the positive real axis.
It should be noted that our current implementation of the other diagrams is limited to center-of-mass energies below the tt threshold, √ s < 2m t , because otherwise a physical threshold would open up in the fermion sub-loop there as well. Nevertheless, an extension to higher center-of-mass energies could be achieved by using the modified dispersion relation for all diagrams, even though it may come at the cost of a slight loss of accuracy and increased integration time.
The numerical precision is primarily limited by the accuracy of the evaluation of the basic one-loop Passarino-Veltman functions B 1 , C 1,2 , C 00,11,12,22 , D 1,2,3 , etc. A high level of numerical precision becomes important (a) for large σ, where the full integrand falls off ∼ σ −2 but individual terms in the integrand decay only ∼ σ −1 , and (b) when the Gram determinant of some Passarino-Veltman functions vanishes at particular points in the integration region. In our current implementation, LoopTools [46] with double-precision floating point arithmetic is used for this purpose. Improvements could be made by using quadruple precision numbers or by performing expansions in the regions where numerical instabilities are encountered. However, a relative precision of 0.1% for the evaluation of two-loop box diagrams is already sufficient for a range of important phenomenological applications, including the e + e − → HZ process at future Higgs factories.
The techniques described in this paper could, in principle, also be applied for the calculation of electroweak corrections to other 2 → 2 process, such as e + e − → W + W − .