arXiv:0704.1271v2  [hep-ph]  23 Jul 2007
Preprint typeset in JHEP style - PAPER VERSION
Numerical Evaluation of Six-Photon Amplitudes
Giovanni Ossola∗, Costas G. Papadopoulos†
Institute of Nuclear Physics, NCSR ”DEMOKRITOS”, 15310 Athens, Greece.
Roberto Pittau‡
Dipartimento di Fisica Teorica, Univ. di Torino and INFN, sez. di Torino, Italy.
Abstract: We apply the recently proposed amplitude reduction at the integrand level
method, to the computation of the scattering process 2γ → 4γ, 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.
Keywords: NLO Computations, Hadronic Colliders, Standard Model, QCD.
∗e-mail: ossola@inp.demokritos.gr
†e-mail: costas.papadopoulos@cern.ch
‡e-mail: roberto.pittau@to.infn.it
http://arxiv.org/abs/0704.1271v2
Contents
1. Introduction 1
2. The method and the computation of the rational terms 2
3. Dealing with numerical instabilities 8
4. Results and comparisons 11
5. Conclusions 14
A. Computing the extra-integrals 15
B. The general basis for the 2-point functions 17
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 an-
alytic 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 re-
duction [4–7]. For multi-particle processes though this method becomes quite cumbersome
because of the large number of terms generated and the appearance of numerical insta-
bilities 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
– 1 –
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–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 imple-
mentation of such a method in the computation of a physical process: namely 2γ → 4γ,
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]
A(q̄) =
D̄0D̄1 · · · D̄m−1
, D̄i = (q̄ + pi)
2 −m2i , p0 6= 0 , (2.1)
where we use a bar to denote objects living in n = 4 + ǫ dimensions, and q̄2 = q2 + q̃2
1. In the previous equation, N(q) is the 4-dimensional part of the numerator function of
1q̃2 is ǫ-dimensional and (q̃ · q) = 0.
– 2 –
the amplitude 2. N(q) depends on the 4-dimensional denominators Di = (q+ pi)
2 −m2i as
follows
N(q) =
i0<i1<i2<i3
d(i0i1i2i3) + d̃(q; i0i1i2i3)
i 6=i0,i1,i2,i3
i0<i1<i2
[c(i0i1i2) + c̃(q; i0i1i2)]
i 6=i0,i1,i2
i0<i1
b(i0i1) + b̃(q; i0i1)
i 6=i0,i1
[a(i0) + ã(q; i0)]
i 6=i0
+ P̃ (q)
Di . (2.2)
Inserted back in Eq. (2.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 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 d̃, c̃, b̃, 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.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 Di 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 the 2 (in general complex)
solutions for which
D0 = D1 = D2 = D3 = 0 , (2.3)
(there are 2 solutions because of the quadratic nature of the propagators) and since the
functional form of d̃(q; 0123) is known, one directly finds the coefficient of the box diagram
containing the above 4 denominators through the two simple equations
N(q±0 ) = [d(0123) + d̃(q
0 ; 0123)]
i 6=0,1,2,3
0 ) . (2.4)
2If needed, the ǫ-dimensional part of the numerator should be treated separately, as explained in [21].
– 3 –
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 d̃, c̃, b̃, ã still depend on q, also the
separation among terms in Eq. (2.2) is somehow arbitrary. Terms containing a different
numbers of denominators can be shifted from one piece to the other in Eq. (2.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 even-
tually lead to integrals that need to be regularized. A typical example are the rank six
6-point functions contributing to the scattering 2γ → 4γ we want to study. Such tensors
are finite, but tensor reduction iteratively leads to rank m m-point tensors with 1 ≤ m ≤ 5,
that are ultraviolet divergent when m ≤ 4. For this reason, we introduced, in Eq. (2.1),
the d-dimensional denominators D̄i, that differs by an amount q̃
2 from their 4-dimensional
counterparts
D̄i = Di + q̃
2 . (2.5)
The result of this is a mismatch in the cancellation of the d-dimensional denominators of
Eq. (2.1) with the 4-dimensional ones of Eq. (2.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. (2.1) as follows
, with Z̄i ≡
. (2.6)
This results in
A(q̄) =
D0D1 · · ·Dm−1
Z̄0Z̄1 · · · Z̄m−1 . (2.7)
Then, by inserting Eq. (2.2) in Eq. (2.7), one obtains
A(q̄) =
i0<i1<i2<i3
d(i0i1i2i3) + d̃(q; i0i1i2i3)
D̄i0D̄i1D̄i2D̄i3
i 6=i0,i1,i2,i3
– 4 –
i0<i1<i2
c(i0i1i2) + c̃(q; i0i1i2)
D̄i0D̄i1D̄i2
i 6=i0,i1,i2
i0<i1
b(i0i1) + b̃(q; i0i1)
D̄i0D̄i1
i 6=i0,i1
a(i0) + ã(q; i0)
i 6=i0
+ P̃ (q)
Z̄i . (2.8)
The rational part of the amplitude is then produced, after integrating over dnq, by the
q̃2 dependence coming from the various Z̄i in Eq. (2.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]
P̃ (q) = 0 ,
ã(q; i0) = ã
µ(i0; 1)(q + pi0)µ ,
b̃(q; i0i1) = b̃
µ(i0i1; 1)(q + pi0)µ + b̃
µν(i0i1; 2)(q + pi0)µ(q + pi0)ν ,
c̃(q; i0i1i2) = c̃
µ(i0i1i2; 1)(q + pi0)µ + c̃
µν(i0i1i1; 2)(q + pi0)µ(q + pi0)ν ,
+ c̃µνρ(i0i1i1; 3)(q + pi0)µ(q + pi0)ν(q + pi0)ρ ,
d̃(q; i0i1i2i3) = d̃
µ(i0i1i2i3; 1)(q + pi0)µ . (2.9)
Eq. (2.9) simply states the fact that ã(q; i0) and d̃(q; i0i1i2i3) are at most linear in (q+pi0),
b̃(q; i0i1) at most quadratic, and c̃(q; i0i1i2) at most cubic. The tensors denoted by (· · · ; 1),
(· · · ; 2) and (· · · ; 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]
c̃µν(i0i1i2; 2) gµν = 0 ,
c̃µνρ(i0i1i2; 3) gµν = c̃
µνρ(i0i1i2; 3) gµρ = c̃
µνρ(i0i1i2; 3) gνρ = 0 and
b̃µν(i0i1; 2) gµν = 0 . (2.10)
The necessary integrals that arise, after a change of variable q → q − pi0 , are of the form
(n;2ℓ)
s;µ1···µr ≡
dnq q̃2ℓ
qµ1 · · · qµr
D̄(k0) · · · D̄(ks)
, with
D̄(ki) ≡ (q̄ + ki)
2 −m2i , ki ≡ pi − p0 (k0 = 0) , (2.11)
where we used a notation introduced in [22] and r ≤ 3. Such integrals (from now on called
extra-integrals) have dimensionality D = 2(1+ ℓ− s)+ r and give a contribution O(1) only
when D ≥ 0, otherwise are of O(ǫ). 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].
– 5 –
We also note that, since all Z̄i are a-dimensional, the dimensionality D of the extra-
integrals generated through Eq. (2.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(i0i1i2i3)
In this case r = 0. All extra-integrals are therefore scalars with D = −4 and do not
contribute.
Contributions proportional to d̃µ(i0i1i2i3; 1)
In this case r = 1. All extra-integrals are therefore rank one tensors with D = −3 and do
not contribute.
Contributions proportional to c(i0i1i2)
In this case r = 0 with D = −2 and no contribution O(1) is developed.
Contributions proportional to c̃µ(i0i1i2; 1)
Here r = 1 and D = −1. Therefore, once again, there is no contribution.
Contributions proportional to c̃µν(i0i1i2; 2)
Now r = 2 with D = 0 and a finite contribution is in principle expected, generated by
extra-integrals of the type
I(n;2(s−2))s;µν . (2.12)
Nevertheless, such contribution is proportional to gµν [22]. Therefore, due to Eq. (2.10), it
vanishes.
Contributions proportional to c̃µνρ(i0i1i2; 3)
Now r = 3 and D = 1. The contributing extra-integrals are of the type
I(n;2(s−2))s;µνρ , (2.13)
and one easily proves that the contributions O(1) are always proportional to gµν or gµρ or
gνρ. Therefore, thanks again to Eq. (2.10), they vanish.
Contributions proportional to b(i0i1)
Those are the first non vanishing contributions. The relevant extra-integrals have r = 0
and D = 0
I(n;2(s−1))s ,
with 2 < s ≤ m − 1. They have been computed, for generic values of s, in [22] (see also
Appendix A)
I(n;2(s−1))s = −iπ
s(s− 1)
+O(ǫ) . (2.14)
– 6 –
Contributions proportional to b̃µ(i0i1; 1)
In this case the relevant extra-integrals are 4-vectors with D = 1
I(n;2(s−1))s;µ with 2 < s ≤ m− 1 .
A computation for generic values of s gives
I(n;2(s−1))s;µ = iπ
(s + 1)s(s − 1)
(kj)µ +O(ǫ) . (2.15)
Contributions proportional to b̃µν(i0i1; 2)
The relevant extra-integrals are now rank two tensor with D = 2
I(n;2(s−1)s;µν with 2 < s ≤ m− 1 .
They read
I(n;2(s−1))s;µν = −2iπ
(s+ 2)(s + 1)s(s− 1)
(kj)µ(kj)ν +
i 6=j
(kj)µ(ki)ν
+ O(gµν) +O(ǫ) . (2.16)
The gµν part is never needed because b̃
µν(i0i1; 2) gµν = 0, according to Eq. (2.10).
Contributions proportional to a(i0)
They involve scalar extra-integrals with D = 2
I(n;2s)s , with 1 < s ≤ m− 1 .
One computes
I(n;2s)s = −2iπ
(s+ 2)(s + 1)s
k2j +
i 6=j
(kj · ki) +
(m2j − k
+ O(ǫ) . (2.17)
Contributions proportional to ãµ(i0; 1)
This last category involves extra-integrals with r = 1 and D = 3
I(n;2s)s;µ , with 1 < s ≤ m− 1 .
One obtains
I(n;2s)s;µ = iπ
(s+ 3)(s + 2)(s + 1)s
k2j (kj)µ + 2
i 6=j
k2j (ki)µ + 2(kj · ki)(kj)µ
– 7 –
i 6=j
ℓ 6=i
(kj · ki)(kℓ)µ + (s+ 3)
(m2j − k
j )(kj)µ
i 6=j
(m2j − k
j )(ki)µ
+O(ǫ) . (2.18)
To conclude, the set of the five formulas in Eqs. (2.14)-(2.18) allows one to compute the
rational part of any one-loop m-point (sub-)amplitude, once all the coefficients of Eq. (2.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 3. 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
A(q̄) =
D̄0D̄1
, (3.1)
in which we supposeN(q) at most quadratic in q. Our purpose is dealing with the situation
in which k21 ≡ (p1−p0)
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.2), we can write an expansion for N(q) as follows:
N(q) = [b(01) + b̃(q; 01)] + [a(0) + ã(q; 0)]D1 + [a(1) + ã(q; 1)]D0 . (3.2)
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 D0 = D1 = 0, needed to determine
the coefficients b and b̃, becomes singular4, in the limit of k21 → 0, when adding the
requirement
dnq b̃(q; 01) = 0 . (3.3)
Then, we must consider two separate cases:
k21 → 0 , but k
1 6= 0 ,
k21 → 0 , because k
1 = 0 . (3.4)
3In this case the Gram-determinant is simply the square of the difference between the momenta of the
two denominators.
4Such a solution goes like 1/k21 .
– 8 –
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 D0 = D1 = 0 by relaxing the further requirement
of Eq. (3.3). Such a solution is given in Appendix B and goes like 1/(k1.v), where v is an
arbitrary massless 4-vector, therefore is never singular in the first case of Eq. (3.4). The
price to pay is that new non zero integrals appear of the type 5
[(q + p0) · v]
D̄0D̄1
with j = 1, 2 and v2 = 0 . (3.5)
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 k21 → 0
appear. The fact that solutions exist to the condition D0 = D1 = 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 k21. Let us
first consider the case j = 1 in Eq. (3.5). The contribution O(1) can be easily obtained
from the observation that6
(q · v)(q · k1)
D̄(k0)D̄(k1)
= O(k21) , (3.6)
as it is evident by performing a tensor decomposition. On the other hand, by reconstructing
denominators, one obtains
(q · k1)
D̄(k1)− D̄(k0)
(q · k1) +
, (3.7)
f = m21 − k
0 . (3.8)
Eq. (3.7), inserted in Eq. (3.6) gives the desired expansion in terms of loop functions with
less points but higher rank, in agreement with well know results [23,24]
(q · v)
D̄(k0)D̄(k1)
dnq (q · v)
D̄(k1)
D̄(k0)
2(q · k1)
+O(k21) .(3.9)
Expansions at arbitrary orders in k21 can be obtained in an analogous way from the two
following equations:
(q · k1)
D̄(k1)− D̄(k0)
i+j=p−1
(q · k1)
(q · v)(q · k1)
D̄(k0)D̄(k1)
= O(k
1 ) . (3.10)
5Since v2 = 0 they still fulfill the third one of Eqs. (2.10), therefore, even in this case, terms O(gµν) can
be neglected in Eq. (2.16).
6From now on, we shift the integration variable: q̄ → q̄ − p0. The definition of the new resulting
denominators is given in Eq. (2.11).
– 9 –
To deal with the case j = 2 in Eq. (3.5) one starts instead from the equation
(q · v)2(q · k1)
D̄(k0)D̄(k1)
= O(k
1 ) . (3.11)
This procedure breaks down when the quantity f vanishes. In this case a double expansion
in k21 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. (3.5). Our starting point is now the
equation
D̄(k0) = D̄(k1)− 2(q · k1) + f . (3.12)
By multiplying and dividing by D̄(k0) one obtains
(q · v)
D̄(k0)D̄(k1)
(q · v)
D̄(k0)2D̄(k1)
D̄(k1)− 2(q · k1) + f
(q · v)
D̄(k0)2
(q · v)(q · k1)
D̄(k0)2D̄(k1)
+O(f) . (3.13)
Applying once more Eq. (3.12) to the last integral gives
(q · v)(q · k1)
D̄(k0)
2D̄(k1)
(q · v)(q · k1)
D̄(k0)
3D̄(k1)
D̄(k1)− 2(q · k1) + f
(q · v)(q · k1)
D̄(k0)3
(q · v)(q · k1)
D̄(k0)3D̄(k1)
+O(f) . (3.14)
Since the last integral in the previous equation is O(k21), the final result reads
(q · v)
D̄(k0)D̄(k1)
(q · v)
D̄(k0)2
(q · v)(q · k1)
D̄(k0)3
+O(k21) +O(f) .
(3.15)
In a similar fashion, expansions at any order can be obtained.
We now turn to the second case of Eq. (3.4), namely k
1 → 0. In this case no solution
can be found to the double cut equation
D(k0) = D(k1) = 0 . (3.16)
The reason is that now D(k1) and D(k0) are no longer independent:
D(k0) = D(k1) + f +O(k1) , (3.17)
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. (3.2). This results is a singularity 1/(k1 · v) in the system given of
Appendix B and we should change our strategy. We than go back to Eq. (3.1) and split
the amplitude from the beginning by multiplying it by
D̄(k0)− D̄(k1)
2(q · k1)
, (3.18)
– 10 –
resulting to
A(q̄) = A(1)(q̄) +A(2)(q̄) +O(k1) , (3.19)
A(1)(q̄) =
D̄(k1)
, A(2)(q̄) = −
D̄(k0)
. (3.20)
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 O(k1) are
perfectly calculable by inserting again Eq. (3.18) in the term O(k1) of Eq. (3.19).
Once again, when f → 0, double expansions in k1 and f can be obtained involving
derivatives of scalar loop functions by using Eq. (3.12). For example, at the zeroth order
in k1 and at the first one in f , one gets
A(q) =
D̄(k0)D̄(k1)
D̄(k0)2D̄(k1)
D̄(k1)− 2(q · k1) + f
D̄(k0)
D̄(k0)
3D̄(k1)
D̄(k1)− 2(q · k1) + f
+O(k1)
D̄(k0)
D̄(k0)
+O(k1) +O(f
2) . (3.21)
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 k1, one introduces the 4-vector
sµ = det
(k2 · k1) (k2 · k2)
, (3.22)
with the properties
s · k2 = 0 , s
2 ∝ ∆(k1, k2) , (k1 · s) ∝ ∆(k1, k2) , (3.23)
where ∆(k1, k2) is the Gram-determinant of the two momenta k1 and k2. Then, instead of
Eq. (3.6) one has, for example,
(q · v)(q · s)2
D̄(k0)D̄(k1)D̄(k2)
= O(∆(k1, k2)) . (3.24)
As before, ∆(k1, k2) can vanish either because s
2 = 0 or sµ = 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γ → 2γ amplitude given
– 11 –
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γ → 4γ 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 contribu-
tion of a fermion loop to the 2γ → 2γ 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γ → 4γ amplitude with zero internal mass7,
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
0 0.5 1 1.5 2 2.5 3
10000
15000
20000
25000
Figure 1: Comparison with Fig. 5 of Ref. [19]. Helicity configurations [+ +−−−−] and
[+ −−++−] for the momenta of Eq. (4.1), represented by black dots and gray diamonds re-
spectively, and comparison with the analytic result of Ref. [18] (continuous line).
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 {~p3, ~p4, ~p5, ~p6}:
~p3 = (33.5, 15.9, 25.0) ,
~p4 = (−12.5, 15.3, 0.3) ,
~p5 = (−10.0,−18.0,−3.3) ,
~p6 = (−11.0,−13.2,−22.0) . (4.1)
After choosing the incoming photons such that they have momenta ~p1 and ~p2 along the
z-axis, we present in the plot the amplitude obtained by rotating the final states of an-
gle θ about the y-axis. This is done for both helicity configurations [+ +−−−−] and
7We thank Andre van Hameren for providing us with his program to compute massless one-loop scalar
integrals.
– 12 –
0 0.5 1 1.5 2 2.5 3
10000
12000
14000
16000
18000
Figure 2: Helicity configurations [+ +−−−−] and [+ +−−+−] for the momenta of Eq. (4.2),
represented by black dots and gray diamonds respectively, and comparison with the analytic result
of Ref. [18] (continuous line).
[+ −−++−]. 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 mo-
menta. Starting from the following choice of {~p3, ~p4, ~p5, ~p6}:
~p3 = (−10.0,−10.0,−10.0) ,
~p4 = (12.0,−15.0,−2.0) ,
~p5 = (10.0, 18.0, 3.0) ,
~p6 = (−12.0, 7.0, 9.0) (4.2)
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.1), we computed the amplitude introducing a non-zero mass mf for the
fermions in the loop 8. The results are plotted in Fig. 3, for the three cases mf = 0.5 GeV,
mf = 4.5 GeV and mf = 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.
8We used here the scalar one-loop functions provided by FF [27].
– 13 –
0 0.5 1 1.5 2
10000
15000
20000
25000
30000
Figure 3: Helicity configuration [+ +−−−−] for the momenta of Eq. (4.1) for different values
of the fermion mass in the loop: mf = 0.5 GeV (diamond), mf = 4.5 GeV (gray box) and mf = 12
GeV (black dots). The continuous line is the result for the massless case.
5. Conclusions
Computing the massless QED amplitude for the reaction 2γ → 4γ, although still unob-
served 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 mass-
less. 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
– 14 –
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 Pier-
paolo 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 re-
search 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
A. Computing the extra-integrals
In this appendix, we compute the extra-integrals listed in Section 2. Since a contribution
O(1) can only develop for non-negative dimensionality D, the integrand in the Feynman
parameter integral is always polynomial. First we decompose the integration as follows
dnq̄ =
d4q dǫµ (q̃2 = −µ2) , (A.1)
then, after using Feynman parametrization and performing first the integral over dǫµ and
then that one over d4q, one derives, for the extra-integrals of Eqs. (2.14)-(2.18)
I(n;2(s−1))s = −iπ
2Γ(s− 1)
[dα]s +O(ǫ) ,
I(n;2(s−1))s;µ = iπ
2Γ(s− 1)
[dα]s (Ps)µ +O(ǫ) ,
I(n;2(s−1))s;µν = −iπ
2Γ(s− 1)
[dα]s (Ps)µ(Ps)ν +O(gµν) +O(ǫ) ,
I(n;2s)s = −iπ
2Γ(s)
[dα]s Xs +O(ǫ) ,
I(n;2s)s;µ = iπ
2Γ(s)
[dα]s Xs(Ps)µ +O(ǫ) , (A.2)
– 15 –
where
[dα]s =
dα0 · · · dαs δ(1 −
αj) , Xs = P
αjkj , M
j − k
j ) , (k0 = 0) . (A.3)
In the following, we compute, as an illustrative example, the first three integrals of Eq. (A.2).
The remaining two can be obtained analogously. We start by changing the integration vari-
ables as follows:
α1 = ρ1ρ2 · · · ρs
α2 = ρ1ρ2 · · · ρs−1(1− ρs)
α3 = ρ1ρ2 · · · ρs−2(1− ρs−1)
αs = ρ1(1− ρ2)
α0 = (1− ρ1) , (A.4)
so that
[dα]s =
dρ2 · · ·
dρs ρ
(s−1)
(s−2)
2 · · · ρs−1 , (A.5)
from which one trivially obtains the first integral
I(n;2(s−1))s = −iπ
2Γ(s− 1)
Γ(s+ 1)
+O(ǫ) . (A.6)
To compute the second integral an integration over (Ps)µ in needed. Since the integrand
is symmetric when interchanging all ki, we concentrate on the coefficient of, say, k1. Since
dρ2 · · ·
dρs ρ
(s−1)
(s−2)
2 · · · ρs−1 α1k1µ
= k1µ
dρ2 · · ·
dρs ρ
(s−1)
2 · · · ρ
s−1ρs
= k1µ
Γ(s+ 2)
, (A.7)
the final result reads
I(n;2(s−1))s;µ = iπ
2Γ(s− 1)
Γ(s+ 2)
(kj)µ +O(ǫ) . (A.8)
To compute the third integral we need to integrate over the product (Ps)µ(Ps)ν . Once again,
given the symmetry of the problem, we can focus on the two contributions proportional to
– 16 –
k1µk1ν and k1µk2ν . The first one gives
dρ2 · · ·
dρs ρ
(s−1)
(s−2)
2 · · · ρs−1 α
1k1µk1ν
= k1µk1ν
dρ2 · · ·
dρs ρ
(s+1)
2 · · · ρ
= k1µk1ν
Γ(s+ 3)
, (A.9)
and the second reads
dρ2 · · ·
dρs ρ
(s−1)
(s−2)
2 · · · ρs−1 α1α2k1µk2ν
= k1µk2ν
dρ2 · · ·
dρs ρ
(s+1)
2 · · · ρ
s−1ρs(1− ρs)
= k1µk2ν
Γ(s+ 3)
. (A.10)
Summing up all of the possibilities one obtains
I(n;2(s−1))s;µν = −2iπ
2Γ(s− 1)
Γ(s+ 3)
(kj)µ(kj)ν +
i 6=j
(kj)µ(ki)ν
+ O(gµν) +O(ǫ) . (A.11)
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
A(q̄) =
D̄0D̄1
, (B.1)
by assuming N(q) at most quadratic in q and k1 ≡ (p1 − p0) 6= 0. In particular also the
case of vanishing k21 is included. First, we introduce a massless arbitrary 4-vector v, such
that (v · k1) 6= 0, that we use to rewrite k1 in terms of two massless 4-vectors (we also take
ℓ2 = 0)
k1 = ℓ+ α v , (B.2)
giving
γ ≡ 2 (k1 · v) = 2 (ℓ · v) and α =
. (B.3)
Then, we introduce two additional independent massless 4-vectors ℓ7,8 defined as
7 = < ℓ|γ
µ|v] , ℓ
8 =< v|γ
µ|ℓ] , (B.4)
– 17 –
for which one finds
(ℓ7 · ℓ8) = −2γ , (B.5)
and we decompose qµ + p
0 in the basis of k1, v, ℓ7 and ℓ8
qµ = −p
0 + yk
1 + yvv
µ + y7ℓ
7 + y8ℓ
8 , (B.6)
so that N(q) takes the form
N(q) = b+ b̂0[(q + p0) · v] + b̂00[(q + p0) · v]
2 + b̃11[(q + p0) · ℓ7] + b̃21[(q + p0) · ℓ8]
+ b̃12[(q + p0) · ℓ7]
2 + b̃22[(q + p0) · ℓ8]
+ b̃01[(q + p0) · ℓ7][(q + p0) · v]
+ b̃02[(q + p0) · ℓ8][(q + p0) · v] +O(D0) +O(D1) . (B.7)
Notice that, because of the identity
2 (q · k1) = D1 −D0 + (d1 − d0) , with di = m
i − p
i , (B.8)
any term proportional to [(q + p0) · k1] either contributes to the constant term b or it is
included in the terms O(D0,1) we are neglecting
9. The same happens for the combination
[(q + p0) · ℓ7][(q + p0) · ℓ8].
To be able to determine all of the coefficients appearing in Eq. (B.7), disentangling
completely the contributions O(D0,1), we look for a q that fulfill the requirement
D0 = D1 = 0 . (B.9)
For a q written as in Eq. (B.6) this implies the system
y7y8 = Fy
d1 − d0 − 2yk
, (B.10)
where
Fy = −
m20 − y (d1 − d0) + y
. (B.11)
It is convenient to introduce two classes of solutions. In the first class, that we call q±yk,
we take y fixed and choose y7 = ±e
iπ/k. In the second class, that we call q′±
, we take y
fixed but choose y8 = ±e
iπ/k. The coefficients b, b̃11, b̃21, b̃12 and b̃22 can be obtained by
evaluating Eq. (B.7) at the values
q±01 , q
02 , q
03 , (B.12)
q′±01 , q
02 , q
03 . (B.13)
9We suppose to determine them at a later stage of the calculation.
– 18 –
In the first case, the coefficients read
b0 = b , b1 = −2γb̃21 , b2 = 4γ
2b̃22 , b−1 = −2γF0b̃11 , b−2 = 4γ
2F 20 b̃12 , (B.14)
b±1 = −
T−(q1)± iT
−(q2)
T+(q1) + T
+(q2)
b±2 =
T+(q1)− T
+(q2)
− e±2iπ/3(T+(q3)− b0)
1− e∓2iπ/3
, (B.15)
and where
T±(qk) ≡
N(q+0k)±N(q
. (B.16)
In the second case, one obtains instead
b′0 = b , b
1 = −2γb̃11 , b
2 = 4γ
2b̃12 , b
−1 = −2γF0b̃21 , b
−2 = 4γ
2F 20 b̃22 , (B.17)
b′±1 = −
T−(q′1)± iT
−(q′2)
b′0 =
T+(q′1) + T
+(q′2)
b′±2 =
T+(q′1)− T
+(q′2)
− e±2iπ/3(T+(q′3)− b
1− e∓2iπ/3
, (B.18)
and where
T±(q′k) ≡
N(q′+0k)±N(q
. (B.19)
The reason why we have chosen two sets of solutions is that, in some special kinematical
configurations, F0 can vanish. Therefore, numerical stable solutions are obtained by taking
b̃21 and b̃22 from Eq. (B.14), and b̃11 and b̃12 from Eq. (B.17), while b is well defined in
both cases.
The coefficients b̂0 and b̂00 can be determined, in terms of additional solutions of the
kind q±
and q±σ1, by defining the combinations
S(q) ≡ N(q)− b− b̃11[(q + p0) · ℓ7]− b̃21[(q + p0) · ℓ8]
− b̃12[(q + p0) · ℓ7]
2 − b̃22[(q + p0) · ℓ8]
U(λ) ≡
S(q+λ1) + S(q
, (B.20)
as the two solutions of the system
. (B.21)
– 19 –
The determinant of the matrix above is always different form zero, for non vanishing λ and
σ, when σ 6= λ, so that numerical inaccuracies never occur.
Finally, the two last coefficients b̃01 and b̃02 are determined, in terms of q
and q′+
as solutions of the system
Z(q+λk)
Z(q′+σk)
−λγ2Fλe
−iπ/k −λγ2eiπ/k
−σγ2eiπ/k −σγ2Fσe
−iπ/k
, (B.22)
where
Z(q) ≡ S(q)− b̂0[(q + p0) · v]− b̂00[(q + p0) · v]
2 . (B.23)
Once again one verifies that when, for example, k = 3 the system never becomes singular.
References
[1] R. K. Ellis, W. T. Giele and G. Zanderighi, JHEP 0605 (2006) 027 [arXiv:hep-ph/0602185];
R. Britto, B. Feng and P. Mastrolia, Phys. Rev. D 73 (2006) 105004 [arXiv:hep-ph/0602178];
C. F. Berger, Z. Bern, L. J. Dixon, D. Forde and D. A. Kosower, Phys. Rev. D 74 (2006)
036009 [arXiv:hep-ph/0604195];
Z. Bern, N. E. J. Bjerrum-Bohr, D. C. Dunbar and H. Ita, JHEP 0511 (2005) 027
[arXiv:hep-ph/0507019] ;
J. Bedford, A. Brandhuber, B. J. Spence and G. Travaglini, Nucl. Phys. B 712 (2005) 59
[arXiv:hep-th/0412108].
[2] See, for example, G. Belanger et al., Phys. Lett. B 576 (2003) 152 [arXiv:hep-ph/0309010];
A. Denner, S. Dittmaier, M. Roth and M. M. Weber, Nucl. Phys. B 680, 85 (2004)
[arXiv:hep-ph/0309274].
[3] See, for example, A. Denner, S. Dittmaier, M. Roth and L. H. Wieders, Nucl. Phys. B 724
(2005) 247 [arXiv:hep-ph/0505042] and Phys. Lett. B 612 (2005) 223 [arXiv:hep-ph/0502063];
K. Kato et al., PoS HEP2005 (2006) 312.
[4] G. ’t Hooft and M. J. G. Veltman, Nucl. Phys. B 153 (1979) 365.
[5] G. Passarino and M. J. G. Veltman, Nucl. Phys. B 160 (1979) 151.
[6] A. Denner and S. Dittmaier, Nucl. Phys. B 734 (2006) 62 [arXiv:hep-ph/0509141].
[7] T. Binoth, J. P. Guillet and G. Heinrich, Nucl. Phys. B 572 (2000) 361
[arXiv:hep-ph/9911342];
[8] A. Ferroglia, M. Passera, G. Passarino and S. Uccirati, Nucl. Phys. B 650 (2003) 162
[arXiv:hep-ph/0209219];
W. T. Giele and E. W. N. Glover, arXiv:hep-ph/0402152;
D. E. Soper, Phys. Rev. D 62 (2000) 014009 [arXiv:hep-ph/9910292] and Phys. Rev. D 64
(2001) 034018 [arXiv:hep-ph/0103262];
Z. Nagy and D. E. Soper, JHEP 0309 (2003) 055 [arXiv:hep-ph/0308127].
[9] Z. Bern, L. J. Dixon, D. C. Dunbar and D. A. Kosower, Nucl. Phys. B 435 (1995) 59
[arXiv:hep-ph/9409265];
Z. Bern, L. J. Dixon, D. C. Dunbar and D. A. Kosower, Nucl. Phys. B 425, 217 (1994).
– 20 –
[10] R. Britto, F. Cachazo and B. Feng, Nucl. Phys. B 725, 275 (2005).
[11] E. Witten, Commun. Math. Phys. 252, 189 (2004);
F. Cachazo, P. Svrcek and E. Witten, JHEP 0409 (2004) 006;
A. Brandhuber, B. J. Spence and G. Travaglini, Nucl. Phys. B 706, 150 (2005);
F. Cachazo, P. Svrcek and E. Witten, JHEP 0410, 074 (2004);
I. Bena, Z. Bern, D. A. Kosower and R. Roiban, Phys. Rev. D 71, 106010 (2005).
[12] C. Anastasiou, R. Britto, B. Feng, Z. Kunszt and P. Mastrolia, JHEP 0703 (2007) 111
[arXiv:hep-ph/0612277];
C. Anastasiou, R. Britto, B. Feng, Z. Kunszt and P. Mastrolia, Phys. Lett. B 645 (2007) 213
[arXiv:hep-ph/0609191];
D. Forde, arXiv:0704.1835 [hep-ph].
[13] T. Binoth, J. P. Guillet, G. Heinrich, E. Pilon and C. Schubert, JHEP 0510 (2005) 015
[arXiv:hep-ph/0504267].
[14] T. Binoth, J. P. Guillet and G. Heinrich, JHEP 0702 (2007) 013 [arXiv:hep-ph/0609054].
[15] Z. G. Xiao, G. Yang and C. J. Zhu, Nucl. Phys. B 758 (2006) 1 [arXiv:hep-ph/0607015];
X. Su, Z. G. Xiao, G. Yang and C. J. Zhu, Nucl. Phys. B 758 (2006) 35
[arXiv:hep-ph/0607016].
[16] Z. Bern, L. J. Dixon and D. A. Kosower, Phys. Rev. D 73 (2006) 065013
[arXiv:hep-ph/0507005];
S. D. Badger, E. W. N. Glover and K. Risager, arXiv:0704.3914 [hep-ph].
[17] G. Ossola, C. G. Papadopoulos and R. Pittau, Nucl. Phys. B 763 (2007) 147
[arXiv:hep-ph/0609007].
[18] G. Mahlon, Phys. Rev. D 49 (1994) 2197 [arXiv:hep-ph/9311213] and arXiv:hep-ph/9412350.
[19] Z. Nagy and D. E. Soper, Phys. Rev. D 74 (2006) 093006 [arXiv:hep-ph/0610028].
[20] T. Binoth, T. Gehrmann, G. Heinrich and P. Mastrolia, arXiv:hep-ph/0703311.
[21] R. Pittau, Comput. Phys. Commun. 104, 23 (1997) [arXiv:hep-ph/9607309] and 111 (1998)
48 [arXiv:hep-ph/9712418].
[22] F. del Aguila and R. Pittau, JHEP 0407 (2004) 017 [arXiv:hep-ph/0404120];
R. Pittau, arXiv:hep-ph/0406105.
[23] G. Devaraj and R. G. Stuart, Nucl. Phys. B 519 (1998) 483 [arXiv:hep-ph/9704308].
[24] A. Denner and S. Dittmaier, Nucl. Phys. Proc. Suppl. 157 (2006) 53 [arXiv:hep-ph/0601085].
[25] G. J. Gounaris, P. I. Porfyriadis and F. M. Renard, Eur. Phys. J. C 9 (1999) 673
[arXiv:hep-ph/9902230].
[26] A. van Hameren, J. Vollinga and S. Weinzierl, Eur. Phys. J. C 41 (2005) 361
[arXiv:hep-ph/0502165].
[27] G. J. van Oldenborgh, Comput. Phys. Commun. 66 (1991) 1.
– 21 –

  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.

Introduction 1
2. The method and the computation of the rational terms 2
3. Dealing with numerical instabilities 8
4. Results and comparisons 11
5. Conclusions 14
A. Computing the extra-integrals 15
B. The general basis for the 2-point functions 17
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 an-
alytic 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 re-
duction [4–7]. For multi-particle processes though this method becomes quite cumbersome
because of the large number of terms generated and the appearance of numerical insta-
bilities 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
– 1 –
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–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 imple-
mentation of such a method in the computation of a physical process: namely 2γ → 4γ,
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]
A(q̄) =
D̄0D̄1 · · · D̄m−1
, D̄i = (q̄ + pi)
2 −m2i , p0 6= 0 , (2.1)
where we use a bar to denote objects living in n = 4 + ǫ dimensions, and q̄2 = q2 + q̃2
1. In the previous equation, N(q) is the 4-dimensional part of the numerator function of
1q̃2 is ǫ-dimensional and (q̃ · q) = 0.
– 2 –
the amplitude 2. N(q) depends on the 4-dimensional denominators Di = (q+ pi)
2 −m2i as
follows
N(q) =
i0<i1<i2<i3
d(i0i1i2i3) + d̃(q; i0i1i2i3)
i 6=i0,i1,i2,i3
i0<i1<i2
[c(i0i1i2) + c̃(q; i0i1i2)]
i 6=i0,i1,i2
i0<i1
b(i0i1) + b̃(q; i0i1)
i 6=i0,i1
[a(i0) + ã(q; i0)]
i 6=i0
+ P̃ (q)
Di . (2.2)
Inserted back in Eq. (2.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 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 d̃, c̃, b̃, 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.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 Di 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 the 2 (in general complex)
solutions for which
D0 = D1 = D2 = D3 = 0 , (2.3)
(there are 2 solutions because of the quadratic nature of the propagators) and since the
functional form of d̃(q; 0123) is known, one directly finds the coefficient of the box diagram
containing the above 4 denominators through the two simple equations
N(q±0 ) = [d(0123) + d̃(q
0 ; 0123)]
i 6=0,1,2,3
0 ) . (2.4)
2If needed, the ǫ-dimensional part of the numerator should be treated separately, as explained in [21].
– 3 –
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 d̃, c̃, b̃, ã still depend on q, also the
separation among terms in Eq. (2.2) is somehow arbitrary. Terms containing a different
numbers of denominators can be shifted from one piece to the other in Eq. (2.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 even-
tually lead to integrals that need to be regularized. A typical example are the rank six
6-point functions contributing to the scattering 2γ → 4γ we want to study. Such tensors
are finite, but tensor reduction iteratively leads to rank m m-point tensors with 1 ≤ m ≤ 5,
that are ultraviolet divergent when m ≤ 4. For this reason, we introduced, in Eq. (2.1),
the d-dimensional denominators D̄i, that differs by an amount q̃
2 from their 4-dimensional
counterparts
D̄i = Di + q̃
2 . (2.5)
The result of this is a mismatch in the cancellation of the d-dimensional denominators of
Eq. (2.1) with the 4-dimensional ones of Eq. (2.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. (2.1) as follows
, with Z̄i ≡
. (2.6)
This results in
A(q̄) =
D0D1 · · ·Dm−1
Z̄0Z̄1 · · · Z̄m−1 . (2.7)
Then, by inserting Eq. (2.2) in Eq. (2.7), one obtains
A(q̄) =
i0<i1<i2<i3
d(i0i1i2i3) + d̃(q; i0i1i2i3)
D̄i0D̄i1D̄i2D̄i3
i 6=i0,i1,i2,i3
– 4 –
i0<i1<i2
c(i0i1i2) + c̃(q; i0i1i2)
D̄i0D̄i1D̄i2
i 6=i0,i1,i2
i0<i1
b(i0i1) + b̃(q; i0i1)
D̄i0D̄i1
i 6=i0,i1
a(i0) + ã(q; i0)
i 6=i0
+ P̃ (q)
Z̄i . (2.8)
The rational part of the amplitude is then produced, after integrating over dnq, by the
q̃2 dependence coming from the various Z̄i in Eq. (2.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]
P̃ (q) = 0 ,
ã(q; i0) = ã
µ(i0; 1)(q + pi0)µ ,
b̃(q; i0i1) = b̃
µ(i0i1; 1)(q + pi0)µ + b̃
µν(i0i1; 2)(q + pi0)µ(q + pi0)ν ,
c̃(q; i0i1i2) = c̃
µ(i0i1i2; 1)(q + pi0)µ + c̃
µν(i0i1i1; 2)(q + pi0)µ(q + pi0)ν ,
+ c̃µνρ(i0i1i1; 3)(q + pi0)µ(q + pi0)ν(q + pi0)ρ ,
d̃(q; i0i1i2i3) = d̃
µ(i0i1i2i3; 1)(q + pi0)µ . (2.9)
Eq. (2.9) simply states the fact that ã(q; i0) and d̃(q; i0i1i2i3) are at most linear in (q+pi0),
b̃(q; i0i1) at most quadratic, and c̃(q; i0i1i2) at most cubic. The tensors denoted by (· · · ; 1),
(· · · ; 2) and (· · · ; 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]
c̃µν(i0i1i2; 2) gµν = 0 ,
c̃µνρ(i0i1i2; 3) gµν = c̃
µνρ(i0i1i2; 3) gµρ = c̃
µνρ(i0i1i2; 3) gνρ = 0 and
b̃µν(i0i1; 2) gµν = 0 . (2.10)
The necessary integrals that arise, after a change of variable q → q − pi0 , are of the form
(n;2ℓ)
s;µ1···µr ≡
dnq q̃2ℓ
qµ1 · · · qµr
D̄(k0) · · · D̄(ks)
, with
D̄(ki) ≡ (q̄ + ki)
2 −m2i , ki ≡ pi − p0 (k0 = 0) , (2.11)
where we used a notation introduced in [22] and r ≤ 3. Such integrals (from now on called
extra-integrals) have dimensionality D = 2(1+ ℓ− s)+ r and give a contribution O(1) only
when D ≥ 0, otherwise are of O(ǫ). 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].
– 5 –
We also note that, since all Z̄i are a-dimensional, the dimensionality D of the extra-
integrals generated through Eq. (2.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(i0i1i2i3)
In this case r = 0. All extra-integrals are therefore scalars with D = −4 and do not
contribute.
Contributions proportional to d̃µ(i0i1i2i3; 1)
In this case r = 1. All extra-integrals are therefore rank one tensors with D = −3 and do
not contribute.
Contributions proportional to c(i0i1i2)
In this case r = 0 with D = −2 and no contribution O(1) is developed.
Contributions proportional to c̃µ(i0i1i2; 1)
Here r = 1 and D = −1. Therefore, once again, there is no contribution.
Contributions proportional to c̃µν(i0i1i2; 2)
Now r = 2 with D = 0 and a finite contribution is in principle expected, generated by
extra-integrals of the type
I(n;2(s−2))s;µν . (2.12)
Nevertheless, such contribution is proportional to gµν [22]. Therefore, due to Eq. (2.10), it
vanishes.
Contributions proportional to c̃µνρ(i0i1i2; 3)
Now r = 3 and D = 1. The contributing extra-integrals are of the type
I(n;2(s−2))s;µνρ , (2.13)
and one easily proves that the contributions O(1) are always proportional to gµν or gµρ or
gνρ. Therefore, thanks again to Eq. (2.10), they vanish.
Contributions proportional to b(i0i1)
Those are the first non vanishing contributions. The relevant extra-integrals have r = 0
and D = 0
I(n;2(s−1))s ,
with 2 < s ≤ m − 1. They have been computed, for generic values of s, in [22] (see also
Appendix A)
I(n;2(s−1))s = −iπ
s(s− 1)
+O(ǫ) . (2.14)
– 6 –
Contributions proportional to b̃µ(i0i1; 1)
In this case the relevant extra-integrals are 4-vectors with D = 1
I(n;2(s−1))s;µ with 2 < s ≤ m− 1 .
A computation for generic values of s gives
I(n;2(s−1))s;µ = iπ
(s + 1)s(s − 1)
(kj)µ +O(ǫ) . (2.15)
Contributions proportional to b̃µν(i0i1; 2)
The relevant extra-integrals are now rank two tensor with D = 2
I(n;2(s−1)s;µν with 2 < s ≤ m− 1 .
They read
I(n;2(s−1))s;µν = −2iπ
(s+ 2)(s + 1)s(s− 1)
(kj)µ(kj)ν +
i 6=j
(kj)µ(ki)ν
+ O(gµν) +O(ǫ) . (2.16)
The gµν part is never needed because b̃
µν(i0i1; 2) gµν = 0, according to Eq. (2.10).
Contributions proportional to a(i0)
They involve scalar extra-integrals with D = 2
I(n;2s)s , with 1 < s ≤ m− 1 .
One computes
I(n;2s)s = −2iπ
(s+ 2)(s + 1)s
k2j +
i 6=j
(kj · ki) +
(m2j − k
+ O(ǫ) . (2.17)
Contributions proportional to ãµ(i0; 1)
This last category involves extra-integrals with r = 1 and D = 3
I(n;2s)s;µ , with 1 < s ≤ m− 1 .
One obtains
I(n;2s)s;µ = iπ
(s+ 3)(s + 2)(s + 1)s
k2j (kj)µ + 2
i 6=j
k2j (ki)µ + 2(kj · ki)(kj)µ
– 7 –
i 6=j
ℓ 6=i
(kj · ki)(kℓ)µ + (s+ 3)
(m2j − k
j )(kj)µ
i 6=j
(m2j − k
j )(ki)µ
+O(ǫ) . (2.18)
To conclude, the set of the five formulas in Eqs. (2.14)-(2.18) allows one to compute the
rational part of any one-loop m-point (sub-)amplitude, once all the coefficients of Eq. (2.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 3. 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
A(q̄) =
D̄0D̄1
, (3.1)
in which we supposeN(q) at most quadratic in q. Our purpose is dealing with the situation
in which k21 ≡ (p1−p0)
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.2), we can write an expansion for N(q) as follows:
N(q) = [b(01) + b̃(q; 01)] + [a(0) + ã(q; 0)]D1 + [a(1) + ã(q; 1)]D0 . (3.2)
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 D0 = D1 = 0, needed to determine
the coefficients b and b̃, becomes singular4, in the limit of k21 → 0, when adding the
requirement
dnq b̃(q; 01) = 0 . (3.3)
Then, we must consider two separate cases:
k21 → 0 , but k
1 6= 0 ,
k21 → 0 , because k
1 = 0 . (3.4)
3In this case the Gram-determinant is simply the square of the difference between the momenta of the
two denominators.
4Such a solution goes like 1/k21 .
– 8 –
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 D0 = D1 = 0 by relaxing the further requirement
of Eq. (3.3). Such a solution is given in Appendix B and goes like 1/(k1.v), where v is an
arbitrary massless 4-vector, therefore is never singular in the first case of Eq. (3.4). The
price to pay is that new non zero integrals appear of the type 5
[(q + p0) · v]
D̄0D̄1
with j = 1, 2 and v2 = 0 . (3.5)
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 k21 → 0
appear. The fact that solutions exist to the condition D0 = D1 = 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 k21. Let us
first consider the case j = 1 in Eq. (3.5). The contribution O(1) can be easily obtained
from the observation that6
(q · v)(q · k1)
D̄(k0)D̄(k1)
= O(k21) , (3.6)
as it is evident by performing a tensor decomposition. On the other hand, by reconstructing
denominators, one obtains
(q · k1)
D̄(k1)− D̄(k0)
(q · k1) +
, (3.7)
f = m21 − k
0 . (3.8)
Eq. (3.7), inserted in Eq. (3.6) gives the desired expansion in terms of loop functions with
less points but higher rank, in agreement with well know results [23,24]
(q · v)
D̄(k0)D̄(k1)
dnq (q · v)
D̄(k1)
D̄(k0)
2(q · k1)
+O(k21) .(3.9)
Expansions at arbitrary orders in k21 can be obtained in an analogous way from the two
following equations:
(q · k1)
D̄(k1)− D̄(k0)
i+j=p−1
(q · k1)
(q · v)(q · k1)
D̄(k0)D̄(k1)
= O(k
1 ) . (3.10)
5Since v2 = 0 they still fulfill the third one of Eqs. (2.10), therefore, even in this case, terms O(gµν) can
be neglected in Eq. (2.16).
6From now on, we shift the integration variable: q̄ → q̄ − p0. The definition of the new resulting
denominators is given in Eq. (2.11).
– 9 –
To deal with the case j = 2 in Eq. (3.5) one starts instead from the equation
(q · v)2(q · k1)
D̄(k0)D̄(k1)
= O(k
1 ) . (3.11)
This procedure breaks down when the quantity f vanishes. In this case a double expansion
in k21 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. (3.5). Our starting point is now the
equation
D̄(k0) = D̄(k1)− 2(q · k1) + f . (3.12)
By multiplying and dividing by D̄(k0) one obtains
(q · v)
D̄(k0)D̄(k1)
(q · v)
D̄(k0)2D̄(k1)
D̄(k1)− 2(q · k1) + f
(q · v)
D̄(k0)2
(q · v)(q · k1)
D̄(k0)2D̄(k1)
+O(f) . (3.13)
Applying once more Eq. (3.12) to the last integral gives
(q · v)(q · k1)
D̄(k0)
2D̄(k1)
(q · v)(q · k1)
D̄(k0)
3D̄(k1)
D̄(k1)− 2(q · k1) + f
(q · v)(q · k1)
D̄(k0)3
(q · v)(q · k1)
D̄(k0)3D̄(k1)
+O(f) . (3.14)
Since the last integral in the previous equation is O(k21), the final result reads
(q · v)
D̄(k0)D̄(k1)
(q · v)
D̄(k0)2
(q · v)(q · k1)
D̄(k0)3
+O(k21) +O(f) .
(3.15)
In a similar fashion, expansions at any order can be obtained.
We now turn to the second case of Eq. (3.4), namely k
1 → 0. In this case no solution
can be found to the double cut equation
D(k0) = D(k1) = 0 . (3.16)
The reason is that now D(k1) and D(k0) are no longer independent:
D(k0) = D(k1) + f +O(k1) , (3.17)
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. (3.2). This results is a singularity 1/(k1 · v) in the system given of
Appendix B and we should change our strategy. We than go back to Eq. (3.1) and split
the amplitude from the beginning by multiplying it by
D̄(k0)− D̄(k1)
2(q · k1)
, (3.18)
– 10 –
resulting to
A(q̄) = A(1)(q̄) +A(2)(q̄) +O(k1) , (3.19)
A(1)(q̄) =
D̄(k1)
, A(2)(q̄) = −
D̄(k0)
. (3.20)
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 O(k1) are
perfectly calculable by inserting again Eq. (3.18) in the term O(k1) of Eq. (3.19).
Once again, when f → 0, double expansions in k1 and f can be obtained involving
derivatives of scalar loop functions by using Eq. (3.12). For example, at the zeroth order
in k1 and at the first one in f , one gets
A(q) =
D̄(k0)D̄(k1)
D̄(k0)2D̄(k1)
D̄(k1)− 2(q · k1) + f
D̄(k0)
D̄(k0)
3D̄(k1)
D̄(k1)− 2(q · k1) + f
+O(k1)
D̄(k0)
D̄(k0)
+O(k1) +O(f
2) . (3.21)
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 k1, one introduces the 4-vector
sµ = det
(k2 · k1) (k2 · k2)
, (3.22)
with the properties
s · k2 = 0 , s
2 ∝ ∆(k1, k2) , (k1 · s) ∝ ∆(k1, k2) , (3.23)
where ∆(k1, k2) is the Gram-determinant of the two momenta k1 and k2. Then, instead of
Eq. (3.6) one has, for example,
(q · v)(q · s)2
D̄(k0)D̄(k1)D̄(k2)
= O(∆(k1, k2)) . (3.24)
As before, ∆(k1, k2) can vanish either because s
2 = 0 or sµ = 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γ → 2γ amplitude given
– 11 –
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γ → 4γ 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 contribu-
tion of a fermion loop to the 2γ → 2γ 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γ → 4γ amplitude with zero internal mass7,
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
0 0.5 1 1.5 2 2.5 3
10000
15000
20000
25000
Figure 1: Comparison with Fig. 5 of Ref. [19]. Helicity configurations [+ +−−−−] and
[+ −−++−] for the momenta of Eq. (4.1), represented by black dots and gray diamonds re-
spectively, and comparison with the analytic result of Ref. [18] (continuous line).
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 {~p3, ~p4, ~p5, ~p6}:
~p3 = (33.5, 15.9, 25.0) ,
~p4 = (−12.5, 15.3, 0.3) ,
~p5 = (−10.0,−18.0,−3.3) ,
~p6 = (−11.0,−13.2,−22.0) . (4.1)
After choosing the incoming photons such that they have momenta ~p1 and ~p2 along the
z-axis, we present in the plot the amplitude obtained by rotating the final states of an-
gle θ about the y-axis. This is done for both helicity configurations [+ +−−−−] and
7We thank Andre van Hameren for providing us with his program to compute massless one-loop scalar
integrals.
– 12 –
0 0.5 1 1.5 2 2.5 3
10000
12000
14000
16000
18000
Figure 2: Helicity configurations [+ +−−−−] and [+ +−−+−] for the momenta of Eq. (4.2),
represented by black dots and gray diamonds respectively, and comparison with the analytic result
of Ref. [18] (continuous line).
[+ −−++−]. 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 mo-
menta. Starting from the following choice of {~p3, ~p4, ~p5, ~p6}:
~p3 = (−10.0,−10.0,−10.0) ,
~p4 = (12.0,−15.0,−2.0) ,
~p5 = (10.0, 18.0, 3.0) ,
~p6 = (−12.0, 7.0, 9.0) (4.2)
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.1), we computed the amplitude introducing a non-zero mass mf for the
fermions in the loop 8. The results are plotted in Fig. 3, for the three cases mf = 0.5 GeV,
mf = 4.5 GeV and mf = 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.
8We used here the scalar one-loop functions provided by FF [27].
– 13 –
0 0.5 1 1.5 2
10000
15000
20000
25000
30000
Figure 3: Helicity configuration [+ +−−−−] for the momenta of Eq. (4.1) for different values
of the fermion mass in the loop: mf = 0.5 GeV (diamond), mf = 4.5 GeV (gray box) and mf = 12
GeV (black dots). The continuous line is the result for the massless case.
5. Conclusions
Computing the massless QED amplitude for the reaction 2γ → 4γ, although still unob-
served 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 mass-
less. 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
– 14 –
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 Pier-
paolo 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 re-
search 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
A. Computing the extra-integrals
In this appendix, we compute the extra-integrals listed in Section 2. Since a contribution
O(1) can only develop for non-negative dimensionality D, the integrand in the Feynman
parameter integral is always polynomial. First we decompose the integration as follows
dnq̄ =
d4q dǫµ (q̃2 = −µ2) , (A.1)
then, after using Feynman parametrization and performing first the integral over dǫµ and
then that one over d4q, one derives, for the extra-integrals of Eqs. (2.14)-(2.18)
I(n;2(s−1))s = −iπ
2Γ(s− 1)
[dα]s +O(ǫ) ,
I(n;2(s−1))s;µ = iπ
2Γ(s− 1)
[dα]s (Ps)µ +O(ǫ) ,
I(n;2(s−1))s;µν = −iπ
2Γ(s− 1)
[dα]s (Ps)µ(Ps)ν +O(gµν) +O(ǫ) ,
I(n;2s)s = −iπ
2Γ(s)
[dα]s Xs +O(ǫ) ,
I(n;2s)s;µ = iπ
2Γ(s)
[dα]s Xs(Ps)µ +O(ǫ) , (A.2)
– 15 –
where
[dα]s =
dα0 · · · dαs δ(1 −
αj) , Xs = P
αjkj , M
j − k
j ) , (k0 = 0) . (A.3)
In the following, we compute, as an illustrative example, the first three integrals of Eq. (A.2).
The remaining two can be obtained analogously. We start by changing the integration vari-
ables as follows:
α1 = ρ1ρ2 · · · ρs
α2 = ρ1ρ2 · · · ρs−1(1− ρs)
α3 = ρ1ρ2 · · · ρs−2(1− ρs−1)
αs = ρ1(1− ρ2)
α0 = (1− ρ1) , (A.4)
so that
[dα]s =
dρ2 · · ·
dρs ρ
(s−1)
(s−2)
2 · · · ρs−1 , (A.5)
from which one trivially obtains the first integral
I(n;2(s−1))s = −iπ
2Γ(s− 1)
Γ(s+ 1)
+O(ǫ) . (A.6)
To compute the second integral an integration over (Ps)µ in needed. Since the integrand
is symmetric when interchanging all ki, we concentrate on the coefficient of, say, k1. Since
dρ2 · · ·
dρs ρ
(s−1)
(s−2)
2 · · · ρs−1 α1k1µ
= k1µ
dρ2 · · ·
dρs ρ
(s−1)
2 · · · ρ
s−1ρs
= k1µ
Γ(s+ 2)
, (A.7)
the final result reads
I(n;2(s−1))s;µ = iπ
2Γ(s− 1)
Γ(s+ 2)
(kj)µ +O(ǫ) . (A.8)
To compute the third integral we need to integrate over the product (Ps)µ(Ps)ν . Once again,
given the symmetry of the problem, we can focus on the two contributions proportional to
– 16 –
k1µk1ν and k1µk2ν . The first one gives
dρ2 · · ·
dρs ρ
(s−1)
(s−2)
2 · · · ρs−1 α
1k1µk1ν
= k1µk1ν
dρ2 · · ·
dρs ρ
(s+1)
2 · · · ρ
= k1µk1ν
Γ(s+ 3)
, (A.9)
and the second reads
dρ2 · · ·
dρs ρ
(s−1)
(s−2)
2 · · · ρs−1 α1α2k1µk2ν
= k1µk2ν
dρ2 · · ·
dρs ρ
(s+1)
2 · · · ρ
s−1ρs(1− ρs)
= k1µk2ν
Γ(s+ 3)
. (A.10)
Summing up all of the possibilities one obtains
I(n;2(s−1))s;µν = −2iπ
2Γ(s− 1)
Γ(s+ 3)
(kj)µ(kj)ν +
i 6=j
(kj)µ(ki)ν
+ O(gµν) +O(ǫ) . (A.11)
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
A(q̄) =
D̄0D̄1
, (B.1)
by assuming N(q) at most quadratic in q and k1 ≡ (p1 − p0) 6= 0. In particular also the
case of vanishing k21 is included. First, we introduce a massless arbitrary 4-vector v, such
that (v · k1) 6= 0, that we use to rewrite k1 in terms of two massless 4-vectors (we also take
ℓ2 = 0)
k1 = ℓ+ α v , (B.2)
giving
γ ≡ 2 (k1 · v) = 2 (ℓ · v) and α =
. (B.3)
Then, we introduce two additional independent massless 4-vectors ℓ7,8 defined as
7 = < ℓ|γ
µ|v] , ℓ
8 =< v|γ
µ|ℓ] , (B.4)
– 17 –
for which one finds
(ℓ7 · ℓ8) = −2γ , (B.5)
and we decompose qµ + p
0 in the basis of k1, v, ℓ7 and ℓ8
qµ = −p
0 + yk
1 + yvv
µ + y7ℓ
7 + y8ℓ
8 , (B.6)
so that N(q) takes the form
N(q) = b+ b̂0[(q + p0) · v] + b̂00[(q + p0) · v]
2 + b̃11[(q + p0) · ℓ7] + b̃21[(q + p0) · ℓ8]
+ b̃12[(q + p0) · ℓ7]
2 + b̃22[(q + p0) · ℓ8]
+ b̃01[(q + p0) · ℓ7][(q + p0) · v]
+ b̃02[(q + p0) · ℓ8][(q + p0) · v] +O(D0) +O(D1) . (B.7)
Notice that, because of the identity
2 (q · k1) = D1 −D0 + (d1 − d0) , with di = m
i − p
i , (B.8)
any term proportional to [(q + p0) · k1] either contributes to the constant term b or it is
included in the terms O(D0,1) we are neglecting
9. The same happens for the combination
[(q + p0) · ℓ7][(q + p0) · ℓ8].
To be able to determine all of the coefficients appearing in Eq. (B.7), disentangling
completely the contributions O(D0,1), we look for a q that fulfill the requirement
D0 = D1 = 0 . (B.9)
For a q written as in Eq. (B.6) this implies the system
y7y8 = Fy
d1 − d0 − 2yk
, (B.10)
where
Fy = −
m20 − y (d1 − d0) + y
. (B.11)
It is convenient to introduce two classes of solutions. In the first class, that we call q±yk,
we take y fixed and choose y7 = ±e
iπ/k. In the second class, that we call q′±
, we take y
fixed but choose y8 = ±e
iπ/k. The coefficients b, b̃11, b̃21, b̃12 and b̃22 can be obtained by
evaluating Eq. (B.7) at the values
q±01 , q
02 , q
03 , (B.12)
q′±01 , q
02 , q
03 . (B.13)
9We suppose to determine them at a later stage of the calculation.
– 18 –
In the first case, the coefficients read
b0 = b , b1 = −2γb̃21 , b2 = 4γ
2b̃22 , b−1 = −2γF0b̃11 , b−2 = 4γ
2F 20 b̃12 , (B.14)
b±1 = −
T−(q1)± iT
−(q2)
T+(q1) + T
+(q2)
b±2 =
T+(q1)− T
+(q2)
− e±2iπ/3(T+(q3)− b0)
1− e∓2iπ/3
, (B.15)
and where
T±(qk) ≡
N(q+0k)±N(q
. (B.16)
In the second case, one obtains instead
b′0 = b , b
1 = −2γb̃11 , b
2 = 4γ
2b̃12 , b
−1 = −2γF0b̃21 , b
−2 = 4γ
2F 20 b̃22 , (B.17)
b′±1 = −
T−(q′1)± iT
−(q′2)
b′0 =
T+(q′1) + T
+(q′2)
b′±2 =
T+(q′1)− T
+(q′2)
− e±2iπ/3(T+(q′3)− b
1− e∓2iπ/3
, (B.18)
and where
T±(q′k) ≡
N(q′+0k)±N(q
. (B.19)
The reason why we have chosen two sets of solutions is that, in some special kinematical
configurations, F0 can vanish. Therefore, numerical stable solutions are obtained by taking
b̃21 and b̃22 from Eq. (B.14), and b̃11 and b̃12 from Eq. (B.17), while b is well defined in
both cases.
The coefficients b̂0 and b̂00 can be determined, in terms of additional solutions of the
kind q±
and q±σ1, by defining the combinations
S(q) ≡ N(q)− b− b̃11[(q + p0) · ℓ7]− b̃21[(q + p0) · ℓ8]
− b̃12[(q + p0) · ℓ7]
2 − b̃22[(q + p0) · ℓ8]
U(λ) ≡
S(q+λ1) + S(q
, (B.20)
as the two solutions of the system
. (B.21)
– 19 –
The determinant of the matrix above is always different form zero, for non vanishing λ and
σ, when σ 6= λ, so that numerical inaccuracies never occur.
Finally, the two last coefficients b̃01 and b̃02 are determined, in terms of q
and q′+
as solutions of the system
Z(q+λk)
Z(q′+σk)
−λγ2Fλe
−iπ/k −λγ2eiπ/k
−σγ2eiπ/k −σγ2Fσe
−iπ/k
, (B.22)
where
Z(q) ≡ S(q)− b̂0[(q + p0) · v]− b̂00[(q + p0) · v]
2 . (B.23)
Once again one verifies that when, for example, k = 3 the system never becomes singular.
References
[1] R. K. Ellis, W. T. Giele and G. Zanderighi, JHEP 0605 (2006) 027 [arXiv:hep-ph/0602185];
R. Britto, B. Feng and P. Mastrolia, Phys. Rev. D 73 (2006) 105004 [arXiv:hep-ph/0602178];
C. F. Berger, Z. Bern, L. J. Dixon, D. Forde and D. A. Kosower, Phys. Rev. D 74 (2006)
036009 [arXiv:hep-ph/0604195];
Z. Bern, N. E. J. Bjerrum-Bohr, D. C. Dunbar and H. Ita, JHEP 0511 (2005) 027
[arXiv:hep-ph/0507019] ;
J. Bedford, A. Brandhuber, B. J. Spence and G. Travaglini, Nucl. Phys. B 712 (2005) 59
[arXiv:hep-th/0412108].
[2] See, for example, G. Belanger et al., Phys. Lett. B 576 (2003) 152 [arXiv:hep-ph/0309010];
A. Denner, S. Dittmaier, M. Roth and M. M. Weber, Nucl. Phys. B 680, 85 (2004)
[arXiv:hep-ph/0309274].
[3] See, for example, A. Denner, S. Dittmaier, M. Roth and L. H. Wieders, Nucl. Phys. B 724
(2005) 247 [arXiv:hep-ph/0505042] and Phys. Lett. B 612 (2005) 223 [arXiv:hep-ph/0502063];
K. Kato et al., PoS HEP2005 (2006) 312.
[4] G. ’t Hooft and M. J. G. Veltman, Nucl. Phys. B 153 (1979) 365.
[5] G. Passarino and M. J. G. Veltman, Nucl. Phys. B 160 (1979) 151.
[6] A. Denner and S. Dittmaier, Nucl. Phys. B 734 (2006) 62 [arXiv:hep-ph/0509141].
[7] T. Binoth, J. P. Guillet and G. Heinrich, Nucl. Phys. B 572 (2000) 361
[arXiv:hep-ph/9911342];
[8] A. Ferroglia, M. Passera, G. Passarino and S. Uccirati, Nucl. Phys. B 650 (2003) 162
[arXiv:hep-ph/0209219];
W. T. Giele and E. W. N. Glover, arXiv:hep-ph/0402152;
D. E. Soper, Phys. Rev. D 62 (2000) 014009 [arXiv:hep-ph/9910292] and Phys. Rev. D 64
(2001) 034018 [arXiv:hep-ph/0103262];
Z. Nagy and D. E. Soper, JHEP 0309 (2003) 055 [arXiv:hep-ph/0308127].
[9] Z. Bern, L. J. Dixon, D. C. Dunbar and D. A. Kosower, Nucl. Phys. B 435 (1995) 59
[arXiv:hep-ph/9409265];
Z. Bern, L. J. Dixon, D. C. Dunbar and D. A. Kosower, Nucl. Phys. B 425, 217 (1994).
– 20 –
[10] R. Britto, F. Cachazo and B. Feng, Nucl. Phys. B 725, 275 (2005).
[11] E. Witten, Commun. Math. Phys. 252, 189 (2004);
F. Cachazo, P. Svrcek and E. Witten, JHEP 0409 (2004) 006;
A. Brandhuber, B. J. Spence and G. Travaglini, Nucl. Phys. B 706, 150 (2005);
F. Cachazo, P. Svrcek and E. Witten, JHEP 0410, 074 (2004);
I. Bena, Z. Bern, D. A. Kosower and R. Roiban, Phys. Rev. D 71, 106010 (2005).
[12] C. Anastasiou, R. Britto, B. Feng, Z. Kunszt and P. Mastrolia, JHEP 0703 (2007) 111
[arXiv:hep-ph/0612277];
C. Anastasiou, R. Britto, B. Feng, Z. Kunszt and P. Mastrolia, Phys. Lett. B 645 (2007) 213
[arXiv:hep-ph/0609191];
D. Forde, arXiv:0704.1835 [hep-ph].
[13] T. Binoth, J. P. Guillet, G. Heinrich, E. Pilon and C. Schubert, JHEP 0510 (2005) 015
[arXiv:hep-ph/0504267].
[14] T. Binoth, J. P. Guillet and G. Heinrich, JHEP 0702 (2007) 013 [arXiv:hep-ph/0609054].
[15] Z. G. Xiao, G. Yang and C. J. Zhu, Nucl. Phys. B 758 (2006) 1 [arXiv:hep-ph/0607015];
X. Su, Z. G. Xiao, G. Yang and C. J. Zhu, Nucl. Phys. B 758 (2006) 35
[arXiv:hep-ph/0607016].
[16] Z. Bern, L. J. Dixon and D. A. Kosower, Phys. Rev. D 73 (2006) 065013
[arXiv:hep-ph/0507005];
S. D. Badger, E. W. N. Glover and K. Risager, arXiv:0704.3914 [hep-ph].
[17] G. Ossola, C. G. Papadopoulos and R. Pittau, Nucl. Phys. B 763 (2007) 147
[arXiv:hep-ph/0609007].
[18] G. Mahlon, Phys. Rev. D 49 (1994) 2197 [arXiv:hep-ph/9311213] and arXiv:hep-ph/9412350.
[19] Z. Nagy and D. E. Soper, Phys. Rev. D 74 (2006) 093006 [arXiv:hep-ph/0610028].
[20] T. Binoth, T. Gehrmann, G. Heinrich and P. Mastrolia, arXiv:hep-ph/0703311.
[21] R. Pittau, Comput. Phys. Commun. 104, 23 (1997) [arXiv:hep-ph/9607309] and 111 (1998)
48 [arXiv:hep-ph/9712418].
[22] F. del Aguila and R. Pittau, JHEP 0407 (2004) 017 [arXiv:hep-ph/0404120];
R. Pittau, arXiv:hep-ph/0406105.
[23] G. Devaraj and R. G. Stuart, Nucl. Phys. B 519 (1998) 483 [arXiv:hep-ph/9704308].
[24] A. Denner and S. Dittmaier, Nucl. Phys. Proc. Suppl. 157 (2006) 53 [arXiv:hep-ph/0601085].
[25] G. J. Gounaris, P. I. Porfyriadis and F. M. Renard, Eur. Phys. J. C 9 (1999) 673
[arXiv:hep-ph/9902230].
[26] A. van Hameren, J. Vollinga and S. Weinzierl, Eur. Phys. J. C 41 (2005) 361
[arXiv:hep-ph/0502165].
[27] G. J. van Oldenborgh, Comput. Phys. Commun. 66 (1991) 1.
– 21 –
