Functional directed acyclical graphs for scattering amplitudes in perturbation theory

I describe a mathematical framework for the efficient processing of the very large sets of Feynman diagrams contributing to the scattering of many particles. I reexpress the established numerical methods for the recursive construction of scattering elements as operations on compact abstract data types. This allows efficient perturbative computations in arbitrary models, as long as they can be described by an effective, not necessarily local, Lagrangian.


Introduction
The efficient and reliable computation of scattering amplitudes for many particles in a large class of models, both on tree level and including higher order corrections, is a central element of all efforts for analyzing the physics at LHC and possible future colliders.
Since the first release of Madgraph [1] about 30 years ago, there has been tremendous progress in the capabilities of the tools that can compute such scattering amplitudes numerically.Replacing sums of Feynman diagrams by recursive numerical evaluation opened the realm of many-legged amplitudes, including loop corrections.In fact, the treatment of QCD corrections has matured so much that tools like Madgraph5 [2] are now employed regularly by endusers for LHC physics.Electroweak radiative corrections are starting to become available in user friendly tools and recursive techniques are being applied to loop calculations.At the same time, de facto standard formats like UFO [3,4] allow the specification of almost any physics model that might be of interest in the near and not so near future.
In this paper, I will elaborate a common mathematical structure behind the recursive calculations.The focus is not on the immediate numerical evaluation, but on the elucidation of an algebraic structure that will later be translated into numerical code.This simplifies supporting more general interactions, because purely numerical codes have to make assumptions that can turn out to be hard to relax later.In addition, algebraic expressions can be used to generate more comprehensive tests of models and implementations.They also simplify the automatic generation of the additional expressions needed for subtractions schemes [5].
Finally, at a time when functional programming and strong type systems are moving more and more from academia into the mainstream, it is a useful exercise to reconstruct the mathematical structures in a way that can easily be translated into efficient programs making use of these paradigms.The mathematical structures presented here have not been developed in a vacuum, but are a distillation of commonalities observed in the concrete data structures implemented for the matrix element generator O'Mega [6] that is part of the Whizard event generator [7].
Nothing in the following discussion will be specific to leading order, tree level matrix elements.Exactly the same structures appear when implementing loops using additional legs [8][9][10][11] or when adding higher order contributions as terms in an effective action using a skeleton expansion.The translation of the algebraic expressions into robust numeric code calling sophisticated external libraries for loop integrals [12] is much more challenging, of course.However, also here the algebraic step offers more options than a purely numerical approach.
The outline of the paper is as follows: in section 2, I briefly review the recursive techniques used for computing scattering amplitudes for processes with many external particles.This section also serves the purpose of establishing the terminology and notation used in the remaining sections.In section 3, I introduce Directed Acyclical Graphs (DAGs), bundles and their relationships.In section 4, I present an algorithm for efficiently constructing the DAGs representing scattering amplitudes.In section 5, I briefly describe how to generate efficient numerical code from DAGs constructed according to the algorithm presented in the previous sections.In appendix A, I sketch the implementation of DAGs and bundles in O'Mega [6,7].

Scattering Amplitudes
It has long been recognized that the textbook representation of scattering amplitudes as a sum of Feynman diagrams becomes very inefficient as the number of external particles rises.Indeed, even though general estimates are hard to derive for realistic models with conserved quantum numbers, analytic formulae for toy models and explicit calculations for specific processes confirm the expectation that the number of tree level Feynman diagrams grows factorially with the number of external particles.If Feynman diagrams with loops are represented by tree diagrams [8][9][10][11], each loop adds two more external particles.In addition to requiring prohibitive computational resources, the destructive interferences inherent in gauge theories lead to a loss of precision if too many terms are added.Starting with 2 → 6 processes at tree level, the need for a more efficient representation became evident.
In order to simplify the notation in this section, I will cross all scattering amplitudes from n in → n out to n = n in +n out → 0. Except for the momentum, I will also suppress all quantum numbers in this introductory section.The treatment of general quantum numbers (spin, flavor, color, etc.) will be the focus of the following sections.

Recursion
The appropriate building blocks to replace Feynman diagrams turned out to be k-particle matrix elements of fields which will be referred to as wavefunctions or of their associated currents as pioneered by Berends and Giele [13].The set of indices I = {i 1 , i 2 , . . ., i k } is a subset of the indices enumerating the external particles or open loop momenta.Since I ∈ 2 {1,2,...,n} , the number of possible different wavefunctions or currents grows only as an exponential 2 n instead of a factorial n! ∼ n n .Furthermore, both can be computed recursively without expanding them into Feynman diagrams, which would reintroduce factorial growth.In (2), P I denotes a propagator and V I,I1,I2 a vertex factor for three legs.The generalization to models containing vertices with more than three legs is obvious.
Note that φ is just j multiplied by a momentum space propagator.Thus the choice between the two is only a matter of convenience.The rest of the paper will mostly refer to wavefunctions (1a), but all constructions can be repeated trivially for the currents (1b).

Topologies
There are many ways in which a scattering amplitude M can be constructed from (1).The first approach observes that the j in (1b) is already amputated.It therefore suffices to set the momentum on the mass shell of particle 1 to obtain the scattering amplitude using the LSZ prescription This is implemented numerically in Helac [14,15] and Recola [10,11].
The second approach glues the φ from (1a) at vertices to obtain the scattering amplitude in the form with obvious generalizations to models containing vertices with more than three legs.The partitions (I 1 , I 2 , I 3 ) of the external particles must be chosen carefully to avoid double counting [6,16] and the keystones K correspond to vertex factors.This approach was pioneered for numerical calculations in the standard model by Alpha [16][17][18] and is implemented as an algebraic algorithm for arbitrary models in O'Mega [6,7].
The third approach combines the DAGs at propagators M({1, 2, . . ., n}) = I∪I ′ ={1,...,n} j(I)P I,I ′ j(I ′ ) (4c) instead of vertices as in (4c).It was pioneered by Comix [19] and OpenLoops [8,9] Algebraically, all expressions (4) will give the same final result, but the number of nodes that need to be evaluated can vary slightly and numerical results will differ due to the different order of evaluation.O'Mega [6,7] allows to compute the amplitude both as (4a) and as (4b) and confirms these expectations.
While it is impossible to give general estimates for the number of wavefunctions that need to be evaluated in realistic models, one can count them for some examples using O'Mega.In the standard model, it appears that (4a) requires some 10% fewer evaluations than (4b) in an optimal implementation.One advantage of (4b) and (4c) is that at most n/2 of the external momenta appear in the φ compared to n − 1 for (4a).Therefore fewer steps with accumulating floating point errors are required in the recursive evaluation of φ(I).While this could in principle be a significant advantage in amplitudes with strong gauge cancellations, the difference appears to be small in practice.
The algorithm adding quantum numbers to (1), ( 2) and (4) described in the following sections is equally applicable for all three variants in (4).

Evaluation
In the case of a fixed physics model with a moderate number of fields and couplings, such as the standard model, the recursive evaluation (2) can be expressed as an iteration of matrix multiplications [10,11,[14][15][16][17][18][19].This approach has the advantage that all scattering amplitudes in the supported model can be computed using the same executable, without the need for recompilation.
However, extending this approach to more complicated models, in particular to models that can be specified by endusers in formats like UFO [3,4], is far from trivial [20].Instead, it is beneficial to represent the recursion relations (2) abstractly by a data structure from which dedicated code can be generated and compiled subsequently, following the pioneering treatment of Feynman diagrams in Madgraph [1,2].This approach has been implemented for the recursive evaluation in [6][7][8][9].
This motivates the search for a data structure that represents the recursion relations (2) concisely and can be constructed efficiently from the Feynman rules of a model.The obvious candidate is a finite Directed Acyclical Graph (DAG), that corresponds to the evaluation of an arithmetical expression in which common subexpressions are evaluated only once and later recalled from memory when needed again.
Additional benefits of algebraic manipulations are that it is easier to prune the computation of wavefunctions that are not needed in the final result, that one can target special hardware or dedicated virtual machines [21] that avoid the need for compilation.Formfactors can be restricted to lightlike momenta at compile time [22][23][24].Finally, one can optionally instrument the code with numerical checks of Ward and Slavnov-Taylor identities for gauge boson wavefunctions, in order to test matrix element generator, numerical libraries and model descriptions.

DAGs and Bundles
In this section, I will focus on universal mathematical constructions, not practical algorithms.The discussion of the latter will follow in section 4.
Given a set N of nodes, a set E of edges and a set C(N ) ⊆ 2 N of children which is typically the set of subsets of nodes with a limited number of elements, any map from N to the powerset of defines a Directed Graph G = (N, E, ∆) in the sense described below.The function ∆ can be specified completely by the set { (n, ∆(n)) | n ∈ N } of ordered pairs.This equivalence will be used below to define transformations on DAGs as set theoretical operations that can be implemented efficiently in computer programs.I will often employ the more intuitive notation In order to avoid excessive nested superscripts, I will sometimes use the notation A → B for set B A of all functions from the set A to the set B.
With wavefunctions as nodes and vertex factors as edges, this definition captures the recursion relations (2) exactly.Note that the map ∆ ( 5) is well defined iff the combination of momenta and other quantum numbers identifies the wavefunctions or currents uniquely.
There are cases where physical quantum numbers are not sufficient to distinguish wavefunctions.For example, if the scattering amplitude is to be expanded in the powers of some coupling constants, these powers can contribute at different levels of the recursive expansion.Therefore a wavefunction can appear more than once with the same momentum and physical quantum numbers.This forces us to add the powers of these coupling constants as unphysical labels that will be combined in the final step (4).Since such labels are later required anyway to disambiguate variable names in the generated numerical code, this adds no additional burden.Such a counting of coupling constants is of course crucial for adding a consistent number of counterterms in calculations involving loops and when adding precomputed loops using a skeleton expansion or effective actions.
The nodes in the preimage of ∅ under ∆ are called leaf nodes and correspond to the external states in scattering amplitudes.Since the elements of C are sets of elements of N , we can derive from ∆ two mutually recursive expansion functions element-by-element and If D = (N, E, ∆) represents an acyclical graph, i.e. a DAG, with a finite number of nodes |N |, the functions ∆ * and ∆ will reach a fixed point after a finite number of steps.This fixed point consists exclusively of mutually nested sets of leaves.If the image of ∆ consists only of singleton sets and ∅, the fixed point reached from any starting node n corresponds to a tree diagram.Otherwise it corresponds to a forest of tree diagrams, if the elements of the sets are distributed recursively.
As an illustration, consider the DAG D with the sets and the map where I have not spelled out the unlabeled edges.
A quick calculation gives This corresponds to the forest consisting of the trees This DAG encodes a stripped down version of the Feynman diagrams for the process e + e − q q → g, that ignores both the details of the couplings and the contributions of Z and Higgs bosons.Note that the common subdiagram e + e − → γ appears only once in the DAG as 5 → {{1, 2}}, but twice in the forest.
A general directed graph can contain cycles and the functions ∆ * and ∆ will not reach a fixed point even if |N | < ∞.As described in section 4.1, it will however always be possible to equip N with a natural order so that the application of ∆ acts strictly decreasing with respect to this order.There can obviously be no cycles and the graph is guaranteed to be a DAG in this case.
If the same node n appears many times in the children, a DAG provides a very efficient encoding of large sets of graphs.The storage and computing time required by typical sets of tree diagrams grows factorially with the number of leaves |L|.In contrast, the space and time required for implementing the DAG scales linearly with |N |, which only grows as an exponentially in |L|.Using persistent functional data structures [25]

Constructing DAGs
Using DAGs as a compact representation has only a marginal benefit if their construction requires the generation of all tree diagrams in intermediate steps or if the applications require a full expansion.Fortunately, the sum of Feynman diagrams encoded in the DAG can be evaluated either using the DAG directly or by generating a dedicated numerical code that evaluates each node n ∈ N only once.As explained in section 4, it turns out that the DAGs representing perturbative scattering amplitudes can be constructed without requiring the construction of the corresponding forest.
For this purpose, I introduce the empty DAG where ∆ = ∅ is the function with empty domain and codomain.I also define a function with the function ω n →(e,c) that adds a node n together with the mapping n → (e, c) where e and (e, c) are shorthands for the sets { e i | i ∈ I } and { (e i , c i ) | i ∈ I } with |I| elements.In particular, they may be empty to allow inserting a leaf node.In order to avoid ambiguities in the definition of ω, I will require that n ∈ (13) where the function applications associate to the right, of course.It is obvious that any finite DAG can be constructed by repeated applications of ω.
For the finite DAGs that are the subject of this paper, the function ω can be implemented easily in programming languages that have efficient support for persistent sets and maps (also known as dictionaries) that can grow without a lot of reallocation.Functional programming languages with garbage collection make such implementations particularly straightforward.The domain and codomain of functions like ω (12) are highly structured sets and static type systems allow to verify already at compile time that only matching functions are being composed.Beyond preventing errors, a strict type discipline helps to uncover mathematical structures, such as the ones described in this section.This paper is based on the implementation in the matrix element generator O'Mega [6] using ocaml [26], as described in appendix A.1.

Lattices of DAGs
For our purposes, DAGs representing scattering amplitudes for the same external states, categories of DAGs that share the same leaf nodes are the most interesting.Since we describe a DAG as a tuple of sets, there is a natural notion of inclusion for pairs of DAGs in It is obvious that this notion of inclusion corresponds to the inclusion of the sets of tree diagrams encoded by the DAGs.
In the same fashion, we can define union and intersection for the DAGs where and in I am careful to avoid adding new leaf nodes to the intersection.
From these definitions, it is obvious that ⊆ turns D L into a partially ordered set and ∪ and ∩ turn it into a lattice.From this point of view, D 1 ∪ D 2 is the least common upper bound of D 1 and D 2 , while D 1 ∩ D 2 is their greatest common lower bound.Finally D L is bounded from below, with as the bottom element.

Mapping and Folding DAGs
The most important functions for manipulating DAGs and extracting the information encoded by them are folds that perform a nested application of a suitable function for all nodes to a starting value x where the elements of ∆ = {δ n1 , δ n2 , . . ., δ n |N | } are arranged in the partial order if the nodes that guarantees acyclicity of the DAG.The only constraint on the function is that the domain and codomain of f δ : X → X must be identical.The computational cost scales with the size of the DAG and not with the size of the forest of tree diagrams described by it.Used with the constructor ω (12) on the empty DAG, the fold performs a complete copy of any Precomposing the first argument of ω in (20) with a function in the first argument using the notation which can encode a different set of tree graphs.
The precomposition (22) can naturally be extended to functions mapping nodes to sets of nodes as with the identity iff the result of f is the empty set ∅.
Finally, I define a function that takes a set S ⊆ N of nodes and a DAG and returns the minimal DAG that contains all the nodes in the set such that the mutually recursive evaluation of the functions ∆ * and ∆ from ( 7) is well defined for the nodes in this set.Intuitively, this corresponds to following all chains of arrows in { n → ∆(n) | n ∈ N } from D that start in S.

Bundles
I am interested in maps between DAGs that respect certain structures.In order to describe these concisely, I borrow the notion of bundles from topology and differential geometry.
A bundle B (X, B, π) is a triple consisting of a set X, called the total set, a set B, called the base, and a projection π : X → B. The preimages π −1 (b) ⊆ X are called fibers.The notation π −1 : B → 2 X must of course not be misunderstood as the inverse of π.The fibers are pairwise disjoint and their union reproduces the set X. A section is a map s : B → X for which π • s : B → B is the identity.It corresponds to choosing one and only one element from each fiber.This definition generalizes the trivial bundle where all fibers are trivially isomorphic to F and a section is the parameterized graph s : B → B × F of a function B → F .Bundles formalize equivalence relations on the set X, with the base B as the set of all equivalence classes and π the canonical projection of an element of X to its equivalence class.The composition π −1 • π : X → 2 X maps each element to the set of the members of its equivalence class.Sections correspond to choosing one element from each equivalence class.An illustrative example is equivalence of nodes up to color quantum numbers, where π corresponds to ignoring color.Flavor, coupling constant and loop expansion order can be treated in the same way.
Bundles can be arranged in a sequence However, since the preimage π −1 i is not the inverse of the projection π i , the preimage of a composition of projections is not the composition of the individual preimages, but instead.
As in the case of DAGs, such structures and the operations on them can be implemented for finite sets X straightforwardly in functional programming languages with static type systems and garbage collection (cf.appendix A.2).In particular, it is efficient to add elements to the set X and update the base B and maps π and π −1 immediately.This allows to grow a bundle simultaneously while building a new DAG in order to maintain the relationships to be introduced in section 3.5.

Projections and Preimages of DAGs
Given a DAG D = (N, E, ∆), where the set of nodes N is also the total set in a bundle B = (N, B, π), it is natural to ask if there is a canonical DAG D ′ = (B, E ′ , ∆ ′ ) with the base of B as its set of nodes.First, we observe that every section s of B and map f : where πf is the distribution of π over the nodes together with the application of f to the edges πf (e, (32b) The formula (32) has to be augmented by the prescription that a b for which s(b) is a leaf node in D and therefore ∆ s,f (b) = ∅ is not added as a leaf node to D s,f , similar to the definition (16d) of the intersection of two DAGs.
In most cases f : E → E ′ will be a simple projection that in our applications will be determined straightforwardly by the two sets of Feynman rules governing the construction of the two DAGs.Therefore we can write D s instead of the more explicit D s,f .
The dependence of this projection on the section s is not satisfactory.However, the DAG where S(B) denotes the set of all sections of the bundle B, is well defined and will be shown to suit our needs.Observe that the union is the correct universal construction for our applications, because the additional quantum numbers in N lead to more selection rules.These selection rules are the reason for the dependency of D s on s.
The DAG corresponding to the more basic set of nodes B should therefore be the combination of all possibilities.As an example consider the scattering of two scalars without and with flavor.Without flavor, there will be s-, t− and u-channel diagrams.With a conserved flavor, only one of them will remain.Note however, that this construction does not guarantee that the set of nodes of the DAG Π(D) is actually the full base B of the bundle B. We must therefore demand in addition compatibility of DAG and bundle, by requiring that the diagram commutes.The function ν in the commuting square (34) just extracts the set of nodes from a DAG ν(N, E, ∆) = N .
(35) The objects in the commuting square (34) can be understood as a combination of a pair of DAGs and a bundle, which I will call a fibered DAG.In programs, nodes can be added to the DAG D and the bundle in concert such that the relationship (34) is maintained.
An immediate benefit of such an universal construction of the projection is that it provides a corresponding preimage Π −1 which maps DAGs with the base B as nodes to all DAGs with the set N as nodes.The maps in the preimage can be written ∆ s,f : where ŝf is to be understood as the distribution of s over the nodes together with the application of f to the edges.Unfortunately, in contrast to (32), there will not be a single function f : E → E ′ .Instead, we must allow that ŝf maps into the powerset 2 E ′ ×C(N ) instead of E ′ × C(N ).In addition, the image of f will depend, via the Feynman rules, on the nodes appearing as children.
Since the resulting notation would be unnecessarily cumbersome, I will refrain from making the nature of f in (36) explicit as a function by specifying its domain and writing out all of its arguments.Nevertheless, the discussion of the example in section 4.2 will demonstrate how a set of Feynman rules defines the maps ∆ s,f unambiguously.
In this picture, the application of Feynman rules amounts to choosing a particular element of the preimage Π −1 .It would however be extremely wasteful to construct the preimage first and to throw away all but one of its elements later.In section 4.2, I will describe an algorithm that can be used to construct the desired element directly.
So far, I have assumed that the DAGs are selected by Feynman rules that are local to each vertex in the case of Feynman diagrams or to each element δ i of the map ∆ in our DAGs individually.There are however important exceptions.The most important is provided by loop expansions.There it is required for consistency that counterterms are inserted a fixed number of times in Feynman diagrams.Such conditions on complete Feynman diagrams do not translate immediately to the DAGs, whose components can enter the scattering amplitudes (4) more than once.Fortunately, this problem can be solved by introducing additional unphysical labels representing loop orders to the physical labels of the nodes and to select the required combinations of wavefunctions in (4) at the end, as will be described in section 4.4.The same applies to selecting fixed orders in the perturbative expansions, as required for comparing to many results from the literature.
We call two DAGs In this case D 1 and D 2 can be viewed as refinements of the same basic DAG D. This notion of equivalence generalizes the notion of topological equivalence for diagrams, where two diagrams are considered equivalent if they agree after stripping off all quantum numbers.With the new notion of equivalence, we can say that the sets of Feynman diagrams encoded in a DAG are equivalent up to flavor or upto color.
Using the basic commuting square (34), we can immediately extend the bundle complex (30) to include the corresponding DAGs In our applications, this complex does not continue further to the left, because for each number of leaf nodes there is a natural leftmost nontrivial DAG D P , described in section 4.1 below.
In the following section 4 I will describe how to use Feynman rules to walk the lower row of (38) to construct a DAG for a scattering amplitude efficiently in stages.

DAGs from Feynman Rules
In principle, it is possible to construct the DAG encoding all Feynman diagrams in a single step.
First one adds leaf nodes for external states, labeled by all quantum numbers (momentum, spin/polarization, flavor, color, . . .).Which states are to be included here depends on the choice of algorithm, as has been discussed in section 2.2.
Then one uses the Feynman rules of the model to add all nodes where the node and its children correspond to an allowed vertex.This proceeds iteratively: in the first step all subsets of the leaf nodes appear as children.In the following steps subsets of all nodes, including the leaf nodes appear as children subject to the constraint that no leaf node appears twice if the DAG is expanded recursively with the functions ∆ * and ∆ from (7).This iteration will terminate after a finite number of steps when all leaf nodes have been combined in all possible ways.While this algorithm inserts nodes that will not appear in the scattering amplitude the function H (27) can be used to harvest the minimal DAG.This is a workable approach, but it is neither the most efficient nor particularly maintainable in actual code.Since the nodes are labeled by all quantum numbers, handling them all at once requires the construction of many nodes that will not appear in the final result.Adding quantum numbers in several stages instead allows us to use the constraints from earlier simpler stages to avoid in later stages the construction of many more complicated nodes that will never be used.While not relevant for the final numerical code, experience with early versions of O'Mega [6,7] revealed that the latter approach requires noticeably less time and memory for constructing the code.
Breaking up the construction of the DAG into several stages also simplifies the implementation of each stage and allows separate testing and swapping of different implementations.Finally, applications often need access to projected DAGs as described in section 3.5 anyway.A prominent example is the construction of phase space parameterizations that only refer to kinematical information, such as propagators and masses.Some of the stages described in the following subsections will be performed in a particular order, while the order of others can be interchanged easily.

Momenta
An element of the set N P of nodes in the first DAG D P = (N P , ∅, ∆ P ) to be constructed is uniquely labeled by a subset of the powerset 2 {1,2,...,n} of labels for the external momenta and the edges are unlabeled.The leaf nodes are the elements n({i}) of N P and the action of the map ∆ P is given by where l is the maximum number of legs of the vertices in the model.Obviously, we can order the nodes n(I) according to the number of elements of I to prove that there are no cycles in D P .
In case of (4a), we only need the elements of 2 {2,...,n} as labels.In the cases (4b) and (4c), only labels with at most n/2 elements are needed.Finally, the function H ( 27) is applied to construct the minimal DAG required for evaluating one of the expressions (4).

Flavors and Lorentz Structures
In the next stage, the momenta of the leaf nodes of D P are combined with the flavor quantum numbers of the corresponding external state.The resulting leaf nodes form the starting point of a new DAG D F = (N F , V F , ∆ F ) and bundle B F = (N F , N P , π F ).The edges V F are vertex factors consisting of coupling constants, Lorentz tensors and Dirac matrices.
Using a fold Φ of D P using the constructor ω of D F with precomposition (23) that maintains the fibration (34) will ensure that the nodes of D P are visited in the correct order of growing label sets.The function f that is precomposed to ω in (23) acts on each element of the map ∆ P as follows: since the n(I i ) ∈ N P have been processed, they are elements of the base of the growing bundle B F .Therefore, the fibers π −1 F (n(I i )) are already complete and we can compute their cartesian product We then use the Feynman rules to select all elements of Γ that can be combined with another flavor to obtain a valid vertex.This defines a function Γ → 2 VF .For each of the resulting flavors, a new node labeled by I and this flavor is added to D F and B F together with the corresponding vertex factors and elements of Γ as edges and children, maintaining the fibration (34).This algorithm has been implemented in O'Mega [6,7] and is completely independent of the kind of Feynman rules.It can accommodate both hardcoded rules and rules derived from a UFO file [3,4].The only potential performance bottleneck is the efficient matching of vertices to the elements of a Γ representing a large number of children.For vertices with few legs, this is not a practical issue, but care has to be taken for vertices with many legs where the factorial growth of the number of permutations might be felt.
Once the flavors have been assigned, it is known which fermion lines contribute in the computation of each node.This information must also be added to the node in order to be able to assign the correct sign to interfering contributions in (4) later.Special care must be taken if the model contains Majorana fermions [27][28][29].
By construction, after the fold is complete, the new DAG D F encodes all the information needed to compute the scattering amplitude for the leaf nodes in a theory without color, using one of the formulae (4).Some nodes in D F might not be needed due to conserved quantum numbers.Therefore the function H (27) from D F is applied again to construct the minimal DAG required to evaluate one of the expressions (4).

Colors
Since the color representation depends on the flavor, the assignment of color quantum numbers in the construction of the DAG D C = (N C , V C , ∆ C ) naturally comes after the construction of D F .
We can now follow the steps of the previous stage, as described in section 4.2, word for word, only replacing the subscripts (F, P) by (C, F).The implementation in O'Mega uses the realization of the color flow basis described in [30], but, except for the labeling of the nodes in N C , the form of the vertices in V C and the Feynman rules to be used, the algorithm is completely independent of the representation of the color algebra.
Having the color information available algebraically allows to compute color factors and color correlators [5] analytically.

Coupling Orders
As already mentioned in section 3.5, there are cases where it is important that the Feynman diagrams encoded by the DAG contain certain coupling constants with fixed powers.The most important examples are the counterterms and the terms of an effective action in a loop expansion.Also the inclusion of self energy type terms will not terminate in at DAG, unless a finite maximum expansion order is prescribed.
For practical purposes it is sometimes also important to compute only a part of a scattering amplitude corresponding to fixed powers of couplings.Such results are often available in the literature from Feynman diagram based calculations and a comparison for the purpose of validation is only possible if the DAG based calculation can select exactly the same contributions.
A priori, this conflicts with the representations (4) of scattering amplitudes as DAGs, since the wavefunctions or currents will have accumulated different powers of couplings that will be mixed by (4).
Fortunately, there is a simple solution.For example, in the case of (4b) we can write to compute the scattering amplitude at the coupling order o.The only change required is that the wavefunctions have to keep track of the coupling orders accumulated in their recursive computation.Since the powers of the couplings are additive, we never have to add the wavefunctions or currents that exceed the requested order to the DAG.This necessitates augmenting the set of labels of the nodes by unphysical "quantum numbers" corresponding to the coupling orders.It can be implemented easily, as long as the number of coupling orders to be tracked remains moderate.

Skeleton Expansion
If we are using DAGs to efficiently implement a skeleton expansion, the remarks in section 4.4 apply word for word by replacing "coupling order" by "loop order".

Multiple Amplitudes
In practical applications [7], it is usually necessary to compute scattering amplitudes for the same external momenta, but more than one combination of flavors and colors at the same time.These flavor and color combinations often overlap pairwise and the sets of leaf nodes will also overlap, i.e.L 1 ∩ L 2 = ∅.In this case, it is efficient to combine the corresponding DAGs D L1 and D L2 into a single DAG and to compute the scattering amplitudes from this DAG in order to reuse nodes from the part of the DAG build on L 1 ∩L 2 .For this purpose, we can generalize the union defined in (16c) to a map in an obvious way.
Identically structured code can be emitted as bytecode for a virtual machine that realizes the operators as basic instructions [21].The improving memory bandwidth for graphical processing units even allows to start targeting GPUs for interesting examples.
As already mentioned in the introduction, the generation of robust numerical code is much more challenging if the DAG encodes diagrams that contain loops.The problem has been solved for the standard model [8][9][10][11].The structures described in the paper will help with the task of extending this approach to general models.

Conclusions
I have described the algebraic structures that organize recursive calculations in perturbative quantum field theory without the need to expand intermediate expressions into Feynman diagrams.In functional programming languages, these algebraic structures translate directly into data structures.In a second step, these data structures are translated to efficient numerical code for any programming language or hardware target required.
This algebraic approach adds flexibility over purely numeric implementations tied to specific models and computing targets.It allows for more extensive consistency checks and paves the way for more challenging applications.this interface is more straightforward and is better tailored to our applications.
The function harvest implements H (27).In particular, harvest dag n dag' finds the subset of the DAG dag that is reachable from the node n and adds it to the DAG dag'.This way, applications can compute a minimal DAG for further processing.
Since the construction of the DAG D P (cf.section 4.1) is very simple, it had been combined with the construction of D F (cf. section 4.2) in O'Mega [6] before the structures described in this paper were elaborated.However, the separation of the remaining stages described in section 4 forms the backbone of the current version of O'Mega.

A.2 Bundles
Here is the signature of the Bundle module in O'Mega [6].Again a functor is applied to the types elt, base and the function pi, corresponding to X, B and π respectively module type Bundle = sig type elt type base val pi : elt -> base type t val empty : t val add : t -> elt -> t val inv_pi : t -> base -> fiber val base : t -> base list end The semantics of the functions is evident from the discussion of bundles in section 3.4.Note that π is universal for all bundles with this type, while π −1 depends on the elements added to the bundle previously.
instead of mutable arrays to implement the function ∆ simplifies the algorithm described below significantly.The additional space and time requirements replace |N | by |N | ln |N | and turn out not to be important for large |N |.