Numerical Evaluation of Six-Photon Amplitudes

TL;DR

This paper demonstrates an improved amplitude reduction method applied to six-photon scattering processes, including massive fermion loops, with strategies for rational part reconstruction and handling vanishing Gram-determinants.

Contribution

It introduces enhancements to the integrand-level amplitude reduction technique for complex scattering processes involving photons and fermion loops.

Findings

Successful computation of 2 photons -> 4 photons scattering amplitudes.

Effective strategies for reconstructing rational parts of one-loop amplitudes.

Methods to address vanishing Gram-determinants in calculations.

Abstract

We apply the recently proposed amplitude reduction at the integrand level method, to the computation of the scattering process 2 photons -> 4 photons, including the case of a massive fermion loop. We also present several improvements of the method, including a general strategy to reconstruct the rational part of any one-loop amplitude and the treatment of vanishing Gram-determinants.

Figures

Equations

$$A(\bar{q})=\frac{N(q)}{{\bar{D}}_0{\bar{D}}_1\cdots {\bar{D}}_{m-1}},{\bar{D}}_i=(\bar{q}+p_i{)}^2-m_i^2,p_0\ne 0,$$

$$N(q)$$

$$D_0=D_1=D_2=D_3=0,$$

$$N(q_0^{\pm })$$

$${\bar{D}}_i=D_i+{\tilde{q}}^2.$$

$$\frac{1}{{\bar{D}}_i}=\frac{{\bar{Z}}_i}{D_i},\mathrm{with}{\bar{Z}}_i\equiv \left(1-\frac{{\tilde{q}}^2}{{\bar{D}}_i}\right).$$

$$A(\bar{q})=\frac{N(q)}{D_0{D}_1\cdots D_{m-1}}{\bar{Z}}_0{\bar{Z}}_1\cdots {\bar{Z}}_{m-1}.$$

$$A(\bar{q})$$

$$\tilde{P}(q)$$

$$\tilde{a}(q;i_0)$$

$$\tilde{b}(q;i_0{i}_1)$$

$$\tilde{c}(q;i_0{i}_1{i}_2)$$

$$\tilde{d}(q;i_0{i}_1{i}_2{i}_3)$$

$${\tilde{c}}^{\mu \nu }(i_0{i}_1{i}_2;2)g_{\mu \nu }$$

$${\tilde{c}}^{\mu \nu \rho }(i_0{i}_1{i}_2;3)g_{\mu \nu }$$

$${\tilde{b}}^{\mu \nu }(i_0{i}_1;2)g_{\mu \nu }$$

$$I_{s;{\mu }_1\cdots {\mu }_r}^{(n;2\ell )}$$

$$\bar{D}(k_i)$$

$$I_{s;\mu \nu }^{(n;2(s-2))}.$$

$$I_{s;\mu \nu \rho }^{(n;2(s-2))},$$

$$I_s^{(n;2(s-1))},$$

$$I_s^{(n;2(s-1))}=-i{\pi }^2\frac{1}{s(s-1)}+\mathcal{O}(ϵ).$$

$$I_{s;\mu }^{(n;2(s-1))}\mathrm{with}2<s\le m-1.$$

$$I_{s;\mu }^{(n;2(s-1))}=i{\pi }^2\frac{1}{(s+1)s(s-1)}\sum _{j=1}^s(k_j{)}_{\mu }+\mathcal{O}(ϵ).$$

$$I_{s;\mu \nu }^{(n;2(s-1)}\mathrm{with}2<s\le m-1.$$

$$I_{s;\mu \nu }^{(n;2(s-1))}$$

$$I_s^{(n;2s)},\mathrm{with}1<s\le m-1.$$

$$I_s^{(n;2s)}$$

$$I_{s;\mu }^{(n;2s)},\mathrm{with}1<s\le m-1.$$

$$I_{s;\mu }^{(n;2s)}$$

$$A(\bar{q})=\frac{N(q)}{{\bar{D}}_0{\bar{D}}_1},$$

$$N(q)$$

$$\int d^{n}q\tilde{b}(q;01)=0.$$

$$k_1^2$$

Full text

NLO Computations, Hadronic Colliders, Standard Model, QCD

1 Introduction

In the last few years a big effort has been devoted by several authors to the problem of an efficient computation of one-loop corrections for multi-particle processes. This is a problem relevant for both LHC and ILC physics. In the case of QCD, the NLO six gluon amplitude has been recently obtained by two different groups [1], and, in the case of $e^{+}e^{-}$ collisions, complete EW calculations, involving 5-point [2] and 6-point [3] loop functions are available at the cross section level. The used techniques range from purely numerical methods to analytic ones, also including semi-numerical approaches. For analytical approaches, the main issue is reducing, using computer algebra, generic one-loop integrals into a minimal set of scalar integrals (and remaining pieces, the so called rational terms), mainly by tensor reduction [4, 5, 6, 7]. For multi-particle processes though this method becomes quite cumbersome because of the large number of terms generated and the appearance of numerical instabilities due to the zeros of Gram-determinants. On the other hand, several numerical or semi-numerical methods aim for a direct numerical computation of the tensor integrals [8]. Although purely numerical methods can in principle deal with any configuration of masses and also allow for a direct computation of both non-rational and rational terms, their applicability remains limited due to the high demand of computational resources and the non-existence of an efficient automation.

In a different approach, the one-loop amplitude rather than individual integrals are evaluated using the unitarity cut method [9], which relies on tree amplitudes and avoids the computation of Feynman diagrams. In another development, the four-dimensional unitarity cut method has been used for the calculation of QCD amplitudes [10], using twistor-based approaches [11]. Moreover, a generalization of the the unitarity cut method in $d$ dimensions, has been pursued recently [12].

Nevertheless, in practice, only the part of the amplitude proportional to the loop scalar functions can be obtained straightforwardly. The remaining piece, the rational part, should then be reconstructed either by using a direct computation based on Feynman diagrams [13, 14, 15] or by using a bootstrap approach [16]. Furthermore the complexity of the calculation increases away from massless theories.

In a recent paper [17], we proposed a reduction technique for arbitrary one-loop sub-amplitudes at the integrand level by exploiting numerically the set of kinematical equations for the integration momentum, that extend the quadruple, triple and double cuts used in the unitarity-cut method. The method requires a minimal information about the form of the one-loop (sub-)amplitude and therefore it is well suited for a numerical implementation. The method works for any set of internal and/or external masses, so that one is able to study the full electroweak model, without being limited to massless theories.

In this paper, we describe our experience with the first practical non-trivial implementation of such a method in the computation of a physical process: namely $2\gamma\to 4\gamma$, including massive fermion loops. For the massless case, there are a few results available in the literature. Analytical expressions were first presented by Mahlon [18] some time ago, however his results do not cover all possible helicity configurations. More recently the complete set of six-photon amplitudes was computed numerically by Nagy and Soper [19]. Very recently the same results were also obtained by Binoth et al. [20], that also provide compact analytical expressions.

In section 2, we recall the basics of our method and, in particular, we show how the knowledge of the rational terms can be inferred, with full generality, once the coefficients of the loop functions have been determined.

In section 3, we outline our solution to cure the numerical inaccuracies related to the appearance of zeros of Gram-determinants. We explicitly illustrate the case of 2-point amplitudes, that we had to implement to deal with the process at hand.

In section 4, we present our numerical results. For massless fermion loops we compare with available results. Moreover, since we are not limited to massless contributions, we also present, for the first time, results with massive fermion loops.

Finally, in the last section, we discuss our conclusions and future applications.

2 The method and the computation of the rational terms

The starting point of the method is the general expression for the integrand of a generic $m$-point one-loop (sub-)amplitude [17]

[TABLE]

where we use a bar to denote objects living in $n=~{}4+\epsilon$ dimensions, and $\bar{q}^{2}=q^{2}+\tilde{q}^{2}$ 444$\tilde{q}^{2}$ is $\epsilon$-dimensional and $(\tilde{q}\cdot q)=0$.. In the previous equation, $N(q)$ is the 4-dimensional part of the numerator function of the amplitude 555If needed, the $\epsilon$-dimensional part of the numerator should be treated separately, as explained in [21].. $N(q)$ depends on the $4$-dimensional denominators $D_{i}=({q}+p_{i})^{2}-m_{i}^{2}$ as follows

[TABLE]

Inserted back in Eq. (1), this expression simply states the multi-pole nature of any $m$-point one-loop amplitude, that, clearly, contains a pole for any propagator in the loop, thus one has terms ranging from 1 to $m$ poles. Notice that the term with no poles, namely that one proportional to $\tilde{P}(q)$ is polynomial and vanishes upon integration in dimensional regularization; therefore does not contribute to the amplitude, as it should be. The coefficients of the poles can be further split in two pieces. A piece that still depend on $q$ (the terms $\tilde{d},\tilde{c},\tilde{b},\tilde{a}$), that vanishes upon integration, and a piece that do not depend on q (the terms $d,c,b,a$). Such a separation is always possible, as shown in Ref. [17], and, with this choice, the latter set of coefficients is therefore immediately interpretable as the ensemble of the coefficients of all possible 4, 3, 2, 1-point one-loop functions contributing to the amplitude.

Once Eq. (2) is established, the task of computing the one-loop amplitude is then reduced to the algebraical problem of determining the coefficients $d,c,b,a$ by evaluating the function $N(q)$ a sufficient number of times, at different values of $q$, and then inverting the system. That can be achieved quite efficiently by singling out particular choices of $q$ such that, systematically, 4, 3, 2 or 1 among all possible denominators $D_{i}$ vanishes. Then the system of equations is solved iteratively. First one determines all possible 4-point functions, then the 3-point functions and so on. For example, calling $q_{0}^{\pm}$ the 2 (in general complex) solutions for which

[TABLE]

(there are 2 solutions because of the quadratic nature of the propagators) and since the functional form of $\tilde{d}(q;0123)$ is known, one directly finds the coefficient of the box diagram containing the above 4 denominators through the two simple equations

[TABLE]

This algorithm also works in the case of complex denominators, namely with complex masses. Notice that the described procedure can be performed at the amplitude level. One does not need to repeat the work for all Feynman diagrams, provided their sum is known: we just suppose to be able to compute $N(q)$ numerically.

As a further point notice that, since the terms $\tilde{d},\tilde{c},\tilde{b},\tilde{a}$ still depend on $q$, also the separation among terms in Eq. (2) is somehow arbitrary. Terms containing a different numbers of denominators can be shifted from one piece to the other in Eq. (2), by relaxing the requirement that the integral over the terms containing $q$ vanishes. This fact provides an handle to cure numerical instabilities occurring at exceptional phase-space points. In Section 3 we will show in detail such a mechanism at work for the 2-point part of the amplitude.

The described procedure works without any modification in 4 dimensions. However, even when starting from a perfectly finite tensor integral, the tensor reduction may eventually lead to integrals that need to be regularized. A typical example are the rank six 6-point functions contributing to the scattering $2\gamma\to 4\gamma$ we want to study. Such tensors are finite, but tensor reduction iteratively leads to rank $m$ $m$-point tensors with $1\leq m\leq 5$, that are ultraviolet divergent when $m\leq 4$. For this reason, we introduced, in Eq. (1), the $d$-dimensional denominators $\bar{D}_{i}$, that differs by an amount $\tilde{q}^{2}$ from their 4-dimensional counterparts

[TABLE]

The result of this is a mismatch in the cancellation of the $d$-dimensional denominators of Eq. (1) with the $4$-dimensional ones of Eq. (2). The rational part of the amplitude comes from such a lack of cancellation.

In [17] the problem of reconstructing this rational piece has been solved by looking at the implicit mass dependence in the coefficients $d,c,b,a$ of the one-loop functions. Such a method is adequate up to 4-point functions; for higher-point functions the dependence becomes too complicated to be used in practice. In addition, it requires the solution of further systems of linear equations, slowing down the whole computation. For those reasons, we suggest here a different method. One starts by rewriting any denominator appearing in Eq. (1) as follows

[TABLE]

This results in

[TABLE]

Then, by inserting Eq. (2) in Eq. (7), one obtains

[TABLE]

The rational part of the amplitude is then produced, after integrating over $d^{n}q$, by the $\tilde{q}^{2}$ dependence coming from the various $\bar{Z}_{i}$ in Eq. (8). It is easy to see what happens, for any value of $m$, by recalling the generic $q$ dependence of the spurious terms. In the renormalizable gauge one has [17]

[TABLE]

Eq. (2) simply states the fact that $\tilde{a}(q;i_{0})$ and $\tilde{d}(q;i_{0}i_{1}i_{2}i_{3})$ are at most linear in $(q+p_{i_{0}})$, $\tilde{b}(q;i_{0}i_{1})$ at most quadratic, and $\tilde{c}(q;i_{0}i_{1}i_{2})$ at most cubic. The tensors denoted by $(\cdots;1)$, $(\cdots;2)$ and $(\cdots;3)$ stand for the respective coefficients. We will also make use of the fact that, due to the explicit form of the spurious terms [17]

[TABLE]

The necessary integrals that arise, after a change of variable $q\to q-p_{i_{0}}$, are of the form

[TABLE]

where we used a notation introduced in [22] and $r\leq 3$. Such integrals (from now on called extra-integrals) have dimensionality ${\cal D}=2(1+\ell-s)+r$ and give a contribution ${\cal}O(1)$ only when ${\cal D}\geq 0$, otherwise are of ${\cal O}(\epsilon)$. This counting remains valid also in the presence of infrared and collinear divergences, as explained, for example, in Appendix B of [22] and in [14].

We also note that, since all $\bar{Z}_{i}$ are a-dimensional, the dimensionality ${\cal D}$ of the extra-integrals generated through Eq. (8) does not depend on $m$. We list, in the following, all possible contributions, collecting the computational details in Appendix A.

Contributions proportional to $d(i_{0}i_{1}i_{2}i_{3})$

In this case $r=0$. All extra-integrals are therefore scalars with ${\cal D}=-4$ and do not contribute.

Contributions proportional to $\tilde{d}^{\mu}(i_{0}i_{1}i_{2}i_{3};1)$

In this case $r=1$. All extra-integrals are therefore rank one tensors with ${\cal D}=-3$ and do not contribute.

Contributions proportional to $c(i_{0}i_{1}i_{2})$

In this case $r=0$ with ${\cal D}=-2$ and no contribution ${\cal O}(1)$ is developed.

Contributions proportional to $\tilde{c}^{\mu}(i_{0}i_{1}i_{2};1)$

Here $r=1$ and ${\cal D}=-1$. Therefore, once again, there is no contribution.

Contributions proportional to $\tilde{c}^{\mu\nu}(i_{0}i_{1}i_{2};2)$

Now $r=2$ with ${\cal D}=0$ and a finite contribution is in principle expected, generated by extra-integrals of the type

[TABLE]

Nevertheless, such contribution is proportional to $g_{\mu\nu}$ [22]. Therefore, due to Eq. (2), it vanishes.

Contributions proportional to $\tilde{c}^{\mu\nu\rho}(i_{0}i_{1}i_{2};3)$

Now $r=3$ and ${\cal D}=1$. The contributing extra-integrals are of the type

[TABLE]

and one easily proves that the contributions ${\cal O}(1)$ are always proportional to $g_{\mu\nu}$ or $g_{\mu\rho}$ or $g_{\nu\rho}$. Therefore, thanks again to Eq. (2), they vanish.

Contributions proportional to $b(i_{0}i_{1})$

Those are the first non vanishing contributions. The relevant extra-integrals have $r=0$ and ${\cal D}=0$

[TABLE]

with $2<s\leq m-1$. They have been computed, for generic values of $s$, in [22] (see also Appendix A)

[TABLE]

Contributions proportional to $\tilde{b}^{\mu}(i_{0}i_{1};1)$

In this case the relevant extra-integrals are 4-vectors with ${\cal D}=1$

[TABLE]

A computation for generic values of $s$ gives

[TABLE]

Contributions proportional to $\tilde{b}^{\mu\nu}(i_{0}i_{1};2)$

The relevant extra-integrals are now rank two tensor with ${\cal D}=2$

[TABLE]

They read

[TABLE]

The $g_{\mu\nu}$ part is never needed because $\tilde{b}^{\mu\nu}(i_{0}i_{1};2)\,g_{\mu\nu}=0$, according to Eq. (2).

Contributions proportional to $a(i_{0})$

They involve scalar extra-integrals with ${\cal D}=2$

[TABLE]

One computes

[TABLE]

Contributions proportional to $\tilde{a}^{\mu}(i_{0};1)$

This last category involves extra-integrals with $r=1$ and ${\cal D}=3$

[TABLE]

One obtains

[TABLE]

To conclude, the set of the five formulas in Eqs. (14)-(18) allows one to compute the rational part of any one-loop $m$-point (sub-)amplitude, once all the coefficients of Eq. (2) have been reconstructed.

3 Dealing with numerical instabilities

In this section we show how to handle, in the framework of the method illustrated in the previous section, the simplest numerical instability appearing in any one-loop calculation, namely that one related to the tensor reduction of 2-point amplitudes in the limit of vanishing Gram-determinant 666In this case the Gram-determinant is simply the square of the difference between the momenta of the two denominators.. This situation is simple enough to allow an easy description, but the outlined strategy is general and not restricted to the 2-point case.

We start from the integrand of a generic 2-point amplitude written in the form

[TABLE]

in which we suppose $N(q)$ at most quadratic in $q$. Our purpose is dealing with the situation in which $k_{1}^{2}\equiv(p_{1}-p_{0})^{2}=0$ exactly (that always occur in processes with massless external particles), as well as to set up an algorithm to write down approximations around this case with arbitrary precision.

According to Eq. (2), we can write an expansion for $N(q)$ as follows:

[TABLE]

If the Gram-determinant of the 2-point function is small, the reduction method introduced in [17] cannot be applied, because the solution for which $D_{0}=D_{1}=0$, needed to determine the coefficients $b$ and $\tilde{b}$, becomes singular777Such a solution goes like $1/k_{1}^{2}$., in the limit of $k_{1}^{2}\to 0$, when adding the requirement

[TABLE]

Then, we must consider two separate cases:

[TABLE]

The former situation may occur because of the Minkowskian metric, while the latter takes place at the edges of the phase-space, where some momenta become collinear. In the first case one can still find a solution for which $D_{0}=D_{1}=0$ by relaxing the further requirement of Eq. (21). Such a solution is given in Appendix B and goes like $1/(k_{1}.v)$, where $v$ is an arbitrary massless 4-vector, therefore is never singular in the first case of Eq. (3). The price to pay is that new non zero integrals appear of the type 888Since $v^{2}=0$ they still fulfill the third one of Eqs. (2), therefore, even in this case, terms ${\cal O}(g_{\mu\nu})$ can be neglected in Eq. (16).

[TABLE]

What has been achieved with this new basis is then moving part of the 1-point functions to the 2-point sector, in such a way that combinations well behaved in the limit $k_{1}^{2}\to 0$ appear. The fact that solutions exist to the condition $D_{0}=D_{1}=0$, still allows one to find the coefficients of such integrals (together with all the others). This solves the first part of the problem, namely reconstructing $N(q)$ without knowing explicitly its analytic structure, but one is left with the problem of computing the new 2-point integrals. In the following, we present our method to determine them at any desired order in $k_{1}^{2}$. Let us first consider the case $j=1$ in Eq. (23). The contribution ${\cal O}(1)$ can be easily obtained from the observation that999From now on, we shift the integration variable: $\bar{q}\to\bar{q}-p_{0}$. The definition of the new resulting denominators is given in Eq. (2).

[TABLE]

as it is evident by performing a tensor decomposition. On the other hand, by reconstructing denominators, one obtains

[TABLE]

with

[TABLE]

Eq. (25), inserted in Eq. (24) gives the desired expansion in terms of loop functions with less points but higher rank, in agreement with well know results [23, 24]

[TABLE]

Expansions at arbitrary orders in $k_{1}^{2}$ can be obtained in an analogous way from the two following equations:

[TABLE]

To deal with the case $j=2$ in Eq. (23) one starts instead from the equation

[TABLE]

This procedure breaks down when the quantity $f$ vanishes. In this case a double expansion in $k_{1}^{2}$ and $f$ can still be found in terms of derivatives of one-loop scalar functions. We illustrate the procedure for the case $j=1$ of Eq. (23). Our starting point is now the equation

[TABLE]

By multiplying and dividing by $\bar{D}(k_{0})$ one obtains

[TABLE]

Applying once more Eq. (30) to the last integral gives

[TABLE]

Since the last integral in the previous equation is ${\cal O}(k_{1}^{2})$, the final result reads

[TABLE]

In a similar fashion, expansions at any order can be obtained.

We now turn to the second case of Eq. (3), namely $k_{1}^{\mu}\to 0$. In this case no solution can be found to the double cut equation

[TABLE]

The reason is that now $D(k_{1})$ and $D(k_{0})$ are no longer independent:

[TABLE]

and clearly no $q$ exists such that the two denominators can be simultaneously zero. Notice that this also implies that one cannot fit separately the coefficients of the 2-point and 1-point functions in Eq. (20). This results is a singularity $1/(k_{1}\cdot v)$ in the system given of Appendix B and we should change our strategy. We than go back to Eq. (19) and split the amplitude from the beginning by multiplying it by

[TABLE]

resulting to

[TABLE]

with

[TABLE]

Now the two amplitudes $A^{(1,2)}$ can be reconstructed separately, without any problem of vanishing Gram-determinant. Notice also that corrections at orders higher than ${\cal O}(k_{1})$ are perfectly calculable by inserting again Eq. (36) in the term ${\cal O}(k_{1})$ of Eq. (37).

Once again, when $f\to 0$, double expansions in $k_{1}$ and $f$ can be obtained involving derivatives of scalar loop functions by using Eq. (30). For example, at the zeroth order in $k_{1}$ and at the first one in $f$, one gets

[TABLE]

This last case exhausts all possibilities.

The same techniques can be applied for higher-point functions. For example, in the case of a 3-point function, instead of $k_{1}$, one introduces the 4-vector

[TABLE]

with the properties

[TABLE]

where $\Delta(k_{1},k_{2})$ is the Gram-determinant of the two momenta $k_{1}$ and $k_{2}$. Then, instead of Eq. (24) one has, for example,

[TABLE]

As before, $\Delta(k_{1},k_{2})$ can vanish either because $s^{2}=0$ or $s^{\mu}=0$ and the two cases should be treated separately.

4 Results and comparisons

We started by checking our implementation of the rational terms. For 4-point functions up to rank four, we reproduced the results obtained with the alternative technique illustrated in [17]. Furthermore, we reproduced the rational part of the full $2\gamma\to 2\gamma$ amplitude given in [25]. We also checked with an independent calculation [26] the rational terms coming from all of the 6-point tensors up to rank six. Finally, we computed the rational piece of the whole $2\gamma\to 4\gamma$ amplitude by summing up all 120 contributing Feynman diagrams and finding zero, as it should be [14].

As a first test on full amplitudes, we checked our method by reproducing the contribution of a fermion loop to the $2\gamma\to 2\gamma$ process. This result is presented in Eqs. (A.18)-(A.20) of Ref. [25], for all possible helicity configurations. We are in perfect agreement with the analytic result, in both massless and massive cases.

The next step was the computation of the $2\gamma\to 4\gamma$ amplitude with zero internal mass101010We thank Andre van Hameren for providing us with his program to compute massless one-loop scalar integrals., finding the results given in Fig. 1 and Fig. 2. It should be mentioned that our results are obtained algebraically, so there is no integration error involved.

In Fig. 1, we reproduce the results presented by Nagy and Soper [19] and very recently also by Binoth et al.[20]. We employ the same values of the external momenta as in Fig. 5 of Ref. [19], namely the following selection of final state three-momenta $\{\vec{p}_{3},\vec{p}_{4},\vec{p}_{5},\vec{p}_{6}\}$:

[TABLE]

After choosing the incoming photons such that they have momenta $\vec{p}_{1}$ and $\vec{p}_{2}$ along the $z$-axis, we present in the plot the amplitude obtained by rotating the final states of angle $\theta$ about the $y$-axis. This is done for both helicity configurations $[++----]$ and $[+--++-]$. In the same plot also appears the analytic results for the configuration $[++----]$ obtained by Mahlon [18]. In Fig. 2, we use a different set of external momenta. Starting from the following choice of $\{\vec{p}_{3},\vec{p}_{4},\vec{p}_{5},\vec{p}_{6}\}$:

[TABLE]

we proceed as in the previous case. The results for the amplitudes are plotted in Fig. 2 for the helicity configurations $[++----]$ and $[++--+-]$. It is known that the six-photons amplitude vanish for the helicity configurations $[++++++]$ and $[+++++-]$, we checked this result for both choices of the external momenta. Finally, using the external momenta of Eq. (4), we computed the amplitude introducing a non-zero mass $m_{f}$ for the fermions in the loop 111111We used here the scalar one-loop functions provided by FF[27].. The results are plotted in Fig. 3, for the three cases $m_{f}=0.5$ GeV, $m_{f}=4.5$ GeV and $m_{f}=12$ GeV.

The code we prepared for producing the results presented in this section is written in FORTRAN 90. Even if we did not spend too much effort in optimizations, it can compute about 3 phase-space points per second, when working in double precision. All figures in this section are actually produced by using double precision, but, to perform a realistic integration, we still need quadruple precision, that slows down the speed by about a factor 60. We are working in implementing the expansions presented in the previous section with the aim of being able to perform a stable integration over the full phase space, that is a “proof of concept” for any method.

5 Conclusions

Computing the massless QED amplitude for the reaction $2\gamma\to 4\gamma$, although still unobserved experimentally, is a very good exercise for checking new methods to calculate one-loop virtual corrections. Such a process posses all complications typical of any multi-leg final state, for example a non trivial tensorial structure, but also keeps, at the same time, enough simplicity such that compact analytical formulas can still be used as a benchmark. However, it is oversimplified in two respects. Firstly, the amplitude it is completely massless. Secondly, the amplitude is cut constructible, namely does not contain any rational part.

In the most general case of one-loop calculations, the presence of both internal and external masses prevents from obtaining compact analytical expressions. Then one has to rely on other computational techniques. For example, it is known that cut-constructible amplitudes can be obtained through recursion relations. But, even then, the presence of rational parts usually requires a separate work.

For such reasons, it would be highly advisable to have a method not restricted to massless theories, in which moreover both cut-constructible and rational parts can be treated at the same time. Such a method has been introduced recently in Ref. [17] and, in this paper, we applied it to the computation of the six-photon amplitude in QED, giving also results for the case with massive fermions in the loop. We also showed in detail how the rational part of any $m$-point one-loop amplitude is intimately connected with the form of the integrand of the amplitude. Once this integrand is numerically computable, cut-constructible and rational terms are easily obtained, at the same time, by solving the same system of linear equations. This is a peculiar property of our method, that we tested in the actual computation of the six-photon amplitude. In practice, we did not use the additional information on its cut-constructibility and verified only a-posteriori that the intermediate rational parts, coming from all pieces separately, drop out in the final sum.

Finally, we presented all relevant formulas needed to infer the rational parts from the integrand of any $m$-point loop functions, in the renormalizable gauges.

In addition, we presented, by analyzing in detail the 2-point case, an idea to cure the numerical instabilities occurring at exceptional phase-space points, outlining a possible way to build up expansions around the zeroes of the Gram-determinants.

Having been able to apply our method to the computation of the massive six-photon amplitude, we are confident that our method can be successfully used for a systematic and efficient computation of one-loop amplitudes relevant at LHC and ILC.

Acknowledgments

We thank Andre van Hameren for numerical comparisons and Zoltan Nagy and Pierpaolo Mastrolia for interesting discussions. G.O. acknowledges the financial support of the ToK Program “ALGOTOOLS” (MTKD-CT-2004-014319). C.G.P’s and R.P.’s research was partially supported by the RTN European Programme MRTN-CT-2006-035505 (HEPTOOLS, Tools and Precision Calculations for Physics Discoveries at Colliders). The research of R.P. was also supported in part by MIUR under contract 2006020509_004.

Appendices

Appendix A Computing the extra-integrals

In this appendix, we compute the extra-integrals listed in Section 2. Since a contribution ${\cal O}(1)$ can only develop for non-negative dimensionality ${\cal D}$, the integrand in the Feynman parameter integral is always polynomial. First we decompose the integration as follows

[TABLE]

then, after using Feynman parametrization and performing first the integral over $d^{\epsilon}\mu$ and then that one over $d^{4}q$, one derives, for the extra-integrals of Eqs. (14)-(18)

[TABLE]

where

[TABLE]

In the following, we compute, as an illustrative example, the first three integrals of Eq. (A). The remaining two can be obtained analogously. We start by changing the integration variables as follows:

[TABLE]

so that

[TABLE]

from which one trivially obtains the first integral

[TABLE]

To compute the second integral an integration over $(P_{s})_{\mu}$ in needed. Since the integrand is symmetric when interchanging all $k_{i}$, we concentrate on the coefficient of, say, $k_{1}$. Since

[TABLE]

the final result reads

[TABLE]

To compute the third integral we need to integrate over the product $(P_{s})_{\mu}(P_{s})_{\nu}$. Once again, given the symmetry of the problem, we can focus on the two contributions proportional to $k_{1\mu}k_{1\nu}$ and $k_{1\mu}k_{2\nu}$. The first one gives

[TABLE]

and the second reads

[TABLE]

Summing up all of the possibilities one obtains

[TABLE]

Appendix B The general basis for the 2-point functions

In this appendix, we solve the problem of reconstructing the coefficients of the 2-point part of the integrand of any amplitude

[TABLE]

by assuming $N(q)$ at most quadratic in $q$ and $k_{1}\equiv(p_{1}-p_{0})\neq 0$. In particular also the case of vanishing $k_{1}^{2}$ is included. First, we introduce a massless arbitrary 4-vector $v$, such that $(v\cdot k_{1})\neq 0$, that we use to rewrite $k_{1}$ in terms of two massless 4-vectors (we also take $\ell^{2}=0$)

[TABLE]

giving

[TABLE]

Then, we introduce two additional independent massless 4-vectors $\ell_{7,8}$ defined as

[TABLE]

for which one finds

[TABLE]

and we decompose $q^{\mu}+p_{0}^{\mu}$ in the basis of $k_{1}$, $v$, $\ell_{7}$ and $\ell_{8}$

[TABLE]

so that $N(q)$ takes the form

[TABLE]

Notice that, because of the identity

[TABLE]

any term proportional to $[(q+p_{0})\cdot k_{1}]$ either contributes to the constant term $b$ or it is included in the terms ${\cal O}(D_{0,1})$ we are neglecting 121212We suppose to determine them at a later stage of the calculation.. The same happens for the combination $[(q+p_{0})\cdot\ell_{7}][(q+p_{0})\cdot\ell_{8}]$.

To be able to determine all of the coefficients appearing in Eq. (64), disentangling completely the contributions ${\cal O}(D_{0,1})$, we look for a $q$ that fulfill the requirement

[TABLE]

For a $q$ written as in Eq. (63) this implies the system

[TABLE]

where

[TABLE]

It is convenient to introduce two classes of solutions. In the first class, that we call $q^{\pm}_{yk}$, we take $y$ fixed and choose $y_{7}=\pm e^{i\pi/k}$. In the second class, that we call $q^{\prime\pm}_{yk}$, we take $y$ fixed but choose $y_{8}=\pm e^{i\pi/k}$. The coefficients $b$, $\tilde{b}_{11}$, $\tilde{b}_{21}$, $\tilde{b}_{12}$ and $\tilde{b}_{22}$ can be obtained by evaluating Eq. (64) at the values

[TABLE]

or

[TABLE]

In the first case, the coefficients read

[TABLE]

with

[TABLE]

and where

[TABLE]

In the second case, one obtains instead

[TABLE]

with

[TABLE]

and where

[TABLE]

The reason why we have chosen two sets of solutions is that, in some special kinematical configurations, $F_{0}$ can vanish. Therefore, numerical stable solutions are obtained by taking $\tilde{b}_{21}$ and $\tilde{b}_{22}$ from Eq. (71), and $\tilde{b}_{11}$ and $\tilde{b}_{12}$ from Eq. (74), while $b$ is well defined in both cases.

The coefficients $\hat{b}_{0}$ and $\hat{b}_{00}$ can be determined, in terms of additional solutions of the kind $q^{\pm}_{\lambda 1}$ and $q^{\pm}_{\sigma 1}$, by defining the combinations

[TABLE]

as the two solutions of the system

[TABLE]

The determinant of the matrix above is always different form zero, for non vanishing $\lambda$ and $\sigma$, when $\sigma\neq\lambda$, so that numerical inaccuracies never occur.

Finally, the two last coefficients $\tilde{b}_{01}$ and $\tilde{b}_{02}$ are determined, in terms of $q^{+}_{\lambda k}$ and $q^{\prime+}_{\sigma k}$, as solutions of the system

[TABLE]

where

[TABLE]

Once again one verifies that when, for example, $k=3$ the system never becomes singular.