
\RequirePackage{ifpdf}

\documentclass[a4paper,11pt]{article}
\pdfoutput=1 % if your are submitting a pdflatex (i.e. if you have
             % images in pdf, png or jpg format)

\usepackage{jheppub} % for details on the use of the package, please
                     % see the JHEP-author-manual

\usepackage[T1]{fontenc} % if needed


\usepackage{ifpdf} % to get it to work with pdflatex and not break for everything else

\usepackage{flexisym}
\usepackage{breqn}
\usepackage{comment}


\usepackage{epsf,epsfig,epstopdf}
\usepackage{graphics}
\usepackage{graphicx}
\usepackage[usenames,dvipsnames]{xcolor}
\definecolor{darkgreen}{RGB}{34,139,34}

%\usepackage[pdftex]{hyperref} %uncomment this line to make the references accessible to firefox %


\usepackage{enumerate}

%\usepackage{url}
\usepackage{amsmath,amsfonts,amssymb}
%\usepackage{cite}

%\usepackage{breqn}
%\usepackage{comment}

%\usepackage{a4wide}
%\numberwithin{equation}{section}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\renewcommand{\topfraction}{1}

\arraycolsep0.1cm

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% additions CS

\DeclareMathOperator{\tr}{tr}
\def\bra#1{\mathinner{\langle{#1}|}}
\def\ket#1{\mathinner{|{#1}\rangle}}
\def\braket#1{\mathinner{\langle{#1}\rangle}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%% %simple case: 2 authors, same institution
%% \author{A. Uthor}
%% \author{and A. Nother Author}
%% \affiliation{Institution,\\Address, Country}

% more complex case: 4 authors, 3 institutions, 2 footnotes
\author{Costas G. Papadopoulos$^{a,b}$, Damiano Tommasini$^{a}$ and Christopher Wever$^{a,c}$}

% The "\note" macro will give a warning: "Ignoring empty anchor..."
% you can safely ignore it.

\preprint{TTP15-042}

\affiliation{$^{a}$Institute of Nuclear and Particle Physics, NCSR `Demokritos', Agia Paraskevi, 15310, Greece}
\affiliation{$^{b}$MTA-DE Particle Physics Group, University of Debrecen, H-4010 Debrecen, Hungary}
\affiliation{$^{c}$Institute for Theoretical Particle Physics (TTP), Karlsruhe Institute of Technology, Engesserstra{\ss}e 7, D-76128 Karlsruhe, Germany $\&$ 
Institute for Nuclear Physics (IKP), Karlsruhe Institute of Technology, Hermann-von-Helmholtz-Platz 1, D-76344 Eggenstein-Leopoldshafen, Germany}

% e-mail addresses: one for each author, in the same order as the authors
\emailAdd{costas.papadopoulos@cern.ch, tommasini@inp.demokritos.gr, wever@inp.demokritos.gr}


\keywords{Feynman integrals, QCD, NLO and NNLO Calculations}





\title{The Pentabox Master Integrals with the Simplified Differential Equations approach}

\abstract{We present the calculation of massless two-loop Master Integrals relevant to five-point amplitudes with one off-shell external leg and derive the complete set of planar Master Integrals with five on-mass-shell legs, that contribute 
to many $2\to 3$ amplitudes of interest  at the LHC, as for instance 
 three jet production, $\gamma, V, H +2$ jets etc., based on the  Simplified Differential 
Equations approach.}





\ifpdf
  \DeclareGraphicsRule{*}{mps}{*}{}
\fi



\begin{document}
\unitlength1cm
\maketitle
\allowdisplaybreaks


\section{Introduction}

With LHC delivering collisions at the highest energy achieved so far, 13 TeV, experiments are analysing data corresponding to an integrated luminosity of $42$ pb$^{-1}$~\cite{Khachatryan:2015uqb} and $85$ pb$^{-1}$~\cite{Atlas}, 
as well as those already collected at an energy of 8 TeV and an integrated luminosity of $20.3$ fb$^{-1}$~\cite{Aad:2015nda} and $19.7$ fb$^{-1}$~\cite{Khachatryan:2015tzo}. 
In order to keep up with the increasing experimental accuracy as more data is collected, more precise theoretical predictions and higher loop calculations are required~\cite{Andersen:2014efa}.

With the better understanding of reduction of one-loop amplitudes to a set of Master Integrals ({\bf MI}) based on unitarity methods~\cite{Bern:1994cg,Bern:1994zx} and at the integrand level via the OPP 
method~\cite{Ossola:2006us,Ossola:2008xq}, one-loop calculations have been fully automated in many numerical tools (some reviews on the topic are~\cite{AlcarazMaestre:2012vp,Ellis:2011cr}).
%\cite{vanHameren:2009dr,Berger:2009zg,Bern:2013gka,Bredenstein:2009aj,Bevilacqua:2009zn,Bevilacqua:2010ve}\footnote{I don't like citations: GoSam,Mad-family-F.Maltoni, Feynart-T.Hahn, OpenLoops-S.Pozzorini...}.
In the recent years, a lot of progress has been made towards the extension of these reduction methods to the two-loop order at the integral~\cite{Gluza:2010ws,Kosower:2011ty,CaronHuot:2012ab} 
as well as the integrand~\cite{Mastrolia:2011pr,Badger:2012dp,Badger:2013gxa,Papadopoulos:2013hra} level. Contrary to the MI at one-loop, which have been known for a long time already~\cite{'tHooft:1978xw}, a complete library of MI at two-loops is still missing.
At the moment this seems to be the main obstacle to obtain a fully automated NNLO calculation framework similar to the one-loop case, that will satisfy the anticipated precision requirements at the LHC~\cite{Butterworth:2014efa}.

Starting from the works of~\cite{Goncharov:1998kja,Remiddi:1999ew,Goncharov:2001iea}, there has been a building consensus that the so-called {\it Goncharov Polylogarithms} ({\bf GPs}) form a functional basis for many MI. A very fruitful method for calculating MI and expressing them in terms of GPs is the differential equations ({\bf DE}) approach~\cite{Kotikov:1990kg,Kotikov:1991pm,Bern:1992em,Remiddi:1997ny,Henn:2013pwa}, which has been used in the past two decades to calculate various MI at two-loops ~\cite{Caffo:1998du,Gehrmann:1999as,Gehrmann:2000zt,Gehrmann:2001ck,Henn:2014lfa,Caola:2014lpa}.
%and more recently a new idea was proposed in~\cite{Henn:2013pwa}. 
In~\cite{Papadopoulos:2014lla} a variant of the traditional DE approach to MI was presented, which was coined the Simplified Differential Equations ({\bf SDE}) approach. 
In this paper we present an application of this method, concerning the calculation of planar massless MI relevant to five-point amplitudes with one off-shell leg, as well as the complete set of planar MI for five-point on-shell amplitudes.
This is an important step towards the calculation of the full set of MI with up to eight internal propagators needed to realise a fully automated reduction scheme, \`a la OPP, for NNLO QCD. 
 
Pentabox integrals are needed in particular in order to compute NNLO QCD corrections to several processes of interest at LHC~\cite{Andersen:2014efa}. The $pp\rightarrow H+2$jets can be used to measure the $HWW$ coupling to a $5\%$ accuracy 
with 300 fb$^{-1}$ data. The $pp\rightarrow 3$jets to study the ratio of $3-$jet to $2-$jet cross sections and measure the running of the strong coupling constant. 
The   $pp\rightarrow V+2$jets for PDF determination and background studies for multi-jet final states.


The paper is organized as follows. In the Section 2 we set the parameterization and notation of the variables describing the two-loop MI of interest. 
In Section 3 we discuss the DE obtained, and the results for the pentabox MI. 
We conclude in Section 4 and provide an overview of the topic and some perspective for future developments. 
In the Appendix we give a few characteristic examples on how the boundary conditions are properly reproduced in our approach by the DE. 
Finally in the ancillary files~\cite{results}, we provide our analytic results for all  two-loop MI in terms of Goncharov polylogarithms together with explicit numerical results.


\section{The pentabox integrals}


The MI in this paper will be calculated with the SDE approach~\cite{Papadopoulos:2014lla}.
Assume that one is interested in calculating an $l-$loop Feynman integral with external momenta $\{p_j\}$, considered incoming, and internal propagators that are massless.
%, since the method has been originally developed for massless internal lines. 
Any $l-$loop  Feynman integral can be then written as
\begin{equation}
G_{a_1\cdots a_n}(\{p_j\},\epsilon)=\int\left(\prod_{r=1}^l \frac{d^dk_r}{i\pi^{d/2}}\right)\frac{1}{D_1^{a_1}\cdots D_n^{a_n}}, \hspace{0.5 cm}
D_i=\left(c_{ij}k_j+d_{ij}p_j\right)^2,\,\, d=4-2\epsilon
\label{eq:loopgen}
\end{equation}
with matrices $\{c_{ij}\}$ and $\{d_{ij}\}$ determined by the topology and the momentum flow of the graph, and
the denominators are defined in such a way that all scalar product invariants can be written as a linear combination of them. The exponents $a_i$ are integers and may be negative in order to accommodate irreducible numerators.

Any integral $G_{a_1\cdots a_n}$ may be written as a linear combination of a finite subset of such integrals, called Master Integrals,  
%$\vec{G}^{MI}(\{ s_{ij} \},\epsilon)$ 
with coefficients depending on the independent scalar products, $s_{ij}=p_i\cdot p_j$, and space-time dimension $d$, by the use of {\it integration by parts} ({\bf IBP}) identities~\cite{Chetyrkin:1981qh,Tkachov:1981wb}.  
In the traditional DE method, the MI $\vec{G}^{MI}(\{s_{ij}\},\epsilon)$, 
are differentiated with respect to $p_i \cdot \frac{\partial}{\partial p_j}$ 
%a Lorentz invariant (or set of invariants) $\textbf{s}$
, and the resulting integrals are reduced by IBP to give a linear system of DE for 
$\vec{G}^{MI}(\{s_{ij}\},\epsilon)$~\cite{Kotikov:1990kg,Remiddi:1997ny}. The  invariants, $s_{ij}$, are then parametrised in terms of dimensionless variables, defined on a case by case basis, so that the resulting DE can be solved in terms of GPs. 
Usually boundary terms corresponding to the appropriate limits of the chosen parameters have to be calculated 
using for instance expansion by regions techniques~\cite{Beneke:1997zp,Smirnov:2002pj}.

SDE approach~\cite{Papadopoulos:2014lla} is an attempt not only to simplify, but also to systematize, as much as possible, the derivation of the appropriate system of DE satisfied by the MI. 
To this end 
the external incoming momenta are {\it parametrized} linearly in terms of $x$ as $p_i(x)=p_i+(1-x)q_i$, where the $q_i$'s are a linear combination of the  momenta $\{p_i\}$ such that $\sum_iq_i=0$. If $p_i^2=0$, the parameter $x$ captures the off-shell-ness of the external legs. 
The class of Feynman integrals in (\ref{eq:loopgen}) are now dependent on $x$ through the external momenta:
\begin{equation}
G_{a_1\cdots a_n}(\{s_{ij}\},\epsilon;x)=\int\left(\prod_{r=1}^l \frac{d^dk_r}{i\pi^{d/2}}\right)\frac{1}{D_1^{a_1}\cdots D_n^{a_n}}, 
\;\;\;
D_i=\left( c_{ij}k_j+d_{ij}p_j(x) \right)^2.
%\hspace{0.5 cm} s=\{p_i.p_j\}|_{i,j}, 
\label{eq:loopgenx}
\end{equation}


By introducing the dimensionless parameter $x$, the vector of MI $\vec{G}^{MI}(\{s_{ij}\},\epsilon;x)$, which now depends on $x$, satisfies
\begin{equation}
\frac{\partial}{\partial x} \vec{G}^{MI}(\{s_{ij}\},\epsilon;x)=\mathbf{H}(\{ s_{ij}\},\epsilon;x)\vec{G}^{MI}(\{s_{ij}\},\epsilon;x)\label{eq:DEx}
\end{equation}
a system of differential equations in one independent variable. 

Experience up to now shows that this simple parametrization can be used universally to deal with up to six kinematical scales involved, as is the case we will present in this paper.
 The expected benefit is that the integration of the DE naturally captures the expressibility of MI in terms of GPs and more importantly makes the problem 
 {\it independent on the number of kinematical scales} (independent invariants) involved. 
Note that as $x\rightarrow 1$, the original configuration of the loop integrals (\ref{eq:loopgen}) is reproduced, which correspond to a simpler one with one scale less. 



\begin{figure}[t!]
\centering
\includegraphics[width=0.20 \linewidth]{Figures/Section2/planar-pentabox-P1_withp-eps-converted-to.pdf} \hspace{0.4 cm}
\includegraphics[width=0.20 \linewidth]{Figures/Section2/planar-pentabox-P2_withp-eps-converted-to.pdf} \hspace{0.4 cm}
\includegraphics[width=0.20 \linewidth]{Figures/Section2/planar-pentabox-P3_withp-eps-converted-to.pdf}
  \caption{The three planar pentaboxes of the families $P_1$ (left), $P_2$ (middle) and $P_3$ (right) with one external massive leg.}
  \label{fig:param-P}
\end{figure}

\begin{figure}[t!]
\centering
\includegraphics[width=0.20 \linewidth]{Figures/Section2/nonplanar-pentabox-N1_withp-eps-converted-to.pdf} \hspace{0.6 cm}
\includegraphics[width=0.20 \linewidth]{Figures/Section2/nonplanar-pentabox-N2_withp-eps-converted-to.pdf} \hspace{0.6 cm}
\includegraphics[width=0.20 \linewidth]{Figures/Section2/nonplanar-pentabox-N3_withp-eps-converted-to.pdf} \\[12pt]
\includegraphics[width=0.20 \linewidth]{Figures/Section2/nonplanar-pentabox-N4_withp-eps-converted-to.pdf} \hspace{0.6 cm}
\includegraphics[width=0.20 \linewidth]{Figures/Section2/nonplanar-pentabox-N5_withp-eps-converted-to.pdf}
  \caption{The five non-planar families with one external massive leg.}
  \label{fig:param-N}
\end{figure}

More specifically we are interested in calculating the MI of two-loop QCD five-point amplitudes. As it is an inherent characteristic of the
SDE method to interpolate among different kinematical configurations of the external momenta the starting point is to compute
five-point amplitudes with one off-shell leg. These amplitudes contribute to the production i.e., of one massive final state $V$, plus two massless final states $j_1,j_2$ at the LHC:
\begin{gather}
p(q_1)p'(q_2)\rightarrow V(q_3)j_1(q_4)j_2(q_5), \ \ q_1^2=q_2^2=0, \ \ q_3^2=M_3^2, \ \ q_4^2=q_5^2=0.\label{eq:ppVV}
\end{gather}
The colliding partons have massless momenta $q_1,q_2$, while the outgoing massive and the two massless particles have momenta $q_3$ and $q_4,q_5$ respectively. Of course, by appropriately taking the limit
$x=1$ the pentabox MI with all external massless momenta on-shell will be obtained, that are relevant for instance to the three-jet production
\begin{gather}
p(q_1)p'(q_2)\rightarrow j_1(q_3)j_2(q_4)j_3(q_5), \ \ q_i^2=0.\label{eq:pp3j}
\end{gather}

%In~\cite{Gehrmann:2015bfy} the planar contribution to the two-loop five-gluon all-plus-helicity gluon amplitudes was calculated. The results in~\cite{Gehrmann:2015bfy} correspond to the massless case $M_3^2\rightarrow 0$, for which there is one planar family of MI. 
For the off-shell case $M_3^2\neq 0$, there are in total three families of planar MI whose members with the maximum amount of denominators, namely eight, are graphically shown in Figure \ref{fig:param-P}. Similarly, there are five non-planar families of MI as given in Figure \ref{fig:param-N}. We have checked that the other five-point integrals with one massive external leg are reducible to MI in one of these eight MI families. The two-loop planar (Fig. \ref{fig:param-P}) and non-planar (Fig. \ref{fig:param-N}) diagrams contributing to (\ref{eq:ppVV}) have not been calculated yet. In fact by taking the limit $x=1$ all planar graphs for massless on-shell external momenta are derived as well~\footnote{During the writing of the present paper, some results related to massless 
planar pentaboxes appeared in~\cite{Gehrmann:2015bfy}.}. 
In this paper we calculate all MI in the family $P_1$ as well as all the the on-shell MI as $M_3^2\rightarrow 0$. We use the c++ implementation of the program {\bf FIRE}~\cite{Smirnov:2014hma} to perform the IBP reduction to the set of MI in $P_1$.

\begin{figure}[t!]
\centering
\includegraphics[width=0.29 \linewidth]{Figures/Section2/fivepoint2-P1_withpx-eps-converted-to.pdf} \hspace{0.4 cm}
\includegraphics[width=0.29 \linewidth]{Figures/Section2/fivepoint3-P1_withpx-eps-converted-to.pdf} \hspace{0.4 cm}
\includegraphics[width=0.22 \linewidth]{Figures/Section2/fivepoint4-P1_withpx-eps-converted-to.pdf} \\
\includegraphics[width=0.29 \linewidth]{Figures/Section2/fivepoint5-P1_withpx-eps-converted-to.pdf} \hspace{0.3 cm}
\includegraphics[width=0.29 \linewidth]{Figures/Section2/fivepoint6-P1_withpx-eps-converted-to.pdf} \hspace{0.3 cm}
\includegraphics[width=0.27 \linewidth]{Figures/Section2/fivepoint7-P1_withpx-eps-converted-to.pdf}
  \caption{The five-point Feynman diagrams, besides the pentabox itself in Figure \ref{fig:param-P}, that are contained in the family $P_1$. All external momenta are incoming.}
  \label{fig:5point}
\end{figure}

The family $P_1$ contains in total 74 MI. Up to five denominators integrals with doubled propagators have been used as MI whereas starting from six denominators integrals with irreducible numerators  are chosen.  Many of the 74 MI already appear in the families of the planar double box integrals discussed in~\cite{Henn:2014lfa,Caola:2014lpa,Papadopoulos:2014hla}. However, there are seventeen new five-point Feynman diagrams that are not contained in the double box integral families. Three of them are pentaboxes, including the scalar and two MI with irreducible numerators. 
There are six seven-denominator, and eight six-denominator ones, the scalar members of which are shown in Figure \ref{fig:5point}. 

%We note that in the massless limit $M_3^2\rightarrow 0$, the family $P_1$ contains the four five-point Feynman diagrams corresponding to MI for the on-shell case~\cite{Gehrmann:%2015bfy}.

\begin{figure}[t!]
\centering
\includegraphics[width=0.35 \linewidth]{Figures/Section2/planar-pentabox-P1_withpx-eps-converted-to.pdf}
  \caption{The parametrization of external momenta in terms of $x$ for the planar pentabox of the family $P_1$. All external momenta are incoming.}
  \label{fig:xparam-P1}
\end{figure}

%In the SDE method~\cite{Papadopoulos:2014lla} we introduce a parameter $x$ in the external momenta. 
For the family of integrals $P_1$ the external momenta are parametrized in $x$ as shown in Figure \ref{fig:xparam-P1}. The parametrization is chosen such that the double box MI with two massive external legs that is contained in the family $P_1$ has exactly the same parametrization as that one chosen in~\cite{Papadopoulos:2014hla}, i.e. two massless external momenta $xp_1$ and $xp_2$ and two massive external momenta $p_{123}-xp_{12}$ and $-p_{123}$. The MI in the family $P_1$ are therefore a function of a parameter $x$ and the following five invariants:
\begin{gather}
s_{12}:=p_{12}^2, \hspace{0.5 cm} s_{23}:=p_{23}^2, \hspace{0.5 cm} s_{34}:=p_{34}^2, \hspace{0.5 cm} s_{45}:=p_{45}^2=p_{123}^2, \hspace{0.5 cm} s_{51}:=p_{15}^2=p_{234}^2, \hspace{0.5 cm} p_i^2=0, \label{eq:qxpar}
\end{gather}
where the notation $p_{i\cdots j}=p_i+\cdots +p_j$ is used and $p_5:=-p_{1234}$. As the parameter $x\rightarrow 1$, the external momentum $q_3$ becomes massless, such that our parametrization (\ref{eq:qxpar}) also captures the on-shell case $M_3^2\rightarrow 0$. 

The MI depend in total on 6 variables, namely the Lorentz products $q_i.q_j$ with $i<j<5$ and the (squared) particle mass $M_3^2=q_3^2$. The $x$-parameterization of the external momenta as in Figure \ref{fig:xparam-P1} results in these variables being related to the parameter $x$ and the five independent scalar products of our choice that are defined in (\ref{eq:qxpar}). The external momenta $q_1,q_2,q_4$ and $q_5$ of the massless external particles can correspond to either of the four massless external legs in Figure \ref{fig:xparam-P1}, while the massive particle $V$ has an external momenta $q_3=p_{123}-xp_{12}$ with a mass:
\begin{gather}
M_3^2=(1 - x) (s_{45} - s_{12} x).\label{eq:Pinv}
\end{gather}

After fixing the $x$-parameterization as in Figure \ref{fig:xparam-P1}, the class of loop integrals describing the planar family%\footnote{In fact the other two planar families $P_2$ and $P_3$ also are part of the same class of integrals in (\ref{eq:P1x}).}
$P_1$ is now explicitly expressed in $x$ as:
\begin{gather}
G^{P_1}_{a_1\cdots a_{11}}(x,s,\epsilon):=e^{2\gamma_E \epsilon} \int \frac{d^dk_1}{i\pi^{d/2}}\frac{d^dk_2}{i\pi^{d/2}}
\frac{1}{k_1^{2a_1} (k_1 + x p_1)^{2a_2} (k_1 + x p_{12})^{2a_3} (k_1 + p_{123})^{2a_4}} \nonumber\\
\times \frac{1}{(k_1 + p_{1234})^{2a_5} k_2^{2a_6} (k_2 - x p_1)^{2a_7} (k_2 - x p_{12})^{2a_8} 
(k_2 - p_{123})^{2a_9} (k_2 - p_{1234})^{2a_{10}} (k_1 + k_2)^{2a_{11}}}, \label{eq:P1x}
\end{gather}
where $\gamma_E$ is the usual Euler-Mascheroni constant.

Using the notation given in Eq. (\ref{eq:P1x}), the indices $a_1\cdots a_{11}$ for the list of MI in the planar family $P_1$ is as follows\footnote{The letter m is used here to indicate the index -1.}:
{\footnotesize
\begin{eqnarray}
P_1: \hspace{-0.1cm} &\{10000000101, 01000000101, 00100000101, 10000001001,01000000011, 00100000011, 10100001100, \nonumber\\
&10100001010,10100101000, 01000101001, 10100100100, 10100000102,10100000101, 10100000011, \nonumber\\
&10000001102, 10000001101, 10000001011, 01000100101, 01000001101, 01000001011,00100100102, \nonumber\\
& 00100100101, 11100000101, 11100000011,11000001102, 11000001101, 11000001012, 11000001011, \nonumber\\
&11000000111, 10100000112, 10000001111, 01100100102, 01100100101, 01100100011, 01100000111, \nonumber\\
&01000101102,01000101101, 01000101011, 01000100111, 01000001111, 00100100111, 10100101100, \nonumber\\
&10100100101, 10100001101,10100001011, 10100000111, 111\text{m}0000111, 110000\text{m}1111, 11000001111, \nonumber\\
&10100101110, 10100100111, 10100001111,011001\text{m}0111, 01100100111, 010\text{m}0101111, 01000101111, \nonumber\\
&11100100101, 11100001101, 11100001011, 11100000111,111\text{m}0101101, 111001\text{m}1101, 11100101101, \nonumber\\
&1110\text{m}101011, 11100101011, 111\text{m}0100111, 11100100111, 111000\text{m}1111,111\text{m}0001111, 11100001111, \nonumber\\
&111001\text{m}0111, 11100101111, 111001\text{m}1111, 111\text{m}0101111\}, 
\label{eq:MIP1}
\end{eqnarray}
}
In the next section we discuss the DE method that we use to calculate the above 74 MI in $P_1$.





\section{Differential equations and their solution}
\label{sec:diff}









% description of pentabox families


%The family P1 is composed by 74 MI. The IBP reduction has been achieved by using {\bf FIRE5}~\cite{Smirnov:2008iw}. 
The resulting differential equation in matrix form can be written as 
\begin{equation}
{\partial _x}{\bf{G}} = {\bf{M}}\left( {\left\{ {{s_{ij}}} \right\},\varepsilon ,x} \right){\bf{G}}
\label{eq:ini}
\end{equation}
where ${\bf{G}}$ stands for the array of the 74 MI given in Eq. (\ref{eq:MIP1}).  
The diagonal part of the matrix at $\varepsilon=0$ defines as usual the integrating factors,  ${\left( {{M_D}} \right)_{IJ}} = {\delta _{IJ}}\,{M_{II}}\left( {\varepsilon  = 0} \right)$, and the equation takes the form
${\partial _x}{\bf{G}} = {\bf{MG}}$ with ${\bf{G}} \to  {{\bf{S}}^{ - 1}}{\bf{G}}$,  ${\bf{S}} = \exp \left( {\int {dx\,{{\bf{M}}_D}} } \right)$ and ${\bf{M}} \to {{\bf{S}}^{{\bf{ - 1}}}}\left( {{\bf{M - }}\,{{\bf{M}}_D}} \right){\bf{S}}$.

We found that, after absorbing the integrating factors, the resulting matrix $\bf{M}$ can be written as 
\begin{equation}
{M_{IJ}} = {N_{IJ}}\left( \varepsilon  \right)\left( {\sum\limits_{i = 1}^{20} {\sum\limits_{j = 1}^2 {\sum\limits_{k = 0}^1 {\frac{{{C_{IJ;ijk}}{\varepsilon ^k}}}{{{{\left( {x - {l_i}} \right)}^j}}}} } }  + \sum\limits_{j = 0}^1 {\sum\limits_{k = 0}^1 {{{\tilde C}_{IJ;jk}}{\varepsilon ^k}{x^j}} } } \right).
\label{eq:mprime}
\end{equation}
The twenty letters $l_i$, are given by
%{\scriptsize
%\begin{equation}
\begin{eqnarray}
%
%\begin{array}{l}
%\begin{align}
& 0,\,\,\,1,\,\,\,\,\,\frac{{{s_{45}}}}{{{s_{45}} - {s_{23}}}},\,\,\,\,\,\frac{{{s_{45}}}}{{{s_{12}}}},\,\,\,\,\,1 - \frac{{{s_{34}}}}{{{s_{12}}}},\,\,\,\,\,1 + \frac{{{s_{23}}}}{{{s_{12}}}}, \nonumber\\
& 1 - \frac{{{s_{34}} - {s_{51}}}}{{{s_{12}}}},\,\,\,\,\,\frac{{{s_{45}} - {s_{23}}}}{{{s_{12}}}},\,\,\,\, - \frac{{{s_{51}}}}{{{s_{12}}}},\,\,\,\,\,\frac{{{s_{45}}}}{{ - {s_{23}} + {s_{45}} + {s_{51}}}},\,\,\,\,\,\frac{{{s_{45}}}}{{{s_{34}} + {s_{45}}}},
\nonumber
\\
& \frac{{{s_{12}}{s_{23}} - 2{s_{12}}{s_{45}} - {s_{12}}{s_{51}} - {s_{23}}{s_{34}} + {s_{34}}{s_{45}} - {s_{45}}{s_{51}} \pm \sqrt {{\Delta _1}} }}{{2{s_{12}}({s_{23}} - {s_{45}} - {s_{51}})}},\,\,\,\,\,\frac{{{s_{12}}{s_{23}} - {s_{12}}{s_{45}} - {s_{12}}{s_{51}} - {s_{23}}{s_{34}} + {s_{34}}{s_{45}} - {s_{45}}{s_{51}} \pm \sqrt {{\Delta _2}} }}{{2{s_{12}}({s_{23}} - {s_{45}} - {s_{51}})}},
\nonumber
\\
& \frac{{{s_{12}}{s_{23}} - {s_{12}}{s_{51}} - {s_{23}}{s_{34}} + {s_{34}}{s_{45}} - {s_{45}}{s_{51}} \pm \sqrt {{\Delta _1}} }}{{2{s_{12}}({s_{23}} + {s_{34}} - {s_{51}})}},\,\,\,\,\,\,\,\,\,\frac{{{s_{12}}{s_{45}} \pm \sqrt {{\Delta _3}} }}{{{s_{12}}{s_{34}} + {s_{12}}{s_{45}}}},\,\,\,\,\,\,\,\,\frac{{{s_{45}}}}{{{s_{12}} + {s_{23}}}},
%\end{align}
%\end{array}
\label{eq:lett}
%
\end{eqnarray}
%\end{equation}
%}
where
%{\scriptsize
\[\begin{array}{l}
{\Delta _1} = {({s_{12}}({s_{51}} - {s_{23}}) + {s_{23}}{s_{34}} + {s_{45}}({s_{51}} - {s_{34}}))^2} + 4{s_{12}}{s_{45}}{s_{51}}({s_{23}} + {s_{34}} - {s_{51}})\\
{\Delta _2} = {({s_{{\rm{12}}}}( - {s_{{\rm{23}}}} + {s_{{\rm{45}}}} + {s_{{\rm{51}}}}) + {s_{{\rm{23}}}}{s_{{\rm{34}}}} + {s_{{\rm{45}}}}({s_{{\rm{51}}}} - {s_{{\rm{34}}}}))^2} - 4{s_{{\rm{12}}}}{s_{{\rm{45}}}}{s_{{\rm{51}}}}( - {s_{{\rm{23}}}} + {s_{{\rm{45}}}} + {s_{{\rm{51}}}})\\
{\Delta _3} =  - ({s_{12}}{s_{34}}{s_{45}}({s_{12}} - {s_{34}} - {s_{45}}))
\end{array}\]
%}
with $\Delta_1$ being the usual Gram determinant. The normalization factors $N_{IJ}\left( \varepsilon  \right)$ can be cast in the factorized form 
$ {N_{IJ}}\left( \varepsilon  \right) = {n_J} \left( \varepsilon  \right) /{n_I}\left( \varepsilon  \right) $ and can be absorbed by redefining ${{G}_I} \to {n_I}\left( \varepsilon  \right){{G}_I}$.

Although the DE can be solved from any form described so far, e.g. (\ref{eq:ini}) and (\ref{eq:mprime}) and the result can be expressed as a sum of GPs with argument $x$ and weights given by the letters in Eq. (\ref{eq:lett}), it is more elegant and easy-to-solve to derive 
a Fuchsian system of equations~\cite{Henn:2014qga}, where only single poles in the variable $x$ will appear. In fact the series of successive transformations   
\[{\bf{G}} \to \left({\bf{I}}-{\bf{K}}_i\right){\bf{G}}, \quad {\bf{M}} \to \left( {{\bf{M}} - {\partial _x}{{\bf{K}}_i} - {{\bf{K}}_i}{\bf{M}}} \right){\left( {{\bf{I}} - {{\bf{K}}_i}} \right)^{ - 1}}\,\,\,\,i = 1,2,3\]
with 
\[{\left( {{{\bf{K}}_1}} \right)_{IJ}} = \left\{ {\begin{array}{*{20}{c}}
{\int {dx} {{\left( {{\bf{M}}\left( {\varepsilon  = 0} \right)} \right)}_{IJ}}}&{\,I,J \ne 69,72,73,74}\\
0&{I,J = 69,72,73,74}
\end{array}} \right.\]
%
\[{\left( {{{\bf{K}}_2}} \right)_{IJ}} = \left\{ {\begin{array}{*{20}{c}}
{\int {dx} {{\left( {{\bf{M}}\left( {\varepsilon  = 0} \right)} \right)}_{IJ}}}&{\,I,J \ne 74}\\
0&{I,J = 74}
\end{array}} \right.\]
and 
\[{\left( {{{\bf{K}}_3}} \right)_{IJ}} = \int {dx} {\left( {{\bf{M}}\left( {\varepsilon  = 0} \right)} \right)_{IJ}}\]
with the enumeration of the MI as given by Eq. (\ref{eq:MIP1}), brings the system in to the form
\[{\partial _x}{\bf{G}} = \left( {\varepsilon \sum\limits_{i = 1}^{19} {\frac{{{{\bf{M}}_i}}}{{\left( {x - {l_i}} \right)}}} } \right){\bf{G}}\]
where the residue matrices ${{{\bf{M}}_i}}$ are independent of $x$ and $\varepsilon$. It should be noticed that the series of the above transformations {\it do not} correspond to the one described by the Moser 
algorithm~\cite{Lee:2014ioa,BarkatouPfluegel2009,BarkatouPfluegel2007,Moser1960}. The result can be straightforwardly given as 
{\scriptsize
%\[\begin{array}{l}
\begin{eqnarray}
\label{eq:solx}
{\bf{G}} &=& {\varepsilon ^{ - 2}}{\bf{b}}_0^{\left( { - 2} \right)}\\
 &+& {\varepsilon ^{ - 1}}\left( {\sum {{G_a}{{\bf{M}}_a}{\bf{b}}_0^{\left( { - 2} \right)}}  + {\bf{b}}_0^{\left( { - 1} \right)}} \right)\nonumber\\
 &+& \left( {\sum {{G_{ab}}{{\bf{M}}_a}{{\bf{M}}_b}} {\bf{b}}_0^{\left( { - 2} \right)} + \sum {{G_a}{{\bf{M}}_a}{\bf{b}}_0^{\left( { - 1} \right)}}  + {\bf{b}}_0^{\left( 0 \right)}} \right)\nonumber\\
 &+& \varepsilon \left( {\sum {{G_{abc}}{{\bf{M}}_a}{{\bf{M}}_b}{{\bf{M}}_c}} b_0^{\left( { - 2} \right)} + \sum {{G_{ab}}{{\bf{M}}_a}{{\bf{M}}_b}{\bf{b}}_0^{\left( { - 1} \right)}}  + \sum {{G_a}{{\bf{M}}_a}{\bf{b}}_0^{\left( 0 \right)} + } {\bf{b}}_0^{\left( 1 \right)}} \right)\nonumber\\
 &+& {\varepsilon ^2}\left( {\sum {{G_{abcd}}{{\bf{M}}_a}{{\bf{M}}_b}{{\bf{M}}_c}{{\bf{M}}_d}} {\bf{b}}_0^{\left( { - 2} \right)} + \sum {{G_{abc}}{{\bf{M}}_a}{{\bf{M}}_b}{{\bf{M}}_c}{\bf{b}}_0^{\left( { - 1} \right)}}  + \sum {{G_{ab}}{{\bf{M}}_a}{{\bf{M}}_b}{\bf{b}}_0^{\left( 0 \right)} + \sum {{G_a}{{\bf{M}}_a}{\bf{b}}_0^{\left( 1 \right)}}  + } {\bf{b}}_0^{\left( 2 \right)}} \right)\nonumber
\end{eqnarray}
%\end{array}\]
}with the arrays ${\bf{b}}_0^{\left( {k} \right)}$, $k=-2,...,2$ representing the $x$-independent boundary terms in the limit $x=0$ at order $\varepsilon^k$. The expression is in terms of Goncharov polylogarithms, 
${G_{{a_1} \cdots {a_n}}} = G\left( {{a_1},{a_2}, \ldots ,{a_n};x} \right)$.
The form of the solution is close to what is usually referred in the literature, as universally transcendental (UT). In fact most of the integrals do 
appear in a form that corresponds to the commonly defined UT, though others exhibit coefficients of GP's depending on the kinematics. In fact working with DE in one independent variable, the concept of UT can accommodate the dependence in the kinematics,
although it might improve numerical efficiency if the form of DE includes matrices independent of the kinematics. In this direction it is also interesting to mention that all the residue matrices $\bf{M}_a$, with a non-trivial dependence on the kinematics, do have eigenvalues that are negative integer numbers independent of the kinematics,  
related to the behaviour ${\left( {x - {l_a}} \right)^{ - {n_a}\varepsilon }}$($n_a$ positive integers) of the integrals at the corresponding limits $x\to l_a$.

The limit $x=1$ represents the solution for all planar pentabox on-shell Feynman integrals. The limit can easily obtained by properly resumming the ${\log ^k}\left( {1 - x} \right)$ terms. Interestingly enough we found a very simple formula 
for this limit given by 
\begin{equation}
{{\bf{G}}_{x = 1}} = \left( {{\bf{I}} + \frac{3}{2}{{\bf{M}}_2} + \frac{1}{2}{\bf{M}}_2^2} \right){{\bf{G}}_{trunc}}
\label{eq:limit}
\end{equation}
with ${\bf{M}}_2$ the residue matrix at $x=1$ and ${{\bf{G}}_{trunc}}$ derived from Eq.(\ref{eq:solx}), by properly removing all divergencies proportional to ${\log ^k}\left( {1 - x} \right)$ and setting $x=1$.


For the majority of the MI in the original basis (\ref{eq:MIP1}), their boundary behaviour are captured by the DE itself as was also the case for the doublebox families~\cite{Papadopoulos:2014hla}. To explain this, we turn to the language of expansion by regions~\cite{Beneke:1997zp,Smirnov:2002pj} which states that all MI can be written as an {\it expansion} in terms of the form $a_{ij}x^{i+j\epsilon}$. Since such functions are linearly independent, linear equations for the coefficients $a_{ij}$ can be found by plugging such expansions for the integrals in the DE. The boundary conditions correspond to the leading terms in the expansion as $x\rightarrow 0$ and thus they are described by the expansion terms $a_{ij}x^{i+j\epsilon}$, where $i=i_0$ is the smallest integer such that $a_{ij}$ is nonvanishing. In this way the DE themselves set constraints on the coefficients $a_{ij}$ and therefore the boundary conditions. 

In general a MI $G$ behaves at the boundary $x\rightarrow 0$ as follows~\cite{Beneke:1997zp,Smirnov:2002pj}:
\begin{equation}
G_{res}=\lim_{x\rightarrow 0}G=\sum_j c_jx^{i_0+j\epsilon}+d_jx^{i_0+1+j\epsilon}+\mathcal{O}(x^{i_0+2}),
\label{eq:exp}
\end{equation}
where in the expression $G_{res}$, the logarithms $\log^k(x)$ have been \textit{resummed} into terms of the form $\sim x^{\alpha\epsilon}$ at the boundary $x=0$. As explained above, by putting the above form (\ref{eq:exp}) for the integrals in the DE and equating the terms $x^{i+j\epsilon}$ with the same exponents on both sides of the DE, linear equations are found for the coefficients $c_i$ and $d_i$. We solved these linear equations for the coefficients in the original basis (\ref{eq:MIP1}) from the {\it bottom-up}. In other words we first solved the linear equations for the non-trivial MI with least amount of denominators and then recursively solved the coefficients for those MI with more denominators. At every step we used the solutions for the coefficients of the MI with less denominators that were found in the previous steps. The expressions for the trivial MI, i.e. those that satisfy a homogeneous DE, 
%of the form $\partial_x G_i=H\, G_i$, 
are plugged in directly and contribute to the {\it inhomogeneous} part of the linear equations for the coefficients. We note here that for the large majority of integrals we did not need to solve for the coefficients $d_i$ that correspond to $x$-suppressed terms. Their calculations were only required for those integrals whose DE had singularities of the form $x^{-2+\alpha\epsilon}$ at the boundary $x=0$ (such singularities were also encountered for the DE of the one-loop pentagon discussed in~\cite{Papadopoulos:2014lla}). Once the resummed terms in equation (\ref{eq:exp}) were calculated, the resulting DE for $G_{fin}:=G-G_{res}$ have no singularities at $x=0$ and can be directly integrated to Goncharov polylogarithms~\cite{Goncharov:1998kja,Remiddi:1999ew,Goncharov:2001iea}.

For the pentabox family $P_1$ specifically, the majority of the coefficients are fixed by these equations, while some others are not. We found in practice that for most of the integrals, the coefficients which are not fixed by the linear equations are zero and we confirmed this by the method of expansion by regions. However, for some integrals we found that the method of expansion by regions predicts that some coefficients that are {\it not} determined by the linear equations, are in fact non-zero and require an explicit calculation. As described in the Appendix \ref{expbyreg}, for those integrals we used other methods to calculate the unknown and nonvanishing coefficients. Once all boundary conditions were found for the integrals in the original basis (\ref{eq:MIP1}), the boundary conditions for the canonical basis followed directly from the relation between the two bases that is described above.

The complete expressions for all MI are available in the ancillary files~\cite{results}. The solution for all 74 MI contains $O(3,000)$ GPs which is approximately six times more than the corresponding double-box with two off-shell legs planar MI.
 We have performed several numerical checks of all our calculations. 
The numerical results, also included in the ancillary files~\cite{results}, have been performed with the {\bf GiNaC} library~\cite{Vollinga:2004sn} and compared with those provided by the numerical code {\bf SecDec}~\cite{Binoth:2000ps,Heinrich:2008si,Borowka:2012yc,Borowka:2013cma,Borowka:2015mxa} in the Euclidean region for all MI and in the physical region whenever possible
(due to CPU time limitations in using {\bf SecDec}) and found perfect agreement. For the physical region we are using the analytic continuation as described in~\cite{Papadopoulos:2014hla}.
At the present stage we are not setting a fully-fledged numerical implementation, which will be done when all families will be computed. Our experience with double-box computations show that using for instance {\tt HyperInt}~\cite{Panzer:2014caa}
 to bring all GPs in their range of convergence, {\it before} evaluating them numerically by {\bf GiNaC}, 
increases efficiency by two orders of magnitude.  Moreover expressing GPs in terms of classical polylogarithms and $Li_{2,2}$, could also reduce substantially the CPU time~\cite{Gehrmann:2015ora}. 
Based on the above we estimate that a target of  $O\left( {{{10}^2-{10}^3}} \right)$ milliseconds can be achieved.





\section{Conclusions}

In this paper we calculated, for the first time, one of the topologies of planar Master Integrals related to massless five-point amplitudes with one off-shell leg as well as the full set of massless planar Master Integrals for on-shell kinematics. 
We have demonstrated that based on the Simplified Differential Equations approach~\cite{Papadopoulos:2014lla} these MI can be expressed in terms of Goncharov polylogarithms. The complexity of the resulting expressions is certainly promising that 
the project of computing all MI relevant to massless QCD, namely all eight-denominator MI with arbitrary configuration of the external momenta, is feasible. Having such a complete library of two-loop MI, the analog of $A_0, B_0, C_0, D_0$ scalar integrals at one loop,  the reduction of an arbitrary two-loop amplitude {\`a la OPP} can pave the road  for a NNLO automation in the near future.

As experience shows, there are several issues that will need to attract our attention in order to accomplish our goal. First of all in order to systematize the whole procedure of reducing an arbitrary Feynman Integral in terms of MI in an efficient way, a deepening of our current understanding of 
IBP identities~\cite{Baikov:1996iu,Smirnov:2003kc}
is necessary~\cite{Ita:2015tya}. Secondly, further standardising the procedure to obtain a canonical form of DE~\cite{Henn:2013pwa}, which drastically simplifies the expression of MI in terms of GPs, is certainly a very desirable feature. Thirdly, the inclusion of MI with massive
internal propagators, at a first stage with one mass scale corresponding to the heavy top quark, will provide the complete basis for NNLO QCD automated computations. Moreover, the calculation of boundary terms for the DE can benefit from further 
developments and exploitations of expansion-by-regions techniques, in conjunction with Mellin-Barnes representation of the resulting integrals. Finally, on the numerical side, a more efficient computation of polylogarithms is also necessary.   




\subsection*{Acknowledgements}

This research was supported by the Research Funding Program ARISTEIA, HOCTools (co-financed by the European Union (European Social Fund ESF) and Greek national funds through the Operational Program "Education and Lifelong Learning" of the National Strategic Reference Framework (NSRF)). 
We gratefully acknowledge fruitful discussions with J.~Henn, A.~V.~Smirnov, V.~A.~Smirnov,  M.~Czakon, S.~Weinzierl, C.~Duhr and E.~Panzer during various stages of this project. We would like also to thank H.~Frellesvig for providing us with useful symbolic code to manipulate
GPs, and A.~Kardos and G.~Somogyi for cross-checking some of the integrals appearing in the expansion-by-region method to compute boundary terms. We would like also to thank S.~Borowka, E.~Panzer and A.~V.~Smirnov for helping with
computer programs {\tt SecDec}, {\tt HyperInt} and {\tt FIRE} respectively. 

\appendix
\section{Methods of calculating the boundary conditions}
\label{expbyreg}

As was mentioned in section \ref{sec:diff}, for some integrals their DE do not completely fix their behaviour as $x\rightarrow x_0=0$. The indices of these integrals are:
{\footnotesize
\begin{gather}
\{(10100000101),(10100000102),(11000001012),(11000001011),(01000101011),(10100100111), \nonumber\\
(10100001111),(111\text{m}0100111),(111000\text{m}1111),(11100001111),(111001\text{m}0111),(10100000011) \nonumber\\
(10000001011),(01000001011),(11100000011),(01100100011),(10100100111),(11100001011), \nonumber\\
(11100101011),(111\text{m}0101011)\}.
\label{eq:boundnot}
\end{gather}
}
In other words, for the above integrals there is some behaviour at $x\rightarrow 0$ which corresponds to coefficients $c_i$ and $d_i$ in equation (\ref{eq:exp}) that are not all determined by the DE itself. For the above integrals we used various methods to compute the undetermined coefficients, which we explain further below.

\paragraph{Expansion by regions}

For the following integrals the method of expansion by regions was used to fix the non-zero coefficients in the boundary behaviour (\ref{eq:exp}) which are not determined by the DE:
{\footnotesize
\begin{gather}
\{(10100000101),(10100000102),(11000001012),(11000001011),(01000101011),(10100100111), \nonumber\\
(10100001111),(111\text{m}0100111),(111000\text{m}1111),(11100001111),(111001\text{m}0111)\}.
\label{eq:boundnotexp}
\end{gather}
}
With the method of expansion by regions~\cite{Beneke:1997zp,Smirnov:2002pj}, the coefficients in equation (\ref{eq:exp}) are expressed as Feynman integrals. For most of the integrals in (\ref{eq:boundnotexp}), these integrals corresponding to the undetermined (by the DE) coefficients are reducible to known single scale integrals. However for a few, the integrals for the coefficients are non-trivial and require further calculation. For example let's consider the MI and its DE expanded around $x=0$:
\begin{gather}
G=G_{0100101011}=e^{2\gamma_E \epsilon}\int \frac{d^dk_1}{i\pi^{d/2}}\frac{d^dk_2}{i\pi^{d/2}}
\frac{1}{(k_1 + x p_1)^2 k_2^2 (k_2 - x p_{12})^2 (k_2 - p_{1234})^2 (k_1 + k_2)^2}, \nonumber\\
\partial_x(x^{1+4\epsilon}G)=B(s,\epsilon)x^{-1+2\epsilon}+\cdots, \label{eq:DEexp}
\end{gather}
where $B(s,\epsilon)$ is independent of the integral parameter $x$ and the dots do not contain any further poles at $x=0$. From expansion by regions it follows that at $x=0$ the integral behaves as $G\sim c_1x^{-1-4\epsilon}+c_2x^{-1-2\epsilon}$ which corresponds to a soft and collinear region. The coefficient $c_2$ is completely determined by the DE (\ref{eq:DEexp}) as explained in section \ref{sec:diff}, while $c_1$ is not. The coefficient $c_1$ is given by the method of regions as a Feynman integral:
\begin{equation}
c_1=e^{2\gamma_E \epsilon}\int \frac{d^dk_1}{i\pi^{d/2}}\frac{d^dk_2}{i\pi^{d/2}}
\frac{1}{(k_1 + p_1)^2 k_2^2 (k_2 - p_{12})^2 (- 2 k_2.p_{1234}) (k_1 + k_2)^2}.
\end{equation}
We can calculate the above integral by the method of DE. As we did for the $x$-parametrization of the $P_1$ family, we introduce a parameter $y$ in the denominators:
\begin{equation}
I(y):=e^{2\gamma_E \epsilon}\int \frac{d^dk_1}{i\pi^{d/2}}\frac{d^dk_2}{i\pi^{d/2}}
\frac{1}{(k_1 + y p_1)^2 k_2^2 (k_2 - p_{12})^2 (- 2 k_2.p_{1234}) (k_1 + k_2)^2}.
\end{equation}
By solving the DE $\partial_y I(y)$ and afterwards setting $y\rightarrow 1$ we could calculate $c_1=I(1)$. We note that the boundary condition for the integral $I(y)$ is completely determined by its DE. If this was not the case, we would use the method of regions again to express the missing coefficient in terms of a Feynman integral and repeat the above step, though this was never necessary.

\paragraph{Shifted boundary point}

For some other integrals we considered their limiting behaviour around another boundary point $x_0$ instead of at $x_0=0$:
{\footnotesize
\begin{equation}
\begin{array}{cc}
x_0=\infty: & \{(10100000011),(10000001011),(11100000011),(01100100011),(10100100111)\} \\
x_0=(s_{12} - s_{34} + s_{51})/s_{12}: & \{(01000001011)\}
\end{array}
\end{equation}
}
In practice, choosing the boundary at $x_0=\infty$ corresponds to performing an inverse transformation $x\rightarrow 1/x$ at the level of the DE and afterwards considering $x_0=0$. This case was already discussed in~\cite{Papadopoulos:2014hla} and therefore here we discuss the case of the MI $G:=G_{01000001011}$. Its DE is of the form:
\begin{gather}
\partial_x\left(\left(1 - \frac{s_{12} x}{s_{12} - s_{34}}\right)^{\epsilon} \left(1 - \frac{s_{12} x}{s_{12} - s_{34} + s_{51}}\right)^{\epsilon}x^{2\epsilon}G\right) \nonumber\\
=C(s,\epsilon)\left(1-\frac{s_{12} x}{s_{12}-s_{34}}\right)^{\epsilon -1} \left(1-\frac{s_{12} x}{s_{12}-s_{34}+s_{51}}\right)^{\epsilon -1},
\label{eq:bound2}
\end{gather}
where $C(s,\epsilon)$ is a factor independent of the integral parameter $x$. By using the method of expansion by regions, it can be checked that $(1 - \frac{s_{12} x}{s_{12} - s_{34}})^{\epsilon} (1 - \frac{s_{12} x}{s_{12} - s_{34} + s_{51}})^{\epsilon}x^{2\epsilon}G\rightarrow 0$ as $x\rightarrow (s_{12} - s_{34} + s_{51})/s_{12}$. The right hand side of equation (\ref{eq:bound2}) has a singularity at $x=(s_{12}-s_{34}+s_{51})/s_{12}$, which exactly captures the behaviour of the integral $G$ around $x=(s_{12}-s_{34}+s_{51})/s_{12}$. Therefore by the transformation $x\rightarrow \frac{s_{12} - s_{34} + s_{51}}{s_{12}}(1-x')$ and then afterwards integrating from $x'=0$ one can compute $G$ as a function of $x'$. Upon tranforming back to $x$, the limiting behaviour of $G$ at $x=0$ and therefore its corresponding boundary condition is found.

\paragraph{Extraction from known integrals}

The last three integrals in (\ref{eq:boundnot}) correspond to taking the $x\rightarrow 1$ limit of other known integrals (that lie in the same family $P_1$) and afterwards redefining the invariants as follows:
{\small
\begin{gather}
G_{11100001011}(x,s_{12},s_{34},s_{51})=G_{11100100101}(x'=1,s_{12}',s_{23}',s_{45}'), \nonumber\\
G_{11100101011}(x,s_{12},s_{34},s_{51})=G_{11100101101}(x'=1,s_{12}',s_{23}',s_{45}'), \nonumber\\
G_{111\text{m}0101011}(x,s_{12},s_{34},s_{51})=G_{111\text{m}0101101}(x'=1,s_{12}',s_{23}',s_{45}'), \nonumber\\
s_{12}'=x^2s_{12}, \hspace{0.5 cm} s_{23}'=x s_{51}, \hspace{0.5 cm} s_{45}'=-xs_{12} + xs_{34} + x^2s_{12}.
\label{eq:limits}
\end{gather}
}
The three integrals on the right hand side of equation (\ref{eq:limits}) are MI that were previously calculated in~\cite{Papadopoulos:2014hla}. From the exact result in $x$ for the three integrals on the left hand side in (\ref{eq:limits}) we could then compute their corresponding non-zero coefficients in (\ref{eq:exp}) that are not determined\footnote{For each of the three integrals, there was exactly one coefficient which was not fixed by the DE as explained in section \ref{sec:diff}.} by the DE.


\bibliographystyle{JHEP}
\bibliography{pentabox_paper}

\end{document}