Master integrals for electroweak corrections to gg→γγ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$gg\rightarrow \gamma \gamma $$\end{document}: light quark contributions

We present a calculation of the master integrals (MI’s) required for the calculation of the Electroweak corrections to gg→γγ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$gg\rightarrow \gamma \gamma $$\end{document} production in which the process contains a light quark loop. The integrals can be broken down into five categories based on the flow of the heavy vector bosons throughout the loop. Three of the families are planar, and two are non-planar. We determine a canonical basis for each family which allows an efficient solution of the resulting differential equations via iterated integrals. We calculate the families in relevant physical kinematics and obtain an efficient numerical evaluation based on an implementation of Chen-iterated integrals.


Introduction
The production of prompt photon pairs has long been a staple process at hadron colliders with many observations ranging back over the last four decades [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19]. A centerpiece of this work and culmination of many years of effort was the discovery of the Higgs boson in this channel in 2012 [20,21]. Going forward into the era of the High Luminosity LHC (HL-LHC) the process will continue to be studied in ever greater experimental precision, mandating a continued theoretical effort to provide increasingly accurate predictions for this channel. At the lowest order in perturbation theory the production of photon pairs proceeds through a qq annihilation channel. For many years the theoretical state of the art corresponded to Next-to-Leading Order (NLO) QCD predictions, implemented into the Monte Carlo code Diphox [22]. However over the last decade significant theoretical progress has been made, including Next-to-Next-to Leading Order (NNLO) QCD predictions [23,24] and q T resummed predictions at order N 3 LL matched to NNLO [25]. Further, the theoretical predictions in association with one or more jets is well under control [26][27][28][29]. These results, in combination with recent computation of the three-loop amplitudes for qq → γγ [30] 1 suggests that it is not unreasonable to assume the availability of N 3 LO predictions over the course of the HL-LHC.
At O(α 2 s ) (NNLO in QCD) the gluon-initiated gg → γγ partonic configuration becomes accessible for the first time. For LHC kinematics there is a large flux of gluons, which somewhat compensates the suppression from the increased powers of the coupling. Thus for a long time [32] the gg initiated production of vector boson pairs has received special attention in the theoretical literature. The NLO QCD corrections to this process (which corresponds to two-loop diagrams) for massless quarks were presented around twenty years ago in refs. [33,34] and more recently massive top loops have also been included at this order [35,36]. The three-loop gg → γγ calculation [37] and two-loop gg → gγγ amplitudes (for light quarks) has further extended the known theoretical structure of this process in QCD.
One of the reasons the gg initiated process is so fascinating is its role in the continued study of the Higgs boson. The gg → γγ amplitude interferes with the amplitude for the production of a Higgs boson through a top loop and subsequent decay to photons. This interference can lead to changes in the measured properties of the Higgs [38][39][40] when compared to a signal only theoretical predictions. These initial studies have generated significant interest, leading to calculations over a wider range of production modes [41][42][43]. In particular, it has been noted that the interference, and hence width extraction, is sensitive to higher order corrections in QCD [44,45]. This motivates the study of the electroweak corrections to the interference pattern, of which the starting point is provided in this paper by the calculation of two-loop integrals required.
Over the last decade signifiant progress has been made in the computation of multiloop Feynman integrals due to the advances related to solution via differential equations. Using differential equations to simplify and solve Feynman integrals is a well-established technique [46][47][48][49][50], however, the field experienced a major advance when, in ref. [51], a refined basis approach was proposed. By modifying the basis from a traditional one, derived for instance via reduction by integration by parts identities, ref. [51] proposed a new canonical basis which obeyed a simpler differential equation which factorized the regulating parameter from the kinematics. If such a basis can be found, solutions to the differential equation can be more readily extracted. Since its introduction, the canonical differential equation technique has led to the calculation of many two-loop process for 2 → 2 and even 2 → 3 kinematics; for an incomplete list see e.g. [52][53][54][55][56][57][58][59][60][61][62][63][64][65][66][67][68].
One of the negative features of calculations beyond one-loop is the lack of common basis. This means that typically each new calculation often requires a dedicated calculation of the relevant MIs for the process, further since different authors use different normalizing factors and conventions, it is often inconvenient or impractical to combine results from different MI calculations. Our aim is to calculate the electroweak corrections to gg → γγ, in this first paper we present the MIs for the simplest set of Feynman diagrams, namely those which proceed through a loop of quarks which can be taken to be massless. Although more complicated MI calculations have been performed in the literature, we were unable to find specific MIs for several of our topologies (relating to diagrams involving three massive W bosons radiating photons). We therefore believe it is worthwhile to discuss these integrals in one place with a consistent framework. We postpone the discussion of the third generation quark loops, which contain many new MIs, to a future publication.
This paper proceeds as follows, in Section 2 we present a discussion of the canonical basis, and the solution of the differential equation via Chen-iterated integrals. Section 3 then presents the calculation of the planar topologies and section 4 presents the calculation of the non-planar topologies, finally we draw our conclusions in section 5. Electronic files containing solutions obtained in this paper are included as ancillary files with this paper.

2.2)
Momentum conservation further constrains the sum of the Mandelstam invariants s+t+u = 0, and for physical scatterings s > 0 and u, t < 0. We employ the 't Hooft-Veltman scheme in which loop momenta are evaluated in d = 4 − 2 dimensions and polarization vectors ( i ) and momenta attributed to the final state bosons are taken in four dimensions. We also note that when computing loop topologies later in this paper we introduce the momenta set p 1 , p 2 and p 3 , where the assigned signs of p i · p j can vary based on which of the k i we assign to each p j . We define the -loop color stripped amplitude as follows A a 1 a 2 = 16π 2 δ a 1 a 2 α 0 where α 0 s represents the bare strong and α 0 denotes the bare electromagnetic coupling, δ a 1 a 2 defines the leading order QCD color factor. In our calculation we will decompose the amplitude into a series of tensor structures following the procedure outlined in refs. [69,70]. In total there are eight relevant tensor structures Here the tensors T i are defined as: Each tensor has a corresponding form-factor F i . In defining the tensor decomposition we make the cyclic gauge choice i · k i+1 = 0, with k 5 = k 1 and impose the transversality condition i ·k i = 0. We use projector operators in order to isolate individual form factors entering into the expansion in eq. 2.4. Each projector is defined through an orthogonality condition with the corresponding tensor such that pol P i T j = δ ij . Since the resulting projectors are the same as those appearing in the literature [37,69,70] we do not list them explicitly here for the sake of brevity. The lowest order gg → γγ amplitude occurs at O(α s α) and proceeds through a closed loop of quarks. Here we are interested in EW corrections to the lowest order amplitude and are therefore calculating two-loop topologies. As a step towards the full O(α) corrections we present in this work the MIs associated with quark loops arising from the first two generations (q = u, d, c, s) for which the quark mass can be set to zero. We postpone analysis of the third generation to a future publication. We work in the t'Hooft-Feynman gauge with ξ = 1. The resulting diagrams can be classified in terms of five-master topologies, three planar and two non-planar for which representative Feynman diagrams are shown in Fig.1. We note that four of the topologies P I-III and N I are also relevant for the mixed QCD-EW corrections to qq → γγ, whereas N II can only occur with a closed loop of fermions. We note that, in addition to the topologies listed above which all contain at least one massive boson, there are contributions from a photon exchange. However, these contributions are well understood from the studies of gg → gg (see e.g. [51] for a canonical basis for these integrals) and we do not consider them explicitly in this paper (although we have checked our results obtained for these diagrams against the literature, finding perfect agreement).
After computing the two-loop diagrams and applying the projectors described above the resulting the form-factors consist of thousands of unique scalar integrals. However, through the application of well-known reduction techniques the amplitude can be written in terms of a much smaller set of Master Integrals (MI's). The choice of MIs is not unique, however, and a judicious choice is a set which satisfy differential equations of the following form [51] here G({x i }) denotes the vector constructed from the MI's which depend on the external scales x i (typically chosen to be dimensionless ratios). A({x i }) is a matrix which depends on these external scales, but crucially not on the dimensional regulating parameter . In this case the kinematics have factorized from the -series expansion and the integrals are said to be in canonical form. If such a basis can be found, the solution to the above differential equation can readily be obtained, and the MIs can be written as the following Chen iterated integral [71,72] Here γ is defined as a smooth map γ : i }) denotes a boundary vector which must be determined by independent means. If such a boundary vector is known, then the solution is obtained at a given weight by successive integrations over the integration kernel dA. The value of the integral does not depend on the chosen path γ, provided no singular points (or branch cuts) are crossed. By writing the vector of MIs as a series expansion in through 4 , and comparing to eq. 2.7, we can quickly define each coefficient: A pleasing feature of the canonical basis is that at a given order in every term has equal transcendental weight (equal to the power of expanded to). It is convenient to write the integration kernel in a dlog form as follows where η k defines a letter, M k is the coefficient matrix relative to said letter and the set of all letters {η k } defines the alphabet of the problem. Following the notation in ref. [53] we see that, when expanded to a given weight each term which enters the solution is of the form defines the integration of the letter along the chosen path. Following the notation of ref. [53] once more, we can more compactly define Our solution is then made up of many terms with different orderings of the letters k j . . . k i . If all of the letters are rational then the term may be written in terms of the well-known Goncharov polylogarithms: G(a n , . . . , a 1 ; x 0 ) = x 0 0 G(a n−1 , . . . , a 1 ; x 0 ) dt t − a n , (2.18) However, we will frequently encounter non-rational letters which must be directly integrated through eq. 2.14. Nominally one would expect to perform k-integrations to obtain a weightk contribution, however Chen-iterated integrals have several pleasing properties which can systematically reduce the complexity of the problem. We refer interested readers to the more complete overviews which can be found, for instance, in ref. [53]. Here we note only the most important identities required in our calculation. Firstly we note that at weight one, path independence renders the integration trivial, one must simply evaluate the integral at its end-points.
Then at weight two one can use this to perform the inner most integration and then perform the following integration by parts identity writing (2.20) such that all weight two contributions can be written as single parameter integrals. Once all weight two functions are obtained one can systematically reduce the number of integrations, by using the following generalized integration by parts identity where [γ t ] is a modified path which integrates over the domain [0, t] (and the complete path is restored as t → 1). Up to weight four order we therefore are able to write weight three contributions as one parameter integrals, and weight four contributions as integrals over two parameters. Finally, situations frequently occur when the inner most letter vanishes at the boundary of integration, which introduces potential singular contributions. In order to extract and remove theses contributions we employ a shuffle algebra to systematically extract singular logarithms (and cancel them with the corresponding terms arising at the boundary) σ(k k ),...σ(k 1 ),σ(j j ),...,σ(j 1 ) , (2.22) where k k k and j j j are the vectors of letters with length k and j and the sum runs over all permutations which preserve the relative order of k k k and j j j.
Our general setup is as follows: starting from the relevant Feynman diagrams for each topology we project out the form factors obtaining scalar integrals. These are then reduced using Kira 2.0 [73,74] to a minimal set of MIs. We then promote the obtained basis to canonical form using the Magnus-Dyson approach outlined in ref. [75]. We then use the physical properties of the Feynman integrals to determine boundaries at special kinematic points. Our final answers are then obtained via path integration away from these boundary vectors. In all instances we are able to check our results numerically using the package AMFlow [76] (we have additionally checked with pySecDec [77] where possible too). During our calculations we make frequent use of the software PolyLogTools [78] and HandyG [79]. The next sections describe these calculations in detail. The first planar topology we encounter involves the exchange of a single massive boson. The associated propagator structure has the following form These propagators form an auxiliary topology which is presented in Fig. 2. We note that, since the external momentum {k i } can be permuted into any ordering, we need to evaluate this topology in three distinct regions defined by the signs of p i · p j . All the integrals that appear in the physical amplitude can be expressed, via integration by parts identities (IBPs), in terms of a minimal basis of Master Integrals (MIs). We used Kira [73] to obtain such set of identities and we found that for the sum over all Feynman diagrams with a single massive boson with a fixed order of the final state momentum we need a total of 25 MIs. Therefore, in order to obtain a basis in canonical form we start with the following 25 ..a 9 ) = (a 1 1 + + a 4 4 + + a 6 6 + + a 8 8 + )(a 2 2 + + a 5 5 + + a 7 7 + + a 9 9 + ) (3.4) Where j + increases the power of the j th propagator by one. In order to write our basis in canonical form, we introduce the following finite and dimensionless integrals where α is defined for J X = f (I P I a 1 ...a 9 ) as ( a i ) − 4, and we utilize the techniques of the Magus-Dyson approach outlined in ref. [75] to obtain the following canonical basis: G P I 10 = yF 10 , G P I 11 = yF 11 , G P I 12 = y 2 F 12 , Here the differential equation is expressed in terms of the following dimensionless ratios, This basis is indeed in canonical form, and satisfies the following differential equations For this topology the MIs depend only on rational functions of x and y. As a result, the solutions of the differential equation can be expressed in terms of Goncharov polylogarithms (GPLs) and efficiently evaluated using HandyG [79] at every point in the phase space. In particular, the kinematic regions of interest in terms of the variables in Eq.3.7 are given by The boundary conditions can be determined for these integrals as follows: • G P I 3 and G P I 4 are fixed since they are finite as x → 0, and in this limit G P I 3 = G P I 1 .
• Similarly G P I 6 and G P I 7 are finite as y → 0, and in this limit G P I 6 = G P I 1 .
• G P I 8 and G P I 9 can be determined in terms of G P I 2 by demanding that these integrals are finite as x → 0 and noting that G P I 8 vanishes in this limit. Further G P I 2 can then be constrained in terms of G P I 1 by demanding that G P I 9 is finite as x → 1. In a similar fashion integrals 10 and 11 and 5 are fixed with x ↔ y in terms of G P I 1 .
• G P I 12 is the product of two trivial one-loop integrals and can be determined quickly via direct integration.
• G P I 13 and {G P I 14 , G P I 22 } are fixed since they must vanish as x → 0 and y → 0 respectively.
• G P I 15−21 , and G P I 23−25 are fixed by demanding finiteness in the limit x → −y.
Our integrals are then completely determined in terms of the overall normalizing integral G P I 1 which is easily computed by direct integration, to our required order we have The evaluation of these integrals for a sample point for each kinematic region is presented in Appendix A.1. We have checked our results against the numerical package AMFlow [76], finding perfect agreement. The second planar topology we will consider arises when one external photon is radiated from a W -boson and the other is radiated from the quark loop. The various permutations of photons and gluons can be written as functions of the following auxiliary Feynman integral where: This topology is shown in Fig. 4. For contributions to gg → γγ considered here, D 8 and D 9 do not appear with positive power. As before, by including the various permutations of physical momenta {k 1 , k 2 , k 3 } we find that we must calculate this auxiliary topology in two distinct regions, depending on the sign of the associated p i · p j . Following the same procedure as before, we reduce all the integrals appearing in the physical amplitude using canonical basis we start with the following set of integrals: We repeat the steps outline in the previous section to rewrite our basis in canonical form by first rendering the integrals finite and dimensionless and ultimately we obtain the following canonical basis: where the differential equation variables are defined to the same as those in the previous topology When considering all possible assignments of physical external momenta, we find we have to consider two kinematic regions, which we define as: The basis in Eq. 3.18 is in canonical form, and satisfies the differential equations where A x and A y have no dependence on . Additionally, since the coefficient matrices depends only on rational functions of the variables x and y, the solution can be written entirely in-terms of Goncharov polylogarithms and can be readily determined, once appropriate boundary conditions are known. The boundary constants can be determined as follows, • The finiteness of G P II 5 at y = 0 fixes G P II 6 and the relation G P II • The finiteness of G P II 9 at x = 0 fixes G P II 8 and the vanishing of G P II 9 at y = 0 fixes G P II 9 , such that G P II 8 and G P II 9 are determined in terms of G P II 2 . G P II 2 can be further constrained by requiring that G P II 8 is finite as x → 1.
• The finiteness of G P II 10 at y = 0 and G P II 11 at y = 1 fixes G P II 11 and G P II 7 . G P II 10 can be further constrained since the integral vanishes at y = 0.
• G P II 13 and G P II 14 vanish as x → 0.
• G P II 15 and G P II 16 , and G P II 22 are finite as y → 0.
• G P II 17 and G P II 18 are finite as x → −y and x → 0.
This fixes all of our integrals in terms of the overall normalizing integral G P II 1 , which can be readily determined by direct integration: We note that our normzaling integral is the same as that of the P I topology. In obtaining the boundary conditions we made use of the PolyLogTools [78] package and it's interface to GiNaC [79] to write combinations of boundary GPLs as simple functions of ζ i . We present results for these integrals in both regions of interest in Appendix A.2, which we obtained using HandyG [79] to numerically evaluate the GPLs. We have checked our results against the numerical package AMFlow [76], finding perfect agreement.   The final, and most complicated, planar topology arises from Feynman diagrams which contain three massive vector bosons. The contributions to the form factors F i arising from the planar topology P III can be written in terms of the following scalar loop integral

Planar topology P III
. (3.24) The propagators entering the loop are defined as follows: These propagators form an auxiliary topology which is shown in Fig. 6. We note that this auxiliary topology is equivalent (up-to re-labeling) as that appearing in P II . For the physical scattering process propagators D 8 and D 9 do not appear with positive powers (a 8 < 0 and a 9 < 0). We have constructed the topology to in the same notation as topologies P II , which minimizes the number of unique MIs which must be computed, however this means that in order to correctly assign the physical momenta k i to the topology momenta p i , mappings such as {p 1 , p 2 , p 3 } → {k 3 , k 1 , k 2 } must occur. For physical scattering processes we must therefore evaluate the auxiliary topology with p 1 · p 2 < 0 and p 2 · p 3 > 0. The set of scalar integrals entering the form factors from this topology can be reduced to 32 MI's, (for which we used Kira [73]). In order to obtain a canonical basis we begin with the following set of We then utilize the same strategy as the previous sections by introducing the following finite and dimensionless integrals where α is defined for J X = f (I P III a 1 ...a 9 ) as ( a i ) − 4. We can introduce the resulting basis: where we have introduced the two dimensionless ratios (which are the same as those used in the previous planar topologies) For physical kinematics, since both photons are radiated from the W bosons, there is only one region to consider x > 0 and y < 0, and in the physical region |x| ≤ |y|. In the differential equations the following roots occur In principle it is possible to rationalize two of these roots simultaneously (r 3 and r 1 ), however the resulting variables are quartic in two parameters, which leads to a proliferation of letters entering the Chen-integrals. In this paper we pursue an approach which rationalizes just r 1 , and results in a smaller alphabet, but with more letters with non-rational components. We find this leads to a more optimal final numerical implementation in our framework. There are two sub-regions of interest, corresponding to whether or not y is above the threshold at y = −4. We define Region A for y ∈ [0, −4] and define we take z A > 0 and therefore x + xz 2 A < 4z 2 A ensures x ≤ −y. Region B is defined as y ∈ (−4, −∞) and here we write with z B ∈ (0, 1). In terms of these variables the differential equations take the following dlog-form where the alphabet {η k } is made from twenty-two letters which may be inspected in the attached ancillary files. The remaining non-rational functions which enter the alphabet are in region A, and in region B. We have written our alphabet in a form which allows for clean extraction in the limit x → 0. In terms of (x, z A/B ) variables, the following 25 MIs contain only rational letters and may be evaluated in terms of GPLs Many of these integrals are recycled from the previous planar topologies, the boundary conditions for the new integrals can be quickly determined via appropriate finiteness of vanishing of the integral in the limits x → 0, y → 0 and x → −y, of which the vanishing of the integrals as y → 0 is particularly useful. This procedure systematically relates MIs to one another (expanded through weight 5 to extract the pertinent limits), up to an overall normalizing integral which we again take be G P III 1 . This integral can be readily computed via direct integration and in our normalization we have which is the same as that used in the previous planar topologies. We evaluate GPLs using the HandyG [79] package, and have also made extensive use of the PolyLogTools [78] package and it's interface to GiNaC [79] to write combinations of boundary GPLs as simple functions of ζ i and log 2. The remaining 7 integrals {G P III 24 , G P III 25 , G P III 28−32 } contain non-rational letters and cannot be evaluated in terms of pure GPLs. As discussed in section 2, a convenient means to evaluate the solution is through Chen-iterated integrals. The pre-requisite for evaluation is the knowledge of the MIs at a boundary point, such that the MIs at any other point can be determined by path integration. For our calculation the remaining task is thus to determine appropriate boundary conditions to evaluate our Chen-iterated integral expressions.
In order to solve for boundary points we note the following limit is regular, and further is free of non-rational functions. We can therefore solve along the x = 0 boundary entirely in terms of GPL solutions, with letters drawn from the following set in region A. The unknown boundary constants attributed to the 7 integrals with nonrational functions in the bulk space can be quickly determined through the vanishing of the integrals as y → 0 or (z A → 0). The equivalent boundary GPLs in Region B can be obtained from these results via analytic continuation (using z A = i/z B ) when y < −4. We have written our alphabet to allow smooth extraction of the x → 0 limit, however in doing so we introduce the spurious letter η 1 = 1 − 2x, which introduces a (spurious) singularity in region A. In order to avoid paths which cross the resulting branch cut, we also solve the following differential equations which again, is entirely made up of GPL solutions (with the same letters as eq. 3.40), and known boundary constants. For region B, the boundary solutions at x = 0 are sufficient.
With our boundary constants in hand we are therefore able to evaluate the remaining 7 integrals via path integration starting at the boundary points described aboves. We present results for all integrals evaluated at specific phase space points in Appendix A.3. In order to demonstrate the flexibility of our code we present the evaluation of the most complicated 7 point integrals along the arbitrary line y = −2.6x − 1 where we show results from x = 1 to x = 3 in Fig.7. We made this choice to compactly demonstrate the lack of issues around x = 1/2 and to evaluate the integrals above and below threshold. Also shown in the figure is a comparison to the result obtained using the numerical package AMFlow [76] with a point taken in each region. We have checked all of our integrals with AMFlow 2 and we find perfect agreement within to the reported uncertainty from the numerical codes.  . Here x ∈ [0, 3] and y is varied along the line y = −2.6x − 1 between Region A and Region B, as defined in Eq.3.32 and Eq.3.33 respectively. Red and blue squares represent the real and imaginary part, respectively, of independent evaluation using AMFlow.   The first non-planar topology which we encounter occurs when a single photon is radiated from the massive line and has the following propagator structure with: These propagators form an auxiliary topology which is presented in Fig. 9. For physical scattering permutations of the external photons and gluons require us to evaluate the integral in all possible kinematic regions, e.g. with p 1 · p 2 and p 2 · p 3 having signs (+−), (−+) and (−−). In order to evaulate the Feynman diagrams occurring in this topology, 40 MI's are required. As a starting point we introduce the following 40 integral combinations (which are presented in Fig. 10): Here the dimension shifting operator is defined as ..a 9 ) = (a 1 1 + + a 5 5 + + a 7 7 + + a 9 9 + )(a 2 2 + + a 6 6 + + a 8 8 + ) As before we begin by rescaling the integrals to become dimensionless where again, α is defined for J X = f (I N I a 1 ...a 9 ) as ( a i ) − 4. Introducing the following dimensionless ratios We can write a canonical basis as follows: (4.7) Here, the 7-propagators family is defined as: (4.8) The MIs above contain the roots of the following functions and satisfy the usual canonical differential equations in x and y. Given the simplicity of the non-rational terms which enter the differential equations, it is possible to simultaneously rationalize r 1 and r 2 to obtain purely GPL solutions. In all regions we find the initial mappings provides a convenient set of variables to express the GPL solutions which only contain rational letters. In total there are 33 such integrals, G 1−31 , G 33 and G 35 . The alphabet for these integrals depends on the following letters Many of the integrals for this topology have been computed already in the previous sections, which allows for quick determination of the boundary constants for the GPL solutions, which includes all three propagator diagrams G 1−8 . The new integrals can be determined by requiring finiteness as v → 0 (G 12 , G 16 , G 18 , G 21 , G 22 , G 26 , G 30 , G 33 , G 35 ) or u → 0 (G 10 ).
The remaining integrals can all be determined by noting that they vanish as u → 0.
After fixing the GPL solutions in the (u, v) variables we are left with seven integrals which contain roots in the initial (x, y) variables. As argued earlier these integrals can be rationalized, puting the entire solution in-terms of GPLs. Since the required variable change for the remaining integrals are sensitive to the sign of x and y we discuss each region separately.
We begin by considering the case where both x and y are less than zero, we introduce the following two variables then our requirements that x < 0 and y < 0 are satisfied by taking j A > 0 and m A ∈ [0, 1]. These variable changes rationalize the roots and result in a solution which be written in terms of GPLs. The remaining boundary vectors for the seven unknown integrals can be quickly determined by noting their vanishing as u → 0 (m A → 0).
Next we consider the case where x > 0 and y < 0, here the roots can be rationalized using the following variables: Here m B > 0 and j B ∈ [0, 1] defines the region where x is below the threshold at 4, and the region 1 < j B < 1 + m 2 B defines the region above threshold. The final region x < 0 and y > 0 can also be expressed in terms of Region B variables, here the domains are taken to be m B > 0 and j B > 1 + m 2 B for y below the threshold at 4 and 1 + m 2 B < j B < 1 + m 2 B defining the region above y > 4.
To summarize, the kinematic regions for this topology are given by

15)
Region C : (p 1 · p 2 < 0, p 1 · p 3 > 0) → x < 0, y > 0, |y|> |x|. The second non-planar topology, and most intricate one, corresponds to diagrams with three massive vector bosons and has the following propagator structure , (4.17) with: These propagators form an auxiliary topology which is presented in Fig. 11. After reduction the Feynman diagrams can be expressed in terms of 32 MI's. Our starting point is the following list of integrals (these are also presented in Fig. 16 Here the dimension shifting operator is defined as D − (I N II a 1 ...a 9 ) = (a 1 1 + + a 5 5 + + a 6 6 + + a 7 7 + )(a 2 2 + + a 8 8 + + a 9 9 + ) +(a 3 3 + + a 4 4 + )( j =3,4 a j j + ) I N II a 1 ...a 9 . (4.20) We proceed as in the earlier sections, first we render each integral dimensionless via the following scaling Where α is defined for J X = f (I N II a 1 ...a 9 ) as ( a i ) − 4. We then write a canonical basis in terms of our original integrals and the following dimensionless ratios Given the restriction that both photons must be radiated from the W bosons all permutations of external photons and gluons result in x > 0 and y > 0. Finally, in terms of these variables we construct the following canonical basis: the 7-propagator family is defined as follows: Six non-rational roots appear in the integrals, defined via The basis integrals and resulting differential equations contain non-rational terms, several of which can be rationalized by an appropriate change of variables. In total there are 6 roots which appear for this topology, shown in Eq. 4.25 and not all of the roots can be simultaneously rationalized. One must therefore make a choice regarding which roots to rationalize to put the set of equations in a form most amenable for a solution in terms of Chen iterated integrals. Since it enters into many of the differential equations √ r 1 is clearly of form which is desirable to rationalize. Since we are interested in solving in the physical region y > 0 and the behavior of the MI's depends on whether one evaluates the integral above threshold (s > 4M 2 W =⇒ y < 1/4), or below (s ≤ 4M 2 W =⇒ y ≥ 1/4), we divide the phase space into relevant regions and introduce the following variable transformations With z A ∈ (0, 1] and z B ∈ [0, 1]. After implementing these transformations the following integrals G N III 1−13 , G N III 16−23 and G N III 27 are entirely rational and may therefore be written in terms of GPLs. We evaluate the remaining MIs as Chen-iterated integrals, for which we need to determine boundary solutions. We do so by following the same procedure as for the planar topology discussed in section 3.3, we solve the differential equation in z j at special points in x, namely  . Here x = 0.25 and and y is varied in the region y > 1/4 (Region B, as defined in Eq. 4.27). Red and blue squares represent the real and imaginary part, respectively, of independent evaluation using AMFlow.  Figure 16: Master integrals for topology N II . Dotted propagators have higher power of the respective denominator, while red propagators are massive.

Conclusions
In this paper we have presented the calculation of the Master Integrals required for the evaluation of electroweak corrections to the amplitude for gg → γγ for the cases in which the closed loop of fermions corresponds to the first two generations. We solved for the MI's using a differential equation method, specifically by choosing a set of MI's which obeyed a differential equation in canonical form (in which the series factorizes from the kinematics in the differential equation) we were able to solve our integrals using iterated integrals.
We are able to classify our integrals based on the flow of the massive vector bosons through the diagrams. The situation in which there is no massive bosons (i.e. QED corrections) is well known and we did not discuss the answer in detail here. Topologies for which there is a single massive boson in the loop was also straightforward, there are only planar topologies with no issues regarding non-rational letters, such that the answer could be written entirely in terms of GPL solutions. For topologies with two massive vector bosons there are planar and non-planar topologies to consider. The planar topologies with two massive bosons are similar to those with one, in that they are readily written entirely in-terms of GPL solutions. For non-planar topologies non-rational letters initially appear in the alphabet, however by an appropriate variable change we were able to rationalize the roots such that a GPL solution was obtained for all two W boson topologies.
Finally the most complicated topologies we had to consider came from diagrams which contain three massive vector bosons. Again here there are both planar and non-planar topologies to consider. The increased number of internal massive lines lead to a proliferation of non-rational letters in the alphabet of the differential equations. In order to solve the differential equations in this instance we used a mixed Goncharov-Chen approach. By writing the non-rational contributions as Chen-iterated integrals we were able to define as solution as a two-parameter integral when expanded up to weight 4.
We evaluated our integrals in all required regions of kinematic phase space for physical scattering. We checked all of our solutions against public numerical packages, finding prefect agreement, in general our analytic solutions were orders of magnitude quicker than the numerical checks, and should be amenable to implementation into a Monte Carlo generator for subsequent phenomenological analysis, which we will pursue in a subsequent publication.