ColorMath—a package for color summed calculations in SU(Nc)

A Mathematica package for color summed calculations in QCD (SU(Nc)) is presented. Color contractions of any color amplitude appearing in QCD may be performed, and the package uses a syntax which is very similar to how color structure is written on paper. It also supports the definition of color vectors and bases, and special functions such as scalar products are defined for such color tensors.


Introduction
With the LHC follows an increased demand of exact calculations in QCD involving many color charged partons.Due to the non-abelian nature of QCD this poses a nontrivial computational problem.This is the issue the Mathematica R package ColorMath1 is set out to tackle.The main feature of ColorMath is thus the ability to automatically perform color summed calculations starting from a QCD color structure which is expressed using a syntax very similar to how the color structure would have been written on paper.
ColorMath allows both partial and full dummy index contractions, and thus deals with color structures having an arbitrary set of free indices.In this sense it is thus more general than [1], on the other hand it is a pure SU(N c ) tool, rather than a high energy physics general purpose tool, such as [2,3].ColorMath works for an arbitrary N c ≥ 2, and for arbitrary trace convention of the squared SU(N c ) generators, T R .It may be used by just giving the color structure in the appropriate form and afterwards running CSimplify to contract all repeated indices.For example, consider q 1 q 2 → q 3 q 4 via gluon exchange in the s-and tchannels.Writing down the amplitude as Amp = S t {g}q1 q2 t {g}q4 q3 + T t {g}q1 q3 t {g}q4 q2 (1) where S and T represent some s-and t-channel kinematics, and t {g}q1 q2 the SU(N c ) generator in the fundamental representation, otherwise typically denoted (t g ) q1 q2 , we may calculate the squared amplitude (defined as in the scalar product eq.( 12)) using immediately giving the answer Each squared amplitude can in principle can be calculated like this.However, ColorMath also facilitates the usage of color vectors and bases and has special functions for the calculations of scalar products and gluon exchanges.This paper, which is intended to be the main reference, is organized as follows: First a general introduction to the basic color building blocks is given in section 2. In section 3 the computational strategy is presented, whereas basic examples are given in section 4, and the usage of vectors and their functions are presented in section 5.In section 6 some remarks concerning validation and scalability are made and in section 7 conclusions are drawn.

QCD building blocks
From a color perspective, the QCD Lagrangian is built out of quark-gluon vertices, and triple-gluon vertices, where we follow the convention of reading of the fully antisymmetric structure constant indices in counterclockwise order.The four gluon vertices can be rewritten in terms of (one gluon) contracted triple-gluon vertices, and thus need no special treatment.Due to confinement, we never observe individual colors and it therefore suffices to calculate color summed/averaged quantities for making predictions in quantum chromodynamics.We may thus constrain ourselves to treat QCD amplitudes carrying a set of external indices with values that need never be specified, as they are always summed over in the end.
In principle, the SU(N c ) generators, along with a delta function for indicating a quark and an antiquark color singlet, δ q1 q2 , constitute a minimal set of objects for treating the color structure in QCD2 .For convenience, and for performance reasons, it is, however, useful to define a larger set of objects.The complete set of basic building blocks for carrying color structure used by the ColorMath package is given in table 1.
Apart from a delta function in quarks indices, δ q1 q2 , a delta function in gluon indices, denoted by ∆ {g1,g2} , is also defined.Note the Mathematica List brackets in {g1, g2}.The gluon delta function can alternatively be expressed using where T R (TR) comes from the normalization of the SU(N c ) generators, and is typically taken to be 1/2 (the Gell-Mann normalization) or 1.A rescaling of the normalization of the SU(N c ) generators can always be absorbed into a normalization of the strong coupling constant.To allow for arbitrary normalization T R is kept as a free parameter, denoted by TR, and may be defined by the user.Also the number of colors, denoted Nc in ColorMath, may be set by the user; by default both Nc and TR are kept as free parameters.Note the ColorMath notation for a trace over two gluons o {g1,g2} in eq. ( 6).Similarly, a general trace over k gluons g1, g2, ..gk is denoted o {g1,g2,...,gk} , and may be thought of as a closed quark-line with k gluons attached.
The totally antisymmetric structure constants, which -along with an extra i -define the triple gluon vertices, are denoted f {g1,g2,g3} .Similarly, the totally symmetric "structure constants" are defined as d {g1,g2,g3} .Recall that, starting from the commutation (anticommutation) relations In ColorMath notation t is used to denote an open quark-line, and the above expression is written similarly, = 1 Pictorially this represents where i is included in the triple gluon vertex.The rationale for putting the gluon index in the SU(N c ) generators t {g2}q1 q2 inside a Mathematica List in eq. ( 9) is to allow for the natural extension of having many gluons attached to an open quark-line, thus t {g1,g2,...,gk}q1 The left hand side has several advantages compared to the right hand side.An open quark-line with an arbitrary number of gluon indices can be written in a compact form with no dummy indices.This is not only more human readable, but also superior from a computational point of view, as it avoids the contraction of unnecessary dummy indices.A third advantage with the above notation is its direct correspondence to the trace type bases [4][5][6][7][8][9][10][11][12].A basis (or spanning set) for the color space for a fixed set of external quarks, anti-quarks and gluons can always be taken to be a sum of products of open and closed quark-lines.
In this context we also remark that the color tensors defined in table 1 are color scalars, i.e., they are invariant under SU(N c ) transformations.This imposes no restriction for our purposes as, for any QCD amplitude, the overall color structure, including both incoming and outgoing particles, always is a color singlet.As the basic building blocks are invariant, each tensor built out of these objects, i.e., each tensor needed for color summed calculations in perturbative QCD is a actually a color scalar.
The scalar product on this vector space is given by summing over all external color indices, i.e.
with a i = 1, ..., N c if parton i is a quark or anti-quark and a i = 1, ..., N 2 c − 1 if parton i is a gluon.As long as the color structures in table 1 are multiplied by real coefficients the scalar product is actually real, which is easy to prove using the computational rules in the next section.

Basic computational strategy
Having defined all the color carrying objects, we turn to describing the basic strategy for carrying out calculations.
Again we note that we need not treat the four gluon vertex as this can be rewritten in terms of three gluon vertices.To treat an arbitrary color structure in QCD we may always compute any squared amplitude in the following way: (i) Rewrite the triple gluon vertices using eq.( 8).This results in a color structure which is a sum of products of open and closed quark-lines, connected to each other via repeated gluon indices.(ii) Contract all internal gluon indices using the Fierz or completeness relation (iv) Use the tracelessness of the SU(N c ) generators, and contract quark and gluon delta functions with repeated indices After this, all dummy indices have been contracted.A squared amplitude is thus represented by a rational function in T R and N c , whereas a general amplitude is represented as a sum of products of closed and open quark-lines, i.e., as sums of products of t {g1,g2,...,gk}q1 q2 , including t {}q1 q2 = δ q1 q2 , and o {g1,g2,...,gk} , which for two gluons may be rewritten as o {g1,g2} = TR ∆ {g1,g2} .
By applying these rules, any amplitude square and any interference term appearing in QCD can be calculated [6,11].To successfully square an arbitrary QCD amplitude using the above set of rules we see that steps (i-ii) have to be performed while keeping the relative order, i.e., by first applying rule (i) and then rule (ii).We also note that these rules increase the number of terms, whereas the rules (iii-iv) decrease the number of terms or keep it fixed.In order not to unnecessarily inflate an expression it may therefore be useful to apply the latter rules at any time during the computation.ColorMath utilizes this and tries to contract indices using (iii-iv) at any time, while rules (i-ii) are used only when needed, i.e., when the non-expanding rules fail.
On top of the rules (iii-iv) there are other cases in which the contraction of a gluon results in at most one term.This will happen if (v) Two neighboring gluons attached to the same (closed or open) quark-line are contracted.In this case the result is simply (vi) Similarly, it is easy to show, using the Fierz identity, eq. ( 13), that the contraction of two next to neighboring gluons results in only one term (17 While the above rules are sufficient for squaring color amplitudes, it is sometimes advantageous to contract indices directly using the structure constants f g1 g2 g3 and d g1 g2 g3 .Therefore rule (i) is not always automatically used for simplifying expressions.Instead, in addition to the above set of rules, a limited set of rules for contraction of repeated gluon indices occurring in the standard and the symmetric structure constants are implemented.More specifically, all rules for contracting two gluons in products of two symmetric or antisymmetric structure constants, for example, and all rules for contracting three gluons in products of three symmetric or antisymmetric structure constants, such as are implemented.For more intricate gluon contractions involving structure constants ColorMath can contract indices by first applying rule (i).

Basic calculations
In this section, basic functionality of ColorMath is explored.As always with Mathematica, the package has to be loaded before it can be used.With where the version number is adjusted to the version in question, this can be done.To start using the package, it is, however, recommended to start from (a version of) the tutorial "ColorMathTutorial.nb"and modify its content according to the desired usage.
The color contractions corresponding to the basic manipulations from section 3 are carried out using Mathematica Rules, i.e., a set of replacement rules based on pattern matching.As always, the rules may be applied using "Expr/.TheRules", and may be applied repeatedly using "Expr//.TheRules".
The rules described in (iii-iv) are, along with eq. ( 6) contained in a set of rules called SimpleRules, as they keep the expression at least as simple (in a term counting sense) as it initially was.Applying these to an expression thus tend to simplify it, for example results in TR(∆ {g1,g2} ) 2 .
To fully utilize the rules we (may) need to apply them repeatedly, giving the fully contracted expression (−1 + Nc 2 )TR.The rules defined in (v-vi), acting on t and o, are similarly contained in OTSimpleRules, and the union of SimpleRules and OTSimpleRules and a few rules for rewriting closed quark-lines with zero to two gluons, are contained in AllSimpleRules.
The special rules for gluon index contraction in structure constants, exemplified in equation (18-19) are defined in FDRules.The complete set of rules are stated in table A.2.
Rather than thinking about how individual rules have to be applied, it is convenient to have a standard procedure for contracting color indices.This is embodied in the function CSimplify, which does what its name suggests; simplifies the color structure as far as possible.For color structures which do not contain structure constants this always implies contracting all repeated indices.
For color structure involving the structure constants CSimplify first attempts simplification using the FDRules.If, after this, the expression still contains structure constants, the structure constants are by default rewritten in terms of traces using eq.( 9), and the indices are fully contracted, resulting in a sum of products of open and closed quark-lines.Sometimes it may, however, be desirable not to rewrite the structure constants, as expressions may be more compact if they are kept.This can be achieved by using the option RemoveFD→ False.
The most useful set of functions are given in table A.3.Apart from CSimplify we especially note the function ReplaceDummyIndices for replacing all repeated indices in a color structure with a new unique set of color indices.
For the purpose of calculating amplitudes square, we need, as in eq. ( 12), the complex conjugated version of the color structures.We note that whereas ∆ {g1,g2} , f {g1,g2,g3} and d {g1,g2,g3} are real.This is implemented in ColorMath via redefinition of the Mathematica's built in function Conjugate.With this in mind, we are ready to perform calculations of the type in eq. ( 2).
Finally we remark that to Mathematica ∆ {g1,g2} = ∆ {g2,g1} , and similarly, o {g2,g3,g4,g1} = o {g1,g2,g3,g4} etc.For this reason a function SortIndices, which writes indices in Mathematica default order is defined, s.t. for example This function is used by CSimplify, and by some of the rules in table A.2. Sometimes more detailed control over the calculation may be desired.For this purpose, and for internal usage, functions manipulating indices and probing the color structure are given in table A.4, in Appendix A.
While each squared amplitude in principle can be calculated as outlined above, ColorMath also offers more efficient tools for dealing with vectors in color space.

Defining and using vectors
For the purpose of studying color space it is often convenient to define a basis for the color space, and sometimes also projection operators.Both of these are examples of color (singlet) tensors, and ColorMath has a set of tools for working directly with such tensors.
As an example, let us consider the color structure for q 1 q 2 → q 3 q 4 g 5 3 .A basis for the color space may be written as [13] Vector181 {q1 ,q2 ,q3 ,q4 ,g5 where the basis vectors are labeled using first the overall multiplet of q 1 q 2 , then the overall multiplet of q 3 q 4 , and finally the overall multiplet of q 3 q 4 g 5 , which -due to color conservation -must equal the multiplet of q 1 q 2 , (implying that the notation is somewhat redundant).
From a Mathematica perspective we note a few things.First, on the left hand side, we see that the indices inside the List in the Subscript are followed by underscore to indicate pattern matching.This is standard in Mathematica and makes it possible to use any symbol to denote the indices in later calculations.Then we note that the last two tensors in eq. ( 25) are defined using Module.This is to ensure that each time the tensor is used, it comes with a fresh set of dummy indices.This is also the reason why set delayed ":=" is used.
This basis is orthogonal since at least one set of partons transform under different representations in the various tensors.It is, however, not normalized.Finding the normalization is easy using ColorMath.To calculate for example Vector181 squared we could enter but it is yet much easier to use the tensor functions for calculating scalar products.Instead we could simply write resulting in Nc −1 + Nc 2 TR.Having a basis we naturally want to calculate all norms (square).ColorMath has special functions for this as well.If we define we may calculate the squares of the basis vectors using This results in a List containing the scalar products between each basis vector and itself Often non-orthogonal bases are used, in particular, this is the case for the trace type bases.In this case, to square a general amplitude, all the scalar products between all the basis vectors are needed.For calculating the scalar product matrix the command CDotMatrix may be used.This returns a List of Lists, i.e. a matrix, where the ij-th element is the scalar product between the i-th and j-th vector in the list of (basis) vectors.For larger vector spaces with complicated scalar products, it may be desirable to get progress information on the calculations.This can be obtained by setting the option Verbose → True for Similarly, it may be useful to be able to simplify potential (normalization) roots assuming a large Nc.This can be done by using the option NcMin.The scalar product functions, along with their options, are listed in table A.5.
Often it is also of interest to investigate the effect of gluon exchange on the color structure expressed in a basis, i.e. starting in a basis vector j, what is the effect on the basis vector of exchanging a gluon between parton1 and parton2.This is useful both for soft gluon resummation and for one-loop corrections via gluon exchange.The result can be expressed in terms of a matrix whose element ij is the i-th component resulting after such an exchange in the initial vector j.This is calculated by the function where Basis is a List of basis vectors, defined using the syntax in eq. ( 25) and parton1 and parton2 are the numbers of the partons in the basis vectors, i.e., in this case numbers between 1 and 5.The sign conventions are such that quark-gluon vertex always comes without additional signs, and the triple-gluon vertex have the indices appearing in the order: external index, internal dummy index and index of the exchanged gluon 4 .For example, we may calculate the effect of gluon exchange between parton 1 and 3 in the above basis using Here we have supplied optional information about the basis type, that the basis is orthogonal, to speed up the calculations.By default CGamma does a few consistency checks.It checks that the vector resulting after gluon exchange, when squared, has the same value as the basis decomposed version.For orthonormal bases (BasisType → OrthonormalBasis) it is also checked that the resulting matrix is symmetric [14].These checks may, however be turned off (MakeChecks → False).The set of options, with default values are listed in table A.5.

Validation and scalability
The computational rules and functions in this package have been used for calculating the three gluon projection operators presented in [15].This imposes highly nontrivial consistency checks, as it is verified that every projector square equals itself, which at intermediate steps often involve several ten thousand terms.Additional consistency checks on the color contraction rules have been made by using the CGamma function which checks that vectors square and basis decomposed vectors square agree, and by changing the order in which the color structure contraction rules are applied.The scalar product matrices have also been compared to the ColorFull code [16] for tree level trace bases with up to six partons out of which one parton is a quark and one an anti-quark.Selected results have been compared against [13], and the package has been tested in Mathematica 7, 8 and 9.
The computational effort needed for exact treatment of the color space grows very quickly with the number of partons.The dimension of the vector space grows roughly as a factorial in the number of gluons plus qq-pairs [15] (strictly speaking an exponential for finite N c in a multiplet basis).The computational effort for ColorMath, or any program operating by direct manipulation of quark-lines, tend to grow roughly as the square of this, as the quark-lines are non-orthogonal.ColorMath (in its current form) is thus rather intended to be an easy to use package for calculations of low and intermediate complexity than a competitive tool for processes with very many colored partons.

Conclusion
In this paper a Mathematica package ColorMath for performing color summed calculations in QCD is presented.This package allows for simple evaluation of QCD color amplitudes which are expressed in a format which very much resembles how the color structure would have been written on paper, see table 1.The idea is that the user -for simple cases -just should give the expression, and then run CSimplify[Expr] rather than Simplify[Expr].The package is based on advanced pattern matching rules, and a list of rules is given in table A.2, whereas functions acting on color structures are given in table A.3.
For calculations of intermediate or high complexity it is often beneficial to use a basis for performing color space calculations.ColorMath allows for definition of color tensors of form C1 {i1 ,i2 ,...,ik } := . .., carrying an arbitrary set of quark, anti-quark and gluon indices.Special functions for calculating scalar products, and investigating the effect of gluon exchange, are given in table A.5. ColorMath is, however, not intended for high speed calculations involving many colored partons.For this purpose a separate C++ package is written [16].
Table A.3: The most useful functions, to be used on any expression carrying color structure of the form given in table 1.
Recall to relabel dummy indices, manually, or by using ReplaceDummyIndices, and to place quark indices that are to be contracted such that one index sits upstairs and one downstairs.Furthermore, check that the correct FullForm is used, cf.table 1.The function WhatIsWrong may be useful for identifying common input mistakes in the basic color objects.

Function with usage Effect
CSimplify [Expr] Options with default value: The most general function for simplifying color structure.If Expr contains structure constants, f {g1,g2,g3} or d {g1,g2,g3} , FDRules are first applied.If, after this, the expression still contains structure constants they are -by default -removed using FDToORules, and all repeated indices are subsequently contracted by repeatedly using AllSimpleRules, then OTThenAllSimpleRules, and finally ExpandThenRules.For color structure containing structure constants, it could happen that it is desired not to replace the structure constants.This can be achieved by setting the option RemoveFD → False.

GluonContract[Expr, Gs]
Contracts a set (List) of gluons Gs = {g1, ..., gn} or a single gluon Gs = g1, in the expression Expr, while leaving other indices uncontracted.This function is intended for quark-lines and will replace structure constants with quark-lines.

ReplaceDummyIndices[Expr]
Replaces the dummy indices in Expr with a new set of unique dummy indices.

SplitConstAndColor[Expr]
Splits an expression Expr into a List of Lists of color structures and corresponding multiplicative factors {{constants 1, color structure 1},{constants 2, color structure 2},...}.This is done by first expanding the expression and then splitting the terms separately.

WhatIsWrong[Expr]
Checks if anything is obviously wrong with an expression, for example if Power is used instead of Superscript or if gluon indices are placed downstairs.The check is performed by first expanding the expression, and then checking each term.Read error messages from above.Describes the effect of gluon exchange between partons p1 and p2, on basis vector Cj as a column vector j in a resulting matrix (technically List of List), i.e., element ij is the resulting color structure's i-th component.The BasisType option should assume one of the values GeneralBasis, Or-thogonalBasis, OrthonormalBasis or TraceBasis.The latter is defined to be any basis where the basis vectors consist of one product of closed and open quark-lines, not a sum.By default it is checked that the color tensor, resulting after gluon exchange, is the same when squared directly, as when basis decomposed, and then squared.This checks that the basis is complete and that the basis decomposition is correct.For orthonormal bases, it is also checked that the resulting matrix is symmetric.Using MakeChecks → False these checks are turned off.By default, progress information is also written out, which can be changed with Verbose → False.

Table 1 :
Along with the number of colors, Nc, and the trace of an SU(Nc) generator squared, TR, the basic building blocks are as below.Note that o {g1,g2,...,gk} represents a trace over gluons with indices g1...gk, tr[t g1 t g2 ...t gk ], and that t {g1,g2,...,gk}q1 q2 represent the q1, q2-component in a trace over generators that has been cut open, (t g1 t g2 ...t gk ) q1 q2 .For convenience, also the totally symmetric "structure constants" d {g1,g2,g3} are defined.Note the Mathematica FullForm.ColorMath is built on pattern matching, and it is therefore essential that the expressions have the correct FullForm.In particular, Power may not be used instead of Superscript.To get the right FullForm it is recommended to use the function form, which is just a function returning the corresponding ColorMath object.
the color structure where the two involved gluons have been contracted,

Table A .
4:The below set of functions are used for probing the color structure, and identifying indices.Concerning the convention of what (lower or upper) quark-indices are referred to as incoming quarks and outgoing anti-quarks, as opposed to incoming anti-quarks and outgoing quarks, ColorMath does not specify any preference, but simply refers to the indices as upper quarks and lower quarks when needed.Table A.5: Special functions for color structures expressed in the form of eq.(25), i.e., the corresponding tensors need to be defined using the pattern matching underscores, potential dummy indices should be hidden inside Modules, and the evaluation should be delayed for expressions containing dummy indices, i.e. ":=" should be used.The vector indices should sit inside a List in the Subscript.Additional options may be supplied to these functions to control for example the verbosity.Thus, CNorm may for example be called using CDot[{C1,. .., Cn}, Verbose → True], to get progress information or simply as CDot[{C1,. . . ,Cn}] to use the default options.Calculates the scalar product between two color tensors C1 and C2.Before using CDot the color tensors C1 {i1 ,i2 ,...,ik } and C2 {i1 ,i2 ,...,ik } should have been defined.For simplifications (in particular of roots) it is by default assumed that Nc is at least 3.However, this may be manually changed by setting the option NcMin, for example CDot[C1, C2, NcMin → 100].CDot[{C1,. . ., Cn}, Options] Options with default value: NcMin → 3 Verbose → False Similar to CDot[C1, C2] but calculates the scalar product of a set of color tensors and themselves and returns the result as a List where the ith element is CDot[Ci, Ci].In addition to assuming a minimal value of Nc for simplifications, it is also possible to get progress information by setting the option Verbose to True.CDotMatrix[{C1,. . ., Cn}, Options] Options with default value: NcMin → 3 Verbose → False Calculates the scalar product matrix given a List of color vectors, i.e., element ij is CDot[Ci, Cj].If the list contains a basis, this function thus returns the scalar product matrix.Calculates the norm of a color tensors using CDot.For simplifying the result, it may be useful to set the option NcMin to a large value.CNorm[{C1,. . .,Cn}, Options] Options with default value: NcMin → 3 Verbose → False Similar to CNorm above, but calculates the norms of a List of color vectors {C1,. . .,Cn}, and returns a List containing the norms.Progress information is available by supplying the option Verbose→ True.CGamma[{C1,. . .,Cn}, p1, p2, Options] Options with default value: NcMin → 3 Verbose → True MakeChecks → True BasisType → GeneralBasis