## MPL - Supplementary User Manual

### Christian Bogner

Version 1.0

### Contents

1 Getting started 2

2 Computing with cubical coordinates 3

2.1 Differential 1-forms and iterated integrals . . . 3

2.2 Differentiation . . . 7

2.3 The symbol map and the unshuffle map . . . 8

2.4 Integration over cubical coordinates . . . 12

3 Computing Feynman integrals 16 3.1 Setting up a simple example . . . 17

3.2 Checking the conditions . . . 18

3.3 Integration over Feynman parameters . . . 22

3.4 Expressing hyperlogarithms in bar-notation . . . 23

4 Automated tests and troubleshooting checklist 24

### 1 Getting started

MPL is a collection of procedures for the computeralgebra system Maple. It was written and tested with Maple 16. It supports a wide range of computations with homotopy invariant iterated integrals on moduli spaces

### M

_{0,n}of curves of genus 0 with n ordered marked points. This class of functions contains the widely used multiple polylogarithms as functions of several variables.

This manual shows how to apply the most important procedures of the program and might be
useful for looking up technical details or examples. It is just meant to supplement the detailed
introduction given in the article [2]. The main algorithms were presented in [4]. They are based
on the mathematical framework elaborated in [6].^{1} When using MPL for a scientific publication,
please cite (some of) these articles.

The latest version of MPL and this manual are available from the webpage http://cbogner.com/software/mpl/

The entire program is obtained by downloading a txt-fileMPLn_m.txt, where the integersnandm indicate the number of the version. For example, the file of MPL version 1.0 is calledMPL1_0.txt and after saving the file in the same directory with your Maple worksheet, the program is started with

1A very useful introduction to the general concept of iterated integrals and generalizations of polylogarithms can be found in [7]. For an overview of the far-reaching applications of such classes of iterated integrals in particle physics, we recommend [13, 14].

>read("MPL1_0.txt"):

Before starting with any computations, it is useful to let Maple read a file by which all multiple zeta values up to a certain weight are automatically expressed in terms of a basis. While writing MPL, we have been using the filemzv-1-12.txtprovided by Bigotte, Jacob, Oussous, Minh and Petitot for this purpose. The file is available from the webpage

http://www.lifl.fr/~petitot/recherche/MZV/mzv12/

After saving the file in the working directory and calling

>read("mzv-1-12.txt"):

every multiple zeta value up to weight 12 appearing from now on will be expressed in terms of a basis, for example:

>zeta(3,2);

−^{11}_{2}ζ(5) +3ζ(2)ζ(3)

This will drastically simplify expressions and it is also necessary for automated tests discussed in section 4 to check that certain parts of the program produce the expected results.

As a supplement to this manual, the above webpage also provides a maple worksheet which contains all examples shown below.

### 2 Computing with cubical coordinates

### 2.1 Differential 1-forms and iterated integrals

In this section, we work with homotopy invariant iterated integrals on

### M

_{0,n}with n>3. Here ho- motopy invariance implies, that these functions only depend on the end-points of some path. We always consider the origin as one of these end-points. The other point has m=n−3 coordinates, which the iterated integral depends on as a function. Due to additional conditions with respect to a basepoint, it is completely determined by the differential 1-forms and the order in which one inte- grates over them. Therefore it is customary to express such an iterated integral by a tensor product of the involved 1-forms a

_{k}⊗a

_{k−1}⊗...⊗a

_{1}which is often replaced by a notation with squared brackets[ak|ak−1|...|a1]called bar-notation.

**bar(...)**

The bar-notation is adapted in MPL. We use the reserved wordbarto express every iterated inte- gral. For example[a3|a2|a1]is expressed by

>bar(a[3],a[2],a[1]);

and stands for an iterated integral where we integrate over some 1-forms a_{1},a_{2},a_{3} in this order-
ing (i.e. from right to left in the argument of bar). More generally, iterated integrals are linear
combinations of such expressions overQ,for example

>3*bar(a[2],a[1])+5*bar(b[8])+9;

Numerical constants are factored out automatically, for example:

>bar(3*a,2*b)-bar(a-2*b,c);

6bar(a,b)−bar(a,c) +2bar(b,c)

>bar(0,c);

0

**BarToLists(expr), ListsToBar(list)**

For some computations, the user may find it convenient to bring an expression of bar-terms in the form of a list. This form is also used internally in some algorithms. For exprbeing a linear combination of words in the bar-notation, the auxiliary procedureBarToLists(expr)returns a list of two minor lists. The first of these lists contains the coefficients of the words and the second list contains the words in the argument ofbar(...), as lists. For example

>f:=bar(a)+bar(b,c)+99*bar(b,c)-y*bar(5*g,c,b):

>BarToLists(f);

[[1,100,−5y],[[a],[b,c],[g,c,b]]]

The commandListsToBar(list)reverses this transformation. It takes such a list of coefficients and words as input and returns the corresponding expression in bar-notation. For example, as continuation of the above:

>ListsToBar(%);

bar(a) +100bar(b,c)−5ybar(g,c,b)
**MPLShuffleProduct(expr1, expr2)**

The multiplication of two iterated integrals expr1 and expr2 is computed by the well-known shuffle-product, implemented inMPLShuffleProduct. For example:

>MPLShuffleProduct(bar(a,b)+3*bar(c)+7, bar(y,z));

3bar(c,y,z)+3bar(y,c,z)+3bar(y,z,c)+bar(a,b,y,z)+bar(a,y,b,z)+bar(a,y,z,b)+

bar(y,a,b,z) +bar(y,a,z,b) +bar(y,z,a,b) +7bar(y,z)

**MPLCoproduct(expr)**

The deconcatenation co-product∆, defined by

∆[ak|ak−1|...|a1] =1⊗[ak|ak−1|...|a1] + [ak]⊗[ak−1|...|a1] +...+ [ak|ak−1|...|a1]⊗1 is implemented in the procedureMPLCoproduct(expr), for example:

>MPLCoproduct(bar(a,b,c))

tens(ONE,bar(a,b,c))+tens(bar(a),bar(b,c))+tens(bar(a,b),bar(c))+tens(bar(a,b,c),ONE) Here the reserved termONEstands for the empty word andtens(...,...) is a formal expression for the tensor-product of two iterated integrals. tens andONEare mostly defined for internal use and will usually not appear in any results below.

**MPLCoordinates(letter, n), MPL_COORDINATES**

As a preparation for many applications, it is necessary to declare a set of variables to be the cubical coordinates we want to work with. For a positive integern, the procedureMPLCoordinates(letter, n) declaresncubical coordinates by use of the first argument letterfor the names of the vari- ables. The list of cubical coordinates defined in this way is stored in the global variableMPL_COORDINATES until the procedure is called again. For example:

>MPLCoordinates(x, 5):

>MPL_COORDINATES;

[x1,x_{2},x_{3},x_{4},x_{5}]

After calling MPLCoordinates, the variables declared in this way are internally recognized to be cubical coordinates and they can be used to construct the corresponding differential 1-forms. For example one has relations liked(x1)∧d(x1) =0 :

>d(x[1]) &^ d(x[1]);

0

**MPLFormsFiber(lifted::boolean), MPLFormsBase(), MPLFormsTotal()**
After declaring x_{1}, ...,x_{m} as cubical coordinates, MPL can compute with iterated integrals with
differential 1-forms of the following sets (as defined in [4]):

Ω_{m} =

(dx_{1}

x_{1} , ...,dx_{m}

x_{m} ,d ∏_{a≤i≤b}xi

∏_{a≤i≤b}x_{i}−1 for 1≤a≤b≤m
)

,

Ω¯^{F}_{m} =

(dx_{m}
xm

, ∏_{a≤i≤n−1}x_{i}
dx_{m}

∏_{a≤i≤m}xi−1 for 1≤a≤m
)

,

Ω^{F}_{m} =

(dx_{m}

x_{m} ,d ∏_{a≤i≤m}x_{i}

∏_{a≤i≤m}x_{i}−1 for 1≤a≤m
)

(1)

where clearlyΩ_{m}=Ω^{F}_{m}∪Ω_{m−1}.In the context discussed in [4, 6] and with respect to a givenn,it
makes sense to speak of these forms as:

• Ω_{m}: forms on the total space,

• Ω¯^{F}_{m}: forms on the fiber,

• Ω^{F}_{m}: forms on the lifted fiber,

• Ω_{m−1}: forms on the base.

Up to signs, these sets can be obtained from the auxiliary proceduresMPLFormsTotal(),

MPLFormsFiber(false),MPLFormsFiber(true),MPLFormsBase()respectively. For example:

>MPLCoordinates(x, 3):

>MPLFormsTotal();

h_{d(x}

1)
x1 ,^{d(x}_{1−x}^{1}^{)}

1,^{d(x}_{x}^{2}^{)}

2 ,^{x}^{2}^{d(x}_{1−x}^{1}^{)+x}^{1}^{d(x}^{2}^{)}

1x2 ,^{d(x}_{1−x}^{2}^{)}

2,^{d(x}_{x}^{3}^{)}

3 ,^{x}^{2}^{x}^{3}^{d(x}^{1}^{)+x}_{1−x}^{1}^{x}^{3}^{d(x}^{2}^{)+x}^{1}^{x}^{2}^{d(x}^{3}^{)}

1x2x3 ,^{x}^{3}^{d(x}_{1−x}^{2}^{)+x}^{2}^{d(x}^{3}^{)}

2x3 ,^{d(x}_{1−x}^{3}^{)}

3

i

>MPLFormsFiber(false);

h_{d(x}

3)

x_{3} ,^{x}_{1−x}^{1}^{x}^{2}^{d(x}^{3}^{)}

1x_{2}x_{3},^{x}_{1−x}^{2}^{d(x}^{3}^{)}

2x_{3},^{d(x}_{1−x}^{3}^{)}

3

i

>MPLFormsFiber(true);

h_{d(x}

3)

x_{3} ,^{x}^{2}^{x}^{3}^{d(x}^{1}^{)+x}_{1−x}^{1}^{x}^{3}^{d(x}^{2}^{)+x}^{1}^{x}^{2}^{d(x}^{3}^{)}

1x_{2}x_{3} ,^{x}^{3}^{d(x}_{1−x}^{2}^{)+x}^{2}^{d(x}^{3}^{)}

2x_{3} ,^{d(x}_{1−x}^{3}^{)}

3

i

>MPLFormsBase();

h_{d(x}

1)

x1 ,^{d(x}_{1−x}^{1}_{1}^{)},^{d(x}_{x}_{2}^{2}^{)},^{x}^{2}^{d(x}_{1−x}^{1}^{)+x}_{1}_{x}^{1}_{2}^{d(x}^{2}^{)},^{d(x}_{1−x}^{2}_{2}^{)}i

In the following we will denote the Q-vectorspaces of differential 1-forms, spanned by the bases
Ω_{m},Ω¯^{F}_{m},Ω^{F}_{m}by

### A

_{m},

### A

¯_{m}

^{F},

### A

_{m}

^{F}respectively. Furthermore, the corresponding Q-vectorspaces of ho- motopy invariant iterated integrals are denoted byV(Ωm),V Ω¯

^{F}

_{m}

,V Ω^{F}_{m}

(see [4, 2]).

**MPLArnoldEquation(form1,form2)**

The above differential 1-forms satisfy the following set of relations due to Arnol’d [1] (also see [4]):

dxm

x_{m} ∧d(xi...xm)

x_{i}...xm−1 = −

m−1

### ∑

k=i

dx_{k}

x_{k} ∧d(xi...xm)
x_{i}...xm−1,
d xj...xm

x_{j}...xm−1∧d(xi...xm)

x_{i}...xm−1 = d xi...xj−1

x_{i}...xj−1−1∧ d(xi...xm)

x_{i}...xm−1−d xj...xm
x_{j}...xm−1

!

−

j−1 k=i

### ∑

dx_{k}

x_{k} ∧d(xi...xm)
x_{i}...xm−1

for 1≤i≤ j≤m.Note that these equations are of the type
ω_{i}∧ω_{j}=

### ∑

k

α_{k}∧ω_{k}

where the 1-formsω_{i}are inΩ^{F}_{m}(the lifted fiber) and the 1-formsα_{i}are inΩ_{m−1}(the base).

Internally, these relations are extensively used in the computations below. More conveniently, for two given forms form1, form2 on the left-hand side of the above equations (in the sense form1∧form2) the right-hand side of the Arnol’d equation can be obtained from

MPLArnoldEquation(form1,form2). For example forn=3, the right-hand side of the equation
dx_{3}

x_{3}−1∧d(x2x_{3})

x_{2}x_{3}−1= dx_{2}
x_{2}−1∧

d(x2x_{3})

x_{2}x_{3}−1− dx_{3}
x_{3}−1

−dx_{2}

x_{2} ∧d(x2x_{3})
x_{2}x_{3}−1
is computed as^{2}:

>MPLCoordinates(x, 3):

>MPLArnoldEquation(d(x[3])/(x[3]-1), d(x[2]*x[3])/(x[2]*x[3]-1));

−^{(d(x}_{(1−x}^{2}^{))&}^{∧}^{(d(x}^{3}^{))}

2)(1−x3) +^{x}^{2}_{(1−x}^{(d(x}^{2}^{))&}^{∧}^{(d(x}^{3}^{))}

2)(1−x2x3) +^{(d(x}^{2}_{(1−x}^{))&}^{∧}^{(d(x}^{3}^{))}

2x3)

### 2.2 Differentiation

**MPLd(expr)**

Differentiation d with respect to the end-point of the path (w.r.t. which the iterated integral is defined) is simply given by the de-concatenation of the leftmost differential 1-form:

d:V(Ωm) →

### A

_{m}⊗V(Ωm),

### ∑

J=(i1, ...,ik)

c_{I}[ωi1|...|ωik] 7→

### ∑

J=(i1, ...,ik)

c_{I}ω_{i}_{1}⊗[ωi2|...|ωik].

For example, one computes d

d(x2x_{3})
1−x2x3

|dx_{3}
x3

−
dx_{2}

x2

|d(x2x_{3})
1−x2x3

= d(x2x_{3})
1−x2x3

⊗
dx_{3}

x3

−dx_{2}
x2

⊗

d(x2x_{3})
1−x2x3

by

2In principle, using the auxiliary procedureGetArnoldRHS(...), these equations can also be obtained in a different output format, where forms of the fiber and the base are explicitely separated. Such a format is used internally.

>MPLCoordinates(x, 3):

>MPLd(bar((x[3]*d(x[2])+x[2]*d(x[3]))/(1-x[2]*x[3]),d(x[3])/x[3]) -bar(d(x[2])/x[2],(x[3]*d(x[2])+x[2]*d(x[3]))/(1-x[2]*x[3])));

tens_{x}

3d(x_{2})+x_{2}d(x_{3})
1−x2x_{3} ,bar

dx_{3}
x_{3}

−tens

dx_{2}

x_{2} ,bar_{x}

3d(x_{2})+x_{2}d(x_{3})
1−x2x_{3}

**MPLTotalConnection(expr)**

In [4] we have defined a connection

∇:V Ω¯^{F}_{m}

→

### A

_{m−1}⊗V Ω¯

^{F}

_{m}by use of the Arnol’d equations and a total connection

∇_{T} :V Ω¯^{F}_{m}

→

### A

_{m}

_{⊗}V Ω¯

^{F}

_{m}by

∇_{T} =d−∇.

For a word in ¯Ω^{F}_{m},the total connection can be computed with the commandMPLTotalConnection.

For example, we obtain

∇_{T}

x_{2}d(x3)
1−x_{2}x_{3}|dx_{3}

x_{3}

= d(x2x_{3})
1−x_{2}x_{3}⊗

dx_{3}
x_{3}

−dx_{2}
x_{2} ⊗

x_{2}d(x3)
1−x_{2}x_{3}

by

>MPLCoordinates(x, 3):

>MPLTotalConnection(bar((x[2]*d(x[3])/(1-x[2]*x[3]),d(x[3])/x[3]));

tens_{x}

3d(x_{2})+x_{2}d(x_{3})
1−x_{2}x_{3} ,bar

dx_{3}
x_{3}

−tens

dx_{2}

x_{2} ,bar_{x}

2d(x_{3})
1−x_{2}x_{3}

### 2.3 The symbol map and the unshuffle map

**MPLSymbolMap(expr)**

It is implied by a theorem of Chen [11] that an iterated integral with differential 1-forms in our framework is homotopy invariant if and only if the corresponding word of 1-forms satisfies a so- called integrability condition. These words are called integrable. Every word in ¯

### A

_{m}

^{F}is integrable, but not every word in

### A

_{m}.The symbol map, as constructed in [3, 4], maps from integrable words in ¯

### A

_{m}

^{F}to integrable words in

### A

_{m}:

Ψ:V Ω¯^{F}_{m}

→V(Ωm).

For example we compute Ψ

x2d(x3)
1−x_{2}x_{3}|dx3

x_{3}

=

d(x2x3)
1−x_{2}x_{3}|dx3

x_{3}

−

d(x2)

x_{2} |d(x2x3)
1−x_{2}x_{3}

(2) by

>MPLCoordinates(x, 3):

>MPLSymbolMap(bar((x[2]*d(x[3]))/(1-x[2]*x[3]),d(x[3])/x[3]));

bar

x3d(x2)+x2d(x3)
1−x2x3 ,^{dx}_{x}^{3}

3

−bar
_{d(x}

2)

x2 ,^{x}^{3}^{d(x}_{1−x}^{2}^{)+x}^{2}^{d(x}^{3}^{)}

2x3

The symbol map is the unique linear map which satisfies
(id⊗Ψ)◦∇_{T} =d◦Ψ.

Let us check this equation for the above example. The left-hand side of the equation is obtained by

>f:=bar(x[2]*d(x[3])/(1-x[2]*x[3]), d(x[3])/x[3]):

>evalindets(MPLTotalConnection(f), ’specfunc(anything, tens)’, proc (h) options operator, arrow; subsop(2 = MPLSymbolMap(op(2,h)),h) end proc);

Note that this command at first computes a total connection of the functionfas a linear com- bination of terms tens(...). Then it replaces the second argument of each tens(...) by the image of this argument under the symbol map.

The right-hand side of the above equation is obtained by

>MPLd(MPLSymbolMap(f));

For both sides we obtain tens

x3d(x2)+x2d(x3) 1−x2x3 ,bar

dx3

x3

−tens

dx2

x2 ,bar

x3d(x2)+x2d(x3) 1−x2x3

Note that this was already the result of the first example in section 2.2, as we used the function in eq. 2 as input there.

**MPLBasis(letter,m,w)**

With the help of the symbol map, the commandMPLBasis(letter,m,w)constructs a basis for the
vectorspace of integrable words (or homotopy invariant iterated integrals ) up to length w in the
alphabetΩ_{m}.

For example the basis of all integrable words up to length 2 in Ω_{2} with variables y_{1},y_{2} is
obtained by

>MPLBasis(y,2,2);

bar

d(y2)
y_{2}

,bar

d(y2)
1−y_{2}

,bar

y_{2}d(y1) +y_{1}d(y2)
1−y_{1}y_{2}

,bar

d(y1)
y_{1}

,bar

d(y1)
1−y_{1}

,

bar

d(y2)

y_{2} ,d(y2)
y_{2}

,bar

d(y2)

y_{2} ,d(y2)
1−y_{2}

,bar

d(y2)

1−y_{2},d(y2)
y_{2}

, bar

d(y2)

y_{2} ,y_{2}d(y1) +y_{1}d(y2)
1−y_{1}y_{2}

+bar

d(y1)

y_{1} ,y_{2}d(y1) +y_{1}d(y2)
1−y_{1}y_{2}

, bar

y_{2}d(y1) +y_{1}d(y2)

1−y_{1}y_{2} ,d(y2)
y_{2}

−bar

d(y1)

y_{1} ,y_{2}d(y1) +y_{1}d(y2)
1−y_{1}y_{2}

,bar

d(y2)

1−y_{2},d(y2)
1−y_{2}

, bar

d(y2)

1−y_{2},y_{2}d(y1) +y_{1}d(y2)
1−y_{1}y_{2}

−bar

d(y1)

1−y_{1},y_{2}d(y1) +y_{1}d(y2)
1−y_{1}y_{2}

+bar

d(y1)

1−y_{1},d(y2)
1−y_{2}

−bar

d(y1)

y_{1} ,y_{2}d(y1) +y_{1}d(y2)
1−y_{1}y_{2}

,bar

y_{2}d(y1) +y_{1}d(y2)

1−y_{1}y_{2} ,d(y2)
1−y_{2}

+bar

d(y1)

1−y_{1},y_{2}d(y1) +y_{1}d(y2)
1−y_{1}y_{2}

−bar

d(y1)

1−y_{1},d(y2)
1−y_{2}

+bar

d(y1)

y_{1} ,y_{2}d(y1) +y_{1}d(y2)
1−y_{1}y_{2}

,bar

y_{2}d(y1) +y_{1}d(y2)

1−y_{1}y_{2} ,y_{2}d(y1) +y_{1}d(y2)
1−y_{1}y_{2}

, bar

d(y1)

y_{1} ,d(y2)
y_{2}

+bar

d(y2)

y_{2} ,d(y1)
y_{1}

,bar

d(y1)

y_{1} ,d(y2)
1−y_{2}

+bar

d(y2)

1−y_{2},d(y1)
y_{1}

,bar

d(y1)

y_{1} ,y_{2}d(y1) +y_{1}d(y2)
1−y_{1}y_{2}

+bar

y_{2}d(y1) +y_{1}d(y2)

1−y_{1}y_{2} ,d(y1)
y_{1}

, bar

d(y1)

1−y_{1},d(y2)
y_{2}

+bar

d(y2)

y_{2} ,d(y1)
1−y_{1}

,bar

d(y1)

1−y_{1},d(y2)
1−y_{2}

+bar

d(y2)

1−y_{2},d(y1)
1−y_{1}

, bar

d(y1)

1−y_{1},y_{2}d(y1) +y_{1}d(y2)
1−y_{1}y_{2}

+bar

y_{2}d(y1) +y_{1}d(y2)

1−y_{1}y_{2} ,d(y1)
1−y_{1}

,bar

d(y1)

y_{1} ,d(y1)
y_{1}

, bar

d(y1)

y_{1} ,d(y1)
1−y_{1}

,bar

d(y1)

1−y_{1},d(y1)
y_{1}

,bar

d(y1)

1−y_{1},d(y1)
1−y_{1}

The output of this command is always a list ofwlists, where fork≤wthek-th list contains the words of lengthk.

**MPLUnshuffle(expr,var)**

The unshuffle map

Φ:V(Ωm)→V(Ωm−1)⊗V Ω¯^{F}_{m}

(3) as defined in [4] is the inverse of the map

µ(id⊗Ψ):V(Ωm−1)⊗V Ω¯^{F}_{m}

→V(Ωm) (4)

and can be used recursively to decompose iterated integrals inV(Ωm)into products of hyperloga-
rithms inV Ω¯^{F}_{k}

with 1≤k≤m.The commandMPLUnshuffle(expr,var)applies the unshuffle mapΦto a functionexprinV(Ωm)where the second argumentvaris the last of the cubical coor- dinatesxm.The hyperlogarithm in the right-hand part of the resulting tensor-product is a function of this variable (ifexprdepends on this variable).

In order to give a non-trivial example, let us at first construct a function inV(Ω3) as follows.

Consider the hyperlogarithms

f_{1} = bar

x_{1}d(x2)
1−x_{1}x_{2}

∈V Ω¯^{F}_{2}
,
f_{2} = bar

x2d(x3)

1−x_{2}x_{3},d(x3)
x_{3}

∈V Ω¯^{F}_{3}
.

We apply the map of eq. 4 twice. At first we compute g1=µ(id⊗Ψ) (1⊗ f1) =Ψ(f1) =bar

d(x1x_{2})
1−x_{1}x_{2}

∈V(Ω2).

Let us furthermore construct the function (or integrable word)g_{12}as
g12 = µ(id⊗Ψ) (g1⊗f2)

= f˜_{1}_{x}Ψ(f_{2})

= bar

x_{3}d(x2) +x_{2}d(x3)
1−x2x3

,d(x3) x3

,x_{2}d(x1) +x_{1}d(x2)
1−x1x2

+bar

x_{3}d(x2) +x_{2}d(x3)
1−x2x3

,x_{2}d(x1) +x_{1}d(x2)
1−x1x2

,d(x3) x3

+bar

x_{2}d(x1) +x_{1}d(x2)
1−x1x2

,x_{3}d(x2) +x_{2}d(x3)
1−x2x3

,d(x3) x3

−bar

d(x2) x2

,x_{3}d(x2) +x_{2}d(x3)
1−x2x3

,x_{2}d(x1) +x_{1}d(x2)
1−x1x2

−bar

d(x2) x2

,x_{2}d(x1) +x_{1}d(x2)
1−x1x2

,x_{3}d(x2) +x_{2}d(x3)
1−x2x3

−bar

x_{2}d(x1) +x_{1}d(x2)
1−x1x2

,d(x2) x2

,x_{3}d(x2) +x_{2}d(x3)
1−x2x3

.

This construction is easily obtained by use of the above procedures MPLShuffleProduct and MPLSymbolMap:

>f1:=bar(x[1]*d(x[2])/(1-x[1]*x[2])):

>f2:=bar(x[2]*d(x[3])/(1-x[2]*x[3]), d(x[3])/x[3]):

>MPLCoordinates(x,2):

>g1:=MPLSymbolMap(f1):

>MPLCoordinates(x,3):

>g12:=MPLShuffleProduct(g1,MPLSymbolMap(f2)):

By construction, the functiong12 is inV(Ω3)and it is easy to check, that it is in the basis of this vectorspace obtained by the commandMPLBasis:

>B3:=MPLBasis(x,3,3):

>has(B3[3],g12);

true

Now let us reverse this construction by use of the unshuffle map. By computing

>MPLUnshuffle(g12,x[3]);

tens

bar
_{d(x}

1x2) 1−x1x2

,bar

x2d(x3)
1−x2x3,^{d(x}_{x}^{3}^{)}

3

we clearly obtain the expected

Φ(g12) =g_{1}⊗f_{2}.
Furthermore computing

>MPLUnshuffle(g1,x[2]);

tens

ONE,bar_{x}

1d(x_{2})
1−x1x_{2}

we confirm the expected

Φ(g1) =1⊗f_{1}.

### 2.4 Integration over cubical coordinates

**MPLPrimitive(a, f, w)**

The procedure MPLPrimitive(a, f, w) computes the primitive^{R}a⊗f ∈V(Ωm) of a function
f ∈V(Ωm)of maximal weight wwith respect to the differential 1-forma∈

### A

¯_{m}

^{F}.For example let us compute

^{R}a⊗ f for

a= x2d(x3)
1−x_{2}x_{3}
and

f =bar

x_{3}d(x_{2}) +x_{2}d(x_{3})
1−x_{2}x_{3} ,dx_{3}

x_{3}

−bar

d(x_{2})

x_{2} ,x_{3}d(x_{2}) +x_{2}d(x_{3})
1−x_{2}x_{3}

+5bar

d(x_{2})
x_{2}

+99.

The words in f are at most of lengthw=2.We compute:

>MPLCoordinates(x,3):

>a:=x[2]*d(x[3])/(1-x[2]*x[3]):

>f:=bar((x[3]*d(x[2])+x[2]*d(x[3]))/(1-x[2]*x[3]), d(x[3])/x[3]) -bar(d(x[2])/x[2], (x[3]*d(x[2])+x[2]*d(x[3]))/(1-x[2]*x[3])) +5*bar(d(x[2])/x[2])+99:

>MPLPrimitive(a, f, 2);

bar

x_{3}d(x2) +x_{2}d(x3)
1−x2x3

,x_{3}d(x2) +x_{2}d(x3)
1−x2x3

,dx_{3}
x3

−bar

x_{3}d(x2) +x_{2}d(x3)
1−x2x3

,d(x2) x2

,x_{3}d(x2) +x_{2}d(x3)
1−x2x3

+5bar

x_{3}d(x2) +x_{2}d(x3)
1−x2x3

,d(x2) x2

+99bar

x_{3}d(x2) +x_{2}d(x3)
1−x2x3

+5bar

d(x2)

x_{2} ,x_{3}d(x2) +x_{2}d(x3)
1−x_{2}x_{3}

−bar

d(x2)

x_{2} ,x_{3}d(x2) +x_{2}d(x3)

1−x_{2}x_{3} ,x_{3}d(x2) +x_{2}d(x3)
1−x_{2}x_{3}

**MPLLimit(f,x,u)**

The procedure MPLLimit(expr,var,limit)returns the limit of a function f inV(Ωm) at x=u where x is one of the m cubical coordinates and u is either 0 or 1. As an example we consider m=3 and

f = Ψ

x1x2d(x3)

1−x_{1}x_{2}x_{3},x2d(x3)
1−x_{2}x_{3}

=

d(x1x2x3)

1−x_{1}x_{2}x_{3}− d(x1)

1−x_{1},d(x2x3)
1−x_{2}x_{3}

+

d(x1)

1−x_{1}+d(x1)

x_{1} ,d(x1x2x3)
1−x_{1}x_{2}x_{3}

(5)
and we compute the limitsg_{1}=lim_{x}_{3}_{→1}f,g_{2}=lim_{x}_{2}_{→1}g_{1},g_{3}=lim_{x}_{1}_{→1}g_{2}as follows:

>f:=MPLSymbolMap(bar(x[1]*x[2]*d(x[3])/(1-x[1]*x[2]*x[3]), x[2]*d(x[3])/(1-x[2]*x[3]))):

>g1:=simplify(MPLLimit(f,x[3],1));

g1 : = bar

x_{2}d(x1) +x_{1}d(x2)

−1+x1x2

, d(x2)

−1+x2

−bar

d(x1)

−1+x1

, d(x2)

−1+x2

+bar

d(x1)

−1+x_{1},x_{2}d(x1) +x_{1}d(x2)

−1+x_{1}x_{2}

−bar

d(x1)

x_{1} ,x_{2}d(x1) +x_{1}d(x2)

−1+x_{1}x_{2}

>g2:=simplify(MPLLimit(g1,x[2],1));

g2 :=−bar

d(x1)

x_{1} , d(x1)

−1+x_{1}

+bar

d(x1)

−1+x_{1}, d(x1)

−1+x_{1}

>g3:=simplify(MPLLimit(g2,x[1],1));

g3 :=ζ(2)

There is a little subtlety we have to take care of, if we compute limits not starting with the last
cubical coordinatexm.For example, let us computeg4=lim_{x}_{2}_{→1}f by

>g4:=simplify(MPLLimit(f,x[2],1));

bar

x3d(x1)+x1d(x3)

−1+x1x3 ,_{−1+x}^{d(x}^{3}^{)}

3

−bar_{d(x}

1)

x1 ,^{x}^{3}^{d(x}_{−1+x}^{1}^{)+x}^{1}^{d(x}^{3}^{)}

1x3

+bar _{d(x}

1)

−1+x1,^{x}^{3}^{d(x}_{−1+x}^{1}^{)+x}^{1}^{d(x}^{3}^{)}

1x3

−

bar
_{d(x}

1)

−1+x1,_{−1+x}^{d(x}^{3}^{)}

3

In this result, the differential 1-form ^{x}^{3}^{d(x}_{−1+x}^{1}^{)+x}^{1}^{d(x}^{3}^{)}

1x3 does not belong to anyΩ_{m} defined in eq. 1,
because it involves the product x1x3 of non-consecutive cobical coordinates. However, we can
simply replace everyx_{3}byx_{2}in this result, for example by

>g5:=subs(x[3]=x[2],g4):

Now ing_{5} every 1-form belongs toΩ_{2} and the result is recognized as a function ofV(Ω2). (It is
easy to check the integrability.) In general, after taking a limit with respect to a variablex_{k} with
1≤k<m,we have to re-name the coordinates by(xk+1, ...,x_{m})7→(xk, ...,x_{m−1}).

**MPLMultipleLimit(expr,list)**

Based on the previous procedure, the procedureMPLMultipleLimit(f,L)computes several limits of a functionfofV(Ωm)successively. HereLis a list of the type

x_{σ(1)}=u_{1},x_{σ(2)}=u_{2}, ...,x_{σ(k)}=u_{k}
for 1≤k≤m, σa permutation on some subset of{1, ...,m}and allu_{i}∈ {0,1}.This list fixes the
limits and also the order in which they are computed. The above command returns

x_{σ(k)→}lim_{uk}... lim

x_{σ(1)→u}_{1} f.

As an example, let us again consider the function f as defined in eq. 5. Then we obtain the above
results forg_{1},g_{2},g_{3}by

>g1:=simplify(MPLMultipleLimit(f,[x[3]=1])):

>g2:=simplify(MPLMultipleLimit(f,[x[3]=1,x[2]=1])):

>g3:=simplify(MPLMultipleLimit(f,[x[3]=1,x[2]=1,x[1]=1])):

It is briefly mentioned in [2] and discussed in more detail in [4] and [6] that the order of such consecutive limits is crucial. For example, we have obtainedg3=ζ(2),but permuting the order of the last two limits, we compute

>MPLMultipleLimit(f,[x[3]=1,x[1]=1,x[2]=1]);

0

Note that in a function of V(Ωm) the terms with weight greater than zero are defined to vanish at the origin of the space of cubical coordinates. Therefore, if we choose the above list L to be

x_{σ(1)}=0,x_{σ(2)}=0, ...,x_{σ(m)}=0

, the procedure MPLMultipleLimit returns the value that remains after by replacing everybar(...) by zero in f ∈V(Ωm).

**MPLCubicalIntegrate(f,var,n)**

MPL can compute multiple integrals of the type I=

Z _{1}

0

dx_{m} q

∏_{i}p^{a}_{i}^{i} f (6)

where where f ∈V(Ωm),qis a polynomial inx_{m}, thea_{i}∈Nand thep_{i}are in{xm,1−x_{m}, ...,1−x_{1}· · ·x_{m}}.
For f being an integrand of this type, the procedure MPLCubicalIntegrate(f,var,n) succes-

sively integrates overx_{m},x_{m−1}, ...,x_{m−n+1} in this order. The first argument of the procedure is the
integrand in bar-notation, the second argument is the variable to be first integrated out (i.e. x_{m} in
the notation chosen here) and n is the number of integrations. As an example, we consider the
integrand

f = u_{1}^{3}(1−u_{1})u_{2}^{4}(1−u_{2})u_{3}^{3}(1−u_{3})u_{4}^{2}(1−u_{4})^{2}
(1−u1u2u3)^{2}(1−u2u3u4)^{2}(1−u1u2u3u4)^{2}(1−u1u2).
We integrate out all four cubical coordinatesu_{4},u_{3},u_{2},u_{1}by

>f:=u[1]^3*(1-u[1])*u[2]^4*(1-u[2])*u[3]^3*(1-u[3])

*u[4]^2*(1-u[4])^2/((1-u[1]*u[2]*u[3])^2*(1-u[4]*u[2]

*u[3])^2*(1-u[4]*u[1]*u[2]*u[3])^2*(1-u[1]*u[2])):

>MPLCubicalIntegrate(f,u[4],4);

(5/3)∗Zeta(3) + (26/9)∗Zeta(2)−(17/5)∗Zeta(2)^{2}+22/9

So we obtain

4

### ∏

i=1_{Z} _{1}

0

dx_{i}

f = 5

3ζ(3) +26

9 ζ(2)−17

5 (ζ(2))^{2}+22
9 .

Intermediate results are too long to be shown here, but can easily be obtained by choosing the last ar- gument ofMPLCubicalIntegratesmaller than 4. Note that here the explicit use ofMPLCoordinates is not necessary (because it is called internally byMPLCubicalIntegrate).

MPLCubicalIntegratealready serves for the analytical computation of a large class of inte- grals and further example applications are shown in [2, 4]. In some cases, the procedure may even be useful for the computation of Feynman integrals. However, in the case of Feynman integrals, the denominators of the integrand usually involve more complicated polynomials than allowed here.

Therefore, appropriate changes of variables are necessary, to express the latter in terms of cubical coordinates. These are subject of the following section.

### 3 Computing Feynman integrals

In this section, we consider integrals of the type
I_{F} =

Z _{∞}

0

dα_{j}∏_{Q}_{i}_{∈}_{Q}Q^{δ}_{i}^{i}L_{w} α_{j}

∏_{P}_{i}_{∈}_{P}_{F}P_{i}^{β}^{i}

(7)

where

### Q

,### P

⊂Qα_{1}, ...,α_{j}, ...,α_{N}

are sets of irreducible polynomials, all δ_{i},β_{i}∈N∪ {0} and
L_{w}(αN)is a hyperlogarithm, given by a wordwin differential 1-forms in

Ω_{F} =

dα_{j}

α_{j} , dα_{j}

α_{j}−ρ_{i} whereρ_{i}=−P_{i}|αj=0

∂Pi

∂αj

,P_{i}∈

### P

.

Such integrals arise in the computation of Feynman integrals, possibly also in other contexts. In [4, 2] we have discussed the conditions under which we can compute such integrals with MPL.

Here let us just recall briefly, that the integrals have to be finiteandunramifiedand the set

### P

has to be linearly reducibleand properly ordered. We assume, that the user has already expressed a divergent Feynman integral in terms of finite integrals of this type, possibly by use of techniques cited and partly demonstrated in [4, 2]. All remaining mentioned conditions can be explicitely checked by procedures of MPL, as demonstrated below.### 3.1 Setting up a simple example

As a simple standard example, let us consider the massless, off-shell one-loop traingle graph throughout this section. The corresponding integral is finite. (More involved examples are dis- cussed in [4, 2]). The two Symanzik polynomials (see e.g. [5]) of this graph are

### U

= x_{1}+x

_{2}+x

_{3},

### F

= −x1x_{2}p

^{2}

_{3}−x

_{2}x

_{3}p

^{2}

_{1}−x

_{1}x

_{3}p

^{2}

_{2},

where the x_{i} are the Feynman (or Schwinger, or alpha) parameters and the p_{i} are the incoming
momenta. In section 4 of [2], we suggested to express the kinematical dependences in terms of
variables, which can be treated as additional Feynman parameters, which are not integrated out. In
this sense, we define the set of generalized Feynman parameters {x1,x_{2},x_{3},x_{4},x_{5}} wherex_{4} and
x_{5}satisfy

p^{2}_{3}

p^{2}_{2} =x_{4}x_{5}and p^{2}_{1}

p^{2}_{2} = (1−x_{4})(1−x_{5}).

(C.f. [10].) We consider the kinematical region where 0≤x_{4},0≤x_{5}.(For a different kinematical
region, we recommend to change the definition ofx_{4}andx_{5}accordingly, such that the region is still
given by 0≤x_{4},0≤x_{5},if possible.) Using these variables, we define the auxiliary polynomial

### F

˜ =−### F

p^{2}_{2} =x_{1}x_{2}x_{4}x_{5}+x_{1}x_{3}+x_{2}x_{3}(1−x_{4}) (1−x_{5}).

The corresponding scalar Feynman integral inD=4−2εdimensions (without any scalar products in the numerator or raised propagator-powers) can be defined by

### I

(x4,x_{5},ε) =Γ(1+ε) −p

^{2}

_{2}−1−ε

I(x4,x_{5},ε)

with

I(x4,x_{5},ε) =

3

### ∏

i=1_{Z} _{∞}

0

dxi

δ(H)

### U

^{−1+2ε}

### F

^{˜}

^{−1−ε}

where H is a hyperplane in the integration domain, which we can choose freely according to the
Cheng-Wu theorem [12]. We chooseδ(1−x_{n})wherex_{n} is our last integration variable. The order
of integration variables still has to be chosen such that the mentioned conditions are satisfied.

We expand the integrand inεand obtain

I(x4,x_{5},ε) =I^{(0)}(x4,x_{5}) +εI^{(0)}(x4,x_{5}) +

### O

ε^{2}

with

I^{(0)}(x4,x_{5}) =

3

### ∏

i=1_{Z} _{∞}

0

dx_{i}

δ(H) 1

### U F

^{˜}

^{,}

^{(8)}

I^{(1)}(x4,x_{5}) =

3

### ∏

i=1_{Z} _{∞}

0

dx_{i}

δ(H)2 ln(

### U

)−ln ˜### F

^{}

### U F

^{˜}

^{.}

^{(9)}

In the following, we compute both of the latter integrals with MPL.

### 3.2 Checking the conditions

Before the computation of any integral over Feynman parameters with MPL, we recommend to check, whether the mentioned conditions are satisfied. We assume, that the user already has checked the finiteness of the integral. The following procedure can be used to check linear re- ducibility.

**MPLPolynomialReduction(S,L[1..n],L), COMPATIBILITY_GRAPH**

The commandMPLPolynomialReduction(S,L[1..n],L)returns the reduction of a set S of poly- nomials with respect tonintegration variables (i.e. Feynman parameters) L[1..n], which are the firstnentries in the listLof all variables of the problem (i.e. generalized Feynman parameters, given

by the Feynman parameters and kinematical invariants). If a global variableCOMPATIBILITY_GRAPH=true, the reduction is computed regarding compatibilities among the polynomials as introduced in [9].

The program uses the version of such compatibilities defined in [14]. IfCOMPATIBILITY_GRAPH=false, the procedure returns a Fubini-reduction as introduced in [8]. See [2] for the details.

In order to interpret the output of this procedure, it is instructive to look at our example of the
triangle graph. For the computation ofI^{(0)}(x4,x_{5})andI^{(1)}(x4,x_{5})we clearly should compute the
reduction of S=

### U

,### F

˜ .The list of generalized Feynman parameters isx_{1},x

_{2},x

_{3},x

_{4},x

_{5}and the integration variables are the first three elements of this list. Therefore we compute

>U:=x[1]+x[2]+x[3]:

>F:=x[1]*x[2]*x[4]*x[5]+x[1]*x[3]+x[2]*x[3]*(1-x[4])*(1-x[5]):

>GenFeyn:=[x[1],x[2],x[3],x[4],x[5]]:

>COMPATIBILITY_GRAPH=true: #This is also the default value.

>MPLPolynomialReduction([U,F],GenFeyn[1..3],GenFeyn);

[[{},[x1+x_{2}+x_{3},x_{1}x_{2}x_{4}x_{5}+x_{1}x_{3}+x_{2}x_{3}(1−x_{4}) (1−x_{5})],[{1,2}]],

[{x3},[x1+x_{2}−x_{2}x_{5}−x_{4}x_{2}+x_{2}x_{4}x_{5},x_{1}+x_{2},x_{1}+x_{2}−x_{2}x_{5},x_{1}+x_{2}−x_{4}x_{2}],

[{3,4},{1,3},{1,4},{2,3},{2,4}]],

[{x2},[x1x_{4}x_{5}+x_{3}−x_{3}x_{5}−x_{3}x_{4}+x_{3}x_{4}x_{5},x_{1}+x_{3},−x3+x_{1}x_{5}+x_{3}x_{5},x_{1}x_{4}−x_{3}+x_{3}x_{4}],
[{3,4},{1,3},{1,4},{2,3},{2,4}]],

[{x1},[x2x_{4}x_{5}+x_{3},x_{2}+x_{3},−1+x_{4},−1+x_{5},x_{4}x_{2}+x_{3},x_{3}+x_{2}x_{5}],

[{3,4},{1,3},{1,4},{2,3},{2,4},{3,5},{4,5},{5,6},{1,5},{1,6},{2,5},{2,6},{3,6}]],
[{x1,x_{2}},[−1+x_{5},−1+x_{4},x_{4}−x_{5}],[{1,2},{1,3},{2,3}]],

[{x1,x_{3}},[−1+x_{5},−1+x_{4},x_{4}−x_{5}],[{1,2},{1,3},{2,3}]],
[{x2,x3},[−1+x4,−1+x5,x4−x5],[{1,2},{1,3},{2,3}]],

[{x1,x_{2},x_{3}},[−1+x_{5},−1+x_{4},x_{4}−x_{5}],[{1,2},{1,3},{2,3}]]]

This reduction is a list of 8 lists, each of which has three entries. The first entry is the set of variables with respect to which the polynomials in S were reduced. The second entry is the list of resulting polynomials at this stage. (In the language of [9, 14], these are the vertices of the compatibility graph.) This list is ordered. With respect to this ordering, the third entry contains all compatible pairs of polynomials (i.e. the edges of the compatibility graph) represented by their numbers with respect to the order in the second entry.

For example let us consider the list in the second line from below:

[{x2,x3},[−1+x4,−1+x_{5},x4−x_{5}],[{1,2},{1,3},{2,3}]].

According to the first entry, reductions were computed with respect tox2andx3.Therefore, in the notation of [2], the second entry is the set

S˜^{{2,3}}={−1+x4,−1+x_{5},x4−x_{5}}

and the third entry tells us, that all three polynomials are compatible with each other. If in this entry
for example{2,3}would be missing,−1+x_{5}would not be compatible withx_{4}−x_{5}.

The compatibilities are usually just of internal importance. If they are considered by the algo- rithm (i.e. ifCOMPATIBILITY_GRAPHis true), the reduction usually gives a better upper bound for

the polynomials appearing in the integration procedure, than the plain Fubini reduction (COMPATIBILITY_GRAPH is false). In the output of latter, the third entry of each list is always the complete graph, i.e. it con-

tains all possible pairs. The plain Fubini reduction obtained by

>COMPATIBILITY_GRAPH=false:

>MPLPolynomialReduction([U,F],GenFeyn[1..3],GenFeyn);

For this simple example, both reductions are exactly the same.

For the user, the important question is, whether the setSis linearly reducible and, if yes, which
are the allowed orders of integrations. This can information is easily obtained from the reduction,
by looking at the first entry of each list. The set Sis linearly reducible with respect to a particu-
lar order x_{σ(1)}, ...,x_{σ(n)} ofn integration variables, given by a permutation σ on{1, ...,n}, if the
reduction contains each of the sets

x_{σ(1)} ,

x_{σ(1)},x_{σ(2)} , ...,

x_{σ(1)}, ...,x_{σ(n)} .
For example, let us ask, whetherS=

### U

,### F

˜ is linearly reducible with respect to the order x_{3},x

_{1},x

_{2}.We see that the above reduction contains lists whose first etries are{x3},{x1,x

_{3}},{x1,x

_{2},x

_{3}}, so the answer is yes. If for example the list with{x3}in the first entry would be missing, neither x

_{3},x

_{1},x

_{2}norx

_{3},x

_{2},x

_{1}would be an allowed order of integration.

**MPLCheckOrder(reduction,L[1..n],L)**

Consider a given polynomial reduction in the form as returned by the previous command (no mat-
ter whether compatibilities are regarded or not), and a given ordered list L of generalized Feyn-
man parameters where the first n entries are the integration variables. For this data, the proce-
dureMPLCheckOrder(reduction,L[1..n],L)checks with respect to the order of variables inL,
whether the integrand is linearly reducible, whether it is properly ordered at a tangential basepoint
(see [2, 4]) and whether all limits which have to be computed in the integration procedure are at 0 or
1.^{3} All three conditions are required for an automatic computation of the integral with MPL. If all
conditions are satisfied, the procedure returns a sequence of messages, all ending with the phrase

“Check OK.” If one of the conditions is violated, the procedure returns a corresponding error mes- sage. In this case, a fully automated computation of the integral in the order of integrations given by L is not possible with MPL. We recommend to try the command again with a different order of variables in L. In some cases, the problem also may be solved by a different choice of kinematical invariants (see remarks in sections 4.2 and 4.3 of [2]).

Let us apply the procedure in our example of the triangle graph whereGenFeynis again the list of generalized Feynman parameters defined above:

>Reduction:=MPLPolynomialReduction([U,F],GenFeyn[1..3],GenFeyn):

>MPLCheckOrder(Reduction,GenFeyn[1..3],GenFeyn);

Checks for integration over x[1]:

"All rho_i are smaller than zero at the tangential basepoint: Check OK."

"Limits in {0, 1}: Check OK."

3Note that the latter property is implied by the condition of the integrand to be unramified as discussed in [2, 4]. The procedure does not check unramifiedness but it checks the condition directly for the required limits at the tangential basepoint.