FERMILAB-PUB-21-263-T

Many-gluon tree amplitudes on modern GPUs
A case study for novel event generators

E. Bothmann1*, W. Giele2, S. Höche2, J. Isaacson2, M. Knobbe1

1 Institut für Theoretische Physik, Göttingen, Germany

2 Fermi National Accelerator Laboratory, Batavia, USA

∗<EMAIL_ADDRESS>
## Abstract

The compute efficiency of Monte-Carlo event generators for the Large Hadron
Collider is expected to become a major bottleneck for simulations in the high-
luminosity phase. Aiming at the development of a full-fledged generator for
modern GPUs, we study the performance of various recursive strategies to
compute multi-gluon tree-level amplitudes. We investigate the scaling of the
algorithms on both CPU and GPU hardware. Finally, we provide practical
recommendations as well as baseline implementations for the development of
future simulation programs. The GPU implementations can be found at:
https://www.gitlab.com/ebothmann/blockgen-archive.

##  1 Introduction

The success of high-energy physics experiments at particle colliders
critically depends on a detailed computer simulation of the collisions. So-
called event generators link theory and experiment through an event-by-event
representation of the observable final-state particles by means of Monte-Carlo
integration methods. A complete event generator consists of modules, which
encapsulate the physics at various energy scales [1, 2]. Interactions at the
highest scales, which are usually regarded as the most relevant to probing our
fundamental understanding of nature, are simulated by matrix element
generators. In a community effort, these programs have been nearly fully
automated. They can be used for computations at tree-level [3, 4, 5, 6, 7, 8,
9, 10], and at one-loop level [11, 12, 13, 14, 15, 16, 17, 18] in perturbation
theory and can handle a large class of fundamental and effective theories [19,
20, 21].

Such flexibility typically comes at the cost of compute efficiency. This is
currently a major problem, and will be prohibitive for the high-luminosity
phase of the LHC, where the projected reach in terms of final-state
multiplicities and large statistics data sets poses new challenges to event
generator codes [22, 23]. Some event generator collaborations have therefore
explored massively parallel computing techniques [24, 25, 26]. The complexity
of existing generators, and the requirement to continuously provide state-of-
the-art simulations for LHC physics have so far prevented the wider adoption
of new computing paradigms, such as GPUs. In the past, several conceptual
studies have shown substantial performance gains, both for tree-level
computations [27, 28, 29, 30, 31, 32, 33, 34, 35, 36] and for one-loop level
computations [37, 38]. In addition, the evaluation of PDFs [39] and the matrix
element method [40, 41, 42] have been explored on GPUs. Nevertheless, a
production-ready GPU-enabled event generator suitable for experimental
applications has not yet become available. The construction of such a
generator must begin with the selection of the computational algorithm which
has the best metrics in terms of flexibility and efficiency. This algorithm
should also be easy to understand and implement, such that portability
questions can be discussed not only by high-energy physicists, but by a
somewhat larger scientific computing and computer science community.

Based on many previous perturbative QCD and Beyond Standard Model
applications, we argue that such an algorithm is given by the Berends–Giele
recursion [43, 44]. Due to its inherently parallel nature, the recursion can
easily be thread-parallelized [6, 32]. In addition, modern computer
architectures with their large memory and efficient memory access allow the
trivial parallelization of the entire event loop in the Monte-Carlo
integration. The computation of color factors needed for the assembly of full
matrix elements can either be carried out in a factorized form, or embedded in
the recursion itself [45, 46]. We refer to the factorized option as a color-
ordered Berends–Giele (COBG) recursion and to the embedded one as a color-
dressed Berends–Giele (CDBG) recursion. From a computational point of view,
both have advantages and disadvantages, and which method is most suitable for
a certain hardware depends critically on the size, bandwidth and access
pattern of fast memory.

In this paper, we therefore perform a comprehensive study of the known
variants of the Berends–Giele recursion and assess their strengths and
weaknesses on various computing platforms. We investigate the scaling of
computation time and memory usage with the number of external gluons in all-
gluon amplitudes, and compare the performance of our new CPU and GPU codes to
the existing general-purpose matrix element generators MadGraph [9], Amegic
[5], and Comix [6]. We start from the thread-scalable Berends–Giele recursion
[46], which uses the leading color approximation, and extend it to full color
and helicity summed amplitudes, by using: naive sums, color sampling in the
color-flow basis [47], color dressing [45] and the continuous color and
helicity representation of [48]. For a realistic performance assessment of the
algorithms, we do not take advantage of symmetries that are specific to all-
gluon scattering, and we only take into account algorithms that generate
strictly positive weights.

This paper is organized as follows. In Sec. 2 we give a brief introduction to
multi-gluon amplitudes, and we review the Berends–Giele recursive relations,
both in their color-ordered and color-dressed form. Section 3 describes our
implementation of the various techniques in practical algorithms, bound by the
constraints of existing computing platforms. The discussion of these
algorithms is continued in Sec. 4, where we compare the efficiency of sampling
versus summing over color and helicity, and consider computational and memory
complexity. We present a comprehensive analysis of our timing studies in Sec.
5 and draw conclusions for future developments of event generators in Sec. 6.

##  2 Efficient computation of multi-gluon tree amplitudes

In this section, we recall various techniques for the efficient computation of
squared $n$-gluon matrix elements. They have a rich structure and are an
excellent proxy for more generic tree-level computations typically encountered
in collider physics. The preferred technique for evaluating matrix elements
will depend on whether it is used for analytic calculations or numerical
evaluation. Furthermore, the optimal numerical evaluation technique depends on
the specific computational hardware.

###  2.1 Summation over unobserved quantum numbers

Due to QCD confinement, it is impossible to detect spin or color of individual
partons. The effect of spin and color correlations at the parton-level only
becomes measurable through intricate correlations in many-particle hadronic
final-states. Multi-parton amplitudes summed over color and spin are therefore
of highest interest for collider phenomenology. At the same time, the
summation over unobserved quantum numbers is one of the major obstacles in the
computation of collider physics observables. We will first discuss how these
computations are performed efficiently.

We define the color and implicitly spin summed squared $n$-gluon tree-level
amplitude $\mathcal{A}_{a_{1}\ldots a_{n}}^{\mu_{1}\cdots\mu_{n}}$ as

$|\mathcal{A}(1,\ldots,n)|^{2}=\sum_{a_{1}\ldots
a_{n}}\mathcal{A}_{a_{1}\ldots
a_{n}}^{\mu_{1}\cdots\mu_{n}}(p_{1},\ldots,p_{n})\left(\mathcal{A}_{a_{1}\ldots
a_{n}}^{\nu_{1}\cdots\nu_{n}}(p_{1},\ldots,p_{n})\right)^{\dagger}\prod_{i=1}^{n}P_{\mu_{i}\nu_{i}}(p_{i},k_{i})\;,$
(1)

where each gluon with label $i$ is characterized by its color $a_{i}$ and its
momentum $p_{i}$. The projection operator $P_{\mu\nu}$ removes the
longitudinal component of the external gluons so no external ghosts particles
need to be included in the calculation. In the lightcone gauge the projection
operator is given by

$P_{\mu\nu}(p,k)=-g_{\mu\nu}+\frac{p_{\mu}k_{\nu}+p_{\nu}k_{\mu}}{p\cdot
k}\;,$ (2)

where $k$ is an arbitrary massless gauge vector. Motivated by the appearance
of high transverse momentum production of hadrons at hadron colliders and the
discovery of 3-jet production at PETRA [49, 50, 51], the first analytic
expressions for 4-gluon [52, 53] and 5-gluon [54, 55] squared matrix elements
were obtained by evaluating Eq. (1) in terms of Feynman diagrams.

To simplify the calculation and to obtain compact expressions, one can use the
dyadic decomposition of the lightcone projection operator $P_{\mu\nu}$ into
helicity eigenstates

$P_{\mu\nu}(p,k)=\sum_{\lambda}\epsilon_{\mu}^{\lambda}(p,k)\epsilon_{\nu}^{\lambda\,\dagger}(p,k)\;,$
(3)

where $\lambda=\pm$ are the helicity labels of the gluon. This was first used
in [56] to obtain a compact expression for the squared 5-gluon matrix elements
by judicious choices of the gauge vectors $k_{i}$ to simplify the calculation.
The squared $n$-gluon amplitude then takes the form

$|\mathcal{A}(1,\ldots,n)|^{2}=\sum_{\lambda_{1}\ldots\lambda_{n}}\sum_{a_{1}\ldots
a_{n}}\mathcal{A}_{a_{1}\ldots
a_{n}}^{\lambda_{1}\ldots\lambda_{n}}(p_{1},\ldots,p_{n})\left(\mathcal{A}_{a_{1}\ldots
a_{n}}^{\lambda_{1}\ldots\lambda_{n}}(p_{1},\ldots,p_{n})\right)^{\dagger}\;,$
(4)

where the $\mathcal{A}^{\lambda}$ are the contraction of the original
$\mathcal{A}^{\mu}$ with the helicity eigenstate $\epsilon_{\mu}^{\lambda}$.

The findings of multiple jet production at CERN (see e.g. [57, 58]) motivated
the calculation of 6-parton helicity amplitudes using supersymmetry relations
and spinor algebra [59, 60, 61]. The obtained expressions were suitable for
numerical evaluations, but making them compact required a refined treatment of
color. The color decomposition introduced in [62, 63] led to the definition of
color-ordered amplitudes and paved the way for more refined methods to compute
helicity amplitudes [64, 65].

To evaluate Eq. (4), we have two choices for the treatment of the helicity
and/or color sums. The traditional method is to explicitly sum over all the
quantum states. The alternative technique is to replace the discrete summation
by a continuous sampling through parametric integrations. For the spin states
of the gluon this is familiar and involves changing helicities to
polarizations. The connection between the two approaches is straightforward
[48, 66]:

$e_{\mu}(p,k,\phi)=e^{i\phi}\epsilon_{\mu}^{+}(p,k)+e^{-i\phi}\epsilon_{\mu}^{-}(p,k)\;,$
(5)

so that

$P_{\mu\nu}(p,k)=\frac{1}{2\pi}\int_{0}^{\,2\pi}{\rm
d}\phi\,e_{\mu}(p,k,\phi)e_{\nu}^{*}(p,k,\phi)\;.$ (6)

Instead of the complex polarization vector definition of Eq. (5) one can also
define a real valued polarization vector (leading to real valued amplitudes
for the gluon-only case). We can choose a unit vector $e_{1}={\rm
Re}(\epsilon^{+})/2$ orthogonal to the gluon momentum $p$ and a subsequent
unit vector $e_{2}=-{\rm Im}(\epsilon^{+})/2$ orthogonal to $p$ and $e_{1}$.
Then the continuous polarization vector is given by

$e^{\mu}(p,\phi)=\cos(\phi)\,e_{1}^{\mu}(p)+\sin(\phi)\,e_{2}^{\mu}(p)\;.$ (7)

In a similar fashion we can replace the discrete sum over external colors by a
continuous color polarization vector which is integrated over. We first
construct continuous color polarizations for the fundamental representation
based on the dyadic decomposition [48, 66]

$\delta_{ij}=\int{\rm d}[z]\,\eta_{i}([z])\eta_{j}([z])\;.$ (8)

Any real valued $\eta_{i}$ can be parametrized in terms of polar and azimuthal
angles, such that

$\eta_{i}([z])\to\eta_{i}(\theta,\phi)=\left(\begin{array}[]{c}\cos\theta\\\
\sin\theta\cos\phi\\\ \sin\theta\sin\phi\\\ \end{array}\right)\;,$ (9)

and

$\int{\rm d}[z]\to\frac{N_{c}}{4\pi}\int_{0}^{\,\pi}{\rm
d}\cos\theta\int_{0}^{\,2\pi}{\rm d}\phi\;.$ (10)

In this way the quark color state is represented by a spherical three
dimensional unit vector.

In order to construct a dyadic decomposition for the adjoint representation,
relevant for gluons, we use the identity

$\begin{split}\delta^{ab}={\rm Tr}(T^{a}T^{b})=&\;\int{\rm d}[z]{\rm
d}[\bar{z}]\,\eta_{i}([z])T^{a}_{ij}\eta_{j}([\bar{z}])\,\eta_{k}([\bar{z}])T^{b}_{kl}\eta_{l}([z])\\\
=&\;\int{\rm d}[z]{\rm
d}[\bar{z}]\,\eta^{a}([z],[\bar{z}])\eta^{b\,\dagger}([z],[\bar{z}])\;,\end{split}$
(11)

where we have defined the gluon color polarization vector

$\eta^{a}([z],[\bar{z}])=\eta_{i}([z])\,T^{a}_{ij}\,\eta_{j}([\bar{z}])\;.$
(12)

Note that this construction differs from [48, 66], in that it uses only four
instead of five integration variables and the gluon color polarization is
represented by two real 3-dimensional spherical unit vectors.

Using the dyadic decompositions above, we can write the summed squared matrix
element in the following form:

$\left|\mathcal{A}(1,\ldots,n)\right|^{2}=\prod_{i}\sum\int\,\frac{{\rm
d}\phi_{i}}{2\pi}\,{\rm d}[z]_{i}\,{\rm
d}[\bar{z}]_{i}\;\left|\mathcal{A}_{[z]_{1}[\bar{z}]_{1}\ldots[z]_{n}[\bar{z}]_{n}}^{\phi_{1}\ldots\phi_{n}}(p_{1},\ldots,p_{n})\right|^{2}\;,$
(13)

where we have defined the color-helicity sub-amplitudes

$\begin{split}\mathcal{A}_{[z]_{1}[\bar{z}]_{1}\ldots[z]_{n}[\bar{z}]_{n}}^{\phi_{1}\ldots\phi_{n}}(p_{1},\ldots,p_{n})=\mathcal{A}_{a_{1}\ldots
a_{n}}^{\mu_{1}\cdots\mu_{n}}(p_{1},\ldots,p_{n})\times\prod_{i}e_{\mu_{i}}(p_{i},k_{i},\phi_{i})\times\prod_{j}\eta_{a_{i}}([z]_{i})\;.\end{split}$
(14)

They can be further simplified analytically, or evaluated numerically with the
help of the Berends–Giele recursion. Furthermore, the generation of continuous
(color) polarizations can easily be included in the Monte-Carlo integration
over phase-space. We will investigate these questions in more detail in Secs.
4 and 5.

###  2.2 Color decomposition and color-ordered amplitudes

The decomposition of summed $n$-gluon amplitudes into helicity amplitudes
$\mathcal{A}^{\lambda_{1}\ldots\lambda_{n}}$ introduced in Eq. (1) is
particularly useful in analytic calculations, because it allows for a
judicious choice of gauge vectors, which in turn can greatly simplify the
calculation [56]. Similarly, color decompositions are useful, because they
allow to factorize the color and kinematics dependent part of $n$-gluon
amplitudes and calculate the color-dependent coefficients once and for all.
This can be beneficial for both analytical and numerical evaluation. In this
subsection, we will discuss color decompositions that are particularly useful
for fast numerical computation in Monte-Carlo programs.

The most intuitive color decomposition of an $n$-gluon amplitude
$\mathcal{A}_{n}$ is based on the adjoint representation of ${\rm SU}(3)$ [62,
67, 68] resulting in

$\mathcal{A}_{a_{1}\ldots
a_{n}}^{\lambda_{1}\ldots\lambda_{n}}(p_{1},\ldots,p_{n})=\sum\limits_{\vec{\sigma}\in
S_{n-2}}(F^{a_{\sigma_{2}}}\ldots
F^{a_{\sigma_{n-1}}})_{a_{1}a_{n}}\;A^{\lambda_{1}\ldots\lambda_{n}}(p_{1},p_{\sigma_{2}},\ldots,p_{\sigma_{n-1}},p_{n})\;,$
(15)

where $F^{a}_{bc}=if^{abc}$. The functions $A$ are called color-ordered or
partial amplitudes. If they carry a helicity label they are often simply
referred to as helicity amplitudes. The multi-index $\vec{\sigma}$ runs over
all permutations $S_{n-2}$ of the ($n-2$) gluon indices $2\ldots n-1$.

Performing an explicit sum over colors in the squared amplitude, Eq. (1),
leads to

$\begin{split}|\mathcal{A}(1,\ldots,n)|^{2}=&\;\sum_{\lambda_{1}\ldots\lambda_{n}}\sum_{a_{1}\ldots
a_{n}}\mathcal{A}_{a_{1}\ldots
a_{n}}^{\lambda_{1}\ldots\lambda_{n}}(p_{1},\ldots,p_{n})\mathcal{A}_{a_{1}\ldots
a_{n}}^{\lambda_{1}\ldots\lambda_{n}}(p_{1},\ldots,p_{n})^{\dagger}\\\
=&\;\sum_{\vec{\sigma},\vec{\sigma}^{\prime}\in
S_{n-2}}\mathcal{C}_{\vec{\sigma}\vec{\sigma}^{\prime}}\sum_{\lambda_{1}\ldots\lambda_{n}}A^{\lambda_{1}\lambda_{\sigma_{2}}\ldots\lambda_{\sigma_{n-1}}\lambda_{n}}(p_{1},p_{\sigma_{2}},\ldots,p_{\sigma_{n-1}},p_{n})\\\
&\;\qquad\times
A^{\lambda_{1}\lambda_{\sigma^{\prime}_{2}}\ldots\lambda_{\sigma^{\prime}_{n-1}}\lambda_{n}}(p_{1},p_{\sigma^{\prime}_{2}},\ldots,p_{\sigma^{\prime}_{n-1}},p_{n})^{\dagger}\;,\end{split}$
(16)

where the $(n-2)!^{\,2}$ color coefficients,
$\mathcal{C}_{\vec{\sigma}\vec{\sigma}^{\prime}}$, are given by

$\mathcal{C}_{\vec{\sigma}\vec{\sigma}^{\prime}}=\sum_{a_{1}\ldots
a_{n}}(F^{a_{\sigma_{2}}}\ldots
F^{a_{\sigma_{n-1}}})_{a_{1}a_{n}}(F^{a_{\sigma^{\prime}_{2}}}\ldots
F^{a_{\sigma^{\prime}_{n-1}}})^{*}_{a_{1}a_{n}}\;.$ (17)

An alternative color decomposition can be obtained by using the definition of
structure constants in terms of ${\rm SU}(3)$ generators [63]

$\mathcal{A}_{a_{1}\ldots
a_{n}}^{\lambda_{1}\ldots\lambda_{n}}(p_{1},\ldots,p_{n})=\sum\limits_{\vec{\sigma}\in
S_{n-1}}{\rm Tr}(T^{a_{1}}T^{a_{\sigma_{2}}}\ldots
T^{a_{\sigma_{n}}})\;A^{\lambda_{1}\ldots\lambda_{n}}(p_{1},p_{\sigma_{2}},\ldots,p_{\sigma_{n}})\;.$
(18)

In this case the sum runs over all permutations $S_{n-1}$ of the ($n-1$) gluon
indices $2\ldots n$, leading to a substantial increase in the number of
partial amplitudes that contribute to the full color-helicity amplitude.
Equation (16) for the explicit sum is correspondingly changed by replacing
$S_{n-2}\to S_{n-1}$, $\lambda_{n}\to\lambda_{\sigma_{n}}$, $p_{n}\to
p_{\sigma_{n}}$ and by making the corresponding replacements in the conjugate
amplitude and the color factors, Eq. (17). Note that the minimal number of
partial amplitudes needed for a complete evaluation is still $(n-2)!$ [69].
The remaining partial amplitudes can be obtained from the calculated partial
amplitudes by a set of linear equations. The larger growth in the number of
color coefficients and the need to construct the remaining $(n-1)$ partial
amplitudes disfavors Eq. (18) for numerical computations and we will therefore
only use Eq. (15).

A third color decomposition, suited especially for Monte-Carlo event
generation, is the color-flow decomposition [47]. In this case, the gluon
vector field is treated as an $N_{c}\times N_{c}$ matrix, $A^{\mu}_{ij}$,
rather than a field with one color index, $A^{\mu}_{a}$. This leads to the
amplitudes

$\begin{split}\mathcal{A}_{i_{1}j_{1}\ldots
i_{n}j_{n}}^{\lambda_{1}\ldots\lambda_{n}}(p_{1},\ldots,p_{n})=&\;\prod_{k=1}^{n}T^{\,a_{k}}_{i_{k}j_{k}}\mathcal{A}_{a_{1}\ldots
a_{n}}^{\lambda_{1}\ldots\lambda_{n}}(p_{1},\ldots,p_{n})\\\
=&\;\sum\limits_{\vec{\sigma}\in
S_{n-1}}\delta^{i_{1}\bar{\jmath}_{\sigma_{2}}}\delta^{i_{\sigma_{2}}\bar{\jmath}_{\sigma_{3}}}\ldots\delta^{i_{\sigma_{n}}\bar{\jmath}_{1}}\;A^{\lambda_{1}\ldots\lambda_{n}}(p_{1},p_{\sigma_{2}},\ldots,p_{\sigma_{n}})\;.\end{split}$
(19)

The maximum number of partial amplitudes to be evaluated for a generic color
configuration, $i_{1}j_{1}\ldots i_{n}j_{n}$, can be as large as in Eq. (18).
However, it can be shown that – when combined with color sampling – the color-
flow decomposition yields the lowest average number of partial amplitudes to
be evaluated per Monte-Carlo event at high parton multiplicity [47].

###  2.3 Recursive computation of (color-)helicity amplitudes

The computation of the color-ordered amplitudes can be accomplished using the
Berends–Giele recursion [43], which corresponds to a dynamic programming
technique that efficiently caches sets of Feynman diagrams with at least one
common propagator. This propagator is dubbed the off-shell current and is
computed recursively as

$\begin{split}&J_{\mu}(1,2,\ldots,n)=\frac{-ig_{\mu\nu}}{p_{1,n}^{2}}\left\\{\;\vphantom{\sum_{k=j+1}^{n-1}}\sum_{k=1}^{n-1}V_{3}^{\nu\kappa\lambda}(p_{1,k},p_{k+1,n})J_{\kappa}(1,\ldots,k)J_{\lambda}(k+1,\ldots,n)\right.\\\
&\qquad\qquad+\left.\sum_{j=1}^{n-2}\sum_{k=j+1}^{n-1}V_{4}^{\nu\rho\kappa\lambda}J_{\rho}(1,\ldots,j)J_{\kappa}(j+1,\ldots,k)J_{\lambda}(k+1,\ldots,n)\right\\}\;.\end{split}$
(20)

Here $p_{i}$ denote the momenta of the gluons, $p_{i,j}=p_{i}+\ldots+p_{j}$
and $V_{3}^{\nu\kappa\lambda}$ and $V_{4}^{\nu\rho\kappa\lambda}$ are the
color-ordered three- and four-gluon vertices defined by

$\begin{split}V_{3}^{\nu\kappa\lambda}(p,q)&=i\,\frac{g_{s}}{\sqrt{2}}\big{(}\,g^{\kappa\lambda}(p-q)^{\nu}+g^{\lambda\nu}(2q+p)^{\kappa}-g^{\nu\kappa}(2p+q)^{\lambda}\,\big{)}\;,\\\
V_{4}^{\nu\rho\kappa\lambda}&=i\,\frac{g_{s}^{2}}{2}\big{(}\,2g^{\nu\kappa}g^{\rho\lambda}-g^{\nu\rho}g^{\kappa\lambda}-g^{\nu\lambda}g^{\rho\kappa}\,\big{)}\;.\end{split}$
(21)

The external particle currents, $J_{\mu}(i)$, are given by the helicity or
polarization vectors of Sec. 2.1. A complete color-ordered $n$-gluon amplitude
$A(1,\ldots,n)$ is obtained by putting the helicity/polarization dependent
($n-1$)-particle off-shell current $J_{\mu}(1,\ldots,n-1)$ on-shell and
contracting it with the external polarization $J_{\mu}(n)$:

$A(1,\ldots,n)=J_{\mu}(n)\,p_{1,n}^{2}\,J^{\mu}(1,\ldots,n-1)\;.$ (22)

The algorithmic complexity scales as $\mathcal{O}(n^{4})$ and is dictated by
the four-gluon vertex. This can be further improved by decomposing the four-
gluon vertex and introducing an auxiliary antisymmetric tensor field with the
“propagator” [48]

$-iD_{\mu\nu}^{\,\kappa\lambda}=-i\big{(}g_{\mu}^{\kappa}g_{\nu}^{\lambda}-g_{\mu}^{\lambda}g_{\nu}^{\kappa}\big{)}\;.$
(23)

Then the recursive relations for the gluon and tensor fields read

$J_{\mu}(1,2,\ldots,n)=\frac{-ig_{\mu\nu}}{p_{1,n}^{2}}\sum_{k=1}^{n-1}\left\\{\vphantom{\sum_{k=1}^{n-1}}\;V_{3}^{\nu\kappa\lambda}(p_{1,k},p_{k+1,n})J_{\kappa}(1,\ldots,k)J_{\lambda}(k+1,\ldots,n)\right.\\\
{}+V_{T}^{\nu\kappa\alpha\beta}J_{\kappa}(1,\ldots,k)J_{\alpha\beta}(k+1,\ldots,n)+\left.V_{T}^{\lambda\nu\alpha\beta}J_{\alpha\beta}(1,\ldots,k)J_{\lambda}(k+1,\ldots,n)\vphantom{\sum_{k=1}^{n-1}}\right\\}\;,$
(24)

and

$J^{\alpha\beta}(1,2,\ldots,n)=-iD^{\,\alpha\beta}_{\gamma\delta}\,\sum_{k=1}^{n-1}V_{T}^{\gamma\delta\kappa\lambda}J_{\kappa}(1,\ldots,k)J_{\lambda}(k+1,\ldots,n)\;,$
(25)

where the tensor-gluon interaction is given by

$V_{T}^{\mu\nu\kappa\lambda}=\frac{i}{2}\frac{g_{s}}{\sqrt{2}}\big{(}g^{\mu\kappa}g^{\nu\lambda}-g^{\mu\lambda}g^{\nu\kappa}\big{)}\;.$
(26)

This will improve the complexity scaling of the algorithmic implementation to
$\mathcal{O}(n^{3})$.

The above recursion can be modified to calculate amplitudes without the need
for an explicit color decomposition by using color-dressed currents instead.
The most convenient form for implementing the color-dressed recursion and the
techniques outlined in Sec. 2.1 is given by the color-flow basis [45], which
works for both the discrete color sampling strategy of [47] and the continuous
sampling of [48, 66]. It is based on the use of the identity $\delta^{ab}={\rm
Tr}(T^{a}T^{b})$ for the gluon propagator. By writing ${\rm
Tr}(T^{a}T^{b})=T^{a}_{ij}T^{b}_{ji}$ and assigning color indices $i$ and $j$
to the intermediate gluon, one can use the generators $T^{a}_{ij}$ to project
the structure constants associated with the elementary interaction vertices
onto fundamental indices as

$F^{a}_{bc}T^{a}_{ij}T^{b}_{kl}T^{c}_{mn}=if^{abc}T^{a}_{ij}T^{b}_{kl}T^{c}_{mn}=\delta_{il}\delta_{kn}\delta_{mj}-\delta_{in}\delta_{kj}\delta_{ml}\;.$
(27)

Note that we have used the conventions of [64] for the normalization of the
${\rm SU}(3)$ generators. Formally, we define the color-dressed gluon and
tensor pseudoparticle currents $\mathcal{J}_{\mu\,ij}$ and
$\mathcal{J}_{\alpha\beta\,ij}$ as

$\begin{split}\mathcal{J}_{\mu\,ij}(1,\ldots,n)&=\sum\limits_{\vec{\sigma}\in
S_{n}}\delta_{ij_{\sigma_{1}}}\delta_{i_{\sigma_{1}}j_{\sigma_{2}}}\ldots\delta_{i_{\sigma_{n}}j}\;J_{\mu}(\sigma_{1},\ldots,\sigma_{n})\;,\\\
\mathcal{J}_{\alpha\beta\,ij}(1,\ldots,n)&=\sum\limits_{\vec{\sigma}\in
S_{n}}\delta_{ij_{\sigma_{1}}}\delta_{i_{\sigma_{1}}j_{\sigma_{2}}}\ldots\delta_{i_{\sigma_{n}}j}\;J_{\alpha\beta}(\sigma_{1},\ldots,\sigma_{n})\;.\end{split}$
(28)

Denoting by $\pi$ the set $(1,\ldots,n)$ of $n$ particle indices, the
following recursive relations for these currents are obtained:

$\begin{split}\mathcal{J}_{\mu\,ij}(\pi)&=D_{\mu\,ij}^{\,\nu\,hg}(\pi)\left\\{\;\sum_{\mathcal{P}_{2}(\pi)}\mathcal{V}_{\nu\,hg}^{\,\kappa\,kl,\,\lambda\,mn}(\pi_{1},\pi_{2})\,\mathcal{J}_{\kappa\,kl}(\pi_{1})\mathcal{J}_{\lambda\,mn}(\pi_{2})\right.\\\
&\text{\hskip
77.49976pt}+\left.\sum_{\mathcal{OP}_{2}(\pi)}\mathcal{V}_{\nu\,hg}^{\,\kappa\,kl,\,\alpha\beta\,mn}\,\mathcal{J}_{\kappa\,kl}(\pi_{1})\mathcal{J}_{\alpha\beta\,mn}(\pi_{2})\;\right\\}\;,\\\\[10.00002pt]
\mathcal{J}_{\alpha\beta\,ij}(\pi)&=D_{\alpha\beta\,ij}^{\,\gamma\delta\,hg}\,\sum_{\mathcal{P}_{2}(\pi)}\mathcal{V}_{\gamma\delta\,hg}^{\,\kappa\,kl,\,\lambda\,mn}\,\mathcal{J}_{\kappa\,kl}(\pi_{1})\mathcal{J}_{\lambda\,mn}(\pi_{2})\;.\end{split}$
(29)

Here we have defined the color-dressed gluon and tensor pseudoparticle
vertices

$\begin{split}\mathcal{V}_{\nu\,hg}^{\,\kappa\,kl,\,\lambda\,mn}(\pi_{1},\pi_{2})=\delta_{lg}\delta_{kn}\delta_{hm}\,V_{3\,\nu}^{\;\;\,\kappa\lambda}(\pi_{1},\pi_{2})+\delta_{hk}\delta_{ml}\delta_{ng}\,V_{3\,\nu}^{\;\;\,\lambda\kappa}(\pi_{2},\pi_{1})\end{split}$
(30)

and

$\begin{split}\mathcal{V}_{\gamma\delta\,hg}^{\,\kappa\,kl,\,\lambda\,mn}=\delta_{lg}\delta_{kn}\delta_{hm}\,V_{T\,\gamma\delta}^{\;\;\,\kappa\lambda}+\delta_{hk}\delta_{ml}\delta_{ng}\,V_{T\,\gamma\delta}^{\;\;\,\lambda\kappa}\;,\end{split}$
(31)

as well as the dressed propagators

${D_{\mu}^{\nu}}_{ij}^{hg}=D_{\mu}^{\nu}\,\delta_{ih}\delta_{jg}\;,\qquad\text{and}\qquad{D_{\alpha\beta}^{\gamma\delta}}_{ij}^{hg}={D_{\alpha\beta}^{\gamma\delta}}\,\delta_{ih}\delta_{jg}\;.$
(32)

The first sum in Eq. (29) runs over all partitions $\mathcal{P}_{2}(\pi)$ of
the set $\pi$ into two subsets $\pi_{1,2}$ while the second sum runs over the
set of ordered partitions $\mathcal{OP}_{2}(\pi)$ into $\pi_{1,2}$.111Our
implementation of the sum over partitions in Eq. (29) uses the equivalence of
Eq. (3.18) and Eq. (3.16) in Ref. [45] to reduce the computational complexity.

##  3 Practical implementation of the algorithms

In this section we discuss different implementations of the methods laid out
in Sec. 2 in actual computer code. In most cases, we provide a CPU and a GPU
version, which are analyzed in terms of computational complexity and memory
requirements in Sec. 4 and compared in terms of practical compute performance
in Sec. 5. We use Rambo [70] to generate the phase-space points for the
different methods. Obtaining a purely real algorithm reduces the memory
requirements by a factor of two compared to a complex valued algorithm. In
memory bound algorithms, reducing the memory of each object allows the CPU and
GPU to fetch more from the cache in a single read, thus improving overall
performance. Therefore, it is ideal to obtain a purely real implementation.

###  3.1 Leading color computation

The first algorithm (Tess) was originally published in [32]:

  * •

Helicity amplitudes are computed using the color-ordered Berends-Giele
recursion, Eqs. (20) and (22). To avoid the larger memory requirements Tess
does not make use of the decomposed four-gluon vertex in Eq. (24). Helicities
are evaluated by sampling real polarization vectors according to Eq. (7). As a
consequence, _all_ numbers in the algorithm are real.

  * •

Parallelization is achieved not only by calculating many events concurrently,
but also by making use of the intrinsically parallel structure of the
recursion, such that up to $n-1$ threads cooperate in the calculation of a
single event.

  * •

Global memory access is minimized by storing only event weights and gluon
momenta. The shared memory per event contains all the currents and momenta
needed to compute the color-ordered amplitudes and scales quadratically with
the number of particles.

  * •

In contrast to the original publication [32], the implementation has been
modified to use double precision numbers, as all of the following algorithms
do.

The second algorithm (BlockGen-LC) is similar to Tess, but differs in the
following details:

  * •

Only event-level parallelization is used, i.e. every threads computes the
matrix element for an independent phase-space point.

  * •

All quantities except for the internal currents of partial amplitudes are
stored in the global memory of the GPU.

  * •

Helicities can be either sampled or summed.

Both of the above algorithms rely on the leading color approximation. The
number of colors $N_{c}$ is considered to be large, such that the color sum in
Eq. (16) can be evaluated analytically. It results in a simple prefactor:

$|\mathcal{A}_{n}(1,\ldots,n)|^{2}=N_{c}^{n-2}\left(N_{c}^{2}-1\right)\left(\sum_{\vec{\sigma}\in
S_{n-1}}\left|A^{\lambda_{1}\ldots\lambda_{n}}(p_{1},p_{\sigma_{2}},\ldots,p_{\sigma_{n}})\right|^{2}+\mathcal{O}\left(\frac{1}{N_{c}^{2}}\right)\right).$
(33)

###  3.2 Color summation and sampling with color-ordered amplitudes

The third and fourth algorithm (BlockGen-COΣ and BlockGen-CO${}_{\text{MC}}$)
use the following approach:

  * •

Helicity amplitudes are computed by utilizing the color-ordered Berends–Giele
recursion, Eqs. (20) and (22). The two BlockGen-CO algorithms do not make use
of the decomposed four-gluon vertex in Eq. (24), to avoid the larger memory
requirements. Helicities are summed using real polarization vectors according
to Eq. (7). As a consequence, _all_ numbers in the algorithms are real.

  * •

In the _color summed_ variant (BlockGen-COΣ), the adjoint representation
decomposition from Eq. (15) is used, which leads to the amplitude assembly
formula, Eq. (16). In the _color sampling_ variant (BlockGen-
CO${}_{\text{MC}}$), the color-flow decomposition of Eq. (19) is used and
color indices are generated event by event using the algorithm outlined in
[6].

  * •

All quantities except for the internal currents of partial amplitudes are
stored in the global memory of the GPU. For the BlockGen-COΣ variant, this
includes the color matrix, Eq. (17), which we compute once using Form [71] and
then read in from storage to perform the matrix element evaluations.

###  3.3 Color sampling with color-dressed amplitudes

The fifth algorithm (BlockGen-CD${}_{\text{MC}}$) uses the following approach:

  * •

Helicity amplitudes are computed using the color-dressed Berends–Giele
recursion, Eq. (29), and in particular making use of the continuous color-
polarizations introduced in Eq. (12). Helicities are sampled using Eq. (5),
which is equivalent to Eq. (7). As a consequence, _all_ numbers in the
algorithm are real.

  * •

All quantities are stored in global memory.

##  4 Computational complexity and memory requirements

In this section, we discuss the characteristics of the various algorithms
introduced in Sec. 3 in terms of computational complexity and memory
requirements for practical implementations of Monte-Carlo programs, which can
be used for parton-level event generation for collider experiments on GPUs. We
will focus on four criteria in particular: i) level of parallelization, ii)
summing vs. sampling colors and helicities, iii) the scaling behavior of
explicit summation and iv) memory requirements. This guides us in our
selection of a set of candidate algorithms that will be studied in more detail
in Sec. 5. Note that unless explicitly mentioned, the timings are measured for
an event, meaning that the time to generate a phase-space point and evaluating
cuts is included. As can be seen in Fig. 1, this is only relevant for the
lower multiplicities and becomes negligible for the higher ones.

###  4.1 Level of parallelization

The GPU implementation of the Tess algorithm is not only parallelized at the
event level, but also at the lower level of the summation over currents in the
recursive amplitude computation [32]. This corresponds to a parallel
computation of the (outer) sums in Eq. (20), such that up to $n-1$ threads can
cooperate. The low-level parallelism can be helpful in a memory bound
algorithm with several threads operating on the same data, which can be stored
on the shared memory of the GPU, a fast cache managed by the program itself.
The downside is a more complex implementation, due to the data sharing and the
combination of different levels of parallelization.

Modern HPC GPU combine shared memory with the chip-controlled L2 cache, and
have a larger throughput for loading missing data from their global memory.
Therefore, the speed gains from the additional level of parallelization might
be diminished. We analyze this effect in Fig. 1 which shows the computational
time per event comparing Tess (orange) to BlockGen-LC (blue). For high gluon
multiplicities, the two implementations perform similarly. The different
scaling of BlockGen-LC at low multiplicities arises from the slightly more
generic handling of phase-space generation and cuts in BlockGen-LC and is
removed when we adopt the same strategy as Tess. This is shown in dashed blue.
For this reason, we only consider pure event-level parallelization in all the
following implementations, avoiding unnecessary complexity. It is interesting
to note the large timing improvement in Fig. 1 when comparing the original
results from Tess (dashed orange) to the current results (orange). This
improvement is purely due to improvements in GPU hardware, and indicates the
increasing importance of modern GPUs for collider physics applications.

Figure 1: A comparison of the runtimes of Tess and BlockGen-LC, measured on a
NVIDIA V100 ($16\text{\,}\mathrm{G}\mathrm{B}$ global memory, 5,120 CUDA
cores, $6144\text{\,}\mathrm{K}\mathrm{B}$ L2 cache). To make a more fair
comparison to Tess, we include the timing for BlockGen-LC only including the
time to evaluate the matrix element (“ME Only”). The results for Tess from the
original publication in 2010 [32] are included as a reference to show the GPU
improvements.

###  4.2 Efficiency of color and helicity sampling

Figure 2: A comparison of the overhead for color, helicity, and combined
color-helicity sampling as a function of gluon multiplicity. The ratio of the
number of points needed to reach a target precision of 1% is shown relative to
a calculation with an explicit sum over both colors and helicities.

Figure 2 shows the overhead for color sampling, helicity sampling, and
combined color and helicity sampling as a function of gluon multiplicity in a
realistic parton-level calculation. For this, we use the recursive phase-space
generator of Comix. Other details of the setup are the same as in Sec. 5. We
require a target precision of 1%. We measure the required number of points
needed to attain this precision and plot the ratio of the number of points for
the various sampling methods to the explicit sum. We observe that the relative
overhead of helicity and color sampling individually is approximately constant
as a function of multiplicity. However, the overhead of combined color and
helicity sampling increases with increasing multiplicity. This can be
understood by looking at the analytic structure of the $n$-gluon amplitudes.
For $2$ and $3$ gluon final states, only MHV amplitudes exist. The numerator
of a mostly plus amplitude is given as the fourth power of the product of the
spinors for the two gluons with negative helicity. Therefore, helicity
configurations can be sampled efficiently by drawing negative helicity pairs
with probability proportional to the magnitude of the two-particle invariants.
For higher multiplicites, the number of MHV configurations is still larger
than the number of NMHV configurations, but the sampling according to two-
particle invariants does not always approximate the true structure of the
matrix element. This effect seems less pronounced in the color summed case, as
can be seen in Fig. 2. However, when colors are sampled, an imperfect
selection of the helicity configuration does create large fluctuations in
$4$\- and $5$-gluon final states (where MHV and NMHV configurations are
present), and yet larger fluctuations in $6$-gluon final states (where MHV,
NMHV and NNMHV configurations exist).222 We were not able to fully confirm the
expected behavior in $7$-gluon final states, because the corresponding color-
summed prediction could not be determined within our computing budget.
However, the ratio between color-helicity summed and color summed only
prediction is indeed the same as in the case of the $6$-gluon final state. The
reason for this is that the numerator and denominator of the partial
amplitudes are uncorrelated, thus an uncorrelated and imperfect sampling of
color and helicity leads to reduced efficiency. Nevertheless, the timing
improvements in the matrix element computation that arise from the use of
sampling algorithms will never be overcompensated by a reduced convergence in
the eventual Monte-Carlo integration. This point is crucial when making a
choice about which algorithm to implement for GPU computing, and in particular
it enables us to choose a memory-lean option.

###  4.3 Scaling of color and helicity sums

When evaluating the sum over pairs of permutations in Eq. (16) (for a given
helicity configuration), the $(n-2)!$ matrix elements needed should first be
calculated and stored. Then the $(n-2)!((n-2)!+1)/2$ independent summands can
be calculated without re-evaluating the same matrix elements over and over.
The time needed to evaluate each summand is then dominated by loading the data
corresponding to $\mathcal{C}_{\vec{\sigma}\vec{\sigma}^{\prime}}$,
$A^{\vec{\sigma}}$ and $A^{\vec{\sigma}^{\prime}}$, and should hence be
individually very small compared to the evaluation of a partial amplitude via
the recursion relations in Eq. (20). However, the different scalings of the
two operations, factorial for precalculating partial amplitudes vs. factorial
squared for evaluating the color sum, have the consequence that there will be
an $n_{\text{out}}$ for which the time required for the summation will
eventually become dominant.

Figure 3 compares the time needed per event for each operation. The crossover
point is found to be between $n_{\text{out}}=6$ and $7$. The inherent
$\mathcal{O}((n-2)!^{\,2})$ scaling of the color-summed algorithm means that a
color-sampled algorithm will eventually become more efficient. However, since
the contribution of the summing itself becomes relevant only beyond
$n_{\text{out}}=6$, this behavior might not affect practically relevant
computations.

Figure 3: A comparison of the time needed to evaluate the required amplitudes
versus the time needed to sum over all the permutations required as a function
of multiplicity (Note that the absolute values are not directly comparable to
Fig. 6, since we have used helicity sampling here.)

###  4.4 Memory requirements of summed vs. sampled algorithms

Figure 4 compares the different memory requirements on the GPU for the three
different full-color algorithms BlockGen-COΣ, and BlockGen-CD${}_{\text{MC}}$.
It also includes leading color results from Tess and BlockGen-LC. The upper
plot displays the global memory usage which is not event-specific, the middle
plot shows the global memory usage per event, and the bottom plot shows the
shared memory usage per event.

Figure 4: Comparison of memory usage in different GPU implementations. The
upper plot displays the global memory required that is allocated independently
of the number of events/threads. The middle plot shows the global memory per
event and the bottom plot shows the shared memory per event. The shared memory
per event for BlockGen-LC and BlockGen-COΣ is identical. Figure 5: Single
threaded CPU heap memory usage for the various algorithms (Note that the
MadGraph executable size is used, since the color matrix, etc. is compiled
into the executable. Therefore, the MadGraph dynamic heap memory allocations
are negligible).

The Tess algorithm uses global memory only to store event weights and gluon
momenta, which gives a small footprint in the per-event global memory plot.
Everything else is stored in shared memory.

In the BlockGen-COΣ algorithm, event-independent global memory is used in
particular to store the color matrix, Eq. (17), leading to a rapid
$\mathcal{O}((n-2)!^{\,2})$ growth. The per-event global memory usage is
dominated by the $\mathcal{O}((n-2)!)$ growth of the number of color-ordered
amplitudes $A(1,{\sigma_{2}},\ldots,{\sigma_{n-1}},n)$. Shared memory is used
only for the intrinsic momenta in the recursion, such that its use grows
linearly with the number of particles.

For BlockGen-CD${}_{\text{MC}}$, the per-event global memory use is governed
by the number of complex valued currents and tensors. Since we only need to
compute the currents for ordered sets of $\pi$ in Eq. (29), the total number
of objects scales as $\mathcal{O}(2^{n})$. Shared memory is not used by the
color-dressed approach.

We first note that the color-summed algorithm is limited to
$n_{\text{out}}\leq 8$. Beyond that the memory required to store the color
matrix will exceed $\mathcal{O}($10\text{\,}\mathrm{GB}$)$, and hence will not
fit any longer into the global memory of the GPU.333 At the expense of
relabeling permutations we could minimize the size of the storage required for
the color matrix and obtain a scaling of $\mathcal{O}((n-2)!)$, down from
$\mathcal{O}((n-2)!^{\,2})$. However, this technique would be beneficial only
for large final-state multiplicities which are generally better to compute in
a color-dressed approach due to the improved scaling. We therefore refrain
from using it.  However, for $n_{\text{out}}\lesssim 7$ the memory use of the
color-summed algorithm is below that of the color-dressed algorithm. Given
that the algorithms are memory bound, we expect the color-summed algorithm to
perform better than the color-dressed one for $n_{\text{out}}\lesssim 7$.

Figure 5 compares the single-thread memory usage of CPU-based implementations
of our algorithms BlockGen-COΣ, BlockGen-CO${}_{\text{MC}}$ and BlockGen-LC,
with the widely used MadGraph, Comix, and Amegic codes. For MadGraph, we use
version 2.9.2, while Comix and Amegic have been run as part of the Sherpa
framework in its version 2.3 series [72]. For MadGraph, we use the standalone
mode, supplemented with a custom Monte-Carlo loop. For Comix, we use its color
summing and its color sampling mode, while for Amegic, we show the default
diagram-based mode, and a mode where analytic CSW rules are employed [73],
which have been derived in [74]. We plot the RAM memory use of all codes
except for MadGraph, where the code generator embeds the color matrices
directly within the source code, such that the memory usage is nearly
exclusively driven by the size of the compiled executable, which we therefore
plot instead in this case. We find that the Comix modes and our custom
algorithms have the smallest memory footprints for multiplicities of
$n_{\text{out}}=5$ and beyond, with Comix in its summed and its sampled mode
using about 10 and 25 times more memory than the respective BlockGen-CO. This
is likely due to Comix storing all contributing sub-currents of Eq. (28), to
reduce evaluation time [6].

##  5 Performance of event generation

We now compare the algorithms described in the previous sections in terms of
their practical performance.444 CUDA Nsight profiling reports can be obtained
from the authors upon request. Of particular interest is the time needed to
generate a single event. Where appropriate, we include commonly used tools in
the comparison, such as MadGraph [9], Amegic [5] and Comix [6] (the versions
used are listed in Sec. 4.4).

We begin by studying single-thread performance on CPU in Sec. 5.1, proceed
with the massively parallelized calculation on GPU in Sec. 5.2, and conclude
with a realistic chip-to-chip comparison of (multi-threaded) CPU and GPU event
generation in Sec. 5.3, in order to find the most promising hardware/algorithm
combination.

For all studies, we use a partonic center-of-mass energy of
$\sqrt{s}=$14\text{\,}\mathrm{TeV}$$ and require the gluon momenta to satisfy
the kinematic constraints

$p_{T}>20~{}{\rm GeV}\;,\qquad|\eta|<2.5\;,\qquad\Delta R>0.4\;.$

A parton density function is not used. For simplicity, we set the strong
coupling to $\alpha_{S}(m_{Z})=0.118$ and use a fixed renormalization scale of
$\mu_{R}=m_{Z}$.555 Even though we only study gluon amplitudes, none of our
algorithms makes use of the process-specific symmetries relating different
helicity amplitudes. We do so in order to reflect the performance that can be
expected from a generic, automated implementation of an event generator for
arbitrary interaction model Lagrangians.

###  5.1 Comparison of single-threaded algorithms

Figure 6: The timings for various CPU-based algorithms run on a single thread
are compared as a function of the gluon multiplicity. The results were all
generated on an Intel${}^{\text{\textregistered}}$
Xeon${}^{\text{\textregistered}}$ E5-2650 v2 8-core CPU
($2.60\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$,
$20\text{\,}\mathrm{M}\mathrm{B}$ cache).

To study the evaluation time per event of the algorithms without
parallelization, we compare in Figure 6 the CPU variants of our algorithms
BlockGen-COΣ, BlockGen-CO${}_{\text{MC}}$ and BlockGen-LC, with the widely
used public codes MadGraph, Comix, and Amegic ME generators, and with Tess. As
described in Sec. 4.4, Comix is used in color summing and sampling mode, while
Amegic is using either its default diagrammatic mode or analytic CSW rules.
Figure 6 shows the evaluation time per event for these different algorithmic
choices, as well as for the public codes. We find that the leading color
implementations, BlockGen-LC and Tess, perform similarly. For the full-color
algorithms, ratio plots show the evaluation times relative to BlockGen-COΣ and
BlockGen-CO${}_{\text{MC}}$. Note that the implementation of BlockGen-COΣ and
BlockGen-CO${}_{\text{MC}}$ is based on purely real numbers, cf. Sec. 3. A
certain overhead compared to the public codes, which work with complex
numbers, is therefore to be expected.

In both cases, Comix has the best scaling behavior due to the reuse of a
maximal number of precomputed sub-currents. This reduces the naively expected
factorial growth in computing time to an exponential growth, but comes at the
cost of a larger memory footprint (see Sec. 4.4). Note that the recycling of
sub-currents is not a practical option for the color-ordered algorithms. While
the color-dressed algorithm only requires to store one current per unordered
set of indices, the color-ordered algorithms (both in the color summed and
color sampled variant) would require to store ordered sets of indices. As the
number of ordered sets grows factorially with the number of members in the
set, the memory required for this technique would become prohibitively large
at large multiplicity. The resulting different scaling behaviors are clearly
visible in Fig. 6. Hence, Comix’ evaluation times are smallest for large
multiplicities, falling below BlockGen-COΣ at ${n_{\text{out}}\approx 7}$ and
below BlockGen-CO${}_{\text{MC}}$ at $n_{\text{out}}\approx 5$, for summing
and sampling, respectively.

Taking Fig. 2 as a ballpark reference for how much a color sampled algorithm
might be penalized in terms of accuracy with respect to a color summed one, we
can also conclude that for single-threaded algorithms, color sampling is much
faster than summing for large multiplicities.

###  5.2 Comparisons of massively parallel execution on GPU

Figure 7 displays the evaluation time per event for our GPU codes as
introduced in Sec. 3.

As opposed to the single thread results in Sec. 5.1, we now observe a similar
performance for the color-summed BlockGen-COΣ and the color-sampled BlockGen-
CD${}_{\text{MC}}$ algorithms. This is caused by the varying number of valid
color configurations for each thread in the color-sampled case: Threads with
few valid configurations need to wait for those with many. However, the
scaling of the color-sampled version still is better, which is mainly due to
the absence of the $\mathcal{O}((n-2)!^{\,2})$ scaling of the color summation
as discussed in Sec. 4.3.

The color-dressed BlockGen-CD${}_{\text{MC}}$ exhibits an excellent scaling,
but is offset by almost two orders of magnitude and only performs better than
the color-ordered codes for ${n_{\mathrm{out}}\gtrsim 6}$. As can be seen in
Eq. (28), the color-dressed approach requires significantly more memory per
current, as discussed in Sec. 4.4, such that the memory handling requires a
significant amount of time during execution.

Figure 7: The timings for various GPU-based algorithms are compared as a
function of gluon multiplicity. All algorithms were run on an NVIDIA V100
($16\text{\,}\mathrm{G}\mathrm{B}$ global memory, 5,120 CUDA cores,
$6144\text{\,}\mathrm{K}\mathrm{B}$ L2 cache).

###  5.3 Comparison of single-threaded and massively parallelized algorithms

To compare the CPU-based with the massively parallelized GPU-based algorithms,
we run the single-threaded Amegic (with CSW rules enabled) and Comix
algorithms parallelized via MPI (Message Passing Interface) over 16 threads.
The resulting performance is very close to dividing the single-thread
evaluation times by 16. By making use of all the threads of our test CPU, we
achieve a realistic chip-to-chip comparison to the performance of our GPU-
based algorithms. However, we note that the CPU is a discontinued
Intel${}^{\text{\textregistered}}$ Xeon${}^{\text{\textregistered}}$ E5-2650
v2 ($2.60\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$,
$20\text{\,}\mathrm{M}\mathrm{B}$ cache), while the GPU is a modern NVIDIA
V100 ($16\text{\,}\mathrm{G}\mathrm{B}$ global memory, 5,120 CUDA cores,
$6144\text{\,}\mathrm{K}\mathrm{B}$ L2 cache)666For additional information
about the NVIDIA V100, see https://images.nvidia.com/content/volta-
architecture/pdf/volta-architecture-whitepaper.pdf. In order to emulate a
realistic chip-to-chip comparison for a similarly capable hardware, we would
need to scale the number of physical CPU cores from 8 to 36, corresponding for
example to an Intel${}^{\text{\textregistered}}$
Xeon${}^{\text{\textregistered}}$ Platinum 8360Y
($2.40\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$,
$54\text{\,}\mathrm{M}\mathrm{B}$ Cache). Moreover, we note again that the
implementation of BlockGen-COΣ and BlockGen-CO${}_{\text{MC}}$ is based on
purely real numbers, cf. Sec. 3. We therefore expect an overhead of about a
factor of two compared to the public codes, which work with complex numbers.

The resulting times are displayed in Fig. 8. We find that for
$n_{\mathrm{out}}\leq 6$, the best-performing GPU algorithm (BlockGen-COΣ) has
evaluation times per event which are at least an order of magnitude smaller
than the best CPU evaluation times. Taking into account hardware differences,
the improvement is reduced by about a factor of four, but remains substantial.
For $n_{\mathrm{out}}=7$, the GPU code performs similarly to the CPU code. For
$n_{\text{out}}\geq 7$, BlockGen-CD${}_{\text{MC}}$ becomes the fastest
algorithm. However, Comix achieves a similar performance on the CPU in that
region, when accounting for hardware differences and the usage of complex
instead of real numbers.

Figure 8: The timings for GPU-based (crosses) and CPU-based (dots) algorithms
are compared against each other as a function of gluon multiplicity. The CPU
numbers are all generated on an Intel${}^{\text{\textregistered}}$
Xeon${}^{\text{\textregistered}}$ E5-2650 v2 8-core CPU,
($2.60\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$,
$20\text{\,}\mathrm{M}\mathrm{B}$ cache), while all the GPU numbers are
generated on a NVIDIA V100 ($16\text{\,}\mathrm{G}\mathrm{B}$ global memory,
5,120 CUDA cores, $6144\text{\,}\mathrm{K}\mathrm{B}$ L2 cache). The MPI
versions are run on 16 threads, and the timing for the color summed algorithm
is divided by a factor of 16 to mimic the improvements that would occur from
MPI. Furthermore, a modified version of Amegic is used in order to perform
helicity sampling.

##  6 Conclusion

As a first step towards GPU-assisted Standard Model Monte-Carlo event
generators, this study explored a variety of methods to calculate leading
order $n$-gluon squared amplitudes with full color dependence. In view of
increased precision requirements for high-multiplicity final-state observables
at the high-luminosity LHC, we focused in particular on the scaling of the
computation time with the number of external gluons. We studied different
algorithms based on the Berends–Giele recursion, which differ in the treatment
of helicity and color sums involved when squaring the scattering amplitude.
The study was performed within a set of constraints relevant for realistic
simulations. In particular, we did not take advantage of special symmetries in
all-gluon scattering amplitudes and only considered algorithms that generate
strictly positive weights. Our results are therefore representative of what
can be expected from a full-fledged matrix-element calculator in the Standard
Model and extensions thereof, which do not involve more than four-particle
interactions in the Lagrangian.

Since the different algorithms vary in their scaling and memory access
patterns, the algorithm with optimal performance depends on the final-state
particle multiplicity. Therefore, there is no one best choice for LHC
phenomenology. At low to medium multiplicity, explicit color summation
combined with color ordered amplitudes, as implemented in BlockGen-COΣ,
provides the best performance. At high multiplicity, color sampling with color
dressed amplitudes, as implemented in BlockGen-CD${}_{\text{MC}}$, is
preferred due to the different scaling pattern. It is interesting to note that
in a direct chip-to-chip comparison the existing parallel CPU-based Comix
generator performs at a similar level as BlockGen-CD${}_{\text{MC}}$ at high
multiplicity. Given that typical event simulations at the LHC will not require
more than seven jets in the final state, we conclude that BlockGen-COΣ is the
optimal choice to implement as a full-fledged low-to-medium multiplicity event
generator on the GPU.

The next steps towards a complete Monte Carlo simulation of Standard Model
events will involve the addition of quark processes and the development of a
GPU-based phase-space generator. Since high-performance computing systems have
both CPUs and GPUs associated with each node, it will be interesting to
explore hybrid computation schemes with load balancing between the CPUs and
GPUs available (see App. A). These hybrid schemes should be able to take
advantage of the computational power available on these systems.

## Acknowledgments

This research was supported by the Fermi National Accelerator Laboratory
(Fermilab), a U.S. Department of Energy, Office of Science, HEP User Facility.
Fermilab is managed by Fermi Research Alliance, LLC (FRA), acting under
Contract No. DE–AC02–07CH11359. MK and EB would like to thank Steffen Schumann
for his support of the project, and the fruitful discussions leading towards
it.

##  Appendix A Proposal of a hybrid approach

In this appendix, we probe the possibility to construct a hybrid event
generator code, where some parts of the computation are performed on the CPU,
and other parts are performed on the GPU. Such a division of tasks can involve
the phase-space generation, the calculation of phase-space weights, and the
calculation of the matrix elements. In Fig. 9, we show timing ratios for these
three tasks in comparison to the color-summed matrix element calculation using
BlockGen-COΣ.

We employ the recursive phase-space generator of Comix [6] to provide the
required time estimates for the phase-space tasks. This corresponds to the
scenario where the CPU provides the phase-space points and weights, while the
GPU only evaluates the matrix elements for those points. In this scenario, a
possible bottleneck might be the communication between the CPU and the GPU. We
therefore also include a time estimate for copying the phase space points to
the GPU. It turns out that the memory copy is not relevant for the overall
timing in this scenario. Here we assume that the current MPI implementation of
Comix’ phase-space generator would be changed to an OpenMP parallelized
algorithm, such that the memory copy to the GPU can take place within a single
process. The results in Fig. 9 show that such a division of tasks would not be
viable for multiplicities $n_{\text{out}}<7$, since the GPU would mostly wait
for the CPU to generate phase space points, weight them, and copy them over.

A more viable option for future development might be a scheme where the CPU
only generates phase space points, while the GPU calculates both the
associated weights and the matrix elements. In this scheme, the CPU would
spend between about $70\text{\,}\mathrm{\char 37\relax}$ at $n_{\text{out}}=3$
and a few percent at $n_{\text{out}}=7$ of the time needed by the GPU and
would thus take over a sizable share of the overall workload without forcing
the GPU to be intermittently idle. Only for $n_{\text{out}}=2$, the CPU time
requirement would exceed the one for the GPU, by about a factor of three. Note
that these estimates assume that the phase-space weight calculation takes the
same time on the GPU as it does for Comix on the CPU. If the GPU
implementation leads to a sizable speed-up, the GPU would need to wait for the
CPU even for $n_{\text{out}}\gtrsim 2$, which might however still be an
acceptable drawback to accelerate the calculation for larger multiplicities.
It is beyond the scope of this work to implement a recursive GPU phase-space
generator needed to further assess this scenario.

Figure 9:  The timing ratio between phase-space point generation and weight
computation using the recursive phase-space generator of Comix on the one
hand, and the matrix element computation based on BlockGen-COΣ on the other
hand. The ratio is also shown for the time required to copy the phase space
points to GPU memory, and for the sum of everything except for BlockGen-COΣ.
Comix timing is measured on an Intel${}^{\text{\textregistered}}$
Xeon${}^{\text{\textregistered}}$ E5-2650 v2 8-core CPU,
($2.60\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$,
$20\text{\,}\mathrm{M}\mathrm{B}$ cache), using MPI over 16 threads, while the
GPU numbers are generated on a NVIDIA V100 ($16\text{\,}\mathrm{G}\mathrm{B}$
global memory, 5,120 CUDA cores, $6144\text{\,}\mathrm{K}\mathrm{B}$ L2
cache).

## References

  * [1] B. Webber,  _Monte Carlo Simulation of Hard Hadronic Processes_ ,  Ann. Rev. Nucl. Part. Sci. 36, 253 (1986).
  * [2] A. Buckley _et al._ ,  _General-purpose event generators for LHC physics_ ,  Phys. Rept. 504, 145 (2011),  http://dx.doi.org/10.1016/j.physrep.2011.03.005,  1101.2599.
  * [3] A. Kanaki and C. G. Papadopoulos,  _HELAC: A package to compute electroweak helicity amplitudes_ ,  Comput. Phys. Commun. 132, 306 (2000),  hep-ph/0002082.
  * [4] C. G. Papadopoulos,  _PHEGAS: A phase-space generator for automatic cross-section computation_ ,  Comput. Phys. Commun. 137, 247 (2001),  hep-ph/0007335.
  * [5] F. Krauss, R. Kuhn and G. Soff,  _AMEGIC++ 1.0: A Matrix Element Generator In C++_ ,  JHEP 02, 044 (2002),  hep-ph/0109036.
  * [6] T. Gleisberg and S. Höche,  _Comix, a new matrix element generator_ ,  JHEP 12, 039 (2008),  10.1088/1126-6708/2008/12/039,  0808.3674.
  * [7] F. Maltoni and T. Stelzer,  _MadEvent: automatic event generation with MadGraph_ ,  JHEP 02, 027 (2003),  hep-ph/0208156.
  * [8] J. Alwall _et al._ ,  _MadGraph/MadEvent v4: The new web generation_ ,  JHEP 09, 028 (2007),  0706.2334.
  * [9] J. Alwall, M. Herquet, F. Maltoni, O. Mattelaer and T. Stelzer,  _MadGraph 5 : Going Beyond_ ,  JHEP 06, 128 (2011),  10.1007/JHEP06(2011)128,  1106.0522.
  * [10] M. L. Mangano, M. Moretti, F. Piccinini, R. Pittau and A. D. Polosa,  _ALPGEN, a generator for hard multiparton processes in hadronic collisions_ ,  JHEP 07, 001 (2003),  hep-ph/0206293.
  * [11] G. Ossola, C. G. Papadopoulos and R. Pittau,  _CutTools: A Program implementing the OPP reduction method to compute one-loop amplitudes_ ,  JHEP 0803, 042 (2008),  10.1088/1126-6708/2008/03/042,  0711.3596.
  * [12] T. Gleisberg and F. Krauss,  _Automating dipole subtraction for QCD NLO calculations_ ,  Eur. Phys. J. C53, 501 (2008),  0709.2881.
  * [13] C. F. Berger, Z. Bern, L. J. Dixon, F. Febres-Cordero, D. Forde, H. Ita, D. A. Kosower and D. Maître,  _Automated implementation of on-shell methods for one-loop amplitudes_ ,  Phys. Rev. D78, 036003 (2008),  10.1103/PhysRevD.78.036003,  0803.4180.
  * [14] G. Bevilacqua, M. Czakon, M. Garzelli, A. van Hameren, A. Kardos _et al._ ,  _HELAC-NLO_ ,  Comput.Phys.Commun. 184, 986 (2013),  10.1016/j.cpc.2012.10.033,  1110.1499.
  * [15] F. Cascioli, P. Maierhöfer and S. Pozzorini,  _Scattering Amplitudes with Open Loops_ ,  Phys.Rev.Lett. 108, 111601 (2012),  10.1103/PhysRevLett.108.111601,  1111.5206.
  * [16] G. Cullen, H. van Deurzen, N. Greiner, G. Heinrich, G. Luisoni _et al._ ,  _GoSam -2.0: a tool for automated one-loop calculations within the Standard Model and beyond_,  Eur.Phys.J. C74(8), 3001 (2014),  10.1140/epjc/s10052-014-3001-5,  1404.7096.
  * [17] J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H.-S. Shao, T. Stelzer, P. Torrielli and M. Zaro,  _The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations_ ,  JHEP 07, 079 (2014),  10.1007/JHEP07(2014)079,  1405.0301.
  * [18] S. Actis, A. Denner, L. Hofer, J.-N. Lang, A. Scharf and S. Uccirati,  _RECOLA: REcursive Computation of One-Loop Amplitudes_ ,  Comput. Phys. Commun. 214, 140 (2017),  10.1016/j.cpc.2017.01.004,  1605.01090.
  * [19] N. D. Christensen, P. de Aquino, C. Degrande, C. Duhr, B. Fuks, M. Herquet, F. Maltoni and S. Schumann,  _A comprehensive approach to new physics simulations_ ,  Eur. Phys. J. C71, 1541 (2011),  10.1140/epjc/s10052-011-1541-5,  0906.2474.
  * [20] C. Degrande, C. Duhr, B. Fuks, D. Grellscheid, O. Mattelaer and T. Reiter,  _UFO - The Universal FeynRules Output_ ,  Comput.Phys.Commun. 183, 1201 (2012),  10.1016/j.cpc.2012.01.022,  1108.2040.
  * [21] F. Staub,  _SARAH 4: A tool for (not only SUSY) model builders_ ,  Comput.Phys.Commun. 185, 1773 (2014),  10.1016/j.cpc.2014.02.018,  1309.7223.
  * [22] P. Calafiura, J. Catmore, D. Costanzo and A. Di Girolamo,  _ATLAS HL-LHC Computing Conceptual Design Report_ ,  Tech. Rep. CERN-LHCC-2020-015. LHCC-G-178, CERN, Geneva (2020).
  * [23] S. Amoroso _et al._ ,  _Challenges in Monte Carlo Event Generator Software for High-Luminosity LHC_ ,  Comput. Softw. Big Sci. 5(1), 12 (2021),  10.1007/s41781-021-00055-1,  2004.13687.
  * [24] S. Höche _et al._ ,  _Working Group Report: Computing for Perturbative QCD_ ,  In _Proceedings, 2013 Community Summer Study on the Future of U.S. Particle Physics: Snowmass on the Mississippi (CSS2013): Minneapolis, MN, USA, July 29-August 6, 2013_ (2013), 1309.3598.
  * [25] J. Albrecht _et al._ ,  _A Roadmap for HEP Software and Computing R &D for the 2020s_,  Comput. Softw. Big Sci. 3(1), 7 (2019),  10.1007/s41781-018-0018-8,  1712.06982.
  * [26] A. Dainese, M. Mangano, A. B. Meyer, A. Nisati, G. Salam and M. A. Vesterinen, eds.,  _Report on the Physics at the HL-LHC,and Perspectives for the HE-LHC_ , vol. 7. 7/2019 of _CERN Yellow Reports: Monographs_ ,  CERN, Geneva, Switzerland,  ISBN 9789290835493,  10.23731/CYRM-2019-007 (2019).
  * [27] K. Hagiwara, J. Kanzaki, N. Okamura, D. Rainwater and T. Stelzer,  _Fast calculation of HELAS amplitudes using graphics processing unit (GPU)_ ,  The European Physical Journal C 66(3-4), 477 (2010),  10.1140/epjc/s10052-010-1276-8.
  * [28] K. Hagiwara, J. Kanzaki, N. Okamura, D. Rainwater and T. Stelzer,  _Calculation of HELAS amplitudes for QCD processes using graphics processing unit (GPU)_ ,  The European Physical Journal C 70(1), 513 (2010),  10.1140/epjc/s10052-010-1465-5.
  * [29] J. Kanzaki,  _Application of graphics processing unit (GPU) to software in elementary particle/high energy physics field_ ,  Procedia Computer Science 4, 869 (2011),  https://doi.org/10.1016/j.procs.2011.04.092,  Proceedings of the International Conference on Computational Science, ICCS 2011.
  * [30] K. Hagiwara, J. Kanzaki, Q. Li, N. Okamura and T. Stelzer,  _Fast computation of MadGraph amplitudes on graphics processing unit (GPU)_ ,  The European Physical Journal C 73(11), 2608 (2013),  10.1140/epjc/s10052-013-2608-2.
  * [31] J. Kanzaki,  _Monte Carlo integration on GPU_ ,  Eur. Phys. J. C71, 1559 (2011),  10.1140/epjc/s10052-011-1559-8,  1010.2107.
  * [32] W. T. Giele, G. C. Stavenga and J. Winter,  _Thread-scalable evaluation of multi-jet observables_ ,  The European Physical Journal C 71(7) (2011),  10.1140/epjc/s10052-011-1703-5.
  * [33] H.-Z. Wu, J.-J. Zhang, L.-G. Pang and Q. Wang,  _ZMCintegral: a Package for Multi-Dimensional Monte Carlo Integration on Multi-GPUs_ ,  Comput. Phys. Commun. 248, 106962 (2020),  10.1016/j.cpc.2019.106962,  1902.07916.
  * [34] G. Grasseau, S. Lisniak and D. Chamont,  _Hybrid implementation of the VEGAS Monte-Carlo algorithm_ ,  In _Proceedings, GPU Computing in High-Energy Physics (GPUHEP2014): Pisa, Italy, September 10-12, 2014_ , pp. 103–108,  10.3204/DESY-PROC-2014-05/19 (2015).
  * [35] S. Carrazza and J. M. Cruz-Martinez,  _VegasFlow: accelerating Monte Carlo simulation across platforms_ ,  arXiv e-prints arXiv:2010.09341 (2020),  2010.09341.
  * [36] S. Carrazza, J. Cruz-Martinez, M. Rossi and M. Zaro,  _Towards the automation of Monte Carlo simulation on GPU for particle physics processes_ ,  In _25th International Conference on Computing in High-Energy and Nuclear Physics_ (2021), 2105.10529.
  * [37] F. Yuasa, T. Ishikawa, N. Hamaguchi, T. Koike and N. Nakasato,  _Acceleration of Feynman loop integrals in high-energy physics on many core GPUs_ ,  J. Phys. Conf. Ser. 454, 012081 (2013),  10.1088/1742-6596/454/1/012081.
  * [38] S. Borowka, G. Heinrich, S. Jahn, S. P. Jones, M. Kerner and J. Schlenk,  _A GPU compatible quasi-Monte Carlo integrator interfaced to pySecDec_ ,  Comput. Phys. Commun. 240, 120 (2019),  10.1016/j.cpc.2019.02.015,  1811.11720.
  * [39] S. Carrazza, J. M. Cruz-Martinez and M. Rossi,  _PDFFlow: Parton distribution functions on GPU_ ,  Comput. Phys. Commun. 264, 107995 (2021),  10.1016/j.cpc.2021.107995,  2009.06635.
  * [40] D. Schouten, A. DeAbreu and B. Stelzer,  _Accelerated matrix element method with parallel computing_ ,  Computer Physics Communications 192, 54 (2015),  https://doi.org/10.1016/j.cpc.2015.02.020.
  * [41] G. Grasseau, D. Chamont, F. Beaudette, L. Bianchini, O. Davignon, L. Mastrolorenzo, C. Ochando, P. Paganini and T. a. Strebler,  _Matrix element method for high performance computing platforms_ ,  Journal of Physics Conference Series 664(9), 092009 (2015),  10.1088/1742-6596/664/9/092009.
  * [42] G. Grasseau, F. Beaudette, C. Martin Perez, A. Zabi, A. Chiron, T. Strebler and G. Hautreux,  _Deployment of a Matrix Element Method code for the ttH channel analysis on GPU’s platform_ ,  EPJ Web Conf. 214, 06028 (2019),  10.1051/epjconf/201921406028.
  * [43] F. A. Berends and W. T. Giele,  _Recursive calculations for processes with $n$ gluons_,  Nucl. Phys. B306, 759 (1988).
  * [44] S. Badger, B. Biedermann, L. Hackl, J. Plefka, T. Schuster and P. Uwer,  _Comparing efficient computation methods for massless QCD tree amplitudes: Closed analytic formulas versus Berends-Giele recursion_ ,  Phys. Rev. D87(3), 034011 (2013),  10.1103/PhysRevD.87.034011,  1206.2381.
  * [45] C. Duhr, S. Höche and F. Maltoni,  _Color-dressed recursive relations for multi-parton amplitudes_ ,  JHEP 08, 062 (2006),  hep-ph/0607057.
  * [46] W. T. Giele, Z. Kunszt and J. Winter,  _Efficient color-dressed calculation of virtual corrections_ ,  Nuclear Physics B 840(1-2), 214 (2010),  10.1016/j.nuclphysb.2010.07.007.
  * [47] F. Maltoni, K. Paul, T. Stelzer and S. Willenbrock,  _Color-flow decomposition of QCD amplitudes_ ,  Phys. Rev. D67, 014026 (2003),  hep-ph/0209271.
  * [48] P. Draggiotis, R. H. Kleiss and C. G. Papadopoulos,  _On the computation of multigluon amplitudes_ ,  Phys. Lett. B 439, 157 (1998),  10.1016/S0370-2693(98)01015-6,  hep-ph/9807207.
  * [49] R. Brandelik _et al._ ,  _Evidence for Planar Events in e $+$ e$-$ Annihilation at High-Energies_,  Phys.Lett. B86, 243 (1979).
  * [50] C. Berger _et al._ ,  _Evidence for Gluon Bremsstrahlung in e $+$ e$-$ Annihilations at High-Energies_,  Phys.Lett. B86, 418 (1979),  10.1016/0370-2693(79)90869-4.
  * [51] D. Barber, U. Becker, H. Benda, A. Boehm, J. Branson _et al._ ,  _Discovery of Three Jet Events and a Test of Quantum Chromodynamics at PETRA Energies_ ,  Phys.Rev.Lett. 43, 830 (1979),  10.1103/PhysRevLett.43.830.
  * [52] B. L. Combridge, J. Kripfganz and J. Ranft,  _Hadron Production at Large Transverse Momentum and QCD_ ,  Phys. Lett. B 70, 234 (1977),  10.1016/0370-2693(77)90528-7.
  * [53] R. Cutler and D. W. Sivers,  _Quantum Chromodynamic Gluon Contributions to Large p(T) Reactions_ ,  Phys. Rev. D 17, 196 (1978),  10.1103/PhysRevD.17.196.
  * [54] T. Gottschalk and D. W. Sivers,  _Basic Processes and Formalism for the Hadronic Production of Three Large p(T) Jets_ ,  Phys. Rev. D 21, 102 (1980),  10.1103/PhysRevD.21.102.
  * [55] Z. Kunszt and E. Pietarinen,  _Production of Three Large p(T) Jets in Hadron-Hadron Collisions_ ,  Nucl. Phys. B 164, 45 (1980),  10.1016/0550-3213(80)90501-5.
  * [56] F. A. Berends, R. Kleiss, P. De Causmaecker, R. Gastmans and T. T. Wu,  _Single Bremsstrahlung Processes in Gauge Theories_ ,  Phys. Lett. B 103, 124 (1981),  10.1016/0370-2693(81)90685-7.
  * [57] Z. Kunszt _et al._ ,  _LA THUILE WORKSHOP ON PHYSICS OF FUTURE ACCELERATORS: LARGE CROSS-SECTION PROCESSES_ ,  In _Workshop on Physics at Future Accelerators_ (1987).
  * [58] J. Alitti _et al._ ,  _A Study of multi - jet events at the CERN anti-p p collider and a search for double parton scattering_ ,  Phys. Lett. B 268, 145 (1991),  10.1016/0370-2693(91)90937-L.
  * [59] S. J. Parke and T. R. Taylor,  _An Amplitude for $n$ Gluon Scattering_,  Phys. Rev. Lett. 56, 2459 (1986),  10.1103/PhysRevLett.56.2459.
  * [60] J. F. Gunion and Z. Kunszt,  _Six Quark Subprocesses in QCD_ ,  Phys. Lett. B 176, 163 (1986),  10.1016/0370-2693(86)90944-5.
  * [61] J. F. Gunion and J. Kalinowski,  _A Direct Calculation of the Six Gluon Subprocess_ ,  Phys. Rev. D 34, 2119 (1986),  10.1103/PhysRevD.34.2119.
  * [62] F. A. Berends and W. Giele,  _The six-gluon process as an example of Weyl-van der Waerden spinor calculus_ ,  Nucl. Phys. B294, 700 (1987).
  * [63] M. L. Mangano, S. J. Parke and Z. Xu,  _Duality and multi-gluon scattering_ ,  Nucl. Phys. B298, 653 (1988).
  * [64] L. J. Dixon,  _Calculating scattering amplitudes efficiently_ ,  In _Theoretical Advanced Study Institute in Elementary Particle Physics (TASI 95): QCD and Beyond_ , pp. 539–584 (1996), hep-ph/9601359.
  * [65] L. J. Dixon,  _A brief introduction to modern amplitude methods_ ,  In _Theoretical Advanced Study Institute in Elementary Particle Physics: Particle Physics: The Higgs Boson and Beyond_ ,  10.5170/CERN-2014-008.31 (2013), 1310.5353.
  * [66] P. D. Draggiotis, R. H. P. Kleiss and C. G. Papadopoulos,  _Multi-jet production in hadron collisions_ ,  Eur. Phys. J. C24, 447 (2002),  hep-ph/0202201.
  * [67] V. Del Duca, L. J. Dixon and F. Maltoni,  _New color decompositions for gauge amplitudes at tree and loop level_ ,  Nucl. Phys. B571, 51 (2000),  hep-ph/9910563.
  * [68] V. del Duca, A. Frizzo and F. Maltoni,  _Factorization of tree QCD amplitudes in the high-energy limit and in the collinear limit_ ,  Nucl. Phys. B568, 211 (2000),  hep-ph/9909464.
  * [69] R. Kleiss and H. Kuijf,  _Multi-gluon cross-sections and five jet production at hadron colliders_ ,  Nucl. Phys. B312, 616 (1989).
  * [70] R. Kleiss, W. J. Stirling and S. D. Ellis,  _A New Monte Carlo Treatment of Multiparticle Phase Space at High-energies_ ,  Comput. Phys. Commun. 40, 359 (1986),  10.1016/0010-4655(86)90119-0.
  * [71] J. Kuipers, T. Ueda, J. Vermaseren and J. Vollinga,  _Form version 4.0_ ,  Computer Physics Communications 184(5), 1453 (2013),  https://doi.org/10.1016/j.cpc.2012.12.028.
  * [72] E. Bothmann _et al._ ,  _Event Generation with Sherpa 2.2_ ,  SciPost Phys. 7(3), 034 (2019),  10.21468/SciPostPhys.7.3.034,  1905.09127.
  * [73] T. Gleisberg, S. Höche, F. Krauss and R. Matyszkiewicz,  _How to calculate colourful cross sections efficiently_ ,  Preprint (2008), 0808.3672.
  * [74] F. Cachazo, P. Svrček and E. Witten,  _MHV vertices and tree amplitudes in gauge theory_ ,  JHEP 09, 006 (2004),  hep-th/0403047.