Top quark mass effects in gg → ZZ at two loops and off-shell Higgs boson interference

We consider top-quark mass effects in the Higgs-interference contribution to Z -boson pair production in gluon fusion. While this production mechanism is formally of next-to-next-to leading order, its contribution is numerically important above the top threshold M 2 ZZ ¼ 4 m 2 t . This region is essential to constrain the width of the Higgs boson, and good control over the top-quark mass dependence is crucial. We determine the form factors that are relevant for the interference contribution at two-loop order using a method based on a conformal mapping and Pad´e approximants constructed from the expansions of the amplitude for large top mass and around the top threshold.


I. INTRODUCTION
A direct measurement of the Higgs boson width Γ H is not possible at the LHC or even the envisioned next generation of collider experiments. However, indirect constraints can be obtained at the LHC by studying the process pp → H → ZZð→ 4lÞ on the Higgs boson peak where the cross section depends on the combination g 2 Hgg g 2 HZZ =Γ H and off the peak where the measurement of the cross section constrains the product g 2 Hgg g 2 HZZ of the effective Higgs boson-gluon coupling g Hgg and the Higgs boson-Z boson coupling g HZZ , as proposed in [1][2][3]. 1 The same strategy can be employed with WW final states [6]. The latest studies from the LHC experiments give an upper limit of 14.4 MeV at 95% C.L. from the ZZ final state at ATLAS [7] and the value 3.2 þ2.8 −2.2 MeV from the combination of VV final states in CMS [8], close to the SM prediction Γ SM H ¼ 4.10 AE 0.06 MeV [9]. Measurements of the Higgs boson signal at large invariant mass can also be used to directly constrain physics beyond the Standard Model in the Higgs sector [10][11][12][13].
Here, we focus on the loop-induced continuum gluon fusion process gg → ZZ and in particular its interference with the off-shell Higgs contribution gg → H Ã → ZZ. Despite the narrow width of the Higgs boson, these interference effects are sizable with 10% of the Higgs signal stemming from the off-shell region where the invariant mass of the two decay products is greater than 2m Z [1] and higher-order corrections are required to control the uncertainties. The Higgs-mediated amplitude only depends on two scales, the mass m q of the quark in the loop and the invariant mass M ZZ of the final state. Next-toleading order (NLO) corrections with the full quark-mass dependence have been known for some time [14][15][16][17], and the top-quark mass dependence at next-to-next-to leading order has been reconstructed very recently [18] (see also [19]). On the other hand, the continuum amplitude depends on four scales m q , m Z , M ZZ , and the transverse momentum p T of one of the Z bosons, and the exact result is only known at leading order (LO) [20] while an analytic NLO calculation appears extremely challenging. In the massless limit m q ¼ 0, the two-loop amplitude has been determined in [21][22][23] and the NLO cross section in [24]. Recently, also the quark-gluon channel has been included [25].
The contribution from top quarks at two-loop order has been computed in a large-mass expansion (LME) [24,26,27] and is known up to 1=m 12 t . While the contribution from massless quarks dominates the interference correction at small invariant masses, the top-quark contribution is of the same size near the top threshold M ZZ ¼ 4m 2 t and dominates in the large invariant-mass regime. Since the LME ceases to provide a reliable description above the top threshold, the authors of [27] have improved their prediction by a conformal mapping and the construction of Padé approximants based on the available number of LME coefficients. In [28], we have extended this method by considering the expansion around the top threshold in addition to the LME and demonstrated that the top-mass effects can be reproduced correctly by comparing results for the two-loop amplitude for gg → HH with the numerical calculation performed in Refs. [29][30][31]. 2 In this work, we consider the form factors of the continuum gg → ZZ amplitude that are relevant for the interference contribution at one and two loops. The nonanalytic terms in the expansion around the top threshold are computed up to at least order ð1 − zÞ 4 , where z ¼ M 2 ZZ =ð4m 2 t Þ þ i0, and used to construct Padé approximants. Together with the exactly known real NLO top quark [27,38] and the massless quark corrections [21][22][23][24][25], this is sufficient to determine the full NLO interference contribution with realistic top-quark mass dependence.

II. FORM FACTORS FOR INTERFERENCE
Up to the two-loop level, the amplitude for the top-mediated nonresonant continuum production process gðμ; A; p 1 Þ þ gðν; B; p 2 Þ → Zðα; p 3 Þ þ Zðβ; p 4 Þ receives contributions from both box and double-triangle diagrams; see Fig. 1. The latter are known for arbitrary quark masses [27,39] and will not be discussed in the following.
The box amplitude jB AB μναβ i has a complicated tensor structure [20][21][22]40]. However, the interference with the Higgs-mediated amplitude is described by a single form factor. Adopting the conventions of [27], it takes the form The form factor can be decomposed into a vector and axialvector part where a t ¼ 1=2 and v t ¼ 1=2 − 4=3 sin 2 θ W denote the axial-vector and vector couplings for an up-type quark.
Mixed v t a t terms are forbidden by charge conjugation symmetry. The order in the strong coupling constant α s is indicated as follows: with i ¼ VV; AA. At order α 2 s , the renormalized form factors contain IR divergences, which cancel in the combination with real corrections, and we define the finite remainder by applying the subtraction [41] The one-loop form-factors jB ð1Þ i i are already finite; we define jF ð1Þ i i ¼ jB ð1Þ i i for the sake of a consistent notation.

A. The amplitude near threshold
Above the top threshold at z ¼ 1, the top quarks in the loop can go on shell which manifests as nonanalytic terms in the expansion of the form factors inz ≡ 1 − z, generating a sizable imaginary part. As shown in [28], the knowledge of these terms alone provides very valuable information for the determination of top-quark mass effects in our approach. The calculation of the nonanalytic terms is significantly simpler than that of the analytic contributions and was described in detail in [28] for the three leading nonanalytic expansion terms of the one-and two-loop form factors for gg → HH. For gg → ZZ, we expand the amplitude up to high orders inz ≡ 1 − z and therefore use the expansion by regions [42,43] to expand the full-theory diagrams instead of an effective field theory approach where a large number of effective vertices is required due to the deep expansion. We use QGRAF [44] to generate the Feynman diagrams which are processed and expanded using private FORM [45] code. The integration-by-parts reduction [46] is performed with FIRE [47] which is based on the Laporta algorithm [48].
Our results are given in Appendix A and an ancillary Mathematica file [49]. They are of the form   2 Recently, an independent numerical calculation [32], several approximations [33][34][35][36] which are consistent with the earlier results and a combined result [37] have appeared. 3 Note that this subtraction differs at order ϵ 0 from the one given in Eq. (2.14) of [27]. wheren 2 is n modulo 2, the coefficients are functions of the dimensionless variables r Z ¼ m 2 We use the symbol ≍ to indicate that terms which are analytic inz and currently unknown have been dropped on the right-hand side.
Threshold logarithms lnz and logarithms lnð−4zÞ related to massless cuts in the amplitude first appear at two-loop order. While we generally compute the expansion coefficients up to n ¼ 8, i.e., expand up toz 4 , we find that for the massless-cut contribution proportional to lnð−4zÞ more input is required to achieve a reliable Padé approximation. We therefore compute the corresponding coefficients b ðn;mÞ i;ln up to n ¼ 9.
As in Higgs pair production, there is no S-wave contribution to the form factors relevant for the interference and the leading nonanalytic terms involve thez-suppressed P-wave Green function [50].

B. Behavior for z → ∞
In addition to the LME and threshold expansions, we can exploit scaling information in the small-mass limit m t → 0 which corresponds to z → ∞. This does not require an additional calculation in this region but relies solely on the symmetries of QCD. The absence of infrared 1=m t power divergences as m t → 0 implies that the form factors can only show logarithmic behavior as z → ∞. Below we show that the difference vanishes as z → ∞. To prove this, we note that chirality is conserved in massless QCD and hence the four-point correlator of two vector currents, a left-handed and a right-handed current, which we denote in short by [V,V,V-A,V+A], vanishes in the limit of zero quark masses. 4 Using that the correlator [V,V,V,A] vanishes due to charge conjugation, we immediately conclude that [V,V,V,V]-½V; V; A; A → 0 as z → ∞. We exploit this below and reconstruct the top-mass dependence of jF ðiÞ VV i and jF ðiÞ AA i − jF ðiÞ VV i where we have one additional condition for the latter.

III. THE METHOD
We approximate the box form factors (2) using our approach from [28]. First, we introduce subtraction functions s The construction of such subtraction functions is detailed in [28], and we give the ones we explicitly need in Appendix B. Note that even after this subtraction the threshold and large mass expansions of the two-loop form factors still receive contributions proportional to a single logarithm L s ≡ lnð−4zÞ from diagrams with massless cuts. We therefore split the subtracted two-loop form factors into a constant and a logarithmic part and construct separate approximants for each part.
The top mass dependence is contained in the variable z and the conformal transformation [51] z is used to map the entire complex z plane onto the unit disc jωj ≤ 1. The physical branch cut for z ≥ 1 corresponds to the perimeter of this disc. Further cuts arise outside the physical region for z ≤ s t and z ≤ s u , where s, t, u are the usual partonic Mandelstam invariants. These cuts are mapped onto negative values of ω.
We then reconstruct the logarithmic and nonlogarithmic parts of the form factors using Padé approximants where the n þ m þ 1 coefficients a i , b j can be fixed by imposing the condition that the expansion of Eq. (8) in the LME and threshold region must reproduce the known coefficients for given, fixed values of r Z andx. The smallmass behavior discussed in Sec. II B is not used to further constrain the Padé coefficients, but is taken into account by a rescaling of the Padé ansatz. Hence, we use approximation functions of the form where P ðjÞ AA−VV is used to approximate the difference between the axial-vector and vector form factors, whereas the vector form factors in isolation are approximated using P ðjÞ VV . The limit z → ∞ corresponds to ω → −1 where the approximants in Eq. (8) approach a constant value. Thus, the rescaling equation (9) enforces the correct asymptotic behavior for z → ∞ discussed in Sec. II B and provides us with free parameters a R;i that can be varied in addition to the polynomial degrees n, m, k, and l to assess the stability of the approximation. We note that these variations are performed independently for all the terms in Eq. (9). Our final ansätze for the form factor approximation are then It should be noted that our ansatz does not account for and cannot reproduce the aforementioned t-and u-channel branch cuts in the unphysical region ω < 0. Even in the physical region not much is known from first principles concerning the convergence behavior of Padé approximation. In the absence of exact results, our approximation could be compared to, we therefore rely on heuristic arguments. The approximation method used in this work was shown to reliably reproduce the exact NLO correction for the very similar scenario of di-Higgs production [28]. The main difference is the structure of the couplings to the produced bosons in the considered box diagrams. However, as will be shown in Sec. IV, this difference appears to have no visible impact when applying the method at LO and the convergence behavior at NLO indeed turns out to be similar to the case of di-Higgs production. We find that including more input terms in the construction of the approximation stabilizes the prediction and reduces the estimated uncertainty.

IV. RESULTS
Before showing our results at NLO for the form factors, we can compare the LO form factors constructed as discussed in the previous sections with the full analytic result. We choose as input for the on-shell Z-boson and top quark masses and show results for two different values ofx in Fig. 2.
The plots contain the maximum information we have available from the LME at LO (see [27]) and our threshold  (11); the solid lines are the full result and the shaded regions are Padé approximants that were constructed using only the information from the LME (cf. text for details). expansion. By construction the Padé ansatz in Eq. (8) may contain poles anywhere in the complex ω plane whereas the functions it approximates are analytic apart from the branch cuts discussed in Sec. III. 5 Since we aim to reconstruct the form factors in the timelike region, we exclude approximants that exhibit poles for We obtain an uncertainty estimate for our results in the following way. For every phase space point, we calculate the mean and standard deviation for each contributing Padé approximant in Eq. (9). To this end, we vary the rescaling parameters a R;i in the region a R;i ∈ ½0.1; 10 ð 14Þ and vary [n=m] within jn − mj ≤ 3, where n þ m þ 1 is the number of available constraints. We construct 100 variants for each Padé approximant. Our final prediction then follows from the sum of the mean values of the Padé approximants, with an uncertainty obtained by adding the individual errors in quadrature. Figure 2 shows the Padé approximants from Eqs. (10) and (11) Fig. 2 with the points and shaded regions corresponding to the Padé approximation constructed from the LME only.
We conclude that the threshold expansion is essential for the reconstruction of the full top mass dependence above the top quark threshold.
We now turn to the NLO form factors. In Fig. 3, we show the results for the virtual corrections to the form factors jF ð2Þ VV i (upper panel) and jF ð2Þ AA i (lower panel) for two values ofx. Note that we do not include the double-triangle contribution to the form factors, as they have been computed analytically in [27]. As at LO, we include only the top quark contributions. The uncertainty associated with the Padé construction increases with M ZZ . Since we input information mainly at low M ZZ , this behavior is expected. With the exception of the vector form factor jF ð2Þ VV i for small transverse momenta (upper left panel in Fig. 3), we find that the Padé approximation based on the LME alone does not yield a realistic reconstruction of the top-quark mass effects of the form factors. In particular, the important axial-vector form factor suffers from very large uncertainties. We remark though that in [27] for the NLO cross section the Padé prediction was improved by a reweighting with the full LO cross section.
We note that jF ð2Þ VV i shows a small oscillation in the region of large M ZZ when the transverse momentum of the Z bosons is small as is evident from the upper left plot in Fig. 3. We trace the appearance of the second peak back to the contribution proportional to L s stemming from diagrams with massless cuts. As shown in Fig. 4, this logarithmic contribution is numerically important in the vector form factor for large M ZZ but significantly smaller than the nonlogarithmic contribution in the axial-vector form factor. In general, we find that the contribution proportional to L s shows worse convergence behavior than the nonlogarithmic terms when including more and more terms in the LME and the threshold expansion. This is shown in Fig. 5 where we compare our results from Fig. 3 to the Padé approximants obtained with the same procedure but only using threshold input up to the orderz 2 andz 3 .
We observe good convergence in the case of the axialvector form factor. On the other hand, the Oðz 2 Þ approximation for the vector form factor does not feature the oscillatory behavior described above and there is no overlap with the full approximation in a significant part of the phase space. However, the Oðz 3 Þ and Oðz 4 Þ results are in good agreement with the full approximation where we have also included the Oðz 5 Þ term in the coefficient of the logarithm L s to verify that this stabilization persists with the addition of higher orders in the threshold expansion. We conclude from this discussion, that the Padé approximation can be improved systematically when including higher orders in the various expansions. Nevertheless, we believe that the prediction for jF ð2Þ VV i should be taken with a grain of salt above M ZZ ≥ 500 GeV because of the slower convergence. We note also that we find very similar convergence behavior as in the case of Higgs boson pair production, discussed in [28]. As we found here for ZZ, the contributions from diagrams proportional to L s and hence massless cuts converge worse than the nonlogarithmic pieces.
In Fig. 6, we show the virtual corrections to the form AA i as it enters in the interference term with the Higgs boson exchange. We expect a large suppression of the vector contribution since v 2 f ≪ a 2 f and the LME of the vector form factor only starts one order higher [27]. Indeed, we find that v 2 f jF ð2Þ VV i has to be amplified by a factor of 300 to be comparable to the axial-vector contribution, cf. the dashed lines in Fig. 6. This clearly demonstrates that the interference term will be dominated by jF ð2Þ AA i, and we therefore choose not to modify the uncertainty estimate for the vector form factor. The fact that jF ð2Þ VV i is negligible compared to jF ð2Þ AA i allows us to make trustworthy predictions for the interference with the Higgs production with subsequent decay to Z bosons up to M ZZ → ∞, even though as stated above we trust our results for jF ð2Þ VV i only for M ZZ ≤ 500 GeV. The numerical implementation of the form factors is available as a FORTRAN routine on request and can be combined with existing computations of the massless loop contributions and the real corrections for the interference of the Higgs exchange with decay to ZZ with the continuum background.

V. CONCLUSIONS AND OUTLOOK
We have considered top-quark mass effects in the continuum process gg → ZZ, focusing on the form factors relevant for the NLO interference with the production of a Higgs boson and its subsequent decay into two Z bosons. We have presented a Padé-based approximation using information from an expansion around a large top quark mass and an expansion around the top quark pair production threshold.
At LO, we have shown that our Padé construction approximates very well the full top mass dependence of the form factors for the whole range of the invariant mass M ZZ of the Z bosons. At NLO, we provide a new prediction with very small uncertainties at small and moderate M ZZ , with an increased uncertainty toward large M ZZ . We expect that adding more information into the Padé construction at large M ZZ would improve the description also in this region.
Our results can be combined both with virtual corrections mediated by massless loops and the real corrections. The latter constitute a one-loop process and can therefore be computed with well-established techniques. The Padé construction can also be applied to the remaining form factors contributing to gg → ZZ, which do not interfere with the Higgs signal.
We note also that while in this work we have applied our method to the production of on-shell Z bosons, there is no obstruction for applying it also to off-shell Z-boson production. Indeed, the LME for off-shell Z-boson production is already known up to the order z 4 [24]. While a calculation of the full top mass dependence for on-shell Z bosons with numerical methods seems to be feasible with current techniques in a reasonable timeframe (see [53,54]), a computation of the off-shell form factors appears to be beyond the current state-of-the-art. where i ∈ fVV; AAg andn 2 is n modulo 2. The coefficients a, b are most conveniently written in terms of the two dimensionless ratios r Z ¼ and use the shorthand notation vanish. Furthermore, coefficients with m ¼ 0 and even n do not contribute to the imaginary part and are therefore not listed here. We have calculated the remaining coefficients a n;0 i ; b n;m i up to n ¼ 8 and the coefficients b ðn;0Þ i;ln up to n ¼ 9, obtaining the following results: AA;ln ¼ 0; ðA10Þ b ð3;1Þ b ð5;0Þ for the expansion of the axial-vector component. The corresponding coefficients in the expansion of the vector part read a ð3;0Þ  VV ¼ 128π 2 315ð1 − 2r Z Þ 6 ½−58 þ 619r Z − 2624r 2 Z þ 5540r 3 Z − 5552r 4 Z þ 1648r 5 Z þ 4ð1 − 2r Z Þ 2 r p T ð−5 þ 6r Z þ 32r 2 Z Þ; ðA33Þ b ð7;0Þ þ 1613260615680 lnð2Þr 7 Z þ 1538421684224r 8 Z þ 694310400π 2 r 8 Z − 902101401600 lnð2Þr 8 Z − 445393100800r 9 Z þ 435456000π 2 r 9 Z þ 220520939520 lnð2Þr 9 Z þ 25144197120r 10 Z − 2ð1 − 2r Z Þð1 − 4r Z Þð1 − r Z Þ 2 r p T ½8486296 − 785295π 2 − 11199720 lnð2Þ − 222273632r Z þ 17777340π 2 r Z þ 267660960 lnð2Þr Z þ 2469324096r 2 Z − 173313000π 2 r 2 Z − 2810969280 lnð2Þr 2 Z − 14926447104r 3 Z þ 941371200π 2 r 3 Z þ 16619420160 lnð2Þr 3 Z þ 52801418112r 4 Z − 3062631600π 2 r 4 Z − 58795390080 lnð2Þr 4 Z − 108303542784r 5 Z þ 5922262080π 2 r 5 Z þ 122121377280 lnð2Þr 5 Z