Two-loop mixed QCD-EW corrections to gg → Hg

We compute the two-loop mixed QCD-Electroweak (QCD-EW) corrections to the production of a Higgs boson and a gluon in gluon fusion through a loop of light quarks. The relevant four-point functions with internal massive propagators are expressed as multiple polylogarithms with algebraic arguments. We perform the calculation by integration over Feynman parameters and, independently, by the method of differential equations. We compute the two independent helicity amplitudes for the process and we find that they are both finite. Moreover, we observe a weight drop when all gluons have the same helicity. We also provide a simplified expression for the all-plus helicity amplitude, which is optimised for fast and reliable numerical evaluation in the physical region.


Introduction
The discovery of the Higgs boson at the LHC [1,2] has marked a turning point in the exploration of the Standard Model of particle physics.Not only is the Higgs boson the only elementary scalar particle in the Standard Model, but it is also related to the Electro-Weak (EW) symmetry breaking mechanism, which is believed to be responsible for the observed values of the masses of all elementary particles.For this reason, the discovery of the Higgs boson and the measurement of its properties allow us to investigate the least studied aspects of the Standard Model.
Theory has to support this program by providing precise predictions for the Higgs production cross sections.A special role here is played by the process gg → H + X, which represents by far the largest Higgs production channel at the LHC.The main contribution to this channel is provided by those Feynman diagrams where the Higgs boson couples to the gluons through a top quark loop.Given its importance, this process started receiving attention already many decades ago, and today it is known up to next-to-leading order (NLO) in QCD [3][4][5][6].While in those papers it was shown that the NLO corrections can be as large as O(100%), an NNLO calculation with full top-mass dependence remains prohibitively complicated still today, primarily due to the complexity of the relevant threeloop massive scattering amplitudes.We note that recently the first numerical results for the relevant three-loop contributions have been obtained in [7].
A surprisingly reliable way to estimate higher-order QCD corrections to the gg → H + X cross-section is provided by studying this process in the limit of infinite top quark mass, where the interaction between gluons and the Higgs boson is shrunk to a pointlike effective vertex.Calculations in this limit are substantially simpler than in the full theory, which made it possible to push the perturbative expansion to NNLO [8][9][10] and more recently up to N 3 LO [11,12] in perturbative QCD.The N 3 LO corrections amount to around ∼ 5% of the total cross section and show a very good convergence of the QCD perturbative series, reducing the scale-uncertainty to ∼ 2% [13].
At this level of precision, other contributions to Higgs production cannot be neglected anymore.One such contribution is given by the class of two-loop diagrams where the gluons couple to a loop of massless quarks, followed by two massive electroweak vector bosons, which finally fuse into a Higgs boson.These EW corrections have been computed at LO and have been shown to contribute up to ∼ 5% to the gluon-fusion cross section [14,15].Given that NLO QCD corrections to gluon induced processes are typically large, it becomes very important to have a reliable estimate of the QCD corrections to this class of diagrams.Unfortunately, the calculation of these mixed NLO QCD-EW corrections is highly nontrivial, as it involves virtual three-loop three-point diagrams and real-emission two-loop four-point diagrams with massive internal propagators.While the former have been recently computed with full dependence on the Higgs and on the vector-boson masses [16], the computation of the latter has remained an outstanding challenge and, before this paper, only the relevant planar master integrals were known analytically [17].
To overcome the complexity of the full calculation, different approximations have been employed to estimate the impact of these corrections.In particular, the mixed QCD-EW corrections have first been computed in the unphysical limit m V m H [18], where they effectively reduce to a Wilson coefficient for the operator O = HG µν G µν and one therefore expects a K-factor similar to the one in the NLO QCD heavy-top approximation.This a priori unphysical approximation has recently been improved in [19], where the exact results for the virtual amplitudes computed in Ref. [16] have been combined with the real radiation computed in the soft-gluon approximation.The soft-gluon approximation is known to work relatively well for Higgs boson production [20][21][22] and the calculation showed that accounting for finite vector boson masses in the virtual corrections provides consistent results with the Wilson-coefficient approximation employed in [18].
The soft-gluon approximation amounts to the factorisation of the QCD and EW corrections in the real corrections.One could therefore wonder if a breaking of this factorisation in the real-radiation pattern could modify the K-factor in a non-trivial way.To estimate how good this approximation is, the mixed QCD-EW corrections have also been considered in the limit m V → 0 [23].This study confirmed that for small vector boson masses, the non-factorisable QCD-EW corrections remain negligible.Clearly, this does not exclude the possibility that keeping full dependence on the masses of the electroweak vector bosons could induce non-negligible modifications to the NLO corrections.It remains therefore very desirable to compute exactly the missing two-loop QCD-EW real amplitudes in order to provide a definite answer to this question.As hinted to above, this calculation is also interesting on a formal level, in particular due to the large number of scales and to the vector boson masses in the internal propagators, which translate into an involved analytic structure of the corresponding Feynman integrals.Specifically, we find that the relevant Feynman integrals can be expressed in terms of multiple polylogarithms [24][25][26][27] with algebraic arguments, involving multiple square roots.While the standard approach to compute such integrals would go through the derivation and solution of differential equations in canonical form [28][29][30][31], the complexity of the alphabet makes this strategy extremely cumbersome in practice.Interestingly, though, we find that all relevant integrals can be computed by integrating over Feynman parameters using the algorithms described in [32,33].The results thus obtained turn out to be very compact, but not extremely efficient for the numerical evaluation of the amplitude in Minkowski kinematics.This provides us with the ground to discuss a general strategy for their simplification and to present alternative results for the amplitude which are of more direct use for phase-space integration.
Finally, we stress that in this paper we only consider the two-loop real scattering amplitudes for the NLO QCD-EW corrections to gg → Hg.While we do not expect them to constitute any additional complexity, we do not consider quark-initiated partonic channels, whose contribution has been shown to be negligible at this precision [34].
The rest of the paper is organised as follows.In section 2 we give our notation and describe how to compute the helicity amplitudes for gg → Hg by decomposing the amplitude into form factors with the help of d-dimensional projection operators.After describing the reduction to master integrals and our choice of basis, we explain in sections 3 and 4 the calculation of the master integrals with two independent approaches, i.e. using differential equations and by parametric integration, respectively.We discuss our final result for the helicity amplitudes in section 5 and finally conclude in section 6.

The scattering amplitudes
We are interested in computing the two-loop mixed QCD-EW corrections to the production of a Higgs boson and a gluon in gluon fusion at the LHC.We begin by considering the process in the decay kinematics where the Higgs couples to the gluons through a pair of massive vector bosons V = Z,W ± and a massless quark loop, see Figure 1.The scattering amplitude for this process depends on the three Mandelstam variables and on the mass of the vector boson that mediates the interaction with the Higgs and which we will generically denote as m V .We use m h to indicate the Higgs mass.Since the QCD-EW contributions to gg → Hg start at two-loop order, the amplitudes computed in this paper are finite, as long as all external gluons are fully resolved.
In order to perform the computation, we begin by decomposing the scattering amplitude for H → ggg into Lorentz-and gauge-invariant tensor structures.We extract the dependence on the SU (3) color structure f abc and write where, for each j = 1, 2, 3, j is the polarisation vector of the gluon of momentum p j , while a j is its color label.A µνρ (s, t, u, m 2 V ) must be a rank-3 tensor under Lorentz transformations and, imposing gauge invariance for each of the external gluons, it can be written as a linear combination of four independent form factors.Following the conventions introduced in [35], we require that the gluons are transverse and make a cyclic choice for their gauge fixing condition With this choice, one easily finds [35] where the F j (s, t, u, m 2 V ) are Lorentz-invariant form factors.We stress that in this decomposition no parity-violating terms appear.This can be justified by noticing that at LO, both for W and Z exchange, the two vector bosons are always attached to the same massless fermion loop and all contributions which have an odd number of γ 5 drop when summing over the two orientations of the loop momentum.Clearly, since we only consider massless quarks, we omit the contribution of the bottom quark in the diagrams where W bosons are exchanged.With this, we can write for each form factor where and O(α s , α) indicates higher order contributions both in the QCD and in the EW coupling.The full QCD-EW corrections can then be obtained by summing the contributions with V = Z or W .
The form factors F j , or equivalently the F (0) j , are not the objects that we are ultimately interested in.Indeed, often substantial simplifications occur when one combines the form factors to compute so-called helicity amplitudes.For the case at hand, each gluon can have two different helicities for a total of eight different combinations.By use of Bose symmetry, parity and charge conjugation, one can easily show that only two of them are independent.We indicate the helicity of the gluon of momentum p j by λ j and write for a generic helicity amplitude and for a given vector boson V where A µνρ (s, t, u, m 2 V ) was defined in eq.(2.3).We proceed by choosing as two independent helicity amplitudes A ++± (s, t, u, m 2 V ).It is straightforward to find compact expressions for these amplitudes in terms of the form factors in eq.(2.5) using the spinor-helicity formalism, see [36] and references therein.We choose for the polarisation vectors of the external gluons where q j is an arbitrary reference vector with q 2 j = 0 and q j • p j = 0.While in principle the vector q j can be chosen freely, the conditions in eq.(2.4) force us to pick q 1 = p 2 , q 2 = p 3 and q 3 = p 1 .With this, the two independent helicity amplitudes become where the Ω ++± are linear combinations of the original form factors Similarly to eq. (2.6), we can explicitly extract the LO EW and QCD couplings from the amplitudes and write for the perturbative expansion of the helicity coefficients such that, again, the full QCD-EW contributions are obtained by summing the corresponding helicity amplitudes with V = Z, W .

The evaluation of the helicity amplitudes
The helicity amplitudes in eq.(2.11) receive contribution from 21 different two-loop Feynman diagrams, see Figure 1 for two representative ones.The contribution of each of these diagrams to the helicity coefficients can be computed by deriving d-dimensional projector operators.The standard approach consists in introducing 4 independent projectors which single out the contribution to each of the form factors defined in eq.(2.5) where, for consistency with eq.(2.4), we must use pol pol We stress at this point that all Lorentz indices in eq.(2.13) have to be understood as d-dimensional.Each projector can be decomposed in terms of the same tensor structures as in eq.(2.5) as follows where j ∈ {1, 2, 3, 4}.By imposing that eq.(2.13) is satisfied we find c d s t u . (2.18) We can either use these projectors to evaluate the four form factors independently, or we can use them, together with the definition of the helicity coefficients in terms of form factors in eq.(2.11), in order to derive new helicity-projectors [37] that directly project on the physical helicity amplitudes such that pol Since the helicity amplitudes are the physical objects that we will be ultimately interested in, we prefer to follow this second approach.
Denominator integral family PL Table 1: Definition of the planar (PL) and non-planar (NP) integral families.The loop momenta are denoted by k and l, while m V indicates the mass of the vector boson.The prescription +i is understood for each propagator and not written explicitly.
In practice, we generate all relevant two-loop diagrams using QGRAF [38] and we use FORM [39] to apply the projectors in eq.(2.19) and write them as linear combinations of scalar two-loop Feynman integrals.We find that all diagrams can be mapped on Feynman integrals of two integral families, one planar (PL) and one non-planar (NP), up to crossings of the external legs.We define these two families as follows: where top ∈ {PL, NP} labels the families and the denominators D 1 , . . ., D 9 are given in table 1.We use dimensional regularization with d = 4 − 2ε, and our convention for the integration measure for each loop is With the definitions given in Table 1, the two diagrams in figure 1 can be described using the first 7 propagators of the two families respectively, and all other diagrams which contribute to the process can be obtained by permutations of the external gluons and by pinching of the propagators. 1Although quite standard, the reduction to a subset of master integrals [40][41][42], is non-trivial due both to the large number of scales and the presence of massive internal propagators.We used Reduze2 [43] to map the diagrams to the relevant integral families and performed a complete reduction of all integrals with KIRA [44]. 2inally, we found it convenient to use FiniteFlow [46] to speed up the substitution of the reduction identities produced by KIRA in the helicity amplitudes of eq.(2.11) and their simplification.
We find that the two independent helicity amplitudes can be expressed in terms of 116 master integrals, counting also the ones obtained through permutations of the external gluons as independent ones.If we limit ourselves to the un-permuted integrals, we find 43 planar and 18 non-planar master integrals, see appendix A for the full list.To construct our initial basis of master integrals, we select integrals whose maximal cuts are defined by integrands with unit leading singularities [31,47,48].Our choice avoids the appearance of irreducible denominator factors that mix the kinematical variables and the dimensional regularization parameter d during the IBP reduction.This reduces the complexity of intermediate expressions, similarly as described in [35], and recently automated in [49,50].
In the next two sections we will describe two different strategies that we used to compute the master integrals in terms of multiple polylogarithms.

Computation of the master integrals with differential equations
The standard approach to compute a complete set of multiloop, multiscale Feynman integrals goes through deriving and solving their system of differential equations with respect to the masses and momenta, as first worked out in full generality in [30].In each of the invariants ξ = (s, t, . . ., m 2 , . . .), a basis of master integrals always fulfils a linear system of differential equations with rational coefficients.By indicating with I the vector of master integrals, we can write this system as where the entries of the matrix B(ε, ξ) are differential one-forms that are rational in the kinematics and in the dimensional regulator ε.One then usually tries to solve these equations as a Laurent series in ε, i.e. for d → 4. The effectiveness of this approach relies on the ability to find a solution of the homogeneous part of the system above, in the limit ε → 0. While this would be a daunting task given a generic system of coupled differential equations, it turns out that an integral representation for the homogeneous solution can always be obtained by analysing the maximal cuts of the corresponding Feynman integrals [51][52][53], whose computation becomes particularly simple using the so-called Baikov representation [54][55][56].
While the approach described above is completely general, it was shown that in many cases the solution of the differential equations can be greatly simplified by the choice of a so-called canonical basis of master integrals [31].If such a basis F can be found, the corresponding system of master integrals becomes where the new matrix A(ξ) does not depend on ε.In addition to the factorisation of ε, an important condition for the basis to be canonical is that the matrix takes a particularly simple, "d log" form where A j are matrices of rational numbers and P j are algebraic functions of ξ, which constitute the alphabet {P 1 , . . ., P J } of the problem.It follows from eq. (3.2), that the master integrals of a canonical d log basis can be expressed, order by order in ε, as iterated integrals of the forms d log(P j ).Furthermore, whenever the alphabet consists entirely of rational functions P j (or if this can be achieved by an algebraic change of variables), then these iterated integrals can be expressed as linear combinations of the functions where the arguments σ i and x will be certain algebraic functions of ξ.The iterated integrals (3.3) are known as multiple polylogarithms [27] and hyperlogarithms [57] of weight k. 3 For most of the Feynman integrals that have been studied so far, finding a canonical d log basis comes along with an expression for the corresponding master integrals in terms of multiple polylogarithms (with potentially complicated algebraic arguments).However, no general method to construct an expression of this kind is known if the alphabet cannot be rationalized 4 such that, in some cases where a canonical form for the differential equations is known, the issue of the existence of a polylogarithmic expression for the master integrals remains matter of discussion, see for example [60]. 5In fact, more recently it was shown that there exist iterated integrals of d log forms which cannot be expressed in terms of multiple polylogarithms [62].In conclusion, whether or not Feynman integrals with a canonical d log form can be expressed through multiple polylogarithms, remains an intricate problem.With these general comments in mind, let us consider now the form of the system of differential equations for the problem at hand.First of all, it is interesting to notice that, in order to evaluate all the master integrals required for the amplitude, we need to introduce two additional master integrals that would otherwise not appear in our problem, namely I PL (2, 2, 0, 0, 0, 1, 0, 0, 0) and I PL (2, 2, 0, 0, 0, 1, 1, 0, 0).These additional master integrals appear in the non-homogeneous part of the differential equations for some of the top-sector master integrals, and it can immediately be seen that they are obtained by pinching some of the internal lines of the diagrams in Fig. 1 (see Appendix A for a complete list of the master integrals).All master integrals are functions of four independent variables, which in this section we choose to be t, u, m 2 h , and m 2 V (V = W, Z).Since Feynman integrals are homogeneous functions in the kinematic invariants and in the masses, it is convenient to introduce the dimensionless variables in order to factorise the dependence of each master integral on m 2 h as a simple power, namely (m 2 h ) d−a 1 −...−a 9 , where the a i are the powers of the propagators, see eq. (2.21).For the remainder of this section, we can hence set m 2 h to 1.To determine the expressions of the master integrals in terms of the remaining variables we derive differential equations in y, z, and ρ for them and cast this system into a canonical form, as in eq.(3.1).This was achieved by starting from a basis of master integrals whose maximal cuts have unit leading singularities (see appendix A), and then applying the algorithm described in [63].
If we limit ourselves to the 48 planar master integrals, see eq. (A.1), then the differential equations take a very simple form and, in particular, all letters are rational functions of y, z, ρ and a single square root, As it is well known, this root can be rationalized by the change of variables and all integrals of the family PL can be expressed in terms of multiple polylogarithms whose arguments are rational functions of x, y, z, as it was shown explicitly in [17], which we refer to for the explicit form of the differential equations and of the canonical basis.Unfortunately, even if there is only a small number of new, non-planar integrals, their differential equations turn out to be substantially more complicated.In this case, the vector of planar and non-planar master integrals F contains 64 entries, and the alphabet (3.7) depends on three additional square roots, defined as We provide both the vector of canonical functions F and the d log forms of eq.(3.2) in the ancillary files of this paper.
As described above, once a canonical form for the differential equations is obtained, the standard procedure consists of constructing a solution as a Dyson series in ε whose coefficients consist of iterated integrals.In case of a three scale problem, the usual strategy consists in solving the partial differential equations sequentially, as described for example in Refs.[30,64].We start with one variable and solve the corresponding linear differential equation up to a function of the other two variables, then we write down a differential equation with respect to a second variable.We check that the right-hand side of this new equation is independent of the first variable and we solve this equation in terms of multiple polylogarithms up to a function of the last variable.After the last differential equation is integrated, a solution is obtained up to constants which are then fixed by choosing appropriately the boundary conditions.While this strategy can be applied straightforwardly when the alphabet is linear in all variables, the presence of several algebraically independent square roots makes it frequently unfeasible in practice, since at a given step, it is not in general possible to find a representation for the result where the corresponding integration variable appears only in the last argument of the various polylogarithms.
Despite this, it turns out that in our problem the four square roots appear in the differential equations in a very structured pattern, which allows us to devise a solution strategy that is always guaranteed to terminate and to produce a result in terms of multiple polylogarithms.First of all, we find it convenient to rationalize the root R 0 , which appears consistently throughout the whole system of equations, by the change of variables of eq.(3.6).We then notice the following crucial structural features of the differential equations: • R 1 appears only in the differential equations for the canonical functions F • all the other equations contain at most the root R 0 ; • when solving the equations for F 61 , F 62 , F 63 , F 64 , at most two square roots are integrated at once and only from weight 3 on: either This separation of square roots allows us to perform different changes of variables depending on the canonical functions we want to evaluate, in particular depending on which roots enter a particular integration.First of all, as customary when dealing with canonical master integrals, we normalise our basis such that all integrals start at order ε 0 with a weight 0 constant (which could of course be zero).We start by solving the equations for the canonical functions F i , i ∈ {1, . . ., 50, 52, 53, 54}, where no square roots appear in (x, y, z).We integrate first in y, then in z and at last in x, obtaining expressions of uniform weight written in terms of multiple polylogarithms up to some constant factors which will be fixed by imposing boundary conditions.We move then to the canonical function F 51 , where also R 1 appears.The relevant letters can be rationalized through the change of variable and then solved first in u and subsequently in x.Only two variables appear here, because F 51 is a three-point function.
The group of canonical functions F {61,62,63,64} , corresponding to the master integrals of the top non-planar sector, is the most difficult one.As observed above, up to order ε 2 no square roots are present in the variables (x, y, z), therefore the integration in terms of multiple polylogarithms is straightforward, and is carried out following the procedure used for F {1,...,50,52,53,54} .
Starting from order ε 3 all square roots R {1,2,3} appear, but always in such a way that a single nested integration contains at most two of them, specifically either R 1 and R 2 , or R 1 and R 3 .We start by removing R 1 via the change of variables of eq.(3.9).This change of variables is sufficient to take care of the nested integrations coming from the homogeneous part of the differential equations, as well as of the one coming from the non-homogeneous terms related to F i , i ∈ {1, . . ., 54}.Despite the fact that only rational functions are now present for this subset of terms (allowing us to represent this part of the solution in terms of multiple polylogarithms), many cumbersome letters arise, as for example Considering now the terms related to F {55,56,57} , a second change of variables to rationalize also R 2 is performed "on the fly", and reads (3.12) An analogous "on the fly" change of variables is performed on the terms related to F {58,59,60} , to get rid of R 3 : Implementing such changes of variables allows us to write the ε 3 coefficients of F {61,62,63,64} again in terms of multiple polylogarithms, at the price of having three different sets of independent variables: (x, u, z), (x, u, v), and (x, u, w).
To integrate one of the above subsystems, we integrate first in z, or v, or w, according to preferred variables just discussed.After the integration in z, we verify that plugging these solutions in one of the remaining differential equations gives a matrix of coefficients which is independent of z, where this condition must be satisfied considering also the hidden dependence through v and w.The expressions that arise are so cumbersome that we do not see any chance to perform this check analytically.On the other hand, we see numerically (we use GiNaC, [65,66] to evaluate multiple polylogarithms) with very high accuracy that our expressions are independent of z.However, we cannot simply substitute z = 0 (which corresponds to v = 0 and w = 0), because individual terms may be singular in this limit.To address this issue, we use shuffle relations to extract carefully all such singular terms as z → 0 explicitly as powers of log(z).We confirm numerically that the sum of the three contributions in terms of different variables as well as its limit at z → 0 is independent of z.Once this is done, we proceed by integrating in u, and we check that the remaining differential equation is independent of u, in a similar way as we did for z.After having solved all differential equations up to integration constants, we fix these constants using boundary conditions in the large-mass limit, x → 0.Here we use well-known prescriptions in a graph-theoretical language for limits typical of Euclidean space -see, e.g., [67].
The procedure described above can be applied to obtain the ε 4 part as well, but the manipulations required are extremely cumbersome and, a posteriori, not needed.Indeed, in the next section we will show how to obtain these integrals in a much simpler way by direct integration over their Feynman/Schwinger parametrisation.In any case, we believe that the approach we used to solve the differential equations presented here can be used also in other situations where many different square roots appear and only subsets of them are rationalizable at once.The key point of this procedure is to check that in each nested integration only one subset of simultaneously rationalizable square roots is present, and then to perform a "local" change of variables "on the fly" to rationalize them.

Computation of the master integrals by parametric integration
While an expression for the master integrals in terms of multiple polylogarithms can in principle be obtained from the differential equations, the procedure was rather cumbersome as explained in the previous section.An entirely different approach, which one might attempt, consists in computing all integrals starting from their Feynman/Schwinger parametrisation.This can be in general quite difficult, in particular in multiloop/multileg problems, where one typically needs to integrate over a large number of Feynman parameters.Nevertheless, it turns out that in our problem all integrands are linearly reducible [33,68], which means that the algorithms described in [32] can be applied rather directly.
In order to make this approach feasible, it is helpful to choose a basis of master integrals that is finite in the limit d → 4. We expect such a change of basis to be particular useful in the case at hand since the two-loop amplitude is, effectively, a leading-order amplitude and therefore expected to be finite.In practice, however, we found it sufficient to replace only the most divergent master integrals with more than 5 propagators by finite counterparts.To achieve this, we generated finite integrals by considering the corresponding six-dimensional integrals, including higher powers of the propagators.For each of the integrals in the families in table 1 we can obtain the corresponding (d + 2)-dimensional integral by [69,70] where G(p 1 , . . ., p n ) is the Gram-determinant of the n momenta p 1 , . . ., p n , and It is pretty easy to see that, at least in the case at hand, as long as we choose UV finite integrals and all powers of the massless propagators equal to unity, the Gram determinant G(k, l, p 1 , p 2 , p 3 ) cures all IR divergences, both in the collinear and in the soft limits.This allows us to easily generate a large number of finite integrals.We stress that this is particularly straightforward here due to the presence of two internal massive propagators.
In fact, even for integrals with fewer propagators (and therefore with poor UV behaviour), we can simply raise the powers of the massive propagators ad libitum in order to obtain UV-finite integrals, without spoiling their IR behaviour.We note that, in a general case, finite integrals can be found algorithmically also in the absence of massive propagators, see for example the algorithm described in [71].We list the finite integrals used in this calculation in Appendix B.

Planar integrals
The parametric representation [72,73] of an integral family such as (2.21) has the form where ω = a 1 + • • • + a 9 − d.The polynomials U = det A and F = U(B A −1 B − C) are determined by the quadratic (A), linear (B) and constant (C) parts of the quadratic form in the two loop momenta = k l , given by the denominators from Table 1.All integrals that we are interested in for this calculation, for both integral families, are chosen such that a 8 = a 9 = 0, which allows us to eliminate the parameters x 8 and x 9 .The remaining denominators D 1 , . . ., D 7 are the inverse scalar propagators of the graphs shown in Figure 1.Concretely, in the planar case we find the Symanzik polynomials to be An analysis by polynomial reduction [68] shows that the set {U, F} is linearly reducible.This means that the integrals (4.2) can be expressed algorithmically in terms of the hyperlogarithms defined in eq.(3.3).In fact, this works to all orders of the ε expansion, and for arbitrary integer values of a 1 , . . ., a 7 .The algorithm described in [33] applies directly only to convergent integrals.As explained above, we therefore adjusted our basis to consist mostly of finite integrals.The remaining divergences in this basis (see Appendix B) occur only in integrals with 4 or fewer propagators, and six further integrals with 5 or 6 propagators, where they can be resolved easily through integration by parts in the parameters x k , following the method of [71,93].In order to perform the polynomial reduction, resolution of divergences, and integration over the Feynman parameters explicitly, we used the code HyperInt [32].Starting from the polynomials (4.3), HyperInt identifies x 1 , x 4 , x 5 , x 6 , x 3 as an admissible order for the first five integrations.They result in expressions with hyperlogarithms whose arguments σ i are rational functions of s, t, u, m 2 V and x 2 , x 7 .We pick j = 7 for the constraint x 7 = 1 in (2.21), leaving the final integral over x 2 .At this stage, the algorithm needs to solve for the roots of the polynomial which we saw also in (3.5).Consequently, the final expressions for the coefficients of the ε-expansion of the integrals I PL are linear combinations of hyperlogarithms, whose coefficients and arguments are rational functions of s, t, u, m 2 V and r.

Non-planar integrals
For the non-planar integral family, the corresponding polynomials are and it was pointed out in [93,Figure 10] that they are linearly reducible too.As an admissible order for the first integrations we use x 1 , x 3 , x 5 , x 2 , x 7 .Setting x 6 = 1, the final integration over x 2 introduces three further square roots in addition to r: which we encountered also in the differential equations, see (3.8).Our results for the integrals I NP from the basis (B.1) therefore consist of linear combinations of hyperlogarithms with coefficients and arguments that are rational functions of s, t, u, m 2 V , r and the three roots in (4.6).In fact, the polynomial reduction shows that the quadratic polynomials responsible for R 2 and R 3 are not compatible [68] with each other.Explicitly, this manifests itself in the fact that our results admit a decomposition into expressions A and B whose hyperlogarithm arguments σ k in (3.3) are rational functions of the listed arguments only, i.e. the roots R 2 and R 3 do not mix.This property corresponds to the structure of the differential equations described in section 3, and makes it possible to rationalize the pieces A and B individually.For the parametric integration, however, such rationalizations provide no advantage.In contrast, the bare expressions with the (unrationalized) roots are much more compact.
Remark.We stress that the hyperlogarithm expressions obtained from HyperInt are valid for all values of the kinematic parameters such that the integral (4.2) converges.In particular, by giving a small positive imaginary part to s, t and u in order to implement the i prescription, these hyperlogarithms can be evaluated directly in the physical region, for example using GiNaC [66].This is a very valuable property, because the analytic continuation of polylogarithms with algebraic arguments is typically much more delicate.
For all ε-expansion coefficients of the integrals (B.1) that contribute to the helicity amplitudes (2.11), we find that only hyperlogarithms of weight k ≤ 4 arise.This weight bound is consistent with other known two-loop amplitudes in four dimensions.In ancillary files to this publication, we provide the explicit expressions thus obtained for all coefficients of the ε-expansions of the integrals in our basis (B.1) that are required for the computation of the helicity amplitudes.The ancillary files also include instructions and code to reproduce these calculations.

The helicity amplitudes
Combining our results for the Feynman integrals, we obtain expressions for the helicity amplitudes Ω (0) ++± .At this step, we see that all poles in ε stemming from individual divergent integrals, as well as from the coefficients in the reduction of the amplitudes to the Feynman integrals, completely cancel each other.As expected, the helicity amplitudes thus turn out to be finite.Furthermore, we notice that: +++ , all hyperlogarithms of weight 4 cancel out, leaving only functions of weight at most 3 in the result.A similar weight drop was found in mixed QCD-EW corrections to gg → H, see [16,74], where the two-and three-loop amplitudes turn out to have maximum weight three and five, respectively.
• In the case of Ω (0) ++− , hyperlogarithms of weight 4 do not cancel completely and persist in the result.
These weights may at first seem surprising, in particular because no such weight drop shows up in the corresponding HEFT amplitudes, see for example [75].But for our mixed QCD-EW corrections, the weight drop can be explained, rather loosely, as follows.If we consider the possible unitarity cuts of the Ω (0) +++ helicity amplitude, we find that angular momentum conservation in the internal quark lines implies that all cuts should be zero in ε = 0, causing the observed weight drop.This argument applies equally well to gg → H, where the only helicity amplitudes different from zero are for equal-helicity gluons.On the other hand, this argument clearly does not apply to Ω (0) ++− and no weight drop can be expected on general grounds.
After some simplification, our result for the helicity amplitude Ω (0) +++ takes the form where the hyperlogarithms H = H 1 + H 2 + H 3 are given in weight 1 and 2 explicitly as The expression for H 3 , the hyperlogarithms of weight 3, is provided in the ancillary files.Their arguments are rational functions of s, t, u, m 2 V and the roots The amplitude Ω ++− is more complicated, not only because it involves hyperlogarithms of weight 4, but also since their arguments require two additional square-roots, These roots arise from R 1 in the crossed versions (s ↔ t or s ↔ u) of the integrals that we computed in section 3 and section 4. Similarly, crossing is responsible for the appearance of the root r s in (5.3).The explicit form of Ω ++− is provided in the ancillary files.The results so obtained can be evaluated rather straightforwardly in any region of phase-space, in particular both in the Euclidean, s, t, u < 0, and in the physical 6 Minkowskian region where Indeed, the hyperlogarithms can be evaluated numerically with GiNaC [66], provided a small imaginary part is given to s, t and u.This is needed, also in the Euclidean region, because individual hyperlogarithms in the expression are not necessarily single-valued, and a consistent determination for all of them must be picked.In the Euclidean region, all choices for the signs of the infinitesimal imaginary parts produce the same, real, result.The correct result in the physical region, however, is obtained by ensuring that both s and m 2 h = s + t + u have a positive imaginary part (according to the i prescription).
point in phase-space For reference, we provide numerical results for the helicity amplitudes in two points in the Euclidean region and two in the physical region.We pick two points in the physical region (5.5) with m 2 h = (125/90) 2 m 2 V such that and two points in the Euclidean region with For these points, the numerical values of the helicity amplitudes are provided in table 2.
It is interesting to notice that, in the bulk of the phase-space, the two helicity amplitudes are numerically similar.This is in part due to the fact that the amplitudes are expected to go to the same value both in the limit m V → ∞ and when the gluon p 3 becomes soft, see section 5.2 for details.

Polylogarithm expressions for Ω (0) +++
While the amplitudes in the form discussed above are guaranteed to produce the correct result, if the Feynman i prescription is applied, the numeric evaluation of the hyperlogarithms is not particularly efficient, especially in the physical region.In order to obtain a fast and stable method to evaluate the helicity amplitudes, we rewrite the hyperlogarithms in terms of simpler functions.In particular, classical polylogarithms [76] of weight k, are readily available for speedy evaluation in many computer algebra systems.It was demonstrated in [77] that every hyperlogarithm of weight 3 can be expressed as a linear combination of Li 3 's with suitable arguments, plus products of Li 2 's and logarithms.However, in deriving such an expression for Ω (0) +++ , great care is required due to the multivaluedness of polylogarithms.The principal branches have discontinuities on the rays (−∞, 0] for log, and [1, ∞) for Li k . (5.9) An expression built out of principal branches of polylogarithms typically develops discontinuities whenever an argument crosses one of these branch-cuts.It is therefore not always possible to find a single expression that captures the desired branches over the entire phasespace.Instead, different expressions must be derived in various sub-regions of phase-space.
In the ancillary files, we therefore provide two different expressions for Ω (0) +++ written in terms of Li 3 , Li 2 and logarithms only: • one expression is valid in the entire Euclidean region defined by s, t, u < 0 < m 2 V .
• one expression is valid in the entire physical region defined in eq.(5.5).
Note that due to the symmetry of Ω +++ under permutations of s, t and u, the latter region completely determines this helicity amplitude in the entire physical region of interest.
In the Euclidean case, the four roots (5.3) are positive and real.The arguments of the polylogarithms Li 2 and Li 3 in our expression are chosen to be real and less than 1, over the entire Euclidean region.Hence, the resulting expression is manifestly real in the entire Euclidean region and efficient to evaluate.
After analytic continuation, in the physical region the roots (5.3) take the values and we ensured that the arguments of all polylogarithms in our corresponding expression stay away from the branch cuts (5.9), throughout the entire region (5.5).Our second polylogarithm expression for Ω +++ , tailored for the physical region and given in the ancillary files, can thus be evaluated in that region efficiently and robustly, without any ambiguities.

Remark.
A priori, it is not guaranteed that such an expression, single-valued throughout the entire physical region, even exists at all.Further subdivisions of phase-space might have been required, see for example [61,63,78].
In order to derive the expressions for Ω (0) +++ discussed above, we followed roughly the approach outlined in [79].First, we computed the symbol of the amplitude, which we find to produce 39 letters, namely The second condition selects different arguments for the Euclidean and physical regions, leading to different final expressions.To check for the factorizations in the first condition, we used integer relation techniques as detailed in [81, section 3].
For the other helicity amplitude Ω (0) ++− , the result includes hyperlogarithms of weight up to and including four, and the corresponding symbol alphabet is more involved due to the presence of the two extra square roots in (5.4).In a similar way as above, it would be possible to rewrite our expressions in terms of simpler polylogarithms, reducing the set of transcendental functions to log, Li 2 , Li 3 , Li 4 and Li 2,2 , as explained for example in [82].We leave this to future work.

Checks on the result
Each master integral, with exception of the weight four piece of the 7-propagator non-planar integrals, has been successfully checked using the Mathematica [83] package PolyLogTools [65,66,84] to numerically compare its expression obtained via differential equations to its expression calculated through integration over Feynman parameters in multiple points inside the Euclidean region.Furthermore, the results from Feynman parameters integration (including weight four for the 7-propagator integrals) have been checked numerically against PySecDec [85][86][87][88][89][90][91] both in the Euclidean and in the Minkowski region, finding excellent agreement in all points.Finally, also the results from the differential equations have been checked in random points in the Euclidean region against FIESTA [92], finding excellent numerical agreement.
To validate our results for the amplitude we considered two different limits for the amplitude: the soft-gluon limit and the limit of a vector boson with infinite mass.
In the soft limit, the gg → Hg amplitude A c 1 c 2 c 3 λ 1 λ 2 λ 3 factorizes into the leading order gg → H amplitude A λ 1 λ 2 times an eikonal factor. 7Using the gauge choice of eq.(2.4), the factorization takes the form which can be rewritten in terms of spinor products as A ++ . (5.13) Using the same normalisation for the EW and QCD couplings, the leading order amplitude for gg → H for gluons of plus helicity can be written schematically as [16] A 7 The color structure of the leading order amplitude has been included in the eikonal factor.
where F is a non-trivial function of the ratio m 2 h /m 2 V .Inserting the expression above into the soft limit we get which correspond to our expressions for the amplitude in eq.(2.10).Indeed, we could check numerically that for t → 0 − , u → 0 − , s → m 2 h , we obtain lim To check the m V m h limit we start by recalling that, in this approximation, the interaction can be encapsulated in a Wilson coefficient for the effective Lagrangian [18,75] where v denotes the vacuum expectation value of the Higgs field.Up to the explicit form of the Wilson coefficient C 1 , this Lagrangian is identical to the heavy-top mass Lagrangian.We can therefore read off the leading order mixed QCD-EW gg → Hg amplitude directly from the corresponding computation in the heavy-top limit, which is presented in [75] as with C 1,EW = −α C W + C Z cos 2 θ W / 16π 2 sin 2 θ W .We then equate this result to our expressions for A c 1 c 2 c 3 eff,+++ and A c 1 c 2 c 3 eff,++− , and then take the limit m W → +∞, numerically finding agreement up to an overall phase.
To check the consistency of this overall phase, we finally perform the same limit for gg → H at leading order (from [16]), finding agreement up to the same global phase.

Conclusions
In this paper we described the first calculation of the two-loop mixed QCD-EW corrections to the production of a Higgs boson and a gluon in gluon fusion, with full dependence on the Higgs and on the vector boson masses.The amplitudes presented here are the last missing building blocks required to compute the NLO mixed QCD-EW corrections to Higgs production in gluon fusion, overcoming the shortcoming of the various approximations that have been used to estimate these corrections in the past.We made use of helicity projector operators to extract the two independent helicity amplitudes from the two-loop Feynman diagrams that contribute to the process in terms of scalar Feynman integrals.We reduced all scalar integrals to master integrals by use of integration by parts identities and computed the master integrals with two independent methods, namely both starting from their differential equations in canonical form and by direct integration over their Feynman/Schwinger parametrisation.In both cases, we find that the result can be expressed in terms of multiple polylogarithms.Achieving this form by integrating the differential equations turned out to be cumbersome in practice, in spite of the fact that a canonical form for the differential equations could be found.In fact, the alphabet of the non-planar master integrals is characterised by the presence of four independent square roots, that we did not manage to rationalize at the same time.For this reason, integrating the equations required us to split the master integrals into different contributions, and to use different changes of variables to rationalize the square roots in each of these pieces.This was doable in practice thanks to the particular structure of the system of differential equations, but it produced rather cumbersome results.

(B.2)
We note that the six integrals in the top two rows of (B.2) have only a single pole as d → 4 and they appear in the amplitude with a factor of (d − 4), so only the pole (leading order) of those integrals contributes to the helicity amplitudes in d = 4.

Figure 1 :
Figure 1: Representative planar (a) and non-planar (b) Feynman diagrams for the LO mixed QCD-EW corrections to gg → Hg.The internal wavy lines represent the massive vector bosons.

Table 2 :
, u, m 2 Numerical values for the two helicity amplitudes in the Euclidean and in the physical region, at the points defined in eqs.(5.6) and (5.7).