Caravel: A C++ Framework for the Computation of Multi-Loop Amplitudes with Numerical Unitarity

We present the first public version of Caravel, a C++17 framework for the computation of multi-loop scattering amplitudes in quantum field theory, based on the numerical unitarity method. Caravel is composed of modules for the $D$-dimensional decomposition of integrands of scattering amplitudes into master and surface terms, the computation of tree-level amplitudes in floating point or finite-field arithmetic, the numerical computation of one- and two-loop amplitudes in QCD and Einstein gravity, and functional reconstruction tools. We provide programs that showcase Caravel's main functionalities and allow to compute selected one- and two-loop amplitudes.


Introduction
The computation of scattering amplitudes in quantum field theory is crucial in our quest to describe high-energy particle interactions. Indeed, these objects allow one to make theoretical predictions which can then be compared with experimental measurements, be it at particle colliders such as the Large Hadron Collider (LHC), which are testing the Standard Model of particle physics, or at experiments such as LIGO or VIRGO, which are testing our understanding of gravity. Computing these amplitudes remains a challenge, in particular when they contain many external particles and when higher-order corrections in their perturbative expansion are necessary. The former implies a dependence on a large number of physical scales, and the latter a number of unconstrained loop momenta which must be integrated over. Combined, these two aspects lead to a considerable level of complexity in the computation of these amplitudes. This is nevertheless a very timely and important problem to tackle. For example, two-loop five-particle amplitudes as well as three-loop four-particle amplitudes are already relevant for phenomenological studies at the LHC, and will become even more so over the next years.
A major obstacle to overcome when evaluating loop amplitudes is the complexity of intermediate computational steps. In order to bypass this is-sue, it has proven fruitful to consider numerical approaches. An outstanding example of this is the current possibility to compute general one-loop amplitudes, which has been powered by the introduction of robust numerical techniques. Among many other developments, these techniques are based on integrand reduction approaches [1,2], on the one-loop numerical unitarity method [3][4][5], and on recursive approaches [6,7]. More recently, there has been great progress in the numerical computation of two-loop amplitudes at high multiplicities. By now, all five-parton amplitudes [8][9][10][11] and the amplitudes for four partons and a W boson [12] have been computed numerically at leading color. Through the use of finite fields and multivariate functional reconstruction techniques [13,14], the frameworks powering these numerical computations have furthermore allowed the calculation of the analytic form of all five-parton leading-color amplitudes [15][16][17], the five-point allplus amplitude at full color [18] and the four-graviton amplitude in Einstein gravity [19]. In related work, the amplitudes for three-photon production at the LHC have also been computed [20].
In this article we present the Caravel C ++ framework. It provides an implementation of many algorithms necessary to perform computations of multi-loop scattering amplitudes within the multi-loop numerical unitarity method. This is the first publicly available code of its kind. It is based on the (generalized) unitarity approach, which was first developed for the analytic computation of one-loop amplitudes [21][22][23][24] and later adapted for numeric calculations [3][4][5]. An extension of the method beyond one loop has been developed recently [25][26][27]. In a nutshell, in this framework the amplitude is computed starting from a parametrization of its integrand. The corresponding free parameters are numerically computed at each phase-space point by constructing systems of linear equations in which the parameters are the unknowns and the numerical entries are associated to products of tree-level amplitudes. With a suitable choice of integrand parametrization [25], this directly gives a decomposition of the amplitude in terms of master integrals. Finally, after inserting the value of the integrals at the required phase-space point we obtain the value of the amplitude.
The current release of Caravel includes a module for computing products of tree-level amplitudes in several theories through off-shell recursion relations [28], and tools that allow the efficient construction and solution of the systems of linear equations that determine the integrand. Whilst these components work for generic multi-loop amplitudes, other components such 5 as the construction of the parametrization are required as input. In this release we showcase the different available tools by providing a series of example programs. Discussions regarding full automation are left to future work. The example programs give a first-hand account of the procedures employed for the calculations presented in refs. [8,10,16,17,19,27].
The rest of this article is organized as follows. Section 2 provides a brief description of our computational methodology, section 3 gives a description of the organization of the internal modules in the Caravel framework, section 4 describes the procedure of installation and setup of the libraries. In section 5 we give details on how to run the example programs that we provide and we conclude in section 6. Appendix A and Appendix B contain technical details about our conventions for color-ordered helicity amplitudes and phase-space parametrizations.

Computational methodology
In this section we briefly review the main features of the multi-loop numerical unitarity method. This framework allows for the numerical evaluation as well as the analytical reconstruction of multi-loop scattering amplitudes. Our approach is generic, in that it facilitates the computation of amplitudes in different quantum field theories. When dealing with QCD amplitudes, we consider the amplitudes to be vectors M in a generic color space. In this section we will mostly be concerned with the components of these vectors, that is with the A(σ) defined as where the C(σ) span the relevant color space. For QCD amplitudes, this is the usual SU(N c ) color space, which can be specialized to N c = 3. For pure gravity processes, no color is present and in eq. (1) there is a single σ with C(σ) = 1. In the remaining of this section, we will always discuss the A(σ) defined in eq. (1), which for brevity we will call amplitudes, and for simplicity we will drop the σ dependence. In later sections, when discussing the calculation of amplitudes in specific theories, we will be more explicit about the C(σ). As is standard in perturbative quantum field theory, we will consider the expansion of A around a small coupling constant, which is associated to an expansion in the number of loops of the contributing Feynman integrals. We refer the reader to Appendix A for more details on the definition of the objects we compute. For more details on the techniques outlined in this section we refer to previous publications [8,10,16,17,19,27].

Integrand parametrization
In full generality, an L-loop amplitude A (L) can be decomposed as a linear combination of master integrals according to Here, Γ defines a propagator structure associated with the amplitude, and we will often refer to it as a diagram (indeed, they are in one-to-one correspondence with scalar Feynman integrals). The set ∆ contains all propagator structures Γ, and can be organized hierarchically according to whether a propagator structure Γ 1 ∈ ∆ can be obtained from another Γ 2 ∈ ∆ by removing some propagators in the latter. In this case we write Γ 1 < Γ 2 . The set M Γ denotes the full set of master integrals associated to Γ. Each master integral I Γ,i is usually expressed as a Laurent series in the dimensional-regularization parameter ǫ = (4 − D)/2. The coefficients in the Laurent expansion involve multi-valued functions with non-trivial branch-cut structures, such as multiple polylogarithms. The coefficients c Γ,i are algebraic functions of the kinematic invariants and rational functions of ǫ.
At its core, the unitarity method is a way to compute the integrand of an amplitude. We therefore start with a parametrization of the integrand A (L) (ℓ l ), where ℓ l represents the set of L loop momenta, of the form The multiset P Γ labels all inverse propagators ρ j in the diagram Γ, and the basis of numerators Q Γ = {m Γ,k (ℓ l )|k ∈ Q Γ } parametrizes every possible integrand up to the maximum allowed power of loop momenta, which is theory specific.
The numerator basis Q Γ is not unique. Let us highlight two natural classes of bases. First, the simplest choice is what we call a tensor basis, denoted by T Γ . It can be built out of independent monomials in a set of variables α j , which parametrize the loop momenta ℓ l ( α) [29][30][31][32]. This type of basis is straightforward to build for generic Γ ∈ ∆. However, with this choice, the relation between eqs. (2) and (3) is not explicit, as it would require solving large systems of integration-by-part (IBP) relations (see e.g. [13,[33][34][35][36][37][38][39][40][41]). A second class of bases is what we call a master-surface basis M Γ , and it is crucial to our multi-loop numerical unitarity method. It was observed in ref. [25] that one can parametrize the integrand of a multi-loop amplitude by functions M Γ = {m Γ,i (ℓ l )|i ∈ M Γ ∪ S Γ } such that the associated integrands either integrate to zero or correspond to the integrand of a master integral: The numerators m Γ,i (ℓ l ) with i ∈ M Γ are called master integrands and the ones with i ∈ S Γ surface terms. A master-surface basis of functions thus makes the relation between eqs. (2) and (3) explicit. As discussed in [25][26][27] the construction of master-surface bases of integrands can be efficiently performed by using unitarity-compatible IBP relations [42,43], employing computational algebraic geometry techniques.

Integrand factorization and cut equations
In order to compute the coefficients c Γ,k in eq. (3) we exploit the factorization properties of multi-loop integrands for on-shell configurations ℓ Γ l of the loop momenta. For a given propagator structure P Γ , these are defined by In most cases, this does not fix the ℓ Γ l completely, and there is some residual degree of freedom. In this limit, the leading pole of eq. (3) is given by where T Γ represents the set of all tree-level amplitudes corresponding to the vertices of the diagram Γ, and the state sum runs over all D s -dimensional particle states that can appear in the loop. On the right-hand side the sum runs over all propagator structures Γ ′ with equal or more propagators than Γ (hence P Γ ⊆ P Γ ′ ). Eq. (6) is a so-called cut equation.
The cut equations allow one to numerically compute the coefficients c Γ,k in eq. (3) by sampling a sufficient set of on-shell momenta ℓ Γ l and solving the resulting system of linear equations. Importantly, some of these coefficients may be identically zero for all phase-space points. To account for this, we identify zero coefficients during a "warm-up run" on a single phase-space point, and remove the corresponding terms from the ansatz for all subsequent evaluations. In order to construct this system of equations we must have an efficient way to evaluate the tree amplitudes on the left-hand side of the cut equations. This is achieved with an implementation of the Berends-Giele off-shell recursion relations [28], which allows one to recursively compute tree amplitudes with an arbitrary number of legs, and where the particles have D s -dimensional states. Being a very general approach, it provides a straightforward way to add new types of fields. Beyond one-loop we also need to consider subleading singular contributions for certain propagator structures for which no generic integrand factorization is known (at two loops, this happens for propagator structures where the same propagator appears twice). To address these cases, we employ cut equations with fewer on-shell constraints and then solve for the corresponding coefficients [26].
Another important aspect in sampling the cut equations is the construction of the on-shell momenta ℓ Γ . These configurations of loop momenta are constructed by solving the quadratic equations in eq. (5), and depending on which number field we use they might not have solutions. To be more precise, we will often do calculations in a field F that is not algebraically closed, such as the field of rational numbers or a finite field (see section 2.5 below), which makes the construction of on-shell momenta with components in F a non-trivial problem. However, it turns out that the square roots originating from the solutions of eq. (5) are only present at intermediate steps of the calculation, and any D-dimensional Lorentz-invariant scalar product of the loop momenta lives in F. In particular, the product of trees in eq. (6) is free of square roots and representable in F. Therefore, we (temporarily) employ an algebraic extension of F for evaluating off-shell currents contributing to the left-hand side of eq. (6). We refer to refs. [10,17,44] for details. In some cases, e.g. for Yang-Mills theory, it is possible to avoid the appearance of square roots altogether [16]. In this case we compute cuts directly in F.
Solving for all coefficients in an amplitude can be efficiently organized in a block-triangular way, using the hierarchical structure of the set ∆. In two-loop five-particle QCD amplitudes each block of equations can have up to a few hundreds of unknowns [8,10], while in Einstein gravity this number is typically an order of magnitude larger. We solve these equations by employing standard linear algebra techniques, such as PLU or QR factorization. Through this procedure, we can compute the coefficients c Γ,i in eq. (2) at a numerical phase-space point, and for numerical values of D and D s .

Special functions
As stated below eq. (2), the master integrals have a Laurent expansion around ǫ = 0, whose coefficients can be written as linear combinations of multivalued special functions. For all the cases currently implemented in Caravel, the special functions are (linear combinations of) multiple polylogarithms. Using modern mathematical techniques [45,46], we can find a basis B for this space of functions, and find an alternative decomposition of the amplitude in terms of the elements h i ∈ B. That is, up to a given order k in the ǫ expansion we can write where the functions r i,j do not depend on ǫ. This decomposition presents a major difference compared to the one of eq. (2): it allows one to write oneand two-loop amplitudes as a linear combination of the same basis of functions. In turn, this then makes it possible to write quantities derived from amplitudes, such as finite remainders, in terms of this basis of functions. This observation was fundamental in reconstructing two-loop five-parton QCD amplitudes [16,17], using the basis of ref. [47], as the coefficients in a decomposition of the form of eq. (7) are much simpler for a two-loop finite remainder than for a two-loop amplitude.

Analytic structure in the dimensional regulators
The cut equation in eq. (6) depends on dimensional regulators, namely on D, the dimension of the loop momenta, and D s , the dimension of the states of loop particles. For convenience, we keep these quantities separate until the final stages of the computation. The coefficient functions c Γ,i in eq. (2) are rational functions in D and polynomials in D s (in some very special cases, this dependence can also be rational, see e.g. [19]).
We evaluate the products of tree-level amplitudes in eq. (6) in integer D s dimensions. To reconstruct the analytic D s dependence we can employ two different approaches. In the first approach, known as dimensional reconstruction [4,8,10,27,48,49], we extract the polynomial dependence by evaluating the tree-level amplitudes for various integer dimensions D s and fit the resulting coefficients of the D s polynomial. This procedure is conceptually straightforward. However, its numerical inefficiency can become a bottleneck for amplitudes with external fermions due to the exponential scaling of the dimension of spinor representations with D s [10,50]. To address this issue, we employ a second approach, which bypasses the dimensional reconstruction method and provides a diagrammatic representation of the coefficients of the D s polynomial [17,44,50]. As an example, this strategy reduces the evaluation time of two-loop five-parton amplitudes with two external fermion lines by about two orders of magnitude.
Regarding the dependence on D, we sample at a sufficient number of values of D in order to reconstruct the rational dependence of each masterintegral coefficient using Thiele's formula [14,51]. This procedure can be computationally intensive. We note nevertheless that the D-dependence in the denominator of the coefficients is rather simple and independent of the phase-space point. When evaluating the same amplitude over several phasespace points, we thus usually perform a warm-up evaluation, dedicated to determining this dependence for each c Γ,k in eq. (3).

Finite fields and functional reconstruction
In ref. [13] it was shown that finite fields can be applied to the Laporta algorithm [34] for the IBP reduction of multi-loop Feynman integrals. The authors were able to not only efficiently perform numerical IBP reductions of integrals, but also to reconstruct the analytic rational dependence of the master integral coefficients in the dimensional regularization parameter ǫ from those numerical reductions. It was later shown in ref. [14] that through a recursive approach one could more generically reconstruct multivariate rational functions. In the same paper, this idea was applied to the computation of scattering amplitudes in generalized unitarity methods. In Caravel, we apply the reconstruction approach to rational coefficient functions in the numerical unitarity method. To do this, we evaluate the amplitudes at rational phase-space points and perform all calculations in a finite field, obtaining exact numerical values for the coefficients. These evaluations can then be used to reconstruct the analytic dependence on the kinematic variables which parametrize the appropriate Lorentz-invariant phase space associated to the amplitude. Caravel contains all the functionalities for the numerical evaluation of scattering amplitudes in a finite field, as well as for the reconstruction of generic rational functions [16,17]. 1 These tools have been fundamental for the computation of the planar two-loop five-parton amplitudes [8,10,16,17], as well as the two-loop four-graviton amplitudes [19].

Internal modules
The Caravel framework is organized in a modular fashion. The structure of the multi-loop numerical unitarity method outlined in the previous section naturally lends itself to this modular implementation. In fig. 1 we show how the different modules of Caravel relate to the main equations of the numerical unitarity method as well as to each other. The red arrows show which module constructs each of the different components of these equations. A black arrow from A to B represents a dependence of B on A. The blue arrows highlight an input that a module receives. The modules highlighted with red boxes receive external input to operate, either in the form of data files or as machine-generated source code.
In the following we list the modules of Caravel, such that the reader gets a general view of the source code of the library. We do not give details on application programming interfaces, as in this release we only include specific example programs for the user, as presented in section 5. The modules are: • Core: This module implements general tools for debugging, arithmetic, kinematics, as well as utilities for linear algebra, rational reconstruction, type traits, and more general operations such as Laurent expansions. Among the arithmetic tools, we include interfaces to the QD [55] library for double-double and quad-double precision floating-point, and the GMP [56] library for arbitrary precision floating-point as well as arbitrary size rational numbers. Furthermore, we use an in-house implementation of 32-bit finite fields based on Barrett reduction [57] which There is also a parser to process headed lists like those employed in Mathematica [60].
• GraphLibrary: A module for the classification and canonicalization of multi-loop graphs. In Caravel, many objects can naturally be associated to graphs, such as the propagator structures Γ in eq. (2). Graph isomorphisms are identified by building a partial order in the representation of the graph (which is ultimately based on the standard C ++ function std::lexicographical compare).
• FunctionalReconstruction: Implementations of algorithms for analytic reconstruction of univariate and multivariate rational functions from exact numerical evaluations (see section 2.5). The reconstruction algorithms are parallelized using either native C ++ threads or using MPI. The latter is useful for use on computer clusters.
• OnShellStrategies: Tools to generate on-shell loop momenta for oneand two-loop diagrams, as required for the construction of the linear system of equations in eq. (6).
• Forest: Implementation of the Berends-Giele off-shell recursion [28] for the computation of general tree-level amplitudes and the products of trees on the left-hand side of eq. (6) in arbitrary D s dimensions. The recursions can be constructed from any given set of Feynman rules 2 and can be evaluated over an arbitrary numerical type (e.g. floating point of different sizes, finite fields, etc.). This release of Caravel includes the Feynman rules for massless QCD and Einstein gravity.
• FunctionSpace: Module for the construction of the integrand ansätze 2 In particular, there is no restriction on the number of particles in vertices.
of eq. (3), both for tensor bases T Γ and master-surface bases M Γ . The former can be constructed for an arbitrary two-loop diagram Γ, while the latter are provided for arbitrary one-loop diagrams and only for the two-loop diagrams required for the calculations in [8,10,16,17,19,27]. Master-surface bases have been produced with private computeralgebra programs, and transformed into C ++ code to be handled by this module.
• Integral Providers: Two separate modules are included to handle analytic expressions of master integrals. The main module is the Inte-gralLibrary, which provides access to the master integrals associated to four-and five-point one-and two-loop planar massless master integrals. In the five-point case, integrals are written in terms of pentagon functions [47]. Internally, all master integrals are normalized as where γ E is the Euler-Mascheroni constant. In the four-point case, master integrals are evaluated with GiNaC [61][62][63] and CLN [64]. In the five-point case, we currently employ a modified version of the library provided in [47]. All master integrals are stored in a format that allows on-the-fly expression of an amplitude in terms of a set of basis functions as described in section 2.3, see eq. (7). In the current version all master integrals are implemented in the Euclidean region. Additionally, Caravel contains the IntegralsBH module, which provides a large collection of one-loop integrals up to O(ǫ 0 ) (including massive and massless propagators) which can be evaluated in up-to quad-double precision using the QD [55] package. This library has been adapted from the BHlibMassive library employed in [65].
• Coefficient Providers: Two different modules perform the hierarchical extraction of master-integrand coefficients via the cut equations (see section 2.2). The AmpEng provider computes general one-loop integral coefficients, building up the hierarchy of diagrams ∆ in an automated fashion. For two-loop calculations, the AmplitudeCoefficients module is employed. For a given amplitude, it requires an input data file, which we call the process library. These process libraries contain all hierarchical relations between the propagator structures in the amplitude, as well as information about the color decomposition [66,67] and the D s -dependence based on the particle content [17,44,50]. In this release, we provide the process libraries employed for the calculations in [8,10,16,17,19,27], which were produced with private computer-algebra code.
• Other modules: further minor modules include miscellaneous functionalities, for example some phase-space generators (including momentum twistor parametrizations [68]) and information on the pole structure of relevant amplitudes.
All modules in Caravel are implemented according to the concept of generic programming, in which algorithms are designed to operate on any data type satisfying certain (minimal) requirements. In particular, our algorithms are well equipped to work with any numerical type, such as floating point number (of fixed or variable size) or numbers in an algebraic field (the rational numbers or numbers in a finite field). This allows us to perform evaluations of amplitudes with different fixed-size floating point numbers and the reconstruction of their analytic form from exact evaluations over finite fields with essentially a single implementation.

Installation and setup
The source code of Caravel can be obtained from a git repository at https://gitlab.com/caravel-public/caravel.git.
Caravel employs the meson [69] build system, which relies on pkg-config for resolving dependencies. For more details on dependence resolution, configuration options and the installation of optional libraries see the README.md and INSTALL.md files in the repository. To build Caravel in the default configuration and install it to the directory <install dir> one can proceed as follows: 1 > git clone https://gitlab.com/caravel-public/caravel.git 2 > cd caravel 3 > mkdir build 4 > cd build All available test suites can be run with the command > ninja test executed in the build directory. Note that, depending on the hardware and build configuration, running all test suites can take a considerable amount of time. This is particularly noticeable for the first time tests are run which produce and store warm-up information.
The default configuration of Caravel provides very limited functionalities. This allows users to customize the installation to suit their particular needs, and additional configuration options and their current values can be queried by running > meson configure in the build directory. An option <option> can be set to a particular value <value> with the command > meson configure -D <option>=<value> These options can be set either at the configuration stage (step 5 of listing 1), or any time before the building is initiated (step 6 of listing 1). It is possible to specify several options at the same time.
Some of the options, listed first when running meson configure, are related to generic C ++ compiler and linker options, which are automatically provided by the meson build system. These options are intended mostly for developers.
Options specific to Caravel can be found in the section Project options, which enable additional features. We first describe the set of options most relevant for a user of the example programs.
• double-inverter: Switch between Eigen and Lapack for solving linear systems in double precision.
• finite-fields: Enable computation using finite fields. Requires the external library GMP. This option can be set to false (default) or true.
• field-ext-fermions: Enable the exact computation of master integral coefficients for amplitudes involving fermions. Off-shell currents are then evaluated in an algebraic extension of the number field in order to handle square roots originating from the solutions of quadratic constraints for on-shell momenta (see section 2.2). The evaluation of cuts in the algebraic extension is slower than in the corresponding field. For this reason, if only amplitudes in Yang-Mills theory are of interest, the option should be left at the default value false.
• gravity-model: Enable/select gravity models. Possible choices are none, Cubic, EH, EH-GB-R3 and all. The last two options increase compilation times considerably. These options give access to the computations of ref. [19], in particular to our implementation of the cubic formulation of the Einstein-Hilbert Feynman rules of [70]. • integrals: Choose whether or not to compile the master integrals of the module IntegralLibrary. The default is none which means that no master integrals are compiled. If goncharovs is selected, GiNaC [63] is required. If pentagons is selected the PentagonLibrary [47] is required. 3 The choice all compiles both representations of the master integrals.
• lapack-path: If necessary, specifies the path of Lapack if meson is unable to find the path to the library.
Beyond these, there are a number of options mainly useful for development: • caravel-debug: Enable a dynamic handling of debugging information from specific source files. To use this feature, simply place a file named debug.dat in the directory where the corresponding binary is run. The file should contain the filenames of the source files (without the full path to it), which should provide additional debugging information. One filename should be listed per line and lines starting with # are ignored. This option can be set to false (default) or true.
• doxygen: Enable the generation of HTML API documentation. This requires Doxygen [72]. This option can be set to false (default) or true.
• ds-strategy: Select the algorithm for the reconstruction of the dependence of the two-loop amplitudes on the dimensional regulator D s (see section 2.4). Possible values are decomposition (default), referring to the decomposition by particle content, and sampling, referring to the reconstruction of the D s polynomial coefficients from the sample values.
The former provides a significant efficiency boost so we recommend not to change its default value. The option decomposition is currently not supported for gravity amplitudes.
• instantiate-rational: Enable selected computations with rational numbers. This option requires finite-fields set to true. This option can be set to false (default) or true.
• precision-arbitrary: Enable computation with arbitrary-precision floating-point types included in the GMP and MPFR libraries. This option can be set to false (default) or true.
• timing: Enable the printout of the time spent in different contributions to amplitude's calculations at the end of each program. This option can be set to false (default) or true.
Enabling some of these features introduces dependencies on third-party libraries, which the user should make available before the start of the building process. Since certain options may significantly increase build times, we suggest to enable only the features necessary for each calculation. Depending on the chosen configuration, building Caravel can take from a few minutes up to half an hour on a modern multi-core processor.

Example programs
In this section we introduce the example programs provided with this release, which demonstrate the main features of Caravel. For each one we specify the required configuration setup, a brief explanation of the computations it performs and instructions for execution. More information about these programs can be found in the file examples/README.md, contained in the Caravel repository. All programs can be found in build/examples, where build is the build directory created in listing 1. Before turning to the description of each example, we first establish some conventions and explain the general structure of the command-line input that the user must provide.

Helicity amplitudes
Each particle q in a multi-particle helicity amplitude is labelled by the particle type f q (here, f q can be a gluon, an (anti)quark, or a graviton), the helicity state h q and the color index c q (for color-charged particles). An n-particle amplitude depends on all this data, that is We assign a momentum index q to particle q. All particles and their momenta are considered outgoing. Note that the color indices c q are not present in pure graviton amplitudes. We consider the perturbative expansion of the (bare) helicity amplitudes and write where a 0 is a generic bare coupling constant and λ denotes the power of the leading-order amplitude. For QCD amplitudes, the expansion is in powers of the strong coupling, i.e. a 0 = g 0 , where g 0 = √ 4πα s . For Einstein gravity amplitudes, instead, a 0 = κ 0 /2, where κ 0 = √ 32πG N and G N is Newton's constant. In this section we will often refer to the index L in M (L) n as the loop order, as for the examples we will consider the two numbers are aligned. L = 0 corresponds to tree-level amplitudes.
As already stated in section 2, see in particular eq. (1), Caravel computes the coefficients A n we must specify several of our conventions and this is done in detail in Appendix A. These conventions are important for defining the output of the example programs described in this section.
The example programs we provide compute tree-level, one-loop and twoloop amplitudes. While the tree-level example program allows one to evaluate amplitudes with an arbitrary number of particles, the one-loop and two-loop example programs evaluate at most five-point amplitudes, up to order ǫ 2 for one loop and ǫ 0 for two loops. The integral normalization for the amplitudes computed numerically in the example programs is defined by equation (10). We stress that this differs from the internal normalization in eq. (8). Indeed, the integrals in the example programs are normalized as Additionally, the QCD loop amplitudes are evaluated in the leading-color limit of QCD. In this limit, we keep the leading term of eq. (1) in the limit of a large number of colors N c , but consider the ratio N f /N c to be fixed, where N f is the number of massless flavors. The leading-color amplitudes have a decomposition in terms of powers of this ratio, specifically

Specifying program input
Many of the programs that we provide can be run for a variety of different scattering amplitudes. To evaluate different amplitudes, we provide a uniform interface by passing (arbitrarily-ordered) command-line arguments. These arguments are formatted similarly to Mathematica lists with heads.
Specifying particles. In the command-line interface, a particle is specified as

Particle[field, index, state]
A (momentum) index needs to be provided for each particle, starting from 1 up to the number of particles in the scattering process. Gluons are specified by setting the field to gluon, and state should be set to either p (+ helicity) or m (− helicity). For massless quarks, Caravel offers multiple possibilities. Generic unflavored quarks and their anti-particles are input as q and qb respectively. Fields of definite flavor can also be specified with u,d,c,s,b and ub,db,cb,sb,bb. We note that scattering processes with identical fermion flavors are currently not supported. However, these can be obtained from anti-symmetrizing distinct flavor amplitudes (see for example [73][74][75]). The associated states for fermions are labeled with qbp (+ helicity) and qm (helicity). The graviton field is labelled by G and its polarization states are given by hmm (--helicity) and hpp (++ helicity). In table 1 we summarize the available field and state options to define particles in Caravel.  where each particle is defined as in the previous section. Note that colorordered amplitudes are invariant under cyclic permutations of the external particles. For gravity the same interface is used, however the ordering of the external gravitons does not matter.
For loop amplitudes, the partial amplitude A (L)[k] of equation (12)    Specifying kinematics. By default, most of the examples we provide evaluate amplitudes on phase-space points which allow one to reproduce the results of ref. [10]. Nevertheless, the user can request evaluations at different phase-space points. Since internally most of the example programs perform calculations in a finite field, one has to be sure that the momenta associated with the chosen phase-space point can be represented in a finite field. Finding such points is in general a non-trivial problem. For all the examples we will be concerned with, however, it can be solved by using a special parametrization of phase-space, called momentum twistor parametrization [68]. In Appendix B, we give more details on this parametrization for four-and five-point massless kinematics. In particular, we give equations that relate our choice of twistor parameters to Mandelstam variables s ij = (p i + p j ) 2 , where p i denotes the momentum of the external particle with particle index i.
To evaluate amplitudes at a chosen phase-space point, the user should provide a list with the head TwistorParameters. For four-point kinematics, the phase space is directly parametrized by the Mandelstam variables s12 and s23 which are passed as For five-point kinematics the parametrization is given in terms of the 5 independent twistor parameters x0, . . ., x4, see Appendix B for more details. These are passed as TwistorParameters[x0, x1, x2, x3, x4]

Numerical amplitude evaluation
In this subsection we present a series of programs which allow the numerical evaluation of a number of QCD scattering amplitudes at tree level, one loop and two loop as well as graviton amplitudes at tree level.

Tree level
The program treeamp evaluates tree-level amplitudes for a variety of processes. It can be executed by specifying the corresponding amplitude with a Particles list (see section 5.2). For example: The program randomly generates a phase-space point and prints the value of the specified tree-level amplitude as well as the point. In Appendix A.2 we provide details about the normalization employed for external helicity states, which are necessary to specify our phase conventions.
If treeamp is to be used to evaluate N-graviton scattering amplitudes, Caravel must be configured with the option -D gravity-model=Cubic Instead of evaluating the amplitudes on randomly-generated phase-space points the user can also provide external momenta by placing a file with the name treeampPSP.dat into the same directory as the executable. The file should list one momentum per line in the format: E PX PY PZ treeamp will find this file, read the momenta and, after performing on-shell and momentum-conservation checks, will compute the corresponding amplitude.
If Caravel has been configured to include high-precision floating-point types (with the option precision-QD set to HP, VHP or all), then highprecision amplitudes can be computed by passing an additional commandline input
Finally, the program can be used to compute tree amplitudes using finitefield arithmetic (-D finite-fields=true required). In this case the cardinality of the finite field has to be passed as a parameter to the program > ./treeamp "Particles[...]" "Cardinality[p]" where p is a prime number smaller than 2 31 and larger than 2 30 , such as p = 2 31 − 1 = 2147483647. For finite-field evaluations, the program does not randomly generate phase-space points and so the user must provide a set of valid external momenta using a file named treeampPSP.dat as explained above. Notice that in the case of finite-field evaluations we use a spacetime metric with alternating signature (+, −, +, −) to enhance performance. Thus, the momentum components should be provided in this signature. For example, using treeampPSP.dat we can pass the rational phase-space point corresponding to s 12 = −3/4 and s 23 = −1/4. The program converts the rational momentum components into their image in the chosen finite field.

One-loop amplitude to O(ǫ 2 )
The program amplitude evaluator 1l numerically evaluates one-loop four-and five-parton helicity amplitudes up to order O(ǫ 2 ) in the dimensional regulator. The master-integral coefficients are rationally reconstructed from finite-field evaluations. The integrals are then computed (in double precision) to obtain the numerical value for the amplitude. If the corresponding tree-level amplitude is non-vanishing the one-loop result is normalized by this tree-level amplitude. Otherwise, the result is normalized by a spinor weight as defined in Appendix A.3, see in particular eq. (A.17). The Laurent expansion of the amplitude is printed to the standard output.
By default, the program runs on the phase-space points defined in ref. [10] so that the user can easily reproduce the results in tables 3 and 4. In the current implementation, this example program runs only in the Euclidean region of phase space as Caravel does not include the analytic continuation of the integrals. We note that the evaluation of five-point amplitudes takes a considerable amount of time because of the evaluation of the one-loop pentagon integral through order ǫ 2 .
In order to enable this example, Caravel has to be configured with the following options: The program can be executed by passing the appropriate amplitude input and kinematic point specifications. evaluates the one-loop color-ordered helicity amplitude A (1)[0] (1 − g , 2 − g , 3 + g , 4 + g ) at s 12 = − 1 3 and s 23 = − 1 5 .

Two-loop amplitude
The program amplitude evaluator 2l numerically evaluates two-loop four-and five-parton helicity amplitudes to O(ǫ 0 ) in the dimensional regulator. It works in the same way as the one-loop program described above. The master-integral coefficients are rationally reconstructed from finite-field evaluations and subsequently combined with the master integrals into a semianalytic object, which is expressed in terms of unevaluated special functions. Next, those special functions are evaluated to produce the amplitude whose Laurent expansion is then printed to the standard output. The computation of the master-integral coefficients is performed in parallel, using all available threads in the CPU, then the special functions are evaluated sequentially in a single thread. As at one-loop, the amplitude is normalized either by the corresponding non-vanishing tree-level amplitude or by the spinor-weight defined in Appendix A.3, see in particular eq. (A.17). On the first evaluation of each amplitude, the program performs a warm-up run, as described in sections 2.2 and section 2.4.
By default, the program runs on phase-space points defined in ref. [10] and can be used to reproduce the results presented in the tables 1 and 2. As at one-loop, this example program is restricted to the Euclidean region of phase space, as master integrals are so far included only for this region.
In order to enable this example the Caravel library has to be configured with the following options: The program can be executed in an analogous way to what was described in the previous example in section 5 evaluates the two-loop color-ordered helicity amplitude A (2) [1] (1 − g , 2 − g , 3 + g , 4 + g ) for s 12 = − 1 3 and s 23 = − 1 5 . The evaluation of two-loop amplitudes is considerably more involved than the corresponding one-loop amplitudes. For example, the runtime of the most complex two-loop five-parton amplitude is of the order of 12 minutes to rationally reconstruct all master-integral coefficients on the default phase-space point (using 22 finite-field evaluations), while the computation of the pentagon functions in double precision by the external library provided with ref. [47] takes about 23 minutes (employing a modern 12-core Intel i7 processor).

Five-point two-loop finite remainder
The program finite remainder 2l numerically computes the finite remainder of planar two-loop five-parton amplitudes. The numerical calculation of finite remainders was instrumental in order to reconstruct the analytic form of the planar two-loop five-parton amplitudes in refs. [16,17]. The program proceeds by building the requested two-loop amplitude and its corresponding infrared subtraction, which also requires the one-loop amplitude to O(ǫ 2 ). The finite remainder obtained in this way is decomposed in terms of the special functions introduced in ref. [47]. As in the two-loop amplitude evaluation example, the computation of the two-loop master integral coefficients is performed in parallel, using all available threads in the CPU and, on the first evaluation, a warm-up run will be performed for the twoloop amplitude. The program automatically verifies that all poles in ǫ cancel exactly, rationally reconstructs the special function coefficients in the finite remainder and subtraction term, and evaluates the special functions. The program then prints the values of the subtraction term, the remainder and the amplitude. We refer to ref. [17] for the precise definition of our subtraction conventions. We note that the run-time of this program is comparable to that of the two-loop amplitude evaluation program.
In order to enable this example Caravel has to be configured with the following options: This example can also be enabled with the integrals option set to -D integrals=all. in order to compute the finite remainder for the color-ordered two-loop amplitude A (2) [2] (1 + g , 2 + g , 3 + g , 4 + g , 5 + g ). As in the two previous examples, the phase-space point can be specified in the command-line input as a list with the head TwistorParameters. The optional argument Verbosity[All] or Verbosity[Remainder] can be passed in the command line to request the printing of additional information.

Analytic reconstruction of amplitudes
An important application of the numerical techniques for amplitude calculation that Caravel provides is the reconstruction of analytic results from numerical evaluations. Here we present two programs which reconstruct either two-loop four-parton or one-loop five-parton amplitudes from numerical evaluations. These examples rely on MPI for parallelization.

Program output
The two analytic reconstruction example programs share common output features, which we describe here. Both programs perform the analytic reconstruction over a single finite field and attempt to rationally reconstruct the result. The computation in the finite field is saved as a text file in the local directory analytics/amplitudes XY/, where XY refers to the relevant cardinality. The master-integral coefficients are then rationally reconstructed employing only the single finite-field evaluation performed. This rational reconstruction is cross checked against a numerical computation at a single phase-space point in a second finite field. In the case of a successful reconstruction, the amplitude is saved in a text file under analytics/amplitudes rational/. The reconstructed amplitudes are normalized either by the corresponding non-vanishing tree-level amplitude or by the spinor weight defined in Appendix A.3, see in particular eq. (A.17).
The files produced by each program are named according to the requested PartialAmplitudeInput. They contain a string that encodes the decomposition of the amplitude as a linear combination of master integral coefficients and master integrals. 4 Each master integral is specified by a string of the form Here, Topology is a human readable name for the topology, e.g. Triangle or DoubleBox. The list of Di is the list of inverse propagators of the master integral, written in terms of the loop momenta (l1, l2) and external momenta (k1, k2, . . .) of the amplitude. Our conventions in Caravel are that all external momenta are outgoing. The label NumeratorLabel denotes the numerator of the master integral. For the example programs provided with this release, there are four possible values for this label, which we list in table 3. For one-loop integrals, these include numerators built from powers of µ 2 , which is the scalar product of the (D − 4) dimensional components of the loop momentum, i.e. In order to aid reading of the output, we provide a collection of Mathematica routines in math/CaravelGraph.m which produce graphical repesentations of the integrals. The output format as described above is not appropriate for these routines. The example programs therefore also writes a textfile in the directory analytics/integral info/, which contains a list of replacement rules that allow one to use these routines.

Univariate amplitude reconstruction
As a simple example of the two-loop analytic reconstruction capabilities of Caravel, the program 4parton 2loop analytics MPI analytically computes the reduction to master integrals. The massless four-point amplitudes depend only on the Mandelstam variables s = (p 1 + p 2 ) 2 and t = (p 2 + p 3 ) 2 . By setting s = 1 and x = t/s the amplitude depends only on a single parameter. Its analytic dependence on x is reconstructed from exact numerical evaluations of the master-integral coefficients over a finite field, which are then fed into Thiele's interpolation formula. The dependence on s is then recovered by dimensional analysis.
To enable this example, Caravel has to be configured with the options -D finite-fields=true -D field-ext-fermions=true

Multivariate amplitude reconstruction
The program computes the analytic form of five-parton one-loop amplitudes, using multivariate functional reconstruction. The five-parton amplitudes depend on five twistor parameters x 0 , . . . , x 4 , see Appendix B. The problem is reduced to a four-dimensional reconstruction by setting x 4 = 1. The dependence on x 4 is recovered from dimensional analysis.
In order to enable this example, Caravel has to be configured with the options -D finite-fields=true -D field-ext-fermions=true To compute the one-loop five gluon all-plus-helicity amplitude, for example, the program should be executed with

Conclusions
We have presented Caravel, a C ++ framework for the computation of multi-loop amplitudes through the multi-loop numerical unitarity method. This is the first publicly available program of its kind. We have provided a series of example programs which showcase the main functionalities of Caravel. In particular, these examples give access to all details of the calculation of the planar two-loop five-parton scattering amplitudes [16,17] and the two-loop four-graviton scattering amplitude in Einstein gravity [19].
In its current form, Caravel is not meant to be able to automatically evaluate arbitrary two-loop multi-leg amplitudes. Rather, it is meant as a framework to do so, which should be complemented with process-specific information such as, for instance, the master integrals. The modular fashion in which Caravel is constructed allows one to easily add these new features. We aim to continue the development of Caravel, to increase the pool of multi-loop amplitude computations that it can perform, while also extending its autonomy.
Tree-level QCD amplitudes. In order to define in a unified way the colorordered tree-level amplitudes that we compute, we shall consider QCD with all fields (that is, both quarks and gluons) in the adjoint representation. In this 'adjoint QCD' theory, amplitudes with any number of partons can be expressed in terms of (n − 2)! color-ordered partial amplitudes using the decomposition of ref. [77], C(a n−1 , a σ 1 , . . . , a σ n−2 , a n ) × A (0) (n−1 h n−1 p n−1 , σ 1 hσ 1 pσ 1 , . . . , σ n−2 hσ n−2 pσ n−2 , n hn pn ). (A.1) Here, S n denotes all permutations of n indices and C is a color structure given by C(a 1 , . . . , a n ) = F a 1 a 2 x 1 F x 1 a 3 x 2 · · · F x n−4 a n−2 x n−3 F x n−3 a n−1 an , (A.2) with the F abc defined above. We stress that we simply use adjoint QCD to define quark and gluon amplitudes in a unified way. It is a well understood procedure to assemble the multi-parton QCD amplitude from these colorordered amplitudes, see for example [67,77,78].
Leading-color QCD loop amplitudes. Beyond tree-level, the example programs compute the color ordered amplitudes relevant for the leading-color limit of QCD. In this limit, we keep the leading term for a large number of colors N c , but consider the ratio N f /N c to be fixed, where N f is the number of massless flavors.
We first discuss the four-point amplitudes and consider amplitudes for the scattering of four gluons, one quark pair and two gluons, and two distinct quark pairs. In the leading-color approximation we write

Four-point kinematics
In order to construct four-point rational momenta we use the momentumtwistor parametrization of ref. [80], which is given by With this parametrization the resulting momenta are rational in s 12 and s 23 .