Two-loop massive QCD and QED helicity amplitudes for light-by-light scattering

: We present the analytic and compact two-loop helicity amplitudes for QCD and QED corrections to the light-by-light scattering process with massive internal fermions. We express the master integrals either in terms of multiple polylogarithms or in terms of iterated integrals with dlog one-forms. We also elaborate on optimizing the analytic results for each phase-space region. This makes the numerical evaluation of the scattering amplitudes fast, stable and suitable for phenomenological applications.


Introduction
Light-by-light (LbL) scattering is a fundamental quantum electrodynamics (QED) process resulting from the quantum effects predicted in the 1930s [1][2][3][4].It is a loop-induced process via one-loop Feynman diagrams with charged internal particles in the loop, as illustrated in the left panel of figure 1.The lowest order of cross section is O(α 4 ), where α is the electromagnetic fine structure constant.In the low energy limit, the process can be effectively described by the Euler-Heisenberg Lagrangian [4].The first complete calculation in QED was carried out by Karplus and Neuman [5] in 1951.The interest in studying LbL revived recently, in particular from the point-of-view of searching for the elusive beyond the Standard Model (BSM) signals.It has been acknowledged that LbL can be used to probe the quartic anomalous gauge couplings [6], the axion-like particles [7], the graviton-like particles [8,9], the nonlinear Born-Infeld extension of QED as well as in the context of the Standard Model (SM) [10], the photon-photon self-interaction from the noncommutative QED [11], large extra dimensions [12,13], and supersymmetric particles [14].LbL is also a background for looking for new particles in the SM, such as ditauonium [15,16] and glueballs [14].
In spite of the early theoretical proposals, the direct experimental observations 1 of the LbL scattering process γγ → γγ were only achieved recently in the ultraperipheral heavyion collisions (UPCs) at the LHC [17][18][19][20].Its experimental feasibility, first suggested in ref. [6], can be mostly attributed to the large coherent photon flux carried by the ultrarelativistic nucleus thanks to its large charge number Z.For instance, the lead (Pb) beam used by the LHC has Z = 82.With UPCs, we have the largest magnitude of electromagnetic fields produced in the laboratory in an experimentally controlled way.It enhances the LbL cross section by Z 4 ≈ 4.5 • 10 7 in Pb-Pb collisions with respect to the counterparts in proton-proton (pp) and electron-positron (e − e + ) collisions.In a UPC exclusive process, the two initial beam hadrons do not break, which can be identified with detectors at very forward angles, Roman Pots and Zero Degree Calorimeters for protonproton and nucleus-nucleus collisions respectively.This usually provides us with distinct signatures compared to the central and peripheral collisions.The UPC exclusive processes have the virtue of low backgrounds due to their much smaller particle multiplicities than in the inclusive scatterings, and thus the experimental triggers can be loosened.This allows us to probe a different kinematic regime with respect to the inclusive reactions that breakup the beam hadrons.
The LHC measurements have been compared with the theoretical calculations by the two Monte Carlo event generators: SuperChic [21] and gamma-UPC [22] 2 .Both of them are based on the leading order (i.e., one-loop) matrix element and yield the consistent LbL cross sections.However, the theoretical cross sections are about 2 standard deviations below the ATLAS measured value.It is therefore interesting to investigate the origin of this discrepancy.On the theory side, it would be a natural step 3 to check the size of the next-to-leading order (NLO) quantum corrections that involve two-loop helicity amplitudes.In the low-energy limit, the two-loop corrections to the Euler-Heisenberg Lagrangian have been computed with the string-inspired approach [26][27][28][29].On the other hand, in the high energy limit, two-loop amplitudes with massless internal fermions in QCD and QED [30] and in supersymmetric QED [31] have also been available for two decades.The aim of this paper is to present the analytic two-loop helicity amplitudes with general massive internal fermions in the theories of QCD and QED, where some representative graphs are shown in figure 1.This will allow us to lift the above two approximations, which will also be relevant for the theory-experiment comparison.We invest efforts in simplifying the initial lengthy expressions into more compact and concise forms.This is not only for the purpose of aesthetics, but also to improve the numerical stability of helicity amplitudes at two 1 LbL has also been indirectly tested by the precise measurements of the anomalous magnetic moments of the electron and the muon and by the elastic scattering of a photon in Coulomb fields of nuclei via the so-called Delbrück scattering.
2 Strictly speaking, gamma-UPC is a library for calculating the coherent photon-photon flux in UPC.It has been integrated into the two event generators HELAC-Onia [23,24] and MadGraph5 aMC@NLO [25] to enable the event generations of exclusive UPC processes. 3Another possibly important missing effect is the higher order Coulomb O(Zα) corrections, which is however beyond the scope of our paper.
loops.This optimisation becomes particularly important when dealing with the substantial cancellations between individual Feynman graphs as a consequence of the gauge invariance and infrared (IR) finiteness.The corresponding two-loop Feynman integrals with massive internal lines have also been studied previously in refs.[32][33][34][35][36], though the results are only partially available in the literature.The rest of the paper is organised as follows.We derive the general structure of scattering amplitude to any loop order in section 2. The two-loop master integrals are discussed in section 3. The analytic expressions of the one-and two-loop helicity amplitudes with full fermion mass dependence are presented in section 4. We draw our conclusion in section 5.The appendices collect the additional information.We define the kinematic regions of our rationalised variables as well as the analytic continuation in appendix A. The analytic boundary conditions for the non-rationalised two-loop master integrals are given in appendix B. We define the one-loop master integrals in appendix C, and present the explicit expressions of our chosen canonical basis of two-loop master integrals in appendix D. Finally, we report the one-loop helicity amplitudes of the W ± boson contribution in appendix E, and collect the one-and two-loop helicity amplitudes in the low-energy limit in appendix F.

General structure of scattering amplitude
In this section, we derive the independent form factors and helicity amplitudes.A similar derivation exists in the literature [31].In addition, an alternative approach of obtaining helicity amplitudes is also described in ref. [37].However, our derivation is neither (completely) identical to ref. [31] nor following ref.[37].Therefore, for completeness, we present our derivation in detail in this section.
Let us consider the process where all the momenta p i of the photons are incoming and their helicities are denoted as λ i .The amplitude is expressed as with ε λ i ,µ i (p i ) being the polarisation vectors of external photons.The Lorentz decomposition of the scattering tensor M µ 1 µ 2 µ 3 µ4 is given as: 4 The coefficients A i , B i jk and C ijkl are functions of the Mandelstam variables s = (p 1 + p 2 ) 2 , t = (p 2 + p 3 ) 2 , u = (p 1 + p 3 ) 2 , as well as the masses of the particles in loops.Therefore, they should be understood as A i (s, t, u), B i jk (s, t, u), C ijkl (s, t, u).The on-shell conditions and the conservation of momenta guarantee that the sum s + t + u = 0. Thanks to the parity invariance in QCD and QED, the term proportional to the Levi-Civita tensor is forbidden.We have 3 + 6 × 3 × 3 + 3 4 = 138 independent functions in (2.3).However, many of the coefficients are related by the following three symmetries: 1. Transversality: ε λ i • p i = 0 for i = 1, 2, 3, 4.This reduces to 57 functions to be considered.

3.
Gauge symmetry: the gauge symmetry guarantees the following four Ward identities They lead to a few more identities to reduce to the 3 independent functions, which are Our derivation and all the relations have been cross checked with explicit calculations.The number of independent form factors (not functions) in M µ 1 µ 2 µ 3 µ 4 is 7, instead of 15 in section 3.1 of ref. [31].They are A 1 (s, t, u), A 1 (t, s, u), A 1 (u, s, t), ∆B 1 11 (s, t, u), ∆B 1 11 (t, s, u), ∆B 1  11 (u, s, t), and ∆C 2111 (s, t, u).
The reduced tensor structure becomes where In the above eqs.(2.7),(2.8),we have used the short-hands for the replacement rules "→" and the exchange rules "↔" with respect to the first terms on the right-hand sides.The tensor In the next step, we need to define the projection operators to obtain the coefficients A 1 , ∆B 1  11 and ∆C 2111 .We can use the following tensor where  is the Gram matrix.P is the projector that projects onto the (d − 3)-dimensional subspace perpendicular to the space spanned by the external momenta (i.e., p i,µ P µν = P µν p i,ν = 0).R j,ν is the dual vector to p ν j relative to this space (i.e., R j,ν p ν i = δ ji ).In addition, we have the following relations P µ µ = d − 3, where we have denoted the space-time dimension by d = 4−2ϵ with ϵ being the dimensional regularisation parameter that ought to be taken infinitesimal at the end.We can obtain the A i tensor coefficients using the following operations In other words, we have and (2.16) In the above, the right hand side appears to introduce the spurious pole 1/(d−4) = −1/(2ϵ).
On the other hand, the sum is free of the spurious pole.As we discuss later, the amplitude only depends on A S .Hence, there are no spurious poles.It is straightforward to define the projectors for the coefficients B 1 kl via Similarly, we have (2.20) With the same spirit, we can define the projectors for the coefficients With the above preparations, the helicity-summed amplitude square takes the form (2.25) As we see from eq.(2.24), at the level of helicity-summed amplitude square, we only need to calculate 5 independent form factors A S (s, t, u), ∆ B1 11 (s, t, u), ∆ B1 11 (u, s, t), ∆ B1 11 (t, u, s), and ∆ Ĉ2111 (s, t, u).Among these, there are no interference terms.The form factors A S (s, t, u), su∆C 2111 (s, t, u), and su∆ Ĉ2111 (s, t, u) are fully symmetric in the Mandelstam variables s, t, u.
Furthermore, the helicity amplitudes can be expressed as a linear combination of the form factors. Like the case of form factors, there are 5 independent helicity amplitudes. 5p to a global phase, we can write them as: (2.30) In addition, if we take into account permutations, the following relations hold between the 3 independent helicity amplitudes: While the helicity amplitudes M ++++ and M −+++ are fully symmetric in the arguments (s, t, u), the remaining ones The amplitudes in other helicity configurations can be related to the above mentioned 5 as follows: (2.32) Summing these helicity amplitudes, it is straightforward to verify that eqs.(2.23) and (2.24) hold.

Two-loop master integrals
It is well known that a form factor or an amplitude at L-loop can be written as a linear combination of a small set of L-loop Feynman integrals, known as master integrals.In this section, we first present the details for obtaining analytic expressions of two-loop master integrals for LbL scattering, before discussing in detail their numerical evaluations.
The one-loop master integrals are provided in appendix C. In dimensional regularisation, the two-loop QED amplitude for a given fermion species with mass m f , generated by Qgraf [38] and FeynArts [39], can be reduced into 103 d-dimensional master integrals across all the permutations by applying integration-by-parts (IBP) identities with the help of Kira [40,41] or FiniteFlow [42]+LiteRed [43].Modulo the permutations in the Mandelstam variables s, t, u, there are 30 master integrals.One among them is a product of a one-loop box integral and a one-loop tadpole integral, which can be solved separately.The remaining 29 integrals form a closed system in the method of differential equations [44][45][46][47] in d dimensions.We now explain how to analytically solve these 29 master integrals.
The two-loop master integrals for LbL in the top-sector have been investigated in ref. [32].In this context, the new ingredients of this article is the presentation of the analytic results in a form that is most suitable for phenomenological applications, facilitating the fast numerical evaluation of the scattering amplitudes across all the physical phase-space regions.In addition, to ensure the comprehensiveness of this article, we elaborate on all the essential steps needed to replicate our results, even at the risk of providing details that may already be well understood.We express the analytic results of the master integrals in terms of Chen iterated integrals [48] with logarithmic one-forms, and cast the results in a form suitable for phenomenological applications.The family of master integrals of interest is expressed as where a i ∈ Z.We represent the loop momenta through ℓ i and Euler-Mascheroni constant by γ E .The inverse propagators are The introduction of i0 + ensures the prescription of the Feynman propagators.We keep the general external momentum dependence with i, j, k ∈ {1, 2, 3, 4} and i ̸ = j, i ̸ = k, j ̸ = k.

The arguments of I
(2) , which are simply a permutation of s, t, u.
As mentioned above, to obtain the analytic results of the master integrals we use the method of differential equations.Within this approach, it is convenient to choose a canonical basis of master integrals that allow us to set up the differential equations into the so-called ϵ-form [49], which enables to solve the system afterwards.Our choice of the canonical basis for the two-loop master integrals is guided by the one outlined in ref. [32], although not completely identical.For completeness of this article, we present our chosen basis ⃗ f (2) = (f 29 ) in eq.(D.1) in appendix D. Each integral in the canonical basis eq.(D.1) is dimensionless and has uniform transcendental (UT) weight properties.Moreover, the system of differential equations of ⃗ f (2) attains the form where the differential equation matrix A (2) does not depend on ϵ.Due to the inclusion of massive loop propagators, the presence of square roots is unavoidable.The square roots present in our case are We can rationalise the first three square roots simultaneously by choosing the same variables as suggested in ref. [32]: However, we still have one square root left, which we are unable to rationalise further.There are in general multiple (complex) solutions of eq.(3.4) for w, z in terms of s ij , s jk .Our choice of the solutions, depending on the real values of s ij and s jk , is given in appendix A. The next step is to integrate the system of differential equations to obtain the analytic results.Our goal is to express the analytic results in terms of iterated integrals with polylogarithmic kernels, with explicit representation in terms of Goncharov polylogarithms (GPLs) [50] [i.e., multiple polylogarithms (MPLs)] wherever possible.We start with rewriting the one-forms as logarithmic one-forms and casting the differential equations in a more compact form.Feynman integrals expressed in terms of iterated integrals can be efficiently described via their symbols [51].In the variables defined in eq.(3.4), the alphabet A consists of the following 16 letters: where the last four depend on the square root ρ.In other words, the matrix A (2) in the differential equation (3.3) can be expressed in terms of these 16 one-forms: Furthermore, this matrix can be rewritten in terms of the following 15 independent linear combinations of dlog one-forms: Here we have dropped the arguments (w, z) in the one-forms.The non-zero matrix elements of A (2) are given in eq.(D.2) in appendix D. In the next step, we take the series expansion over the dimensional regulator ϵ in the master integrals ⃗ f (2) = ω=0 ϵ ω ⃗ f (2,ω) and the differential equation system eq.(3.3) becomes Since our choice of ⃗ f (2) is a UT basis, all terms contributing to ⃗ f (2,ω) should have the same transcendental weight ω.With a proper boundary condition, we can solve eq.(3.8) iteratively from lower to higher weights.The solution can be written in terms of iterated integrals with logarithmic one-forms.Up to weight ω = 4 that we need, most of the master integrals f (2,ω) n (x ij , x jk , x ik ) can be expressed through Goncharov polylogarithm G's, while only the iterated integrals for the following master integrals depend on the kernels with the square root ρ: (3.9) In our two-loop amplitudes, we actually only need 5 integrals f (2,4)

, and f
(2,4) 29 which cannot be written in G's.For the integrals with G's, 6 their numerical evaluations are rather straightforward with the algorithm outlined in ref. [52].In this paper, we use FastGPL [53] and take HandyG [54] as a cross check.Now, as a new ingredient from this section, we outline the steps taken to prepare all the master integrals present in the scattering amplitudes for phenomenological applications.We explain how to numerically evaluate iterated integrals with algebraic kernels using an idea that has been recently successfully utilised for other processes (see, e.g., ref. [55]).We also find different analytic expressions for all the master integrals valid in different regions of phase space in order to avoid performing analytic continuations.Moreover, we identify relations among the Feynman integrals in the scattering amplitudes, which reduces the number of iterated integrals to be computed numerically at each point of the phase space.For the numerical evaluations of our master integrals with square roots, we choose a technique elucidated in refs.[32,55].The kernels in the first two-fold integration can be rationalised and therefore can be expressed through logarithms and classical polylogarithms Li n in terms of x ij , x jk , x ik . 7For the weight-3 integral f (2,3) 18 , we are only left with a single integration to be carried out numerically, while there are two-fold integrations left for the weight-4 integrals.For the latter, the basic idea is to convert the last two-fold integration into a one-fold integration, and to perform the last one-fold integration numerically (cf.an example given in section 4.4 of ref. [32]).For numerical integration in all phase-space regions in appendix A, we also use an approach that makes the numerical evaluation of all the master integrals faster as we avoid crossing any physical singularity.Instead of integrating our results within a single region and then relying on analytic continuation to cross the boundaries into other regions, we obtain distinct analytic results that are specifically valid within each of these regions.Hence, our numerical evaluations of the integrals are valid within specific phase-space regions.We also obtain analytic boundary constants f (2,ω) n (x ij , 0, −x ij ) and f (2,ω) n (0, x jk , −x jk ) valid in all these regions.For these analytic boundary constants for the 5 iterated integrals with roots, refer to appendix B. For the last one-fold integration, the convergence of the numerical integration strongly depends on the path.We opt for the following path parameterisation from the starting point (s ij,0 , s jk,0 ) to the end point (s ij,1 , s jk,1 ) where r ∈ [0, 1].It is easy to verify that (s ij (0), s jk (0)) = (s ij,0 , s jk,0 ) and (s ij (1), s jk (1)) = of the former are much faster than the latter.This is possible thanks to the methods outlined in ref. [56] by matching symbols in both the (w, z) as well as (sij, s jk ) coordinate systems.However, the Riemann sheets for the multi-valued functions should be carefully picked up.
(s ij,1 , s jk,1 ).In this article, we choose the following boundary conditions: The integrals have been cross checked numerically for a few phase space points in each region against AMFlow [57].
As the final step in making the integrals pheno-ready and in order to simplify the final helicity amplitudes, it is necessary to find all possible relations among the master integrals needed for our amplitudes.Otherwise, not only can we not check pole cancellations at the analytic level, as a similar finding in ref. [58] for another process, but large numerical cancellations will also happen among different pieces, which prevents us from taking any interesting kinematic limit.In this work, we use the symbol techniques [59] and the shuffle algebra that the iterated integrals satisfy to find the following relations, by including the one-loop master integrals f (1,ω) n defined in appendix C. The relations can be categorised into three parts: • The following integrals are identical over permutations of external momenta: n (x ij , x ik , x jk ) for n ∈ {2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 24} , (3.12) • Some integrals with same transcendental weight can be related: • Some higher transcendental weight two-loop master integrals can be expressed into lower weight integrals: We remind readers that x ij = s ij /m 2 f , x jk = s jk /m 2 f , and x ik = s ik /m 2 f .Using the above relations, we manage to cancel the poles analytically.Furthermore, these relations also considerably reduce the number of master integrals that we need to compute at each transcendental weight from around 300 to 84 in the two-loop helicity amplitudes.

Analytic expressions of helicity amplitudes
In this section, we present our main results, i.e., the analytic expressions of helicity amplitudes for LbL.To obtain the helicity amplitudes we use Qgraf and FeynArts to generate the Feynman diagrams, and further manipulations such as performing d-dimensional Lorentz and Dirac algebra are done in FORM [60,61] or Mathematica.
We first present the one-loop form factors and helicity amplitudes as a warm-up exercise.In this section, we will only consider a given fermion loop with the fermion mass m f , while the one-loop expressions for the W ± loop can be found in appendix E. Throughout this section, let us define the following dimensionless variables with x s + x t + x u = 0. Let us assume the fermion has charge Q f and colour charge N c,f .8 In addition, we also denote the LO (one-loop) helicity amplitudes with the fermion f as M (0,0,f ) λ 1 λ 2 λ 3 λ 4 , while the two-loop QCD and QED amplitudes with the fermion are M (1,0,f ) and M (0,1,f ) λ 1 λ 2 λ 3 λ 4 respectively.

One-loop results
The one-loop form factors are 4 (x t , x u , x s ) x s x t (x s x t + 4x u ) + r (1) where the basis of rational coefficients is as follows: 3 5 6 3) The one-loop master integrals, defined in appendix C, can be expressed through GPLs: The prefactors (−1) θ(xs−4) and (−1) θ(xs)+θ(xt) ensure the correct analytic continuation when using the convention in section A, so that the multi-valued square-root function in eq.(4.2) takes the values in the first Riemann sheet.
Alternatively, the one-loop helicity amplitudes with a given massive fermion f loop M (0,0,f ) where 10 11 It is interesting to note that, while the conversion between the helicity amplitudes and the form factors is merely a linear transformation (see eq.(2.30)), the expressions for the amplitudes are shorter than those for the form factors.Such an observation motivates us to simplify at the level of the helicity amplitude rather than at the form factor.Note that the one-loop amplitudes have been given in ref. [62] in terms of one-loop scalar integrals, rather than as presented here in terms of a UT basis.The UT basis has the advantage of rendering shorter expressions.Our one-loop results have been numerically checked against MadLoop [25,63].

Two-loop helicity amplitudes
We are now ready to present the UV-renormalised two-loop helicity amplitudes 9 with the full fermion mass dependence.Following the IBP reduction, we obtain large rational coefficients in front of the master integrals.We put dedicated efforts to simplify these coefficients into more compact and well-organised forms as given below.This is achieved with the help of FiniteFlow [42] and MultivariateApart [64].Also, the relations between the master integrals at different weights as given in eqs.(3.12),(3.13)and (3.14) were one of the crucial ingredients for simplifying the final expressions.The size of the amplitude file is reduced from ∼ 300 MB just after IBP to around 180 KB counting with Mathematica.
After simplification, the two-loop QCD helicity amplitudes are given as follows, with ).The QED counterparts are obtained with a simple rescaling It is implicitly understood that the two-loop QCD amplitudes only apply to quarks, while the QED amplitudes are non-zero for all charged fermions, i.e., quarks and charged leptons in the SM.
The most-lengthy amplitudes are the two-minus-two-plus helicity amplitudes.Thanks to the relation in eq.(2.31), we only need to present M (1,0,f ) −−++ , which is: x j (x j − 4)f (1,1) 3 (2) 16 (x i , x j , x k )f (1,2) 5 15 (x i , x j , x k ) x j (x j − 4) 14 Here the additional rational coefficients are: 13 14 15 16 17 18 (x i , x j , x k ) = x j (x j − 4) 19 20 21 22 23 24 25 26 27 28 29 30 31 The amplitude M (1,0,f ) −−++ is symmetric in x t and x u .The expressions of all these twoloop amplitudes are also provided in the ancillary files.We have also implemented them into a Fortran code, that will be published elsewhere.The timing of evaluating the 2loop amplitude square after summing all helicity configurations is around 0.3 second for a random phase space point on a single core with the 2.6 GHz Intel Core i7 CPU.
We have also carried out a few checks.First, the crossing symmetry relations of the amplitudes have been explicitly verified by plugging in the concrete analytic expressions.Our high-energy limit approaches the massless results in ref. [30] numerically, and our lowenergy limit approaches the results present in ref. [29] (also see appendix F).Finally, at the cross section level, our results agree perfectly with the pure numerical approach based on the local unitarity construction [65,66].

Conclusion
In this paper, we have presented the compact and analytic two-loop helicity amplitudes of the LbL process γγ → γγ keeping the full dependence of the internal fermion mass.The results remain simple enough even at two loops to be written in a few pages thanks to the numerous symmetries of the process.They could serve as a benchmark for developing novel multi-loop techniques, such as the recent local unitarity approach [65,66].The connection between the simplicity of our amplitudes and the string-derived Bern-Kosower formalism [67][68][69] needs to be explored in the future.The phenomenological applications and the comparisons with the local unitarity method will be reported in a companion paper [70].As an outlook, it is intriguing in the future to take into account the full electroweak corrections and next-to-next-to-leading order QCD and QED corrections.Their phenomenological relevance may show up in UPCs at the high-luminosity LHC by tagging proton(s) with forward proton spectrometers and future e + e − colliders.

A Kinematic regions and analytic continuation
In this section, we determine the kinematic regions that will be needed in two-loop amplitude calculations.In particular, we use the following parametrisation for the Mandelstam variables for all the polylogarithmic integrals: In the above parametrisation, there are multiple solutions of w and z in terms of x ij and x jk .We choose a unique solution that covers the whole 2 → 2 scattering phase space.First of all, since there is a symmetry in w and z, we assume ℜz ≥ ℜw if w, z ∈ R. we derive the regions of w and z that cover the required three kinematic regions: • Region II: s ij > 0, s jk < 0, s ik < 0 −→ this region can be further divided into 2 sub-regions. - can also be rewritten as 0 < ℜw < 1 and max(0, 1 • Region III: s ij < 0, s jk > 0, s ik < 0 −→ this region can also be further divided into 2 sub-regions. - The regions of (w, z) in the above cover the entire phase space only once.For w, z ∈ R, we need to know the sign of their imaginary parts in order to choose the correct Riemann sheet in GPLs.Since is sufficient for our purpose to take w → w − i0.For the square roots in the prefactors of the canonical basis eq.(D.1), we always evaluate them in the first Riemann sheet for simplicity.However, a minus sign stemming from the square roots should be introduced in the canonical master integrals f (2,ω) n (x ij , x jk , x ik ) with the corresponding square roots when the following condition is satisfied For instance, for the two square roots s ij (s ij − 4m 2 f ) and 27 , we should multiply the prefactor (−1) θ(x ij −4) (−1) θ(x ij )+θ(x jk ) in front of the expression consisting of G's as what is similarly done in eq.(4.4).We have taken care of these prefactors in the relations given in eqs.(3.13) and (3.14), and the remaining square roots in the amplitudes can be evaluated on the first Riemann sheet.The same applies to the one-loop master integrals f (1,ω) n (x ij , x jk , x ik ), which are defined in appendix C.

B Analytic boundary conditions
In this section, we present the analytic boundary conditions in the regions II and III defined in appendix A. For the region I, we simply choose x s = x t = x u = 0 as the boundary point, where all the master integrals are zero except the first one f (2,ω) 1 . In the regions II and III, we use the parametrisation x s = s/m 2 f = −(y s − 1) 2 /y s > 0 to rationalise the square roots in the boundary conditions.The solution is Given the analytic continuation of x s → x s + i0, we have y s → y s − i0 when x s > 4.

Figure 1 :
Figure 1: Representative Feynman diagrams of γγ → γγ at one loop (left) and two loops (rest).The internal wavy line can be either a gluon or a photon.