Universal four-dimensional representation of $H \to \gamma \gamma$ at two loops through the Loop-Tree Duality

We extend useful properties of the $H\to\gamma\gamma$ unintegrated dual amplitudes from one- to two-loop level, using the Loop-Tree Duality formalism. In particular, we show that the universality of the functional form -- regardless of the nature of the internal particle -- still holds at this order. We also present an algorithmic way to renormalise two-loop amplitudes, by locally cancelling the ultraviolet singularities at integrand level, thus allowing a full four-dimensional numerical implementation of the method. Our results are compared with analytic expressions already available in the literature, finding a perfect numerical agreement. The success of this computation plays a crucial role for the development of a fully local four-dimensional framework to compute physical observables at Next-to-Next-to Leading order and beyond.


Introduction
The calculation of observables the physics the LHC delivers has nowadays been improved with several techniques to make predictions at the Next-to-Leading order (NLO) accuracy. In general, it is aimed at reducing the scale uncertainties from O(10 %) to O(1 %), and to even less for some specific processes such as Drell-Yan. In order to provide these observables, we rely our prediction on the calculations of scattering amplitudes through perturbation theory. For the calculation of the latter at NLO, apart from evaluating Feynman integrals, we also need to deal with the evaluation of integrals in the momentum space. At one-loop level, the basis of integrals is known and their evaluation has been implemented in several codes (for instance, in refs. [1,2]). Nevertheless, the evaluation of multi-loop integrals still remains a work in progress [3,4]. On top of it, depending on the process under consideration, there may appear infrared (IR) and ultraviolet (UV) singularities that are canceled out by adding real corrections and proper counter-terms.
The two-loop QCD corrections to the decay process H → γγ have been first evaluated in the heavytop limit [5][6][7][8] and with the full top-mass dependence [9,10]. The two-loop electroweak corrections have been investigated in refs. [11][12][13][14]. Combining the two-loop QCD and electroweak corrections, it is possible to observe a nearly complete cancellation between these two contributions for M H = 126 GeV [15]. At Next-to-Next-to-Leading order (NNLO) the non-singlet [16] and singlet QCD contributions [15] have been calculated in the heavy top quark limit.
In view of the enormous success of the Standard Model (SM) of particle physics with the detection of the Higgs boson, new directions have been taken to discuss in more details the consequences of this discovery. In particular, from the phenomenological point of view, the background of the experiment has to be removed. Hence, QCD predictions up to the Next-to-Next-to-Next-to-Leading order (N 3 LO) have been provided in an effective theory [17]. Also, it has been shown that the mixed effects of QCD-electroweak contribution to the amplitude are relevant [18].
Nevertheless, hidden mathematical properties of the amplitudes gg → H and H → γγ have been extensively studied in the full theory at Leading order (LO) in ref. [19], in which, it was showed that these amplitudes exhibited remarkable properties when computed using the Loop-Tree Duality (LTD) theorem [20][21][22]. The dual contributions we have obtained for different internal particles -charged electroweak gauge bosons, massive fermions and charged scalars -featured the very same functional forms, and could be written in a universal way using scalar parameters depending only on the space-time dimension (d), and the mass of the particles involved in the process. We also obtained a pure four-dimensional (d = 4) representation of the renormalised amplitude and recovered the well-known results found in the literature [23][24][25][26][27][28].
In this manuscript we push the computation further by considering the H → γγ process at two-loop level, and show that the above-mentioned properties are still present. To this end, we make use of the LTD theorem, which converts loop integrals into phase-space ones. In order to provide the renormalised amplitude, we perform a local UV renormalisation that leads to a finite integrand in four space-time dimensions. This algorithm is based on the refinement of the expansion around the UV propagator [29][30][31] to account for the different singular behaviours of the internal loop momenta in the UV region. Furthermore, since this amplitude is IR safe, we can directly treat the virtual integrand in four dimensions. We remark that the calculation of this amplitude is the first two-loop application to a physical process done through LTD, and it is computed below the mass threshold limit in the MS renormalisation scheme. We note that for individual diagrams, unphysical threshold singularities appear but they cancel among themselves when the full amplitude is considered.
In the same spirit of the universality that these amplitudes exhibit at LO, we consider as internal particles charged scalars and top quarks. While we only consider QED corrections, they can be straightforwardly promoted to QCD ones by replacing the couplings accordingly. We compare our results with known results [9,10] finding full agreement.
With this paper we verify that the LTD approach holds also at multi-loop level and, therefore, N k LO predictions involving virtual amplitudes can be achieved by means of the former. Additionally, we remark that the traditional approach based on the use of integration-by-parts identities [32,33] is not needed to evaluate the actual amplitude. In fact, we overcome the calculation of the latter making our procedure much lighter as we shall describe here.
The paper is organised as follows. In section 2, we recall the basics of the LTD formalism, at oneand two-loop level. We introduce there the notations used in this paper, and provide the master formula for obtaining the dual representation of a two-loop Feynman integral. In section 3, we sketch the algorithm to algebraically reduce the integrand-level expressions of two-loop dual amplitudes, and rewrite every scalar product involved in terms of denominators. We provide the tensor structure of the H → γγ amplitude in section 4, and briefly recall the one-loop result obtained in ref. [19]. Then, in section 5, we collect and write the universal coefficients involved in the universal structure of the two-loop dual expressions. We discuss in section 6 the cancellations of unphysical threshold singularities that appear among the dual contributions, and we explicitly show how they occur. In section 7, we discuss an algorithmic approach to locally renormalise two-loop amplitudes within the LTD formalism. In particular, we focus on the determination of the schemefixing parameters in the MS scheme. In section 8, we present our numerical results and show there is a complete agreement with the analytical expressions. We draw our conclusions and discuss future directions of this work in section 9.
Algebraic manipulations have been carried out by using an in-house implementation of LTD which is based on the MATHEMATICA packages FEYNARTS [34] and FEYNCALC [35,36].

The Loop-Tree Duality at two loops
The Loop-Tree Duality (LTD) theorem [20][21][22] transforms any loop integral or loop scattering amplitude into a sum of tree-level like objects that are constructed by setting on shell a number of internal loop propagators equal to the number of loops. Explicitly, LTD is realised by modifying the ı0 prescription of the Feynman propagators that remain off shell with k ji = q j − q i , and η µ an arbitrary future-like vector. The most convenient choice is η µ = (1, 0), which is equivalent to integrate out the loop energy components of the loop momenta through the Cauchy residue theorem. The left-over integration is then restricted to the Euclidean space of the loop three-momenta. The dual prescription can hence be either −ı0 ηk ji = −ı0 for some dual propagators or −ı0 ηk ji = +ı0 for the others, since indeed only the sign matters. In fact, this prescription encodes in a compact and elegant way the contribution of the multiple cuts that are introduced by the Feynman tree theorem [37]. The on-shell condition is given byδ , and determines that the loop integration is restricted to the positive energy modes, q i,0 > 0, of the on-shell hyperboloids (light-cones for massless particles) of the internal propagators. We also introduce the short-hand notation for the loop integration measure in d dimensions, where µ is an arbitrary mass scale to compensate the extra dimensions generated by d-dimensional integration measure. In the following, we use d = d s = 4 − 2 according to the convention of ref. [38]. In order to generalise LTD to higher orders, we introduce the following functions [21] where α k labels all the propagators, Feynman or dual, of a given subset. An interesting identity fulfilled by these functions is the following involving the union of two subsets α i and α j . These are all the ingredients necessary to iteratively extend LTD to multi-loop level. For example, at one loop, the Feynman and the dual representations of a N -leg scattering amplitude are is the numerator that depends on the loop momentum 1 and the fourmomenta of the N external partons {p i } N . In the absence of multiple powers of the Feynman propagators, the numerator is not altered by the application of the Cauchy theorem. However, the calculation of the residues of multiple poles to obtain the corresponding LTD representation requires the participation of the numerator. This is represented in eq. (2.5) by the symbol ⊗. At two-loop level, all the internal propagators can be classified into three different subsets (e.g. those depending on 1 , 2 and their sum 12 = 1 + 2 , as shown in figure 1). Starting from the Feynman representation of a two-loop scattering amplitude by applying LTD to one of the loops (eq. (2.5)), we obtain in a first step: Before applying LTD to the second loop, it is necessary to use eq. (2.4) to express the dual function G D (α 2 ∪ α 3 ) in a suitable form. The identity in eq. (2.4) splits the dual integrand into a first term that contains two dual functions -and therefore two internal lines on shell -and two more terms with a single dual function and Feynman propagators involving the other two sets of propagators, to which we can recursively apply LTD. The final dual representation of the two-loop amplitude in eq. (2.6) is In eq. (2.8), it is necessary to take into account that the momentum flow in the loop formed by the union of α 1 and α 2 occurs in opposite directions. Therefore, it is compulsory to change the direction of the momentum flow in one of the two sets. This is represented by adding a sign in front of e.g. α 2 , explicitly we have Changing the momentum flow is equivalent to select the negative energy modes. For the internal momenta in the set α 2 , this meansδ (2.10) The dual representation gets its simplest form if the Feynman representation contains only single powers of the Feynman propagators. This restriction cannot be avoided anymore at two-loops where, for example, self-energy insertions in internal lines lead automatically to double powers of one propagator. However, all the double poles can be included with a clever labelling of the internal momenta in the set α 1 , exclusively, which is not integrated in the first instance. Therefore, we have assumed that the numerator in eq. (2.7) is not affected by the application of LTD. The final dual representation in eq. (2.8) depends, in general, on the explicit form of the numerator. Again, this is represented by the symbol ⊗.
The number of independent double cuts in eq. (2.8) per Feynman diagram is 11) where N counts the number of propagators in each set. Therefore, it is convenient to have α 2 as the set with the smallest number of propagators. For planar diagrams, the set α 2 contains one single propagator. It is interesting to note that although the integration over the loop three-momenta is unrestricted, after analysing the singular behaviour of the loop integrand one realises that thanks to a partial cancellation of singularities among different dual components, all the physical threshold and IR singularities remain confined to a compact region of the loop three-momentum [39,40]. This relevant fact allows to construct mappings between the virtual and real kinematics, which are based on the factorisation properties of QCD, to implement the summation over degenerate soft and collinear states for physical observables in the Four-Dimensional Unsubtraction (FDU) formalism [30,31,41]. This framework, however, is not going to be needed in the following. This is because, as we shall see in section 6 where we analyse the cancellation of singularities at two loops, no infrared singularities remain when considering the full amplitude.

Algebraic reduction of two-loop dual amplitudes
In order to make the two-loop expressions more compact, we perform an algebraic reduction of the dual amplitudes to dual integrals that involve both positive and negative powers of dual propagator denominators. We analyse only the case of planar diagrams that are those that appear in the practical example that we present. Let us first consider scattering amplitude with N external legs and ordered external momenta {p 1 , p 2 , . . . , p N }. At one-loop, we have N different propagators and N − 1 independent scalar products 1 · p i (indeed, because of momentum conservation, 1 1 · p i , with 1 the loop fourmomentum). In the Feynman representation, the propagators are quadratic in 1 , while in LTD the dual propagators are linear. In both formalisms, however, it is possible to write numerators in terms of propagators diagram by diagram. Now, we consider the set of two-loop planar Feynman diagrams constructed from the ordered one-loop seed diagram -that is, all the two-loop diagrams that have one loop line involving one single propagator (see figure 2 for the assignment of momenta). These planar two-loop Feynman diagrams can be constructed from the sets of propagators In LTD at two loops, two internal particles are set on shell, which means that only 2N − 1 dual propagators remain for a given double cut. Moreover, dual propagators are linear in each of the loop momenta, and the dual numerators do not involve squared loop-momenta. Therefore, for each double cut, and considering all the Feynman diagrams with the same ordering of the external particles, it is possible to rewrite all the scalar products involved (and thus the numerators) in terms of dual propagators, and this in a unique way. Due to the assignment of the loop momenta, all the required irreducible scalar products (ISP) are automatically introduced. This is because the set of Feynman diagrams contains all the necessary propagators to perform the complete algebraic reduction.
The algebraic reduction of a planar two-loop dual amplitude with N external legs, and at most squared propagators in one single loop line, leads to The scalar coefficients c a0;a1,...,a 2N −1 depend only on the external momenta, and are not necessarily independent. Our purpose is to rearrange the expressions for the dual amplitudes in order to obtain the minimal set of independent coefficients c a0;a1,...,a 2N −1 . Another relevant issue to obtain the most compact integrand expressions is to label the internal momenta in the most symmetric way. In figure 2, we show the assignments that we use in the most general case for planar two-loop diagrams. Computer algebra programs for the automatic generation of two-loop amplitudes, like FEYNARTS [34], might use a different criteria which require a relabelling of the internal propagators to achieve the most suitable assignment. In the next sections, we shall illustrate the full procedure with the benchmark amplitude H → γγ.

Tensor projection and H → γγ at one loop
The scattering amplitude describing the Higgs boson decay to two photons is given by with e the electromagnetic coupling, e f (respectively N f C ) the electric charge (respectively the number of colours) of the virtual particle f , and ε(p i ) the polarisations vectors of the external photons. Here, we assume that the two photons are coupled to the same flavour. We are interested in the QED corrections at two-loop where A (1,f ) µν is the one-loop amplitude, and A (2,f ) µν is the two-loop QED correction. The tensor amplitudes A (L,f ) µν can be decomposed through Lorentz and gauge invariance as in terms of the tensor basis The tensor structure T µν 6 may appear for the first time at two-loop, because of the potential presence of a γ 5 , but its interference with the one-loop amplitude vanishes. As in ref. [19], we use the projectors to extract the scalar amplitudes A , though, as it can be used to simplify expressions. 1 At this point is is not necessary to make the distinction, because...

The dual representations of the one-loop amplitudes A
(1,f ) 1 and A (1,f ) 2 were calculated in ref. [19] in terms of the global factor g f and scalar coefficients c For the three different virtual particles that we considered, these coefficients are given by The vanishing one-loop amplitude A (1,f ) 2 is given by where The on-shell loop energies are given by As stated above, we know that A (L,f ) 2 = 0 after integration due to gauge invariance. Therefore, we can use this feature to simplify the expression for A µν . Already at one loop (L = 1), we showed that the following transformation [19] A with c reduces the number of necessary independent scalar coefficients c from three to two. Notice that while they are labelled differently because they do not have the same integrand-level expressions, we have A (L,f ) = A (L,f ) 1 after integration 2 . For the complete expression of the amplitude, we obtain It is remarkable to note that the dependence on the nature of the internal particle appears in eq. (4.7) and eq. (4.11) only through the scalar coefficients c (f ) i , defined in eq. (4.6). Although A (1,f ) is UV finite, because there is no direct coupling of the Higgs boson to photons at tree level, its integrand expression still exhibits a local UV behaviour. For that reason, we defined in ref. [19] the UV counter-term with q (+) The UV counter-term in eq. (4.12) integrates to zero in d-dimensions, though, it is used to cancel the local UV behaviour of eq. (4.11) in such a way that the locally renormalised one-loop amplitude A (1,f ) R can be obtained without altering the dimensions of the space-time (4.13)  figure 3 where all those that share the same global topology are superimposed in the so-called mandala diagrams. In this paper, we only consider QED corrections, and therefore photons as the extra internal particle, and do not take into account "mixed" diagrams where different massive particles may appear. In the traditional approach, massless snail diagrams are usually ignored, since they integrate to zero. However, within our approach, we need them to preserve the universal structure of the integrands. All the diagrams are planar and can be constructed from the following internal momenta: Only q 4 (and q 4 ) is massless (it labels the photon), while all the other internal momenta have mass M f . If the Higgs boson is on shell, the loop amplitude is below threshold and is therefore purely real. In that kinematical regime the dual prescriptions become irrelevant, and the dual functions fulfil the identity Hence, the LTD representation in eq. (2.8) adopts the simpler form Following the algebraic reduction defined in section 3, For a given dual or Feynman propagator, we have defined the dimensionless denominator For example, in terms of these dimensionless denominators, the one-loop amplitude in eq. (4.11) takes the form From now on, we use a different global factor depending on whether the expressions it multiplies has been algebraically reduced or not. The usual g f will be used for unreduced expressions, and g (L) f for reduced ones.

The two-loop amplitude
µν . Due to gauge invariance, it has to vanish after integration. Still, it is interesting to obtain an explicit expression because it establishes a useful integrand relation that can be used afterwards. Remarkably, it can be written in a very compact form for f = {φ, t}: and where the coefficients c  f and g (2) f . This is due to the presence of the second loop measure and the additionalδ, whose product has dimension of mass squared. The following contribution Consequently, the sum over integrals in eq. (5.10) vanishes.

5.2
The two-loop amplitude A (2,f ) In this section, we apply the same transformation in eq. (4.9) at two-loop (L = 2) to simplify the expressions for the amplitude A For the top quark and the scalar, the two independent coefficients that appeared at one loop are (see eq. (4.6)) At two loops, they are still present, along with the following extra coefficients: 14) Notice that these coefficients can be highly simplified in the particular case d = 4, motivating further the study of a full four-dimensional representation of this process at two-loop level.

Cancellation of integrand singularities at two loops
In addition to be infrared safe, the scattering amplitude for H → γγ is purely real if the Higgs boson is on shell. It is therefore completely free of soft and physical threshold singularities. Still, some of the dual propagators might go on shell inside the integration domain, leading to singularities of the integrand. As we demonstrated in refs. [39,40] at one-loop level, one of the advantages of LTD is the partial cancellation of potential singularities of the integrand. In this section, we extend the analysis of the integrand singularities of H → γγ at two-loop level. If a dual propagator G D (q i ; q j ) becomes on shell where both internal momenta q i and q j belong to the same loop line α k , then the corresponding singular behaviour is equivalent to the one-loop case. For example, consider the dual propagator G D (q 12 ; q 3 ). A physical threshold (forward-backward singularity 3 ) occurs if (see ref. [39]) H . An integrand singularity (forward-forward singularity) occurs if This would be the case of, e.g., G D (q 1 ; q 12 ) because k ji = p 2 , but this integrand singularity cancels in the sum of the two dual contributions involvingδ(q 1 ) andδ(q 12 ). The analysis can be extended easily to the other dual cuts, and similar conclusions are found. Now, we consider the genuine two-loop case where a dual propagator becomes on shell in the double cut of two propagators that do not belong to the same loop line. For the rest of this section, we consider q i ∈ α 1 , q j ∈ α 2 and q k ∈ α 3 . (6.3) We study study the quantity 4 where j indicates that we reverse the momentum flow of q j , as explained in section 2, namely We therefore have only depends on external momenta. The quantity S ijk becomes singular in the limits summarised by the condition λ ±±± → 0, with Accorded to eq. (6.6), only four independent solutions have to be considered, namely λ +++ , λ ++− , λ +−− and λ −−− . For two of these limits, there is a perfect cancellation of the integrand singularities, as indeed while for the remaining two limits, Although these singularities remain in the general case, it is possible to show that λ +++ (resp. λ −−− ) can only cancel if, in the particular case where m i = m k = M f and m j = 0, Because of the fact k does not depend on the loop momenta, the highest possible value of k 2 is (p 1 + p 2 ) 2 = s 12 , reached for instance when considering G D (q 1 , q 4 ; q 12 ). This means the second condition of eq. (6.11) reduces to 4M 2 f s 12 < 1 , k,0 . In that case λ +− and λ −+ can vanish, then enhancing the integrand singularity in 2 , but the overall singularity cancels between dual contributions. The other solution, with λ ++ = 0, is not possible because q If the three propagators do not interact in the same vertex, then k k(ij),0 = 0. This configuration includes the cases where q i and q k belong to the same or to different loop lines. There are solutions to eq. (6.13), but again either they cancel among dual contributions or eq. (6.1) is not fulfilled. In all the cases, the soft singularities of the integrand in 2 do not translate into soft singularities of the amplitude because of the integration measure.

Algorithmic approach to two-loop local UV renormalisation
Let us consider a generic unrenormalised two-loop amplitude A (2) , written as which we will assume to be completely free of any infrared singularity. Thus, only UV singularities may appear when either or both of | 1 | and | 2 | go to infinity. The local two-loop UV counter-terms are built recursively by first fixing one of the two loop momenta, say j , and expanding the integrand I( 1 , 2 ) up to logarithmic order around the UV propagator [29] G F (q i,UV ) = 1 where the arbitrary scale µ UV represents the renormalisation scale, and q i,UV = i + k i,UV . For simplicity, we take k i,UV = 0. The quantity where A (2) i,UV denotes the two-loop amplitude A (2) in the limit | i | → ∞, is not necessarily UV finite when both loop momenta are simultaneously large. It is necessary to subtract also the double UV behaviour of eq. (7.3). With this contribution, which is represented by A and is UV safe in all the limits. As an example, we consider This integrand produces a UV singularity when | 1 | → ∞, but it is superficially regular in 2 , meaning A 2,UV = 0. Computing the remaining counter-term gives which effectively removes the UV behaviour in 1 , but at the same time also introduces a singularity in 2 . It is therefore necessary to introduce the additional counter-term A UV 2 to fix the UV behaviour when both loop momenta go to infinity. This is done by expanding eq. (7.3) for very high values of 1 and 2 , while never neglecting one compared to the other. In this example, we get (7.7) Then, the renormalised amplitude, as defined in eq. (7.4), is finite in the UV. It is still necessary, though, to introduce subleading contributions to fix the renormalisation scheme. This is better explained in the following for the H → γγ two-loop amplitude.
With the labelling of the internal momenta that we have adopted for the H → γγ amplitude, it is more convenient to express the UV behaviour at two loops in terms of q 1,UV = 1 and q 12,UV = 12 , with 2 = 12 − 1 . Explicitly, the single and double UV behaviours are implemented by making use of the following transformations then expanding for λ → ∞ and truncating the corresponding series in λ up to logarithmic degree. This last operation is represented by the function L λ . In particular, the UV counter-terms are defined as A (2,f )

Higgs boson vertex renormalisation
In the Feynman gauge, the one-loop QED correction to the Higgs boson vertex exhibits the UV behaviour where Γ From the expression of the vertex counter-term in eq. (7.11) we can construct the UV counter-term of the two-loop scattering amplitude in the limit | 1 | → ∞ with | 12 | fixed. It reads instead of A (1,f ) in eq. (7.13) would therefore alter the UV behaviour of the single counter-term and not properly remove the corresponding infinities. It is only when considering the double UV counter-term (section 7.4) that the one-loop amplitude implicitly gets renormalised. (7.14) where q (+) 1,UV = 2 1 + µ 2 UV . Since the diagrams (2 for the top quark, 3 for the charged scalar) that contribute to the Higgs vertex correction are the only ones that are divergent when | 1 | → ∞, we directly have A

Photon vertex renormalisation
The one-loop correction to the photon interaction vertex to top quarks in the UV is given in the Feynman gauge by In eq. (7.15), the term proportional to ∆ (1,t) γ,UV integrates to zero in d space-time dimensions. Similarly to eq. (7.12), the coefficient of the integrated vertex counter-term is given by Although the integrated UV vertex correction is proportional to the tree-level vertex Γ (0,t) γ = ı e e t γ µi , thanks to the replacement q µ1 12,UV q µ2 12,UV → q 2 12,UV g µ1µ2 /d, we cannot use this replacement in the unintegrated form because it would alter the local UV behaviour. We must keep the full expression in eq. (7.15), including especially the term proportional to ∆ (1,t) γ,UV , to construct the local UV counter-term of the two-loop scattering amplitude.
For charged scalars as internal particles, we need to consider both the three-point and the four-point interaction vertices. For the three-point vertex there are three contributing diagrams, and the corresponding counter-term reads where q i−1 and q i are the outgoing and incoming internal momenta, respectively. For example, q i−1 = q 3 +p 1 and q i = q 3 +p 12 for the vertex correction with emission of a photon with momentum p 2 in the lower corner of the two-loop Feynman diagram.
For the four-point interaction vertex, there are nine contributing diagrams 6 , and we have γ,UV is the same as for the three-point interaction vertex. Again, in eq. (7.18) and eq. (7.20), we cannot apply the replacement q µ1 12,UV q µ2 12,UV → q 2 12,UV g µ1µ2 /d at integrand level even though the terms ∆ (1,φ) γ in eq. (7.18) and ∆ (1,φ) γγ in eq. (7.20) integrate to zero. Also notice that it was not necessary to introduce subleading terms for the scalar vertices because the finite part of the corresponding integrated counter-term is already 0.
The integrated counter-term reads where A (1,f ) γ is the sum of the two one-loop amplitudes involving triangle diagrams and A

Self-energies renormalisation
With our labelling of the momenta the self-energy insertions are defined in terms of the internal momenta q 4 = 2 and q i = 12 + k i , with i = {1, 2, 12, 3}. Explicitly, in the Feynman gauge we have (notice the relative sign in q i with respect to [31] because of the fact the momentum flows in the opposite direction) The UV expansion of eq. (7.26) reads where the coefficients d k,UV are subleading contributions to be fixed through the renormalisation scheme. It is remarkable that it has been possible to write the quark self-energy in terms of the same coefficients that appear in the Higgs boson and photon vertices. Notice that the expression in eq. (7.27) is simpler than the corresponding expression provided in ref. [31] and only differs at O( ), which does not have any consequence at the considered order.
The scalar self-energy corrections, which also include the snail diagrams, is written 4,u ). One possibility would be to subtract the contribution which is proportional to c (φ) T,UV before expanding in the UV, which would be equivalent to subtract a zero (actually this would not only work for the term generated by the snail diagrams, but also for any term that exclusively contains propagators depending only on 1 ). However, it would modify only the double cutsδ(q i , q 4 ). The UV expansion of eq. (7.29) reads The terms ∆ T,UV integrate to zero independently in d dimensions. In the latter, i.e. eq. (7.32), it was necessary to include a subleading contribution, proportional to µ 4 UV . We will not provide the integrated and dual expressions for A

The double UV counter-term
While it is entirely possible to compute A (2,f ) UV 2 by directly taking the sum of all contributions, it is more interesting to consider well-chosen subsets of diagrams -it also lightens intermediate expressions. For the top as the internal particle, it is logical to consider all the contributions to the Higgs boson vertex corrections, all the contributions to the photon vertex corrections and all the contributions to the self-energy corrections, as there is no ambiguity or cross-contributions for the H → γγ process at two-loop. They will be written A (2,t) H,UV 2 , A (2,t) γ,UV 2 , A (2,t) Σ,UV 2 , and account for 2, 4 and 6 diagrams, respectively. For the charged scalar as an internal particle, there is a subtlety. The three diagrams that contribute to the Higgs vertex correction also contribute to the γγφφ † vertex correction. For this reason, if we want to split A (2,φ) UV 2 , we have to consider both corrections together. For the photon and self-energy corrections, though, there is no ambiguity whatsoever. Thus, we define A (2,φ) Σ,UV 2 , that account for 9, 12 and 16 diagrams, respectively.
According to eq. (7.10), the unintegrated double UV counter-terms have the form and with n i being positive integers. Note that even though N (f ) {H,γ,Σ} (q 1,UV , q 12,UV ) should be expected to also depend on the external momenta p 1 and p 2 , it is a remarkable feature that, thanks to welcome cancellations, it does not when considering the sum of all contributing diagrams. The expression in eq. (7.34) is therefore free of irreducible scalar products and can very easily be reduced through integrations by parts to the form A where Cl 2 is the Clausen function of order 2, is the sunrise scalar integral with missing external momenta and all the internal masses equal, and is the massive tadpole. Because of the presence of the double pole inside eq. (7.37), it is, in the general case, necessary to keep track of the different normalisations as choosing one over another could lead to a shift in the finite part. As we are working in the MS scheme, we should rather factorise the usual S MS = (4π) e − γ E , but doing so, a global factor equal to S /S MS = 1+δS = 1+π 2 2 /12+O( 3 ) would appear.
We will see that in the end considering one normalisation over the other does not introduce any mismatch, so we chose to keep S for simplicity. These two master integrals (MI) are the only ones needed to evaluate A H+γγ,UV 2 µ 4 , (7.48) which after IBP reduction gives We obtain H+γγ,UV 2 + O( ) .
(7.52) After integration, we find and leading to The coefficients d It is remarkable that for both particles, the full double UV counter-terms does not exhibit any -poles, justifying in the mean time that using S instead of S MS does not introduce any discrepancy in the final result. While for the top quark each intermediate double UV counter-terms are all finite, for the charged scalar the cancellation of divergences only occurs when taking into account the sum of all 37 diagrams. This is due to the fact there are more subtle interplays between different contributing topologies, because of the more complex gauge structure. In both cases, the absence of divergences after integration means that in the traditional approach, these double UV counter-terms should not be needed. In our formalism however, they are essential in order to cancel local divergent behaviours appearing inside the amplitude and the single UV counter-terms. Note also that the unintegrated counter-terms exhibit terms proportional to d − 4 that vanish when taking the four-dimensional limit at integrand level, even though they still lead to finite parts when keeping the d dependence, because they multiply quantities that exhibit -poles. This means that while the counter-terms A (2,f ) UV 2 are finite, they will not lead to the same result if computed in d or 4 dimensions. It is only when considering the renormalised amplitude A (2,f ) R that we have a strict and rigorous commutativity between integrating and taking the limit → 0. The same happens at one-loop level [19], where the counterterm A (1,f ) UV would not be needed (it integrates to 0 in the MS) in d dimensions. One can wonder if this property holds at three-loop order and beyond, if it is exclusive to the H → γγ process and why, and if there is an underlying reason behind it.

Numerical integration
The renormalised unintegrated amplitude A ξ 1 (sin(θ 1 ) sin(ϕ 1 ), sin(θ 1 ) cos(ϕ 1 ), cos(θ 1 )) , Note that we can assume 1 , for instance, to belong to the (x, z) plane, as there is a global rotational symmetry; this means we can trivially perform the integration over one azimuthal angle -leading to a factor 2π -and set ϕ 1 to 0. Furthermore, it is better to use 12 instead of 2 as integration variable, because of the way the UV counter-terms have been defined. We therefore integrate over the two three-momenta 1 = √ s 12 2 ξ 1 (sin(θ 1 ), 0, cos(θ 1 )) , 12 = √ s 12 2 ξ 12 (sin(θ 12 ) cos(ϕ 12 ), sin(θ 12 ) sin(ϕ 12 ), cos(θ 12 )) , (8.2) which leads to five integration variables, namely ξ 1 , ξ 12 ∈ [0, ∞), θ 1 , θ 12 ∈ [0, π] and ϕ 12 ∈ [0, 2π]. In addition, the usual compactification of the integration domain is performed, where the domains of ξ 1 and ξ 12 are mapped from [0, ∞) onto to [0, 1], thanks to the change of variables with x i ∈ [0, 1], and i ∈ {1, 12}. This increases the stability of the numerical integration for very high energies, because it restricts the local cancellation of UV singularities to a compact region. The integration measure, after applying this change of variables and in four dimensions, reads Directly writing all the dual cuts explicitly in terms of the integration variables would not lead to a reasonable computational time, as the integration would be too heavy numerically speaking. Moreover, it would be far from being optimal, since many identical terms would have to be evaluated more than once. Instead, we take advantage of the fact it was possible to write all dual cuts in terms of the reduced denominators D i , as explained in section 3 and shown in appendix A. The first step, which only needs to be done once, is to express all the reduced denominators in terms of the 5 integration variables, and this for each dual cut 7 . We can 7 Recall that the expression of a given D i differs from dual cut to dual cut. Figure 4. Integrated renormalised amplitude of the two-loop corrections to the H → γγ process, as a function of the inverse mass square of the particle running inside the loop. On the left we show the top quark contribution and on the right, the charged scalar contribution. The solid blue lines represent the analytical result using DREG, while the red dots have been obtained numerically with the LTD formalism. For each particle, two values of the renormalisation scale have been considered, namely µUV = MH /2 = √ s12/2, and µUV = M f . We have rescaled our results as we used different normalisations. Explicitly, F then compute their numerical values for a given point in the integration domain, which allows us to quickly evaluate the integrand at this very point by appropriately replacing each D i . Thus, regardless of the complexity of a given integrand, only a limited amount of objects have to be numerically evaluated. Although quite simple to implement, this strategy helped decreasing the integration time by more than one order of magnitude. To perform the actual integration, we used the in-built MATHEMATICA function NIntegrate. Our results are shown in figure 4, where they are compared with the analytical results given in appendix B. The integration time is of a few minutes for each point. The biggest source of numerical error comes from the cancellation, between the amplitude and the double UV counter-term, of the non-decoupling term going as O(M 2 f /s 12 ). For M 2 t = 2s 12 for example, the part being removed from the amplitude by the counterterm is two orders of magnitude higher than the actual result, effectively multiplying the relative numerical error by a factor 100, roughly. Note that this error is twice as big for the top quark as for the charged scalar, because of the relative factor -2 between their respective O(M 2 f /s 12 ) terms. Nevertheless, the agreement with the analytical result is excellent for all values of the internal and renormalisation masses considered. Numerical instabilities may however appear when considering r very close to 4 (or equivalently M f very close to √ s 12 /2), as we approach the mass threshold. Beyond this limit (r > 4), contour deformation would be needed. Within the LTD formalism, numerical implementations of contour deformation have already been successfully implemented in [40].

Conclusion
In this paper, we have studied a purely four-dimensional representation of the renormalised Higgs boson decay amplitudes at two-loop order, so that it can directly be evaluated numerically. For this purpose, we have applied the Loop-Tree Duality (LTD) theorem to the H → γγ process, where internal particles are charged scalars and fermions. Working at the two-loop level within this formalism involves performing double-cuts and several algebraic manipulations of the expressions. This has been done through a fully automatised MATHEMATICA code which can be adapted to deal with any two-loop scattering amplitude.
Once the full sets of double-cuts were generated, we have taken advantage of the previously known one-loop results [19], in order to infer the universal integrand-level structure of the two-loop expressions.
Surprisingly, we have found many similarities between both cases, which has allowed us to write the full two-loop amplitude at integrand level using the same functional forms independently of the nature of the particles circulating inside the loop. As in the one-loop case, the explicit process dependence is codified inside specific scalar coefficients.
We have kept the d dependence in all intermediate steps, although a noticeable simplification of the universal coefficients takes place in the limit d = 4. We have studied the singular structure of the two-loop amplitudes to understand how to achieve a purely four-dimensional representation of the finite parts. In the first place, we have studied the cancellation of spurious threshold singularities that appear as a consequence of splitting the different cuts. When considering H → γγ below the physical threshold, the amplitude is guaranteed to be infrared safe as well as free of any physical threshold singularity. In fact, we have managed to prove that all the spurious singularities vanish in this configuration after putting together all possible double cuts.
Then, we have developed a fully local framework to remove ultraviolet singularities. We have implemented an algorithmic approach to renormalise, at integrand level, the two-loop amplitudes. It is based on the refinement of the expansion around the UV propagator strategy [19,30,31], including an iteration to remove all the possible UV divergences. A careful study of the UV structure of vertices and self-energies has been performed, which allowed to impose constraints on the finite remainder containing the specific renormalisation scheme dependence. Furthermore, we have also shown that the UV scale we used in the derivation of the local counter-terms corresponds to the renormalisation scale used in the traditional approach. These last two points allowed to build the required counter-terms at integrand level to reproduce the MS results.
Finally, we have proceeded to combine the universal dual representations of the H → γγ amplitudes together with the local UV counter-terms, achieving a fully renormalised integrand in four space-time dimensions. This has allowed a purely numerical implementation which fully agrees with the available results in the literature. On top of that, intermediate checks with the scalar sunrise diagram were implemented, in order to test the reliability of the codes.
The results presented in this paper constitutes a major advance into the extension of the Four-Dimensional Unsubtraction (FDU) formalism [30,31] to NNLO accuracy, paving the way for a more efficient and purely numerical implementation of currently relevant physical processes.