#  HepLib: A C++ Library for High Energy Physics

Feng Feng<EMAIL_ADDRESS>Yi-Fan Xie  Qiu-Chen Zhou  Shan-Rong Tang
China University of Mining and Technology, Beijing 100083, China

###### Abstract

HepLib is a C++ Library for computations in High Energy Physics, it works on
top of GiNaC, a well-established C++ library used to perform symbolic
computations. HepLib combines serval well-known packages to get high
efficiency, including qgraf to generate Feynman aptitudes, FORM to perform
Dirac/Color matrix related computations, and FIRE or KIRA for integration-by-
parts (IBP) reduction. Another core feature of HepLib lies in the numerical
evaluation of master integrals using sector decomposition, which is a general
method widely used in high-order numerical computation and has been
implemented in many public packages in many different languages, and we
present another implementation in the language of C++ with many new features.
We use GiNaC to handle the symbolic operations, and export the corresponding
integrand into an optimized C++ code, that will be compiled internally and
linked dynamically, a customizable numerical integrator is selected to perform
the numerical integration, while the integrand can be evaluated in different
float precisions, including the arbitrary precision supported by MPFR.

###### keywords:

Loop Corrections , Feynman Integral , Sector Decomposition

††journal: Computer Physics Communications

Program Summary

Program title: HepLib

GitHub pages: https://heplib.github.io/

Licensing provisions: GPLv3

Program language: C++

Operating System: Linux, Mac OS

External libraries: GiNaC, qgraf, MPFR, QHull, MinUit2, CUBA

External programs: FORM, Fermat, FIRE, KIRA

Nature of the problem: There are many independent programs or packages in high
energy physics, which are written or developed in different programming
languages, here we are trying to provide an integrated interface in the C++
language, for generating Feynman diagrams/amplitudes, performing algebraic
simplifications on Dirac/Color matrix objects, and reducing scalar integrals
to master integrals. The evaluation of master integrals is in the core part
during the high-order calculations, one needs to resolve ultraviolet/infrared
divergence within dimensional regularization and the singularities inside the
integration domain while performing multiple dimensional numerical
integrations.

Method of solution: HepLib uses qgraf to generate the Feynman
diagrams/aptitudes, FORM to perform Dirac/Color matrix related computations,
FIRE or KIRA for the IBP reduction. A C++ implementation of sector
decomposition method is used to extract both ultraviolet and infrared
singularity of the integral, the contour deformation is adopted to avoid the
singularity inside the integration domain, a parallelized numerical
integration is chosen to achieve high efficiency or performance.

Restrictions: Depending on the complexity of the problem itself.

##  1 Introduction

The higher-order computation in perturbative quantum field theory or Standard
Model is indispensable, as the precision or accuracy increases at the Large
Hadron Collider, many scattering processes have been calculated up to next-to-
next-to-leading order to meet the precision demand, some processes to even
higher order [1]. The traditional basic procedures in the higher-order
computations in high energy physics involve generating Feynman
diagrams/amplitudes, performing algebraic simplifications or traces on
Dirac/Color matrix objects, applying partial fraction decomposition and
integration-by-parts (IBP) reduction on the scalar integrals to result in the
master integrals, finally one needs a way to evaluate those master integrals
analytically or numerically. There are many public packages available to
complete those procedures, here we just mention some (not complete) commonly
used packages in each step. FeynArts [2, 3] is a Mathematica package used to
generate Feynman diagrams/amplitudes, qgraf [4, 5] is another widely used
program in amplitudes generation, which is written in Fortran language to get
much higher efficiency. FeynCalc [6, 7, 8] developed in Mathematica language
is commonly used to perform Dirac/Color matrix trace or simplifications, FORM
[9, 10] is another highly efficient program developed in C language to handle
the operations involving Dirac-$\gamma$ or indexed objects, FeynCalc/FormLink
[11, 12] is written as an interface between FeynCalc and FORM. The Mathematica
package Apart [13, 14] can be used to perform partial fractions decomposition.
One can use AIR [15, 16] written in Maple language, FIRE [18, 19, 20, 21]
written in Mathematica and C++ language, LiteRed [22, 23, 24] a powerful
Mathematica package, Reduze [25, 26, 27] or KIRA [28, 29, 30], both in C++
language, to perform the IBP reductions on scalar integrals to get master
integrals. As for the evaluation of master integrals, Mellin-Barnes (MB) [31,
32] or sector decomposition (SD) [34, 35] methods are widely used in numerical
calculations. There are many useful subroutines on the MBTools website [33] to
facilitate the usage of MB method, one can also consult the AMBRE [36, 37, 38]
packages for the automatic Mellin-Barnes representation. FIESTA [39, 40, 41,
42, 43] and SecDec [44, 45, 46, 47] are two widely used packages developed in
Mathematica language, both provide implementations of SD method, pySecDec [48,
47] is a new version of the program SecDec, it uses Python as a framework to
organize different parts of the evaluation while the main work is done using
FORM and CUBA [49, 50, 51]/QMC [52], sector_decomposition [53, 54] is an
implementation in C++ language with the help of GiNaC [55, 56], CSector [57,
58] is a Mathematica interface to sector_decomposition. The differential
equation (DE) [59, 60] method becomes a powerful tool to evaluate the master
integrals, especially when the equations can be transformed into a canonical
form or $\epsilon$-form [61, 62], there are also many public packages to find
the possible canonical form automatically: Fuchsia [63, 64], epsilon [65, 66],
CANONICA [67, 68] and so on. The numerical or semi-numerical DE method with
powers series expansions can also be found in [69, 70], and DESS [69, 71] is a
Mathematica implementation of such a kind of method near singular points. Last
but not least, there are also many useful self-integrated programs for the
perturbative calculation of cross sections mainly up to next-to-leading order:
FormCalc [72, 73], CompHEP [74, 75, 76], FDC [77], CalcHEP [78, 79], GoSam
[80, 81, 82], DIANA [83], MadGraph5_aMC@NLO [84], Herwig++ [85], SHERPA [86],
WHIZARD [87] and many others, we refer the reader to hepforge projects [88]
for more related programs.

Mathematica is a widely used language for generic symbolic or algebraic
computations, Python and other languages can also be used to achieve the same
purpose, while for the numerical evaluation, Fortran or C++/C is still most
favored. Usually, one needs to jump from one language to another during the
high order calculations through an interface, here we want to present the
library HepLib with those interfaces implemented internally, so one does not
need to know the internal conversions between different languages. HepLib
tries to provide a general C++ framework or interface involving both symbolic
computation and numerical evaluation in the high order calculations. The core
symbolic operations in C++ language are built on GiNaC (we refer the reader to
the official tutorial [89] for a basic introduction to GiNaC). HepLib makes
small extensions to GiNaC by adding some basic classes or objects which are
familiar (especially to the users of FeynCalc) in the higher-order
computations, it provides a few interfaces to other well-known sophisticated
programs including Fermat [90] to do fast simplifications on multivariate
rational polynomials, qgraf to generate Feynman diagrams and amplitudes, FORM
to handle Dirac/Color object related operations, FIRE or KIRA to perform the
IBP reductions. As for numerical evaluation, HepLib provides another native
GiNaC implementation of the sector decomposition method with many new
features. We will try to add other methods including DE or MB in the near
future versions.

##  2 Description of HepLib

###  2.1 Extension to GiNaC

HepLib works on top of GiNaC which is a well-established C++ library for
symbolic calculations, in this section we want to give a short introduction to
GiNaC and make a few extensions to it, while for a general description of
GiNaC, we refer the reader to its official tutorial [89].

GiNaC::ex is one of the most important and commonly used class in GiNaC, it
represents a general symbolic expression which contains variables, numbers,
functions and so on. The GiNaC expressions work like handles, i.e., a pointer
to the internal actual object, but one needs not to care about the memory
allocation and deallocation, since GiNaC uses an intrusive reference-counting
mechanism to handle the garbage collection, i.e., when an instance of
GiNaC::ex runs out of its scope, the destructor will check whether any other
expression points to such internal actual object too, and deletes the object
from memory if that turns out not to be the case. Due to such a reference-
counting mechanism, GiNaC expressions are not thread-safe, and it would be
problematic if one tries to use GiNaC from different threads simultaneously.

A GiNaC::ex object can contain complex variables instanced from GiNaC::symbol
class, numbers instanced from GiNaC::numeric class, functions instanced from
GiNaC::function class and so on. Expressions can also be put together to form
a new expression or passed as arguments to a function, and so on. Here is a
little collection of GiNaC expressions (we omit the namespace GiNaC for each
class for simplicity):

⬇

symbol x("x"), y("y"), z("z"); //complex variables

numeric r(2,3); //exact fraction 2/3

ex e1 = r*x+2*y+pow(y,10); //polynomial in x and y

ex e2 = (x+1)/(x-1); //rational expression

ex e3 = sin(x+2*y)+3*z+41; //containing a function

ex e4 = e3+e2/exp(e1); //expression formed from other expressions

where the pow(y,10) refers to $y^{10}$, and the function sin and exp are
objects of GiNaC::function, they are actually pseudo-functions, i.e., their
evaluations may be halted if the arguments are not numerical. There are many
other mathematical functions defined in GiNaC, we refer the reader to the
tutorial [89] for a complete list.

It should be noted that GiNaC::symbol refers to a complex variable, and the
string name in its constructor symbol::symbol(string name) is not the
identification to the object, so if we write an expression as ex e =
symbol("x")-symbol("x"), the two symbol("x") are distinct in this expression,
the expression e is not zero, and one will get something like x-x in the
output. Furthermore, most of the variables encountered in a physical process
are real, so we introduce a class HepLib::Symbol, its usage is almost the same
as GiNaC::symbol, except that it represents a real variable, and the name in
its constructor Symbol::Symbol(string name) is indeed the identification to
the object, i.e., the expression Symbol("x")-Symbol("x") will produce the
result 0 as expected. Finally, we use the following segment of the code
0.cpp111The file xxx.cpp in the current version can be found at codes/xxx.cpp
in HepLib archive. to show the difference of GiNaC::symbol and HepLib::Symbol
in calling the GiNaC::conjugate function which returns the complex conjugation
of an expression (we also omit the namespace HepLib for the corresponding
class for simplicity).

⬇

symbol x("x"); //a complex variable

Symbol y("y"); //a real variable

ex e1 = conjugate(x); //complex conjugation of x \--> conjugate(x)

ex e2 = conjugate(y); //complex conjugation of y \--> y

where one can see conjugate(y) will be evaluated to y, while conjugate(x) will
still return as conjugate(x).

HepLib::Symbol class is very convenient when one converts a string text into a
GiNaC expression, and such a conversion usually happens in the data exchange
between HepLib and other packages (e.g., Fermat and FORM). The conversion from
a string text to a GiNaC or HepLib expression is encoded in the class
HepLib::Parser which has been introduced as follows:

⬇

class Parser {

public:

Parser(); Parser(symtab st); symtab STable;

ex Read(string instr,bool s2S=true); ex ReadFile(string filename,bool
s2S=true);

};

where symtab is defined in GiNaC itself as typedef std::map<string,ex> symtab,
so the member STable is used to map a variable appearing in the string text to
the corresponding ex object, and an unknown variable (one does not provide the
mapping explicitly in STable) will be converted to a symbol or Symbol object.
The method Read or ReadFile will parse the expression from a string text or a
file content respectively, and when s2S is true, an unknown variable is
converted to a Symbol object, which is also the default case. Here is a snip
of the code 1.cpp to exemplify the usage of Parser,

⬇

Parser parser;

string expr_str = "WF(1)+x(1)^2+sin(5)+pow(a,n)";

ex e1 = parser.Read(expr_str);

Symbol a("a"), n("n");

ex e2 = WF(1)+pow(x(1),2)+sin(5)+pow(a,n);

cout << e1-e2 << endl; // \--> 0

e1 is now a GiNaC expression x(1)^2+WF(1)+a^n+sin(5), the unknown variables a
and n are parsed to Symbol("a") and Symbol("n") respectively, and the
functions sin, WF and x are also parsed correctly, since sin is defined in
GiNaC, WF and x are pseudo-functions introduced in HepLib, the both functions
can take any ex object as the argument and will do nothing on it. We also
construct the expression e2 explicitly to see e1 and e2 are indeed identical.

GiNaC::lst is another commonly used class that serves for holding a list of
arbitrary expressions. One can use the nops method to get the number of
expressions in a list, the op method to access individual elements, and let_op
method to change an element in the list. The op and nops methods are actually
available to any general GiNaC expression too, and one can use them to access
the internal subexpressions. Here is a snip of the code 2.cpp to demonstrate
the basic usage of those methods:

⬇

Symbol x("x"), y("y"), z("z"), n("n");

lst l1 = lst{ x, y, x+z };

ex e1 = pow(sin(x),n);

int tot = l1.nops(); // tot is 3

ex item1 = l1.op(0); // 1st item in l1, \--> x

l1.let_op(2) = e1; // now l1 is {x, y, sin(x)^n}

tot = e1.nops(); // now tot is 2

ex item2 = e1.op(1); // 2nd item in e1, \--> n

Note that one can not modify an item in a general expression by using the
let_op method, to achieve this purpose one can use one of the following two
methods, one is the expression substitution to make some simple replacements
with the help of GiNaC function GiNaC::subs, the other is to use the class
HepLib::MapFunction to perform more complicated operations on parts of the
expression. For the first method, one usually use the class GiNaC::wildcard to
do the pattern matching, a wildcard is a special kind of GiNaC object which
represents an arbitrary expression, and to match different expressions one can
use different label in its constructor ex wildcard::wildcard(unsigned
label=0), one can also use w and w0 in HepLib for the abbreviation of
wildcard(0) object, w1 for wildcard(1), and so are the w2 up to w9. While for
the second method, the HepLib::MapFunction class is introduced as follows:

⬇

class MapFunction {

public:

ex operator()(const ex &e);

MapFunction(std::function<ex(const ex &, MapFunction &)>);

};

We exemplify both methods by making two substitutions on a simple polynomial
$x^{4}+x^{3}+x^{2}+x$, one substitution is to replace any term which
explicitly matches the pattern $x^{n}$ by $y^{n+2}$, i.e., $x^{n}\to y^{2+n}$,
the other one is to do the replacement only for the terms with an even
exponent, i.e., $x^{2n}\to y^{2+2n}$, using the following snip of the code
3.cpp:

⬇

Symbol x("x"), y("y");

//HepLib::str2ex: convert string to expression, using Parser internally

ex e0 = str2ex("x^4+x^3+x^2+x");

ex e1 = e0.subs( lst{ pow(x,w)==pow(y,w+2) } ); //e1 is $y^{6}+y^{5}+y^{4}+x$

MapFunction map([&](const ex &e, MapFunction self)->ex{

if(e.match(pow(x,w)) && e.op(1).info(info_flags::even)) return
pow(y,e.op(1)+2);

else return e.map(self);

});

ex e2 = map(expr); //e2 is $y^{6}+x^{3}+y^{4}+x$

where we have used the C++ lambda expression as the std::function object, and
note that the statement e.map(self) will apply the specified function on all
subexpressions (in the sense of op()), non-recursively [89].

Last but not least, one usually gets a long list of expressions (e.g.,
hundreds or thousands of Feynman amplitudes) in high order computations, and
performs same or similar operators on each item, since the operations on each
item are independent, one can do those operations simultaneously, i.e., a
single operation executes simultaneously on multiple elements of data (SIMD).
As we have mentioned before, one can not use GiNaC from different threads
simultaneously due to the internal reference-counting mechanism, so one can
not run the operations from different threads in one process simultaneously,
and we need to turn to a process-based parallel method. We extend GiNaC to a
parallel version by using the fork() system call, fork() will create a new
process (child process) by duplicating the calling process (parent process),
so the input data to the child process is an exact duplicate of the parent
process, while we still need a way to get data back from each child process to
the parent process.

GiNaC allows creating archives of expressions which can be stored to or
retrieved from files, we wrap these operations into HepLib functions garWrite
(to write an expression into an archive file) and garRead (to read an
expression from an archive file), for example, garWrite("e.gar", e) will
export the expression e into the archive file e.gar, and the file e.gar
contains all information which is needed to reconstruct the expressions e, and
garRead("e.gar") will return the expression identical to the original e. So to
get data back from a child process to the parent process, the child process
exports the data to an intermediate archive file by using garWrite, and the
parent process can get the data by using garRead. All those steps are encoded
in the HepLib function GiNaC_Parallel which is introduced as follows:

⬇

vector<ex> GiNaC_Parallel(int ntotal, std::function<ex(int)> f,

const string & key="", bool rm=true, const string &pre=" ");

where ntotal is the number of items, f is a function taking an int as input
and ex in return, the function f will be called on each i (0$\leq$i$<$ntotal)
in parallel, key is used to name the directory where the intermediate files
store (the process id will be used if key is empty, which is also the default
value), rm=true is to remove the intermediate files at the end of the function
call, and the last parameter pre refers to some indent space appended to the
output message (when the global variable Verbose is large than 1).
GiNaC_Parallel is equivalent to a parallel version of the following for
statement:

⬇

vector<ex> ret;

for(int i=0; i<ntotal; i++) ret.push_back(f(i));

return ret;

note that the number of parallel processes can be changed by setting the
variable GiNaC_Parallel_Process, when GiNaC_Parallel_Process is negative, the
number of parallel processes will be set to the number of cpu cores minus 1,
which is also the default case, and when GiNaC_Parallel_Process is 0,
GiNaC_Parallel is equivalent to the for statement above. Finally, we end this
section by a snip of the code 4.cpp to illustrate the basic usage of
GiNaC_Parallel:

⬇

Symbol x("x"), y("y");

vector<ex> in_vec;

int total = 1000;

for(int i=0; i<total; i++) in_vec.push_back(sin(exp(x+y*i)));

Verbose = 100;

vector<ex> ret = GiNaC_Parallel(total, [&](int idx)->ex {

ex data = in_vec[idx];

data = series_ex(data,x,5); // Taylor expansion around x=0

return data;

});

where we have filled the in_vec with some simple expressions, and taken the
Taylor expansion on each item parallelly by using HepLib::ex mma_series(const
ex & e, const symbol &s, int n) which returns the Taylor or Laurent series
expansion of e up to the n-th order with respect to s at the point s=0.

Before ending this section, it should be noted that the terms order in the
output may be different in each run, as stated in GiNaC official tutorial
[89], for example, the code cout << x-y; can produce the output x-y or -y+x,
and to get a fixed or term-ordered output, one needs to implement one’s own
output format, and the class HepLib::HepFormat is introduced for such a
purpose. The global object hout is instanced from HepFormat class, its usage
is very similar to the standard C++ object cout, and one can use the code hout
<< x-y; to get the fixed or term-ordered output.

###  2.2 Basic objects in HepLib

Object |  Example |  Description
---|---|---
Symbol |  Symbol("s") |  a real variable s.
Index |  Index("mu",Type::VD) |  a Lorentz index mu with dimension D.
|  Index("a",Type::CA) |  a color index a with dimension NA.
|  Index("i",Type::CF) |  a color index i with dimension NF.
Vector |  Vector("p") |  a vector/momentum p.
Pair |  Pair(mu,nu) |  a Kronecker delta $\delta_{\mu\nu}$ with Index mu and nu.
|  Pair(p,mu) |  a Vector p with Lorentz Index mu, i.e., ${p}^{\mu}$.
|  Pair(p,q) |  a scalar product $p\cdot q$ between Vector p and q.
SUNT |  SUNT(a,i,j) |  a $T$-matrix element $T^{a}_{ij}$ for SU(N) group.
|  SUNT(lst{a,b,c},i,j) |  a matrix element of a product of $T$, i.e., $(T^{a}T^{b}T^{c})_{ij}$
SUNF |  SUNF(a,b,c) |  a structure constant $f^{abc}$ of SU(N) group.
SUNF4 |  SUNF4(a,b,c,d) |  a contract of two SUNF, i.e., $f^{abe}f^{ecd}$.
Eps |  Eps(mu1,mu2,mu3,mu4) |  a Levi-Civita tensor $\varepsilon_{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$.
|  Eps(p1,p2,mu1,mu2) |  a partially contracted Levi-Civita tensor $\varepsilon_{p_{1}p_{2}\mu_{1}\mu_{2}}$.
|  Eps(p1,p2,p3,p4) |  a fully contracted Levi-Civita tensor $\varepsilon_{p_{1}p_{2}p_{3}p_{4}}$.
DGamma |  DGamma(mu,l) |  a Dirac-$\gamma$ matrix $\gamma_{\mu}$ for a fermion line labeled by l.
|  DGamma(p,l) |  a Dirac slash ${p}\\!\\!\\!/={p^{\mu}}\gamma_{\mu}$ for a fermion line labeled by l.
|  DGamma(1|5|6|7,l) |  a unit matrix, $\gamma_{5}$,$\gamma_{6}$,$\gamma_{7}$ for a fermion line labeled by l.
SP |  SP(mu,nu) |  evaluated to $\delta_{\mu\nu}$.
|  SP(p+s*q,mu) |  evaluated to ${p}^{\mu}+{s}{q}^{\mu}$.
|  SP(2*p+q,p+s*q) |  evaluated to $2{p}^{2}+(2{s}+1){p}\cdot{q}+{s}{q}^{2}$.
GAS |  GAS(mu) |  evaluated to $\gamma_{\mu}$.
|  GAS(3*p+s*q) |  evaluated to $3\,{p}\\!\\!\\!/+{s}{q}\\!\\!\\!/$.
|  GAS(1|5|6|7) |  evaluated to a unit matrix, $\gamma_{5}$,$\gamma_{6}$,$\gamma_{7}$, respectively.
LC |  LC(p,mu,p+s*q,k) |  evaluated to ${s}\varepsilon_{kpq\mu}$
TR |  TR(expr) |  a wrapper for the Dirac trace of expression expr.
TTR |  TTR(lst{a,b,c,d}) |  a wrapper for the SU(N) trace of $T^{a}T^{b}T^{c}T^{d}$.
form |  form(expr) |  evaluate the expression expr using FORM program.
Table 1: The basic objects defined in HepLib. Here, mu, nu and mu1 up to mu4
are objects of Lorentz Index class with Type::VD, a, b, c, d and e are objects
of Index with Type::CA, i and j are objects of Index with Type::CF, p, q, k
and p1 up to p4 are objects of Vector, and s is an object of Symbol.

The basic objects introduced in HepLib (under namespace HepLib) for high order
computations of high energy physics are listed in Table 1, including the
following basic classes: Index (a Lorentz index or color index of SU(N)),
Vector (a vector or momentum), Pair (a metric tensor or scalar product between
vectors), SUNT (a T matrix element of SU(N)), SUNF (a structure constant of
SU(N)), Eps (a Levi-Civita tensor) and DGamma (a Dirac-$\gamma$ matrix in a
specific fermion line). To facilitate the usage of these classes, the
functions SP/GAS/LC can be used to construct the Pair/DGamma/Eps objects
directly from the expressions involving Index or Vector objects respectively,
the linear expressions provided to the arguments of those functions will be
expanded automatically. The pseudo-functions TR and TTR are provided as
function wrappers for the Dirac trace and SU(N) trace respectively, and all
the actual evaluations (both trace and contraction) will be performed by FORM
program internally through the interface form. These objects introduced above
are similar to the ones in FORM program and FeynCalc package with the
correspondence in Table 2.

HepLib | Index | Vector | Pair | Eps | SUNT | SUNF | DGamma
---|---|---|---|---|---|---|---
FORM | Index | Vector | d_ | e_ | T | f | g_
FeynCalc | LorentzIndex | Momentum | MTD/FVD/SPD | LC | SUNT | SUNF | DiracGamma
Table 2: The correspondence between HepLib and FORM or FeynCalc.

There are a few things that need to be noted:

  * 1.

Similar to the Symbol class introduced in the previous section, two objects of
class Index, Vector, SUNT and so on are equal or identical if the inputs
provided to the constructor are equal, e.g., if we defined two Vector objects
with Vector p1("p"), p2("p");, then p1 and p2 are identical.

  * 2.

There is no need to introduce an extra operator for the product of two non-
commutative DGamma objects, the operator * between DGamma objects will be
treated as non-commutative multiplication automatically.

  * 3.

Sometimes, it may be useful to use the commutative matrix elements, instead of
the non-commutative matrix itself, for example in generating Feynman
amplitudes, the pseudo-functions Matrix(M,i,j) (the matrix element $M_{ij}$)
has been introduced for this purpose, the multiplication of Matrix objects is
commutative, and MatrixContract can contract the matrix indices and result in
the properly ordered non-commutative objects.

  * 4.

All actual operations involving Lorentz index contraction or trace of
Dirac/Color matrix will be handled by FORM program internally through the
function call form(expr). The internal details behind the form function is
almost the same as FormLink: two pipes will be created between HepLib and FORM
program for data exchange, the input expr will be translated into a FORM
script and sent to FORM for execution through one of the pipes, and the output
will be transferred back to HepLib through the other pipe and convert back to
HepLib objects at the finish of FORM execution. To avoid the expanse in the
creation and destruction of processes again and again, we keep the FORM
process alive after the first call to form during the entire program life.
Since the FORM program will be invoked internally, we need the shell command
form can be found in the environment variable PATH in order to use form
function from HepLib.

  * 5.

The default dimension for a Lorentz Index is D (represented by Symbol("D")),
and the default conversion from TR to FORM is tracen, so an Error will be
thrown if one uses $\gamma_{5}$ in TR, while one can change the default
dimension to 4 by the statement form_using_dim4=true; before calling the
function form, and the trace4 will be used for the Dirac trace in that case.
Furthermore, HepLib does not distinguish the upper and lower index objects,
when the contraction occurs between two Index objects, one index is supposed
to be an upper index and the other is supposed to be a lower index.

  * 6.

The FORM procedure SUn.prc from the color package [91, 92, 93] has been used
for the treatment of SUNT and SUNF color objects.

  * 7.

The function conjugate is introduced in GiNaC itself, this function has been
overridden for the objects introduced in Table 1, it can be used for complex
conjugation involving Dirac-$\gamma$ and SUNT objects, it is similar to the
function ComplexConjugate in FeynCalc.

  * 8.

To facilitate the conversion to FeynCalc, the FCFormat class is provided in
HepLib, the global object fcout is an instance of FCFormat, which is similar
to hout or the standard C++ cout object, one can use fcout to print an
expression expr in FeynCalc format with fcout << expr << endl;.

We end this section by a physical expression: the squared matrix element
$|{\cal M}|^{2}$ for the process
$e^{-}(p)+e^{+}(p^{\prime})\to\gamma^{*}(q)\to\mu^{+}(k)+\mu^{-}({k^{\prime}})$,
which is taken from (5.4) on page 132 of the textbook of quantum field theory
[94]:

$\frac{1}{4}\sum_{\rm spins}|{\cal M}|^{2}=\frac{e^{4}}{4q^{4}}{\rm
tr}\Big{[}(p\\!\\!\\!/^{\prime}-m_{e})\gamma^{\mu}(p\\!\\!\\!/+m_{e})\gamma^{\nu}\Big{]}\,{\rm
tr}\Big{[}(k\\!\\!\\!/+m_{\mu})\gamma_{\mu}(k\\!\\!\\!/^{\prime}-m_{\mu})\gamma_{\nu}\Big{]}$
(1)

We use the following snip of the code 5.cpp to illustrate how to input the
expression above in HepLib:

⬇

Symbol me("me"), mm("mm"), e("e");

Index mu("mu"), nu("nu");

Vector p("p"), P("P"), k("k"), K("K"), q("q");

letSP(p)=me*me; letSP(P)=me*me;

letSP(k)=mm*mm; letSP(K)=mm*mm;

#define gpm(p,m) (GAS(p)+m*GAS(1))

ex tr1 = TR( gpm(P,-me)*GAS(mu)*gpm(p,me)*GAS(nu) );

ex tr2 = TR( gpm(k,mm)*GAS(mu)*gpm(K,-mm)*GAS(nu) );

ex res = pow(e,4) / (4*pow(SP(q),2)) * tr1 * tr2;

form_using_dim4 = true; //using 4 dimension

res = form(res); //using form(res,n) with n>1 to print FORM script

res = exfactor(res);

hout << res.subs(me==0) << endl; // note hout is used

where the Vector objects p, P, k, K and q refer to the momenta $p$,
$p^{\prime}$, $k$, $k^{\prime}$ and $q$ respectively, the Index objects mu and
nu refer to the Lorentz indices $\mu$ and $\nu$ respectively, the symbols me,
mm and e refer to $m_{e}$, $m_{\mu}$ and $e$ respectively. We also use letSP
to assign the scalar product between momenta (to clear the SP assignment, one
can use clearSP). After the function call to form, the expression will be
translated to FORM script (to view the FORM script, one can replace the
function call form(res) with form(res,10)) and evaluated by FORM program, and
finally one gets the output as 8*e^4*q.q^(-2)*(k.P*p.K+mm^2*p.P+P.K*p.k),
which corresponds to the final result in (5.10) on page 135 of textbook [94].

$\frac{1}{4}\sum_{\rm spins}|{\cal M}|^{2}=\frac{8e^{4}}{q^{4}}\Big{[}(p\cdot
k)(p^{\prime}\cdot k^{\prime})+(p\cdot k^{\prime})(p^{\prime}\cdot
k)+m_{\mu}^{2}(p\cdot p^{\prime})\Big{]}$ (2)

One needs the SU(N) group factor ${\rm tr}(t^{a}t^{b})$ when one replaces
$\mu^{+}\mu^{-}$ by $q\bar{q}$, i.e., the process $e^{+}e^{-}\to q\bar{q}$
(see (17.6) on page 549 of [94]), and one can use the pseudo-function TTR for
the trace of a product of SU(N) T-matrix objects:

⬇

Index a("a",Index::Type::CA);

ex tr = TTR(lst{a,a});

tr = form(tr);

Another way to achieve the same purpose is to use the SU(N) T-matrix element,
i.e., ${\rm tr}(t^{a}t^{a})=t^{a}_{ij}t^{a}_{ji}$, as follows:

⬇

Index a("a",Index::Type::CA), i("i",Index::Type::CF), j("j",Index::Type::CF);

ex tr = SUNT(a,i,j) * SUNT(a,j,i); // or SUNT(lst{a,a},i,i)

tr = form(tr);

###  2.3 Generating Feynman diagrams/amplitudes

The namespace HepLib::QGRAF is introduced for the generation of Feynman
diagrams/amplitudes for the underlying scattering process. The main interface
to qgraf program is encoded in the class Process introduced as follows:

⬇

class Process {

public:

string Model; string In; string Out; string LoopPrefix = "q";

int Loops; string Options; vector<string> Others;

lst Amplitudes(symtab st, bool debug=false);

}

Most of the members are the same as those in the input file qgraf.dat of qgraf
program, here we exemplify the basic usage by the process
$e^{-}(p)+e^{+}(P)\to\gamma^{*}(q)\to\mu^{+}(k)+\mu^{-}(K)$ introduced in the
previous section, using the following snip of the code 6.cpp:

⬇

Vector p("p"), P("P"), k("k"), K("K");

Process proc;

proc.Model = R"EOF(

[e, ebar, -]

[mu, mubar, -]

[A, A, +]

[ebar, e, A]

[mubar, mu, A]

)EOF";

proc.In = "e[p],ebar[P]";

proc.Out = "mubar[k],mu[K]";

proc.Options = "onshell";

proc.Loops = 0;

symtab st;

st["p"] = p; st["P"] = P; st["k"] = k; st["K"] = K;

auto amps = proc.Amplitudes(st);

where we define the QED model involving three particles e ($e^{-}$) and mu
($\mu^{-}$) and A ($\gamma$), and their anti-particles ebar ($e^{+}$) and
mubar ($\mu^{+}$) as well, and define the In and Out states for our process
respectively. The function Amplitudes will invoke qgraf program internally,
and use the following style222One can modify the style through the static
member Process::Style, but if one introduces new functions other than
Propagator, Vertex and so on, one needs to define these new functions as GiNaC
pseudo-functions. to control and generate the output,

⬇

<diagram>

((<sign><symmetry_factor>)*

<in_loop>InField(<field>,<field_index>,<momentum>)*

<end>

<back><out_loop>OutField(<field>,<field_index>,<momentum>)*

<end>

<back><propagator_loop>Propagator(

<back>Field(<field>,<field_index>),

<back>Field(<dual-field>,<dual-field_index>),

<back><momentum>)*

<end>

<back><vertex_loop>Vertex(

<back><ray_loop>Field(<field>,<field_index>,<momentum>),

<back><end><back>)*

<end><back><back>)

,

<epilogue>

<exit>

and finally parse the output to HepLib objects wrapped in the pseudo-functions
(Field, Propagator, Vertex, etc.). Note that one needs to provide a symtab
object as the argument in calling Amplitudes, since we want the objects p, k
and so on to be parsed as Vector instead of Symbol. All amplitudes are
returned as a lst (amps above), and the typical amps looks as follows:

⬇

{
-InField(e,-1,p)*InField(ebar,-3,P)*OutField(mu,-4,K)*OutField(mubar,-2,k)*Propagator(Field(A,1),Field(A,2),-P-p)*Vertex(Field(ebar,-3,P),Field(e,-1,p),Field(A,1,-P-p))*Vertex(Field(mubar,-4,-K),Field(mu,-2,-k),Field(A,2,P+p))
}

where one can see the field_index appearing in the second argument of function
Field, and note that qgraf program labels the incoming particles as -1, -3,
$\cdots$, and outgoing particles as -2, -4, $\cdots$, any field_index
associated with an internal line is positive.

Currently, HepLib does not provide other kinds of diagram filtering options
than what is available in qgraf program itself. One can use the function
TopoLines to generate the topological lines from an amplitude, and draw the
corresponding Feynman diagram with the help of the LaTeX package Tikz-Feynman
[95, 96], all the internals are encoded in the function call DrawPDF, as
exemplified by the snip of the code 6.cpp:  ⬇ InOutTeX[-1]="$e^-(p)$";
InOutTeX[-3]="$e^+(P)$"; InOutTeX[-2]="$\\\mu^+(k)$";
InOutTeX[-4]="$\\\mu^-(K)$"; LineTeX[e] = "fermion, edge label=$e$";
LineTeX[ebar] = "anti fermion, edge label=$e$"; LineTeX[mu] = "fermion, edge
label=$\\\mu$"; LineTeX[mubar] = "anti fermion, edge label=$\\\mu$";
LineTeX[A] = "photon, edge label=$\\\gamma$"; DrawPDF(amp, "amps.pdf"); note
that, since the Tikz-Feynman package is used, it is required that the shell
command lualatex and pdflatex should be located in the environment variable
PATH, and Tikz-Feynman package should be installed as well.

As for the returned amplitudes from Amplitudes, one still needs to provide the
corresponding Feynman rules explicitly, as illustrated in the following snip
of the code 6.cpp:

⬇

auto amps_feyn_rule = MapFunction([&](const ex &e, MapFunction &self)->ex {

if(isFunction(e,"OutField") || isFunction(e,"InField")) return 1;

else if(isFunction(e, "Propagator")) {

auto fi1 = e.op(0).op(1);

auto fi2 = e.op(1).op(1);

auto mom = e.op(2);

if(e.op(0).op(0)==A) {

return (-I) * SP(LI(fi1),LI(fi2))/SP(mom); //Feynman Gauge

} else if(e.op(0).op(0)==ebar) {

return I * Matrix(GAS(mom)+GAS(1)*me,DI(fi1),DI(fi2))/(SP(mom)-me*me);

} else if(e.op(0).op(0)==mubar) {

return I * Matrix(GAS(mom)+GAS(1)*mm,DI(fi1),DI(fi2))/(SP(mom)-mm*mm);

}

} else if(isFunction(e, "Vertex")) {

auto fi1 = e.op(0).op(1);

auto fi2 = e.op(1).op(1);

auto fi3 = e.op(2).op(1);

if(e.op(0).op(0)==ebar) {

return I*Symbol("e")*Matrix(GAS(LI(fi3)),DI(fi1),DI(fi2));

} else if(e.op(0).op(0)==mubar) {

return I*Symbol("e")*Matrix(GAS(LI(fi3)),DI(fi1),DI(fi2));

}

}

return e.map(self);

})(amps);

where the uniqueness of field_index from qgraf program is used to label the
Lorentz index, Dirac-$\gamma$ or Color index (through the functions
LI/DI/CI333Each function accepts an integer as input, e.g., LI(1) refers to
Index("li1"), LI(-1) refers to Index("lim1"), more details can be found in
HepLib document [97].), and the final expression of amps_feyn_rule looks like:

⬇

{-I*e^2*Matrix(($\gamma$.li1),dim3,dim1)*Matrix(($\gamma$.li2),dim4,dim2)*li2.li1*(P.P+p.p+2*p.P)^(-1)}

and to get the squared matrix element in (1), one still needs to multiply it
by its conjugation and the spin polarization summation as follows:

⬇

ex ampL = amps_feyn_rule.op(0);

ex ampR = IndexL2R(conjugate(ampL));

#define SS1(p,m,i) Matrix(GAS(p)+m*GAS(1),DI(i),RDI(i))

#define SS2(p,m,i) Matrix(GAS(p)+m*GAS(1),RDI(i),DI(i))

ex M2 = ampL*ampR*SS1(p,me,-1)*SS2(P,-me,-3)*SS1(k,-mm,-2)*SS2(K,mm,-4);

M2 = MatrixContract(M2);

note that the function IndexL2R (usually used in the case that the dummy index
needs to be remapped) will change, and only change, the objects generated from
LI/DI/CI functions to their conjugated parts, e.g., LI(1) (equivalent to
Index("li1")) will be mapped to RLI(1) (equivalent to Index("rli1")), and so
on.

Last but not least in this section, to illustrate the operations on more
complicated indexed objects, or IBP reduction for loop corrections in the
following sections, let’s consider another physical process
$\psi\to\gamma^{*}\to e^{+}e^{-}$ [98, 99]. The form factor for this process
can be extracted from the amplitude of the corresponding partonic process
$Q(p)+\bar{Q}(p)\to\gamma^{*}(2p)$ ($Q$ refers to the heavy quark), and one
can use the following Process object to get the related amplitudes up to next-
to-next-leading order:444All code can be found at codes/nnlo.cpp in HepLib
archive.

⬇

Vector p("p"), q1("q1"), q2("q2");

ex p1 = p, p2 = p;

QGRAF::Process proc;

proc.Model = model;

proc.In = "Q[p1],Qbar[p2]";

proc.Out = "A[pA]";

proc.Options = "onshell";

proc.LoopPrefix = "q";

symtab st; st["pA"] = p1+p2;

st["p1"] = p1; st["p2"] = p2; st["q1"] = q1; st["q2"] = q2;

proc.Loops = 0; auto amp0 = proc.Amplitudes(st);

proc.Loops = 1; auto amp1 = proc.Amplitudes(st);

proc.Loops = 2; auto amp2 = proc.Amplitudes(st);

and the underlying model is defined as follows:

⬇

std::string model = R"EOF(

[q, qbar, -]

[Q, Qbar, -]

[gh, ghbar, -]

[g, g, +, notadpole]

[A, A, +, external]

[qbar, q, g; QCD=’+1’]

[Qbar, Q, g; QCD=’+1’]

[g, g, g, g; QCD=’+2’]

[g, g, g; QCD=’+1’]

[ghbar, gh, g; QCD=’+1’]

[qbar, q, A; QCD=’+0’]

[Qbar, Q, A; QCD=’+0’]

)EOF";

where q refers to the light quark, Q to heavy quark, g to gluon, and gh to the
ghost of gluon, furthermore we define the photon field as external to avoid
QED loop corrections. One can use the following MapFunction object to assign
the corresponding Feynman rules to the amplitudes amp0 (leading order), amp1
(next-to-leading order) and amp2 (next-to-next-to-leading order):

⬇

auto map = MapFunction([&](const ex &e, MapFunction &self)->ex{

if(isFunction(e,"OutField")||isFunction(e,"InField")) return 1;

else if(isFunction(e,"Propagator")){

if(e.op(0).op(0)==q) return QuarkPropagator(e,0);

else if(e.op(0).op(0)==Q) return QuarkPropagator(e,m);

else if(e.op(0).op(0)==g) return GluonPropagator(e);

else if(e.op(0).op(0)==gh) return GhostPropagator(e);

}else if(isFunction(e,"Vertex")){

if(e.nops()==3&&e.op(0).op(0)==ghbar && e.op(1).op(0)==gh){

return gh2gVertex(e); // ghbar-gh-g

}else if(e.nops()==3&&e.op(0).op(0)==g&&e.op(1).op(0)==g&&e.op(2).op(0)==g){

return g3Vertex(e); // g^3

}else if(e.nops()==4&&e.op(0).op(0)==g&&e.op(1).op(0)==g&&e.op(2).op(0)==g){

return g4Vertex(e); // g^4

}else if(e.nops()==3 &&((e.op(0).op(0)==qbar)||(e.op(0).op(0)==Qbar))){

if(e.op(2).op(0)==g) return q2gVertex(e); // qbar-q-g or Qbar-Q-g

else{ // qbar-q-A or Qbar-Q-A

ex fi1 = e.op(0).op(1), fi2 = e.op(1).op(1), fi3 = e.op(2).op(1);

return Matrix(GAS(LI(fi3)),DI(fi1),DI(fi2)) * SP(TI(fi1),TI(fi2));

}

}

}

return e.map(self);

});

amp0 = ex_to<lst>(map(amp0)); amp1 = ex_to<lst>(map(amp1)); amp2 =
ex_to<lst>(map(amp2));

where the functions xxxPropagator/xxxVertex under namespace HepLib::QGRAF
correspond to the Feynman rules for the propagator/vertex in QCD, and the
first argument for the xxxPropagator/xxxVertex are always Propagator/Vertex
objects directly generated from Amplitudes, we refer the reader to HepLib
document [97] for detailed information.

Finally, to save the amplitudes for later processing, one can export the whole
list to a gar file with the function garWrite, and use garRead to get the data
back as follows:

⬇

ex amps = lst { amp0, amp1, amp2 };

garWrite("amps.gar", amps);

amps = garRead("amps.gar"); // read back to amps

###  2.4 More operations on indexed objects

Generally, one gets the amplitude with non-commutative objects wrapped in the
pseudo-function Matrix with matrix index, and one can use MatrixContract to
perform the matrix multiplication, while for the Lorentz or Color indexed
objects, the related operations are handled by FORM program internally. In
this section, we take amp1 (the amplitude for $Q\bar{Q}\to\gamma^{*}$ at next-
to-leading order) in last section as an example to illustrate more operations
on the indexed objects. Note that all the C++ outputs in this section are
generated by the code codes/nlo.cpp in HepLib archive.

One can use garRead to get amp1 back from the gar file as follows:

⬇

(-I)*gs^2*Matrix(($\gamma$.li1),di4,dim1)*Matrix(($\gamma$.li2),dim3,di5)*Matrix(($\gamma$.mu),di6,di3)*Matrix(($\gamma$.p)+m*$\mathbbm{1}$+(-1)*($\gamma$.q1),di3,di4)*Matrix(m*$\mathbbm{1}$+(-1)*($\gamma$.p)+(-1)*($\gamma$.q1),di5,di6)*ci2.ci1*li2.li1*ti4.ti3*ti6.ti3*ti6.ti5*q1.q1^(-1)*T(ci1,ti4,tim1)*T(ci2,tim3,ti5)*(p.p+q1.q1+(-1)*m^2+(-2)*q1.p)^(-1)*(p.p+q1.q1+(-1)*m^2+2*q1.p)^(-1)

where q1 is the loop momentum, and mu is the Lorentz index associated with the
virtual photon $\gamma^{*}$, and one can get a more compact expression as
follows, with the help of function call MatrixContract:

⬇

(-I)*gs^2*Matrix(($\gamma$.li2)*(m*$\mathbbm{1}$+(-1)*($\gamma$.p)+(-1)*($\gamma$.q1))*($\gamma$.lim2)*(($\gamma$.p)+m*$\mathbbm{1}$+(-1)*($\gamma$.q1))*($\gamma$.li1),dim3,dim1)*ci2.ci1*li2.li1*ti4.ti3*ti6.ti3*ti6.ti5*q1.q1^(-1)*T(ci1,ti4,tim1)*T(ci2,tim3,ti5)*(p.p+q1.q1+(-1)*m^2+(-2)*q1.p)^(-1)*(p.p+q1.q1+(-1)*m^2+2*q1.p)^(-1)

as one can see there is only one Matrix object left with external field index
dim3 and dim1, which needs to be contracted with the Dirac spinors
$\bar{v}_{\tt dim3}$ (from heavy anti-quark) and $u_{\tt dim1}$ (from heavy
quark) respectively. We use the the spin projector technique [100, 101] to
make the following replacement: u_dim1 ¯v_dim3 →[
-142m2(p/+m)ε/(p/-m) ]_dim1 dim3  →-142m2 ε^ν
Matrix((p/+m)γ_ν(p/-m),dim1,dim3) note that one can pull/drop out the
overall polarization vector $\varepsilon^{\nu}$ in the actual computation, and
get the result involving the Dirac trace TR after using MatrixContract again:

⬇

(-I)*NF^(-1/2)*gs^2*TR((1/32)*sqrt(2)*m^(-3)*(($\gamma$.li2)*(m*$\mathbbm{1}$+(-1)*($\gamma$.p)+(-1)*($\gamma$.q1))*($\gamma$.mu)*(($\gamma$.p)+m*$\mathbbm{1}$+(-1)*($\gamma$.q1))*($\gamma$.li1)*(($\gamma$.p)+m*$\mathbbm{1}$)*(2*m*$\mathbbm{1}$+2*($\gamma$.p))*($\gamma$.nu)*(($\gamma$.p)+(-1)*m*$\mathbbm{1}$)))*ci2.ci1*li2.li1*ti4.ti3*ti6.ti3*ti6.ti5*tim3.tim1*q1.q1^(-1)*T(ci1,ti4,tim1)*T(ci2,tim3,ti5)*(p.p+q1.q1+(-1)*m^2+(-2)*q1.p)^(-1)*(p.p+q1.q1+(-1)*m^2+2*q1.p)^(-1)

one can send the expression to FORM program through form function, and get the
evaluated result as follows:

⬇

(8/3*I)*sqrt(2)*sqrt(3)*gs^2*q1.q1^(-1)*(q1.q1+2*q1.p)^(-1)*(2*q1.p+(-1)*q1.q1)^(-1)*(2*m^2*nu.mu+2*q1.mu*q1.nu+(-2)*ep*q1.mu*q1.nu+(-1)*nu.mu*q1.q1+ep*nu.mu*q1.q1)

furthermore, one can write down the expression above in a more readable form:

$\sqrt{\frac{2}{3}}8ig_{s}^{2}\frac{1}{q_{1}^{2}}\frac{1}{2p\cdot
q_{1}-q_{1}^{2}}\frac{1}{2p\cdot
q_{1}+q_{1}^{2}}\left\\{\left[(\epsilon-1)q_{1}^{2}+2m^{2}\right]g^{\mu\nu}-2(\epsilon-1)q_{1}^{\mu}q_{1}^{\nu}\right\\}$
(3)

with the help of the Mathematica package HepLib.m using the following
Mathematica code:

⬇

Get["https://heplib.github.io/HepLib.m"]; "output from HepLib"//C2M//Simplify

One can see the appearance of tensor integrals (the integrals with free
indices) involving $q_{1}^{\mu}$ or $q_{1}^{\nu}$ in (3), those can not be
handled directly by IBP reduction tools, so one needs to transfer the Lorentz
indices from loop momenta to external momenta in the loop integration, we call
it Tensor Index Reduction (TIR). We illustrate the basic ideas behind TIR by
the tensor integral above, by writing down the following substitution
according to the Lorentz invariance:

$q_{1}^{\mu}q_{1}^{\nu}\to c_{00}(p^{2},p\cdot
q_{1},q_{1}^{2})\,g^{\mu\nu}+c_{11}(p^{2},p\cdot
q_{1},q_{1}^{2})\,p^{\mu}p^{\nu}$ (4)

and multiple both sides of (4) by $g^{\mu\nu}$ and $p^{\mu}p^{\nu}$, and get
coefficients $c_{00}$ and $c_{11}$ by solving the resulting linear equations.
It should be noted that TIR is only valid when we perform loop integration on
both sides, i.e.,

$\int\frac{d^{D}q_{1}}{(2\pi)^{D}}\,{\cal F}(p^{2},p\cdot
q_{1},q_{1}^{2})\,q_{1}^{\mu}q_{1}^{\nu}=g^{\mu\nu}\int\frac{d^{D}q_{1}}{(2\pi)^{D}}\,c_{00}{\cal
F}(p^{2},p\cdot
q_{1},q_{1}^{2})+p^{\mu}p^{\nu}\int\frac{d^{D}q_{1}}{(2\pi)^{D}}\,c_{11}{\cal
F}(p^{2},p\cdot q_{1},q_{1}^{2})$

where ${\cal F}$ is any function depending on the scalar products $p^{2}$,
$p\cdot q_{1}$ and $q_{1}^{2}$. The Tensor Index Reduction in HepLib is
encoded in the function TIR555The related functions in FeynCalc are known as
TID and TIDL. introduced as follows:

⬇

ex TIR(const ex &expr_in, const lst &loop_ps, const lst &ext_ps);

where loop_ps refers to the list of loop momenta, ext_ps to the list of
external momenta, and one can get the result as follows, after invoking
function call TIR(amp1,lst{q1},lst{p}):

$\sqrt{\frac{2}{3}}8ig_{s}^{2}\,g^{\mu\nu}\frac{1}{q_{1}^{2}}\frac{1}{2p\cdot
q_{1}-q_{1}^{2}}\frac{1}{2p\cdot
q_{1}+q_{1}^{2}}\frac{(2\epsilon^{2}-3\epsilon+1)m^{2}q_{1}^{2}+2(2\epsilon-3)m^{4}-2(\epsilon-1)(p\cdot
q_{1})^{2}}{(2\epsilon-3)m^{2}}$ (5)

as one can see the tensor integrals are reduced to the scalar integrals, and
it is ready to perform the IBP reduction on those integrals which will be
introduced in the next section.

###  2.5 Partial fraction decomposition and IBP reduction

Usually, one gets the scalar loop integrals with propagators that form an
overdetermined basis, for example, one of the scalar integrals appearing in
(5) is as follows:

${\cal I}=\int\frac{d^{D}q_{1}}{(2\pi)^{D}}(q_{1}^{2})^{-1}(2p\cdot
q_{1}-q_{1}^{2})^{-1}(2p\cdot q_{1}+q_{1}^{2})^{-1}$ (6)

if we take the scalar products involving loop momentum (i.e., $q_{1}^{2}$ and
$p\cdot q_{1}$) as independent vectors, then the three propagators $q_{1}^{2}$
and $(2p\cdot q_{1}\pm q_{1}^{2})$ are linear dependent. The original
Mathematica package Apart can be used to carry out the partial fraction
decomposition automatically on such integrals, these operations have been
reimplemented in HepLib too, through the function HepLib::Apart introduced as
follows:

⬇

ex Apart(const ex &expr_in, const lst &loop_ps, const lst & ext_ps);

where loop_ps and ext_ps are the same as those in TIR function introduced in
last section. One can carry out the partial fraction decomposition on (6)
using the following snip of the code 7.cpp:

⬇

Vector p("p"), q1("q1");

ex expr = 1/SP(q1) * 1/(2*SP(p,q1)-SP(q1)) * 1/(2*SP(p,q1)+SP(q1));

ex r = Apart(expr,lst{q1},lst{p});

ex r1 = ApartIR2ex(r); ex r2 = ApartIR2F(r);

where the returned expression r from Apart involves the pseudo-function
ApartIR (an internal matrix representation) which can be transformed into a
normal ex expression or F function (as in the FIRE program) with the help of
ApartIR2ex and ApartIR2F respectively, for example, the results r1 and r2 from
ApartIR2ex and ApartIR2F look as follows:

⬇

r1:
(-1/2)*q1.q1^(-2)*(q1.q1+2*q1.p)^(-1)+(1/2)*q1.q1^(-2)*(2*q1.p+(-1)*q1.q1)^(-1)

r2:
(-1/2)*F({q1.q1,q1.q1+(-2)*q1.p},{2,1})+(-1/2)*F({q1.q1,q1.q1+2*q1.p},{2,1})

and in more readable form as follows:

$-\frac{1}{2}(q_{1}^{2})^{-2}(q_{1}^{2}+2p\cdot
q_{1})^{-1}-\frac{1}{2}(q_{1}^{2})^{-2}(q_{1}^{2}-2p\cdot q_{1})^{-1}$ (7)

where one can see that there are only two independent propagators in each term
after the decomposition. Furthermore, there is a more general version of Apart
function introduced as follows:

⬇

ex Apart(const ex &expr_in, const lst &vars, exmap sign_map={});

where vars refers to a list of independent vectors, and the internal details
can be found in [13], here one can also get the same result as r with the
following function call:

⬇

rr = Apart(expr, lst{SP(q1),SP(p,q1)}); // rr is the same as r

Now it is ready to perform the actual IBP reduction by the functions or
classes introduced under the namespace HepLib::IBP. The classes FIRE and KIRA
are introduced in the current version, both are derived from the Base class:

⬇

class Base {

public:

lst Internal; lst External; lst Replacements; lst Propagators;

lst Integrals; lst Cuts; lst Rules; lst MIntegrals;

virtual void Export(); virtual void Run(); virtual void Import();

void Reduce();

};

where the function Reduce is just sequential calls to the virtual functions
Export, Run and Import. The classes FIRE and KIRA have implemented the virtual
function Export to prepare the input information for calling FIRE and KIRA
programs, the virtual function Run to invoke the shell command FIRE5/FIRE6 or
kira respectively, and the virtual function Import to get the reduced results
back to HepLib.

The basic usage of each class is similar to FIRE program by providing the
necessary input parameters, including Internal (the loop momenta), External
(the external momenta), Propagators (a list of complete and linear-independent
propagators), Replacements (the scalar product substitution) and Integrals
(the integrals needed to be reduced), then calling Reduce subroutine, HepLib
will generate the IBP identities automatically with respect to provided
Internal and External, and export the corresponding files which are needed for
the FIRE or KIRA program. After Reduce, the members Rules and MIntegrals will
get updated accordingly if everything works fine.

We exemplify the basic usage of FIRE on the first integral in (7), by the
following snip of the code 8.cpp:

⬇

Symbol q1("q1"), p("p"), m("m");

FIRE fire; fire.Internal=lst{ q1 }; fire.External=lst{ p };

fire.Replacements=lst{ p*p == m*m };

fire.Propagators=lst{ q1*q1, 2*p*q1-q1*q1 };

fire.Integrals.append(lst{2,1}); // one can add more integrals

fire.WorkingDir="IBPdir";

fire.Reduce(); // fire.Rules and fire.MIntegrals got updated

where the C++ version of FIRE program is used, both FIRE6 and FIRE5 are
supported, the default one is to use FIRE6 (to use FIRE5, one needs to set the
static member FIRE::Version to 5). Note that the .start files are generated by
HepLib itself, which are different from the original ones generated from the
Mathematica package FIRE.m. After the function call fire.Reduce(), one can get
the reduction result from the member fire.Rules as follows:

⬇

{ F(0,{2,1})==(1/4)*m^(-4)*F(0,{0,1})*((-5)+d)^(-1)*((-2)+d) }

which represents the replacing rule $\\{F(0,\\{2,1\\})\to(d-2)/(d-5)\times
F(0,\\{0,1\\})/(4m^{4})\\}$, and the resulting master integrals from
fire.MIntegrals, and there is only one master integral F(0,{0,1}) in this
case. Usually, one will get a number of master integrals that are not
independent to each other, the function IBP::FindRules is provided to find the
relations between those master integrals, the basic idea is almost the same as
the one in FIRE package, which uses permuted UF function (a routine to
calculate Symanzik polynomials) to find identical integrals.

All those operations (from partial fraction decomposition to IBP reduction)
are wrapped in a single function call HepLib::ApartIBP introduced as follows:

⬇

void ApartIBP(int IBPmethod, exvector &io_vec, const lst & loops_exts=lst{});

with IBPmethod=1 for using FIRE and 2 for KIRA, the function will call Apart
for partial fraction decomposition, Reduce for IBP reduction and FindRules to
eliminate identical master integrals. Note that the input argument io_vec will
get updated at the end of the function call.

Now, let’s come back to the amp1 (the amplitude for $Q\bar{Q}\to\gamma^{*}$ at
next-to-leading order), we can summarize all the necessary function calls in
the following snip of the code nlo.cpp:

⬇

auto amp1 = garRead("amps.gar").op(1).op(0);

Index mu("mu"), nu("nu");

letSP(p)=m*m; letSP(p,nu)=0;

lst extps = lst{ p };

lst loops = lst{ q1 };

ex proj = SpinProj(IO::In,1,p1,p2,m,m,nu,-1,-3) * ColorProj(-1,-3);

amp1 = amp1 * proj; amp1 = amp1.subs(LI(-2)==mu);

amp1 = MatrixContract(amp1);

amp1 = form(amp1);

amp1 = TIR(amp1, loops, extps);

amp1 = amp1.subs(lst{NF==3,D==4-2*ep});

exvector res_vec;

res_vec.push_back(amp1);

ApartIBP(1, res_vec, lst{loops, extps});

ex res = 0;

for(auto item : res_vec) res += item.subs(lst{d==4-2*ep,D==4-2*ep});

ex mi = str2ex("F({q1^2+(-2)*p*q1},{1})"); // the master integral

ex miv = -tgamma(-1+ep)*pow(m,2*(1-ep))*I*pow(Pi,2-ep)*pow(2*Pi,2*ep-4);

res = res.subs(lst{ mi==miv }); // substitutions here

res = series_ex(res,ep,0); // taylor expand to ep^0

where we use the replacement subs to substitute the master integral with its
explicit expression, i.e.,

$\mbox{\tt
F({q1\^{}2+(-2)*p*q1},{1})}=\int\frac{dq_{1}^{D}}{(2\pi)^{D}}\frac{1}{q_{1}^{2}-2p\cdot
q_{1}}=-\frac{i\pi^{D/2}}{(2\pi)^{D}}\Gamma(-1+\epsilon)m^{2(1-\epsilon)}$ (8)

and we end with the following result for the amplitude at next-to-leading
order (the bare result, not including the renormalization yet):

⬇

(-1/4)*sqrt(2)*sqrt(3)*ep^(-1)*gs^2*Pi^(-2)*nu.mu+(1/12)*sqrt(2)*sqrt(3)*gs^2*Pi^(-2)*nu.mu*(4+6*log(m)+(-6)*log(2*Pi)+3*Euler+3*log(Pi))

and in more readable form as follows:

$\sqrt{6}\frac{g_{s}^{2}}{4\pi^{2}}g^{\mu\nu}\left[-\frac{1}{\epsilon}+\gamma_{E}-\ln(4\pi)+\ln
m^{2}+\frac{4}{3}\right]$ (9)

where the coupling constant $g_{s}$ is still bare, not renormalized, and each
bare $g_{s}$ will be associated with a $\mu^{\epsilon}$, so the logarithm $\ln
m^{2}$ will be transformed into $\ln(m^{2}/\mu^{2})$ in the final renormalized
result.

One can do similar operations on the amplitudes at next-to-next-to-leading
order, while one still needs the renormalization to get the final physical
result. One can adopt the multiplicative renormalization, e.g., in the
amplitudes generated from QGRAF, the quark mass in each propagator is treated
as a bare mass $m_{0}$, its relation to the renormalized mass $m$ is
$m_{0}=Z_{m}m$, where the $Z_{m}$ is the renormalization constant, the QCD
related renormalization constants up to two-loop order are taken from [102]
and encoded in the namespace HepLib::QGRAF::RC.

Finally, all the related operations for the process $Q\bar{Q}\to\gamma^{*}$ up
to next-to-next-to-leading order (from amplitude/diagram generation, amplitude
evaluation, IBP reduction to the substitution of master integrals, as well as
the renormalization) can be found at codes/nnlo.cpp in the HepLib archive, and
one can get the following final expression from the output of nnlo.cpp:

$\frac{{\cal A}_{\rm NNLO}}{{\cal A}_{\rm
LO}}=1-\frac{8}{3}\frac{\alpha_{s}}{\pi}+\frac{\alpha_{s}^{2}}{\pi^{2}}\left[\frac{35\pi^{2}}{54}\left(2\ln\frac{m^{2}}{\mu^{2}}-\frac{1}{\epsilon}\right)+\frac{50}{9}\ln\frac{m^{2}}{\mu^{2}}-\frac{125\zeta(3)}{9}-\frac{511\pi^{2}}{324}-\frac{23}{54}-\frac{14}{9}\pi^{2}\ln
2\right]$ (10)

where we have substituted $C_{A}=3$, $C_{F}=4/3$ and $n_{f}=4$ (4 active light
quarks) for simplicity, and the remaining single pole in (10) still needs to
be factorized/absorbed into the NRQCD long-distance matrix element,
furthermore the $\mu$ in the first logarithm corresponds to the factorization
scale $\mu_{F}$, and $\mu$ in the other logarithm refers to the
renormalization scale $\mu_{R}$, so one arrives at the final expression after
these manipulations:

$1-\frac{8}{3}\frac{\alpha_{s}}{\pi}+\frac{\alpha_{s}^{2}}{\pi^{2}}\left[\frac{35\pi^{2}}{27}\ln\frac{m^{2}}{\mu_{F}^{2}}+\frac{50}{9}\ln\frac{m^{2}}{\mu_{R}^{2}}-\frac{125\zeta(3)}{9}-\frac{511\pi^{2}}{324}-\frac{23}{54}-\frac{14}{9}\pi^{2}\ln
2\right]$ (11)

and one can check the correctness with [98, 99], e.g., the $c_{2}(m_{Q}/\mu)$
in [99] as follows:

$\displaystyle c_{2}(m_{Q},\mu)$ $\displaystyle=$ $\displaystyle
C_{F}^{2}c_{\rm 2,A}+C_{F}C_{A}c_{\rm 2,NA}+C_{F}T_{F}n_{f}c_{\rm
2,L}+C_{F}T_{F}c_{\rm 2,H}$ (12) $\displaystyle=$
$\displaystyle\frac{35}{27}\pi^{2}\ln\frac{m_{Q}^{2}}{\mu^{2}}-\frac{125\zeta(3)}{9}-\frac{511\pi^{2}}{324}-\frac{23}{54}-\frac{14}{9}\pi^{2}\ln
2$

where the coefficients $c_{\rm 2,A}$, $c_{\rm 2,NA}$, $c_{\rm 2,L}$ and
$c_{\rm 2,H}$ can be found in [99], and note that $\mu_{R}=m_{Q}$ is taken
explicitly in those expressions.

###  2.6 Numerical evaluation by sector decomposition method

Sector Decomposition is a general numerical method used to extract the
singularity of Feynman integral within dimensional regularization, we refer
the reader to [35] for a good review on the method of sector decomposition.

Principally, sector decomposition can only resolve the singularity located in
the boundary of the integration domain, to handle the divergence from the
internal of the domain, one usually resorts to the integration contour
deformation to avoid crossing the internal singularity. A practical
programable contour deformation is to use the following transformation for the
corresponding $\cal{F}$-term:

$\displaystyle x_{i}\to
z_{i}=x_{i}-i\lambda_{i}x_{i}(1-x_{i})\frac{\partial\cal F}{\partial x_{i}}$
(13)

where we have adopted the notation as FIESTA, one can always find some
{$\lambda_{i}$} for the transformation above, such that the imaginary part of
${\cal F}$ is always negative for the whole integral domain, which is
consistent with the $-i\epsilon$-prescription in the original Feynman
integral. One also notes that the efficiency for the numerical integration is
highly dependent on the chosen {$\lambda_{i}$} in some extreme cases, and
there is no priori method to select an optimized {$\lambda_{i}$} yet. The
public programs FIESTA and SecDec have already implemented such a
transformation.

We present another implementation in C++ language with the help of GiNaC,
different from the package sector_decomposition which also chooses the GiNaC
to handle the symbolic operations, we implement the sector decomposition by
the geometric strategy with the help of qhull library (it is almost a GiNaC
conversion of STRATEGY_KU2 in FIESTA), we also develop some new features,
including contour deformation as well.

All related classes or functions are introduced under the namespace
HepLib::SD, the basic usage is very similar to FIESTA program by calling the
Evaluate method of SecDec class after providing the necessary information for
the integral. To facilitate the parameter input, two structures are introduced
in HepLib, one is FeynmanParameter which is used to wrap the information of
loop integration:

⬇

struct FeynmanParameter {

lst LoopMomenta; lst tLoopMomenta; lst Propagators; lst Exponents;

exmap lReplacements; exmap tReplacements; exmap nReplacements; ex Prefactor =
1;

};

and the other is XIntegrand which is introduced for the generic parameter
integration:

⬇

struct XIntegrand {

lst Functions; lst Exponents; exmap nReplacements; lst Deltas; bool isAsy =
false;

};

The basic usage involving sector decomposition method can be illustrated by
numerically evaluating the integral ${\cal I}_{1}$ in dimensional
regularization with dimension $D=4-2\varepsilon$:

$\displaystyle{\cal
I}_{1}=\int\frac{d^{D}k\,d^{D}r\,d^{D}q}{(i\pi^{2-\varepsilon})^{3}\Gamma^{3}(1-\varepsilon)}\frac{1}{(-k^{2})^{1+3\varepsilon}}\frac{1}{-(k+p_{1}+p_{2})^{2}}\frac{1}{-(-k+r)^{2}}\frac{1}{-(p_{1}+r)}\frac{1}{-(k-q)^{2}}\frac{1}{-(p_{1}+q)^{2}}$
(14)

where $p_{1}^{2}=p_{2}^{2}=0$ and $s=2p_{1}\cdot p_{2}=-1$, using the
following snip of the code 9.cpp:

⬇

Symbol k("k"),r("r"),q("q"),p1("p1"),p2("p2"),s("s");

FeynmanParameter fp;

fp.LoopMomenta = lst{k,r,q};

fp.Propagators= lst { -pow(k,2),-pow(k+p1+p2,2),-pow(-k+r,2),-pow(p1+r,2),

-pow(k-q,2),-pow(p1+q,2)};

fp.Exponents = lst{1+3*ep,1,1,1,1,1};

fp.lReplacements[p1*p1] = 0;

fp.lReplacements[p2*p2] = 0;

fp.lReplacements[p2*p1] = s/2;

fp.lReplacements[s] = -1;

fp.Prefactor = pow(I*pow(Pi,2-ep),-3) * pow(tgamma(1-ep),3);

SecDec work; work.epN = 0; Verbose = 2; work.Evaluate(fp);

hout << work.VEResult() << endl;

where Verbose=2 is used for some detailed output, and epN=0 refers to
evaluating the integral to ${\cal O}(\varepsilon^{0})$, the hout statement
will result in the output as follows:

⬇

1.2532095(6)E2+ep^(-3)*1.666666667(4)E-1+ep^(-2)*1.83333333(6)+ep^(-1)*1.812319(2)E1

which corresponds to

$\displaystyle{\cal I}_{1}=125.32095(6)+\varepsilon^{-3}\times
0.1666666667(4)+\varepsilon^{-2}\times 1.83333333(6)+\varepsilon^{-1}\times
18.12319(2)+{\cal O}(\varepsilon)$ (15)

One can see that all operations are encoded in the function call to Evaluate
method of SecDec, which actually contains the following basic procedures while
using sector decomposition method:

⬇

1. Initialize(fp); 2. Scalelesses(); 3. KillPowers();

4. RemoveDeltas(); 5. KillPowers(); 6. SDPrepares();

7. EpsEpExpands(); 8. CIPrepares(); 9. Contours(); 10.Integrates();

They are the member functions of class SecDec parts of which are listed below:

⬇

class SecDec {

public:

int epN = 0; int epsN = 0; bool CheckEnd = false;

vector<ex> FunExp; vector<ex> Integrands; vector<ex> expResult;

SecDecBase *SecDec = NULL; IntegratorBase *Integrator = NULL;

MinimizeBase *Minimizer = NULL;

ex ResultError;

void Initialize(FeynmanParameter fpi); void Initialize(XIntegrand xint);

void Normalizes(); void Scalelesses();

void SDPrepares(); void EpsEpExpands(); void RemoveDeltas();

void CIPrepares(const string & key = "");

void Contours(const string & key = "", const string & pkey = "");

void Integrates(const string & key="", const string & pkey="", int kid=0);

void Evaluate(FeynmanParameter fpi, const string & key = "");

void Evaluate(XIntegrand xint, const string & key = "");

void Evaluate(vector<ex> FunExp, const string & key = "");

void MB(); void XEnd(); void ChengWu();

}

The detailed descriptions of the members or functions in SecDec, we refer the
reader to the document [97].

The key operations involving sector decomposition methods are sequential calls
from SDPrepares to Integrates: SDPrepares is used to perform the sector
decomposition currently with the geometric strategy, EpsEpExpands is used to
expand the integrand around eps=0 and ep=0 (eps is another infinitesimal
regulator similar to the dimensional regulator ep, and note that the expansion
around eps=0 will be performed before the expansion of ep), CIPrepares is used
to generate C++ source code for the integrand, the C++ code will be compiled
to a shared .so library, Contours is used to perform contour deformation if
the ${\cal F}$-term has negative terms, Integrates is used to perform the
actually numerical integrations, with the help of dl library, the integrands
are dynamically loaded from the compiled .so files generated in CIPrepares.

HepLib currently supports two numerical integrators CUBA and HCubature [103],
a tiny modified version666A PrintHooker has been added to the original version
of HCubature to retrieve the intermediate integral result and estimated error.
of HCubature is shipped along with the library, the interface to each
integrator is derived from the base class IntegratorBase [97] which wraps a
few basic parameters to control the numerical integration. Since the quadruple
precision is adopted in the base class IntegratorBase, the library libcubaq is
actually used. Although the default precision of the numerical integrator is
quadruple, one can still use other numerical precision to evaluate the
integrand at a specific input point and cast the computed values to quadruple
precision at the end. There are three internal types of numeric precision in
the HepLib: long double, quadruple and arbitrary precision by MPFR, one of the
precisions will be chosen according to the input arguments supplied to the
integrand function call, the controlling parameters can be found in the
IntegratorBase.

For negative ${\cal F}$-terms, HepLib adopts contour deformation method to
avoid the singularity inside the integration domain, the $\lambda$-tuning is
performed dynamically during numerical integration, this is achieved by
exporting related ${\cal F}$-terms and their derivatives to the generated C++
source, and the variable substitution defined in (13) actually occurs at
evaluation running time, instead of the static compilation, generally it will
results in a smaller size of source code and a faster running time by avoiding
evaluating the same variable substitution multiple times.

In addition, we want to address some other methods than contour deformation,
including WickRotation and ChengWu method, by using the following integral
${\cal I}_{2}$:

$\displaystyle{\cal
I}_{2}=\int_{0}^{1}dx_{1}\int_{0}^{1}dx_{2}\int_{0}^{1}dx_{4}\;x_{1}^{-\frac{5}{2}+3\varepsilon}x_{2}^{-\frac{3}{2}+\varepsilon}\left(-4x_{4}^{2}-12x_{1}x_{4}+x_{1}x_{2}\right)^{\frac{1}{2}-2\varepsilon}\delta(x_{1}+x_{2}+x_{4}-1)$
(16)

The WickRotation method is illustrated by the following snip of the code
wr.cpp:

⬇

ex wra = WRA(Pi/5);

XIntegrand xint;

xint.Functions = lst{1,-4*x(4)*x(4)-12*x(1)*x(4)+x(1)*x(2),x(1),x(2)};

xint.Exponents = lst{1,1/ex(2)-2*ep,3*ep-5/ex(2),ep-3/ex(2)};

xint.Deltas = lst{lst{x(1),x(2),x(4)}};

xint.Functions = ex_to<lst>(subs(xint.Functions,{x(4)==x(4)*exp(I*wra)}));

xint.Functions.let_op(0) = exp(I*wra); // Jacobi

SecDec work;

work.Initialize(xint);

... ...

where we have used XIntegrand for the generic parameter integration in sector
decomposition method, note that HepLib adopts the notation used in FIESTA: the
${\cal F}$-term is always supplied at 2nd position of Functions and compatible
with $-i\epsilon$ convention. In WickRoation evaluation of the integral ${\cal
I}_{2}$, we can first pull out $x_{4}$ from the $\delta$-function with the
help of Cheng-Wu theorem [104, 105], then make the transformation $x_{4}\to
x_{4}e^{i\pi/5}$ with $\pi/5$ as the rotation angle, finally add $x_{4}$ back
to the $\delta$-function, it is ready to see the imaginary of the ${\cal
F}$-term is compatible with $-i\epsilon$ convention when the rotation angle
changes from $0$ to $\pi/5$. The calculated numerical value of ${\cal I}_{2}$
up to ${\cal O}(\varepsilon)$ (produced by the codes/wr.cpp) is as follows:

$\displaystyle{\cal
I}_{2}=-\frac{1.5}{\varepsilon}+(3.29583687(5)-4.71238898(5)i)+(14.716159(1)+10.354177(1)i)\varepsilon+{\cal
O}(\varepsilon^{2})$ (17)

The integral ${\cal I}_{2}$ can also be evaluated by the ChengWu method by the
snip code cw.cpp

⬇

XIntegrand xint;

xint.Functions = lst{1,-4*x(4)*x(4)-12*x(1)*x(4)+x(1)*x(2),x(1),x(2)};

xint.Exponents = lst{1,1/ex(2)-2*ep,3*ep-5/ex(2),ep-3/ex(2)};

xint.Deltas = lst{lst{x(1),x(2),x(4)}};

SecDec work;

work.Initialize(xint);

work.ChengWu(); // call ChengWu method

... ...

note the function call to ChengWu method, by Cheng-Wu theorem, $x_{1}$ and
$x_{2}$ can be pulled out of the $\delta$-function, sequential transformations
are then performed automatically in HepLib: (1) $x_{1}\to x_{1}^{\prime}/3$,
(2) $x_{2}\to
4x_{4}x_{2}^{\prime}/x_{1}=12x_{4}x_{2}^{\prime}/x_{1}^{\prime}$, with Jacobi
terms properly added, the corresponding transformed ${\cal F}$-term now
becomes $(x_{2}^{\prime}-x_{1}^{\prime}-x_{4})$ (generally, the transformed
${\cal F}$-term will be $(x_{1}+\cdots+x_{m}-x_{m+1}-\cdots-x_{n})$), again
the transformation $x_{2}^{\prime}\to y_{2}(x_{1}^{\prime}+x_{4})/y_{1}$ is
made to result in the final ${\cal F}$-term $(y_{1}-y_{2})$, note that the
positive coefficient of ${\cal F}$-term is always factorized out. Both $y_{1}$
and $y_{1}$ are integrated from $0$ to $\infty$, so we can divide the domain
into $y_{1}\geq y_{2}$ and $y_{1}\leq y_{2}$, make transformation $y_{1}\to
y_{1}+y_{2}$ and $y_{2}\to y_{1}+y_{2}$ for those two regions respectively,
the resulted ${\cal F}$-term will be proportional to $y_{1}$ and $(-y_{2})$,
and the negative one can be readily changed to positive with an overall
constant factor using
$(-y_{2}-i\epsilon)^{a+b\varepsilon}=e^{-i\pi(a+b\varepsilon)}y_{2}^{a+b\varepsilon}$,
finally $y_{1}$ or $y_{2}$ is changed back to $x_{1}$ or $x_{2}$ respectively,
which will be added back to $\delta$-function with the help of Cheng-Wu
theorem. All those internal operations are actually performed by the functions
from a separate class ChengWu [97]. Furthermore, the integral ${\cal I}_{2}$
can also be cross-checked by contour deformation by using the code
codes/ct.cpp in HepLib archive.

Note that, when the integrand does not encounter any singularity in the whole
domain with the help of contour deformation or the WickRoation method, the
numerical integration itself does not pose any problem and always results in a
finite result, however we need pay attention to the branch cut of multiple-
valued functions (including the log or power). HepLib treats the branch cut in
a simple way, taking a log term $\ln(g(z_{i}(\vec{\lambda}))$ in the method of
contour deformation as an example, we multiple $\vec{\lambda}$ by a scaling
parameter $s$, when the parameter $s$ varies from $0$ to $1$, we can check
whether the argument of log function $g(z_{i}(\vec{\lambda}))$ goes through
the branch cut or not. Practically, we let $s$ take $n$ values $s_{j}=j/n$
($0<j\leq n$), then compare the corresponding $g(z_{i}(s_{k}\vec{\lambda}))$
and $g(z_{i}(s_{k+1}\vec{\lambda}))$, if the line segment between the two
points goes through the branch cut, a $(2\,i\,\pi)$ is added or subtracted
accordingly.

Finally, we want to demonstrate the automatic expansion of loop/parameter
integrals with respect to a small parameter by combining the sector
decompositions and Mellin–Barnes representations (initially suggested in
[106]), by using the example from [40] in the following snip of the code
exp.cpp:

⬇

auto t = vs; // the global Symbol vs

FeynmanParameter fp;

fp.LoopMomenta = lst{ k };

fp.Propagators = lst{
-pow(k,2),-pow(k+p1,2),-pow(k+p1+p2,2),-pow(k+p1+p2+p4,2) };

fp.Exponents = lst{ 1, 1, 1, 1 };

fp.lReplacements[p2*p4] = -t/2;

fp.lReplacements[p1*p4] = (s+t)/2;

... ...

SecDec work;

work.Evaluate(fp);

the only difference from the unexpanded case is to use the global Symbol vs to
indicate the small parameter, SecDec class will handle the expansion
automatically and result in the following expansion:

$\displaystyle\int\frac{d^{4-2\varepsilon}k}{i\pi^{2-\varepsilon}e^{-\varepsilon\gamma_{E}}}\frac{1}{k^{2}}\frac{1}{(k+p_{1})^{2}}\frac{1}{(k+p_{1}+p_{2})^{2}}\frac{1}{(k+p_{1}+p_{2}+p_{4})^{2}}=\frac{4}{\varepsilon^{2}t}-\frac{2\ln
t}{\varepsilon t}-\frac{13.159472535(0)}{t}+{\cal O}(\varepsilon,t)$ (18)

with the kinematics $p_{1}^{2}=p_{2}^{2}=p_{4}^{2}=0$, $p_{1}\cdot
p_{2}=-s/2$, $p_{2}\cdot p_{4}=-t/2$, $p_{1}\cdot p_{4}=(s+t)/2$ and $s=1$.

##  3 Installation and Usage

Before we go into the details of the installation, one can use the script
install.sh or makefile to install the required external libraries and programs
automatically, it will install HepLib as well, by invoking one of the
following commands in the terminal:

  * 1.

using install.sh

⬇

wget https://heplib.github.io/install.sh

chmod +x install.sh

INSTALL_PATH=<INSTALL PATH> jn=<jn> ./install.sh

  * 2.

using makefile

⬇

wget https://heplib.github.io/makefile

make INSTALL_PATH=<INSTALL PATH> jn=<jn>

where the <INSTALL PATH> is the path for the libraries to be installed to, and
<jn> is the number of jobs used in invoking make -j $jn in the script.

Since GiNaC or HepLib will get updated from time to time, one can use the
update script (located in the directory <INSTALL PATH>/bin/) to only update
GiNaC or HepLib to their latest versions:

⬇

jn=<jn> <INSTALL PATH>/bin/update ginac # update both GiNaC and HepLib

jn=<jn> <INSTALL PATH>/bin/update heplib # only update HepLib

###  3.1 Prerequisites

Since HepLib uses several external routines or libraries, one needs to install
these required libraries before the installation of HepLib, one also needs
cmake and GNU compiler system for the compilation as well. Those prerequisites
can be divided into two categories: external libraries and external programs,
only external libraries are required for the compilation or installation of
HepLib, while to run the user program properly which uses HepLib, it’s still
required for the external programs to be found in the environment variable
PATH. Note that there is no need to adjust the variable PATH if one installs
HepLib using the install script install.sh or makefile.

Before the installation, one needs to provide the installation destination
through the environment variable $INSTALL_PATH by typing the commands in the
terminal:

⬇

export INSTALL_PATH="<INSTALL PATH>"

  * 1.

External libraries:

    1. (a)

GMP is required by MPFR and GiNaC. The typical instructions to download and
install GMP (v6.2.1) is as follows:

⬇

curl -L -O https://gmplib.org/download/gmp/gmp-6.2.1.tar.bz2

tar jxf gmp-6.2.0.tar.bz2;cd gmp-6.2.1

./configure \--prefix=$INSTALL_PATH

make -j 16 && make install

    2. (b)

MPFR is used to handle the multiple precision in the numerical integration
when large number cancelation occurs, since we choose quadruple precision as
the default float precision type in the numerical integrator, MPFR needs to be
compiled with the option \--enable-float128. Note that the quadruple precision
type __float128 has been changed to _Float128 since MPFR 4.1.0, so we prefer
the version MPFR 4.0.2 for the moment, furthermore the MPFR C++ [107] wrapper
is included in HepLib archive. The typical instructions to download and
install MPFR (v4.0.2) is as follows:

⬇

curl -L -O https://heplib.github.io/download/mpfr-4.0.2.tar.gz

tar zxf mpfr-4.0.2.tar.gz;cd mpfr-4.0.2

./configure \--prefix=$INSTALL_PATH \--enable-float128 \--enable-thread-safe

make -j 16 && make install

    3. (c)

CLN is required by GiNaC. The typical instructions to download and install CLN
(v1.3.6) is as follows:

⬇

curl -L -O https://www.ginac.de/CLN/cln-1.3.6.tar.bz2

tar jxf cln-1.3.6.tar.bz2;cd cln-1.3.6

./configure \--prefix=$INSTALL_PATH \--with-gmp=$INSTALL_PATH

make -j 16 && make install

    4. (d)

GiNaC is actually the underlying language of HepLib, all symbolic functions
are built on top of it. The typical instructions to download and install GiNaC
(v1.8.0) is as follows:

⬇

curl -L -O https://www.ginac.de/ginac-1.8.0.tar.bz2

tar jxf ginac-1.8.0.tar.bz2;cd ginac-1.8.0

./configure \--prefix=$INSTALL_PATH
PKG_CONFIG_PATH=$INSTALL_PATH/lib/pkgconfig

make -j 16 && make install

    5. (e)

QHull is used in the sector decomposition when using the geometric strategy.
The typical instructions to download and install QHull (v2020.2) is as
follows:

⬇

curl -L -O http://www.qhull.org/download/qhull-2020.2.zip

unzip -q qhull-2020.2.zip;cd qhull-2020.2

cp Makefile Makefile.bak

cat Makefile.bak | sed "s/\/usr\/local/\$\$INSTALL_PATH/g" > Makefile

make && make install

    6. (f)

MinUit2 is used to find the numerical minimum of a function. The typical
instructions to download and install MinUit2 (v5.34.14) is as follows:

⬇

curl -L -O http://project-mathlibs.web.cern.ch/project-
mathlibs/sw/5_34_14/Minuit2/Minuit2-5.34.14.tar.gz

tar zxf Minuit2-5.34.14.tar.gz;cd Minuit2-5.34.14

./configure \--prefix=$INSTALL_PATH

make -j 16 && make install

    7. (g)

CUBA is one of the numerical integrators with many different numerical
integration routines, note that libcubaq is actually used by providing the
option \--with-real=16 to the configure script. The typical instructions to
download and install CUBA (v4.2) is as follows:

⬇

curl -L -O http://www.feynarts.de/cuba/Cuba-4.2.tar.gz

tar zxf Cuba-4.2.tar.gz;cd Cuba-4.2

./configure \--prefix=$INSTALL_PATH \--with-real=16 CFLAGS="-fPIC -fcommon"

make && make install

  * 2.

External programs:

    1. (a)

Fermat is used for high-performance operations on matrix or rational
polynomial. The typical installation instructions for Linux are as follows:

⬇

curl -L -O http://home.bway.net/lewis/fermat64/ferl6.tar.gz

tar zxf ferl6.tar.gz;mv ferl6 $INSTALL_PATH;

cd $INSTALL_PATH/bin;ln -s -f ../ferl6/fer64 .

    2. (b)

FORM is used for Dirac/Color matrix trace, Lorentz index contraction and so
on. The typical installation instructions for Linux are as follows:

⬇

curl -L -O
https://github.com/vermaseren/form/releases/download/v4.2.1/form-4.2.1-x86_64-linux.tar.gz

tar zxf form-4.2.1-x86_64-linux.tar.gz

cp -rf form-4.2.1-x86_64-linux/form $INSTALL_PATH/bin/

cp -rf form-4.2.1-x86_64-linux/tform $INSTALL_PATH/bin/

    3. (c)

FIRE is required for IBP reduction using FIRE class. The typical installation
instructions for Linux are as follows:

⬇

git clone https://bitbucket.org/feynmanIntegrals/fire.git

mv fire/FIRE6 $INSTALL_PATH/FIRE6

cd $INSTALL_PATH/FIRE6

./configure && make -j 16 && make

    4. (d)

KIRA is used for IBP reduction using KIRA class. The typical installation
instructions for Linux are as follows:

⬇

curl -L -o kira https://kira.hepforge.org/downloads?f=binaries/kira-2.0

chmod +x kira;mv -f kira $INSTALL_PATH/bin/kira

###  3.2 Installation of HepLib

As the required external libraries have been installed, the most recent
version of HepLib is available for download from
https://heplib.github.io/HepLib.tar.gz as a compressed archive. One can
uncompress the archive and change the current directory into the extracted
directory HepLib/src by typing the commands in the terminal:

⬇

wget https://heplib.github.io/HepLib.tar.gz

tar zxfv HepLib.tar.gz;cd HepLib/src

and create a directory build for cmake to build the library as follows:

⬇

mkdir build && cd build

cmake -DCMAKE_INSTALL_PREFIX=<INSTALL PATH> ..

make -j 4 && make install

where the standard cmake variable CMAKE_INSTALL_PREFIX refers to the directory
to which HepLib will be installed, i.e., the library libHepLib.so (the file
name may be system dependent) will be installed to <INSTALL PATH>/lib, the
related C++ header files will be installed to <INSTALL PATH>/include, and the
binary programs, including heplib++, garview, etc., will be installed to
<INSTALL PATH>/bin.

If GiNaC or other dependent external library is not installed to
CMAKE_INSTALL_PREFIX, one needs to specify the locations by supplying the
variables INC_PATH and LIB_PATH in the cmake arguments as follows:

⬇

cmake -DCMAKE_INSTALL_PREFIX=path -DINC_PATH="inc1;inc2"
-DLIB_PATH="lib1;lib2" ..

The most recent version of the detailed document is available on the website
https://heplib.github.io/doxygen/, and to view the html document on the local
machine, one can run the shell command doxygen in the directory HepLib/src/doc
and open html/index.html in a browser.

It is not required but will be convenient for later use to adjust the related
environment variable by appending the following command to one’s .bashrc:

⬇

export PATH=<INSTALL PATH>/bin:$PATH

###  3.3 Basic usage within C++

The basic usage of HepLib is similar to other C++ library by including the
proper header files in the C++ source code, compiling the program and linking
with HepLib and other necessary libraries. The following code shows how to
perform the $D$-dimensional trace on a Dirac-$\gamma$ chain ${\rm
Tr}[p\\!\\!\\!/_{1}\gamma^{\mu}(p\\!\\!\\!/_{2}+m)\gamma_{\mu}]$ and print the
result at the end:

⬇

#include "HepLib.h"

using namespace HepLib;

int main(int argc, char** argv) {

Index mu("mu"), nu("nu");

Vector p1("p1"), p2("p2");

Symbol m("m");

ex gline = GAS(p1)*GAS(mu)*(GAS(p2)+m*GAS(1))*GAS(mu);

ex trace = form(TR(gline));

hout << trace << endl;

return 0;

}

Assuming the file is named trace.cpp777trace.cpp can be found in the directory
codes in HepLib archive, one can compile and run it as follows:

⬇

$ <INSTALL PATH>/bin/heplib++ -o trace trace.cpp

$ ./trace

8*p2.p1+(-4)*D*p2.p1

where heplib++ is a shell script installed to <INSTALL PATH>/bin, heplib++
encodes the g++ flags including -I, -L and -l flags. One can also compile the
program with pkg-config as follows:

⬇

export PKG_CONFIG_PATH=<INSTALL PATH>/lib/pkgconfig:$PKG_CONFIG_PATH

g++ $(pkg-config \--cflags \--libs HepLib) -o trace trace.cpp

or provide the g++ flags explicitly as follows:

⬇

g++ -I <INSTALL PATH>/include -L <INSTALL PATH>/lib -Wl,-rpath,<INSTALL
PATH>/lib -lHepLib -lginac -o trace trace.cpp

###  3.4 Examples

To facilitate the usage of HepLib library, a list of examples are also shipped
along with HepLib.tar.gz, here we provide a few descriptions for helping to
use the examples for the current version.

  * 1.

The files in examples/Basic provide some basic usages of the functions
extended to GiNaC itself, including the interfaces to Fermat and FORM
programs.

  * 2.

The files in examples/Hep demonstrate the basic objects in High Energy
Physics, including a complete example using qgraf, usage of form to evaluate
Dirac-$\gamma$ or SUNT traces, TIR for tensor index reduction, IBP reductions
with FIRE and so on.

  * 3.

The files in examples/SD illustrate basic usage of sector decomposition
method, including numerical evaluation of loop integrals or general parametric
integration using SecDec class.

More examples to be added in the future will be available at
https://heplib.github.io/example.html.

##  4 Summary and outlook

A C++ library HepLib is presented for computations in High Energy Physics, it
uses GiNaC for the symbolic operations. The library tries to provide an
integrated C++ interface to several well-established programs, including qgraf
for generating Feynman diagrams and amplitudes, Fermat for fast
simplifications on rational polynomials, FORM for highly efficient
calculations of Dirac-$\gamma$ and Color matrix. The geometric sector
decomposition method has also been implemented in the native GiNaC language
with the help of qhull library with many new features. In the future version,
we are going to introduce the namespace HepLib::DE for the related operations
using DE method, the namespace HepLib::MB for the related operations using MB
method, other IBP reduction programs (like Reduze) to HepLib::IBP, interfaces
to other programming languages (including Mathematica, Python, etc.) and many
other features.

## Acknowledgment

The authors are grateful to the developers of GiNaC for answering questions
and providing supports on the GiNaC mailing list [108]. This work is supported
by the National Natural Science Foundation of China under Grants No. 11875318,
and by the Yue Qi Young Scholar Project in China University of Mining and
Technology (Beijing).

## References

  * [1] S. Amoroso, et al. [arXiv:2003.01700 [hep-ph]].
  * [2] T. Hahn, Comput. Phys. Commun. 140, 418-431 (2001) doi:10.1016/S0010-4655(01)00290-9 [arXiv:hep-ph/0012260 [hep-ph]].
  * [3] FeynArts website: http://www.feynarts.de.
  * [4] P. Nogueira, J. Comput. Phys. 105, 279-289 (1993) doi:10.1006/jcph.1993.1074
  * [5] qgraf website: http://cfif.ist.utl.pt/~paulo/qgraf.html.
  * [6] R. Mertig, M. Bohm and A. Denner, Comput. Phys. Commun. 64, 345-359 (1991) doi:10.1016/0010-4655(91)90130-D
  * [7] V. Shtabovenko, R. Mertig and F. Orellana, Comput. Phys. Commun. 207, 432-444 (2016) doi:10.1016/j.cpc.2016.06.008 [arXiv:1601.01167 [hep-ph]].
  * [8] FeynCalc website: https://feyncalc.github.io.
  * [9] J. A. M. Vermaseren, [arXiv:math-ph/0010025 [math-ph]].
  * [10] FORM website: https://www.nikhef.nl/~form/.
  * [11] F. Feng and R. Mertig, [arXiv:1212.3522 [hep-ph]].
  * [12] FeynCalc/FormLink website: https://github.com/FormLink.
  * [13] F. Feng, Comput. Phys. Commun. 183, 2158-2164 (2012) doi:10.1016/j.cpc.2012.03.025 [arXiv:1204.2314 [hep-ph]].
  * [14] APart website: https://github.com/F-Feng/APart.
  * [15] C. Anastasiou and A. Lazopoulos, JHEP 07, 046 (2004) doi:10.1088/1126-6708/2004/07/046 [arXiv:hep-ph/0404258 [hep-ph]].
  * [16] AIR website: https://people.phys.ethz.ch/~pheno/air/.
  * [17] A. V. Smirnov, JHEP 10, 107 (2008) doi:10.1088/1126-6708/2008/10/107 [arXiv:0807.3243 [hep-ph]].
  * [18] A. V. Smirnov and V. A. Smirnov, Comput. Phys. Commun. 184, 2820-2827 (2013) doi:10.1016/j.cpc.2013.06.016 [arXiv:1302.5885 [hep-ph]].
  * [19] A. V. Smirnov, Comput. Phys. Commun. 189, 182-191 (2015) doi:10.1016/j.cpc.2014.11.024 [arXiv:1408.2372 [hep-ph]].
  * [20] A. V. Smirnov and F. S. Chuharev, doi:10.1016/j.cpc.2019.106877 [arXiv:1901.07808 [hep-ph]].
  * [21] FIRE website: https://bitbucket.org/feynmanIntegrals/fire/.
  * [22] R. N. Lee, [arXiv:1212.2685 [hep-ph]].
  * [23] R. N. Lee, J. Phys. Conf. Ser. 523, 012059 (2014) doi:10.1088/1742-6596/523/1/012059 [arXiv:1310.1145 [hep-ph]].
  * [24] LiteRed website: http://www.inp.nsk.su/~lee/programs/LiteRed/.
  * [25] C. Studerus, Comput. Phys. Commun. 181, 1293-1300 (2010) doi:10.1016/j.cpc.2010.03.012 [arXiv:0912.2546 [physics.comp-ph]].
  * [26] A. von Manteuffel and C. Studerus, [arXiv:1201.4330 [hep-ph]].
  * [27] Reduze website: https://reduze.hepforge.org.
  * [28] P. Maierhöfer, J. Usovitsch and P. Uwer, Comput. Phys. Commun. 230, 99-112 (2018) doi:10.1016/j.cpc.2018.04.012 [arXiv:1705.05610 [hep-ph]].
  * [29] J. Klappert, F. Lange, P. Maierhöfer and J. Usovitsch, [arXiv:2008.06494 [hep-ph]].
  * [30] KIRA website: https://kira.hepforge.org.
  * [31] V. A. Smirnov, Phys. Lett. B 460, 397-404 (1999) doi:10.1016/S0370-2693(99)00777-7 [arXiv:hep-ph/9905323 [hep-ph]].
  * [32] J. B. Tausk, Phys. Lett. B 469, 225-234 (1999) doi:10.1016/S0370-2693(99)01277-0 [arXiv:hep-ph/9909506 [hep-ph]].
  * [33] MBTools website: https://mbtools.hepforge.org.
  * [34] T. Binoth and G. Heinrich, Nucl. Phys. B 680, 375-388 (2004) doi:10.1016/j.nuclphysb.2003.12.023 [arXiv:hep-ph/0305234 [hep-ph]].
  * [35] G. Heinrich, Int. J. Mod. Phys. A 23, 1457 (2008) doi:10.1142/S0217751X08040263 [arXiv:0803.4177 [hep-ph]].
  * [36] I. Dubovyk, J. Gluza, T. Riemann and J. Usovitsch, PoS LL2016, 034 (2016) doi:10.22323/1.260.0034 [arXiv:1607.07538 [hep-ph]].
  * [37] J. Blümlein, I. Dubovyk, J. Gluza, M. Ochman, C. G. Raab, T. Riemann and C. Schneider, PoS LL2014, 052 (2014) doi:10.22323/1.211.0052 [arXiv:1407.7832 [hep-ph]].
  * [38] AMBRE website: https://prac.us.edu.pl/~gluza/ambre/.
  * [39] A. V. Smirnov and M. N. Tentyukov, Comput. Phys. Commun. 180, 735 (2009) doi:10.1016/j.cpc.2008.11.006 [arXiv:0807.4129 [hep-ph]].
  * [40] A. V. Smirnov, V. A. Smirnov and M. Tentyukov, Comput. Phys. Commun. 182, 790 (2011) doi:10.1016/j.cpc.2010.11.025 [arXiv:0912.0158 [hep-ph]].
  * [41] A. V. Smirnov, Comput. Phys. Commun. 185, 2090 (2014) doi:10.1016/j.cpc.2014.03.015 [arXiv:1312.3186 [hep-ph]].
  * [42] A. V. Smirnov, Comput. Phys. Commun. 204, 189 (2016) doi:10.1016/j.cpc.2016.03.013 [arXiv:1511.03614 [hep-ph]].
  * [43] FIESTA website: https://bitbucket.org/feynmanIntegrals/fiesta/.
  * [44] S. Borowka, J. Carter and G. Heinrich, Comput. Phys. Commun. 184, 396-408 (2013) doi:10.1016/j.cpc.2012.09.020 [arXiv:1204.4152 [hep-ph]].
  * [45] S. Borowka and G. Heinrich, Comput. Phys. Commun. 184, 2552-2561 (2013) doi:10.1016/j.cpc.2013.05.022 [arXiv:1303.1157 [hep-ph]].
  * [46] S. Borowka, G. Heinrich, S. P. Jones, M. Kerner, J. Schlenk and T. Zirke, Comput. Phys. Commun. 196, 470-491 (2015) doi:10.1016/j.cpc.2015.05.022 [arXiv:1502.06595 [hep-ph]].
  * [47] SecDec website: https://secdec.hepforge.org.
  * [48] S. Borowka, G. Heinrich, S. Jahn, S. P. Jones, M. Kerner, J. Schlenk and T. Zirke, Comput. Phys. Commun. 222, 313-326 (2018) doi:10.1016/j.cpc.2017.09.015 [arXiv:1703.09692 [hep-ph]].
  * [49] T. Hahn, Comput. Phys. Commun. 168, 78-95 (2005) doi:10.1016/j.cpc.2005.01.010 [arXiv:hep-ph/0404043 [hep-ph]].
  * [50] T. Hahn, J. Phys. Conf. Ser. 608, no.1, 012066 (2015) doi:10.1088/1742-6596/608/1/012066 [arXiv:1408.6373 [physics.comp-ph]].
  * [51] CUBA website: http://www.feynarts.de/cuba/.
  * [52] Z. Li, J. Wang, Q. S. Yan and X. Zhao, Chin. Phys. C 40, no.3, 033103 (2016) doi:10.1088/1674-1137/40/3/033103 [arXiv:1508.02512 [hep-ph]].
  * [53] C. Bogner and S. Weinzierl, Comput. Phys. Commun. 178, 596-610 (2008) doi:10.1016/j.cpc.2007.11.012 [arXiv:0709.4092 [hep-ph]].
  * [54] sector_decomposition website: http://wwwthep.physik.uni-mainz.de/stefanw/~sector_decomposition/.
  * [55] C. W. Bauer, A. Frink and R. Kreckel, J. Symb. Comput. 33, 1-12 (2002) doi:10.1006/jsco.2001.0494 [arXiv:cs/0004015 [cs.SC]].
  * [56] GiNaC website: https://ginac.de/.
  * [57] J. Gluza, K. Kajda, T. Riemann and V. Yundin, Eur. Phys. J. C 71, 1516 (2011) doi:10.1140/epjc/s10052-010-1516-y [arXiv:1010.1667 [hep-ph]].
  * [58] CSector website: http://prac.us.edu.pl/~gluza/csectors/.
  * [59] A. V. Kotikov, Phys. Lett. B 267, 123-127 (1991) [erratum: Phys. Lett. B 295, 409-409 (1992)] doi:10.1016/0370-2693(91)90536-Y
  * [60] E. Remiddi, Nuovo Cim. A 110, 1435-1452 (1997) [arXiv:hep-th/9711188 [hep-th]].
  * [61] J. M. Henn, Phys. Rev. Lett. 110, 251601 (2013) doi:10.1103/PhysRevLett.110.251601 [arXiv:1304.1806 [hep-th]].
  * [62] J. M. Henn, J. Phys. A 48, 153001 (2015) doi:10.1088/1751-8113/48/15/153001 [arXiv:1412.2296 [hep-ph]].
  * [63] O. Gituliar and V. Magerya, Comput. Phys. Commun. 219, 329-338 (2017) doi:10.1016/j.cpc.2017.05.004 [arXiv:1701.04269 [hep-ph]].
  * [64] Fuchsia website: https://github.com/gituliar/fuchsia.
  * [65] M. Prausa, Comput. Phys. Commun. 219, 361-376 (2017) doi:10.1016/j.cpc.2017.05.026 [arXiv:1701.00725 [hep-ph]].
  * [66] epsilon website: https://github.com/mprausa/epsilon.
  * [67] C. Meyer, Comput. Phys. Commun. 222, 295-312 (2018) doi:10.1016/j.cpc.2017.09.014 [arXiv:1705.06252 [hep-ph]].
  * [68] CANONICA website: https://github.com/christophmeyer/CANONICA.
  * [69] R. N. Lee, A. V. Smirnov and V. A. Smirnov, JHEP 03, 008 (2018) doi:10.1007/JHEP03(2018)008 [arXiv:1709.07525 [hep-ph]].
  * [70] X. Liu, Y. Q. Ma and C. Y. Wang, Phys. Lett. B 779, 353-357 (2018) doi:10.1016/j.physletb.2018.02.026 [arXiv:1711.09572 [hep-ph]].
  * [71] DESS website: https://bitbucket.org/feynmanIntegrals/fire/.
  * [72] T. Hahn and M. Perez-Victoria, Comput. Phys. Commun. 118, 153-165 (1999) doi:10.1016/S0010-4655(98)00173-8 [arXiv:hep-ph/9807565 [hep-ph]].
  * [73] FormCalc website: http://www.feynarts.de/formcalc/.
  * [74] A. Pukhov, E. Boos, M. Dubinin, V. Edneral, V. Ilyin, D. Kovalenko, A. Kryukov, V. Savrin, S. Shichanin and A. Semenov, [arXiv:hep-ph/9908288 [hep-ph]].
  * [75] E. Boos et al. [CompHEP], Nucl. Instrum. Meth. A 534, 250-259 (2004) doi:10.1016/j.nima.2004.07.096 [arXiv:hep-ph/0403113 [hep-ph]].
  * [76] CompHEP website: http://comphep.sinp.msu.ru.
  * [77] J. X. Wang, Nucl. Instrum. Meth. A 534, 241-245 (2004) doi:10.1016/j.nima.2004.07.094 [arXiv:hep-ph/0407058 [hep-ph]].
  * [78] A. Belyaev, N. D. Christensen and A. Pukhov, Comput. Phys. Commun. 184, 1729-1769 (2013) doi:10.1016/j.cpc.2013.01.014 [arXiv:1207.6082 [hep-ph]].
  * [79] CalcHEP website: https://theory.sinp.msu.ru/~pukhov/calchep.html.
  * [80] G. Cullen, N. Greiner, G. Heinrich, G. Luisoni, P. Mastrolia, G. Ossola, T. Reiter and F. Tramontano, Eur. Phys. J. C 72, 1889 (2012) doi:10.1140/epjc/s10052-012-1889-1 [arXiv:1111.2034 [hep-ph]].
  * [81] G. Cullen, H. van Deurzen, N. Greiner, G. Heinrich, G. Luisoni, P. Mastrolia, E. Mirabella, G. Ossola, T. Peraro, J. Schlenk, J. F. von Soden-Fraunhofen and F. Tramontano, Eur. Phys. J. C 74, no.8, 3001 (2014) doi:10.1140/epjc/s10052-014-3001-5 [arXiv:1404.7096 [hep-ph]].
  * [82] GoSam website: https://gosam.hepforge.org.
  * [83] DIANA website: https://github.com/apik/diana.
  * [84] J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H. S. Shao, T. Stelzer, P. Torrielli and M. Zaro, JHEP 07, 079 (2014) doi:10.1007/JHEP07(2014)079 [arXiv:1405.0301 [hep-ph]].
  * [85] M. Bahr, S. Gieseke, M. A. Gigg, D. Grellscheid, K. Hamilton, O. Latunde-Dada, S. Platzer, P. Richardson, M. H. Seymour, A. Sherstnev and B. R. Webber, Eur. Phys. J. C 58, 639-707 (2008) doi:10.1140/epjc/s10052-008-0798-9 [arXiv:0803.0883 [hep-ph]].
  * [86] T. Gleisberg, S. Hoeche, F. Krauss, M. Schonherr, S. Schumann, F. Siegert and J. Winter, JHEP 02, 007 (2009) doi:10.1088/1126-6708/2009/02/007 [arXiv:0811.4622 [hep-ph]].
  * [87] W. Kilian, T. Ohl and J. Reuter, Eur. Phys. J. C 71, 1742 (2011) doi:10.1140/epjc/s10052-011-1742-y [arXiv:0708.4233 [hep-ph]].
  * [88] hepforge projects website: https://www.hepforge.org/projects.
  * [89] GiNaC official tutorial: https://ginac.de/tutorial/.
  * [90] Fermat website: http://home.bway.net/lewis/.
  * [91] T. van Ritbergen, A. N. Schellekens and J. A. M. Vermaseren: ”Group theory factors for Feynman diagrams”, Int. J. Mod. Phys. vol. 14 (1999) 41.
  * [92] FORM color package website: https://www.nikhef.nl/~form/maindir/packages/color/color.html.
  * [93] P. Cvitanovic, Phys. Rev. D 14, 1536-1553 (1976) doi:10.1103/PhysRevD.14.1536
  * [94] Michael E. Peskin and Daniel V. Schroeder. _An introduction to Quantum Field Theory_ , Published 2018 by CRC Press.
  * [95] J. Ellis, Comput. Phys. Commun. 210, 103-123 (2017) doi:10.1016/j.cpc.2016.08.019 [arXiv:1601.05437 [hep-ph]].
  * [96] Tikz-Feynman website: https://jpellis.me/projects/tikz-feynman/.
  * [97] HepLib document: https://heplib.github.io/doxygen/.
  * [98] A. Czarnecki and K. Melnikov, Phys. Rev. Lett. 80, 2531-2534 (1998) doi:10.1103/PhysRevLett.80.2531 [arXiv:hep-ph/9712222 [hep-ph]].
  * [99] M. Beneke, A. Signer and V. A. Smirnov, Phys. Rev. Lett. 80, 2535-2538 (1998) doi:10.1103/PhysRevLett.80.2535 [arXiv:hep-ph/9712302 [hep-ph]].
  * [100] G. T. Bodwin and A. Petrelli, Phys. Rev. D 66, 094011 (2002) [erratum: Phys. Rev. D 87, no.3, 039902 (2013)] doi:10.1103/PhysRevD.66.094011 [arXiv:hep-ph/0205210 [hep-ph]].
  * [101] G. T. Bodwin, X. Garcia Tormo, i and J. Lee, Phys. Rev. D 81, 114014 (2010) doi:10.1103/PhysRevD.81.114014 [arXiv:1003.0061 [hep-ph]].
  * [102] P. Bärnreuther, M. Czakon and P. Fiedler, JHEP 02, 078 (2014) doi:10.1007/JHEP02(2014)078 [arXiv:1312.6279 [hep-ph]].
  * [103] Cubature website: https://github.com/stevengj/cubature.
  * [104] H. Cheng, T. Wu, Expanding Protons: Scattering at High Energies. The MIT Press (1987).
  * [105] V. A. Smirnov, Feynman integral calculus. Springer, 2006.
  * [106] V. Pilipp, JHEP 09, 135 (2008) doi:10.1088/1126-6708/2008/09/135 [arXiv:0808.2555 [hep-ph]].
  * [107] MPFR C++ website: http://www.holoborodko.com/pavel/mpfr/.
  * [108] The GiNaC-list Archives: https://www.ginac.de/pipermail/ginac-list/.