text:
\section{Master Integrals for the two-loop QED vertices}
\label{results:vertex2L}
In this Appendix we collect   the  $17$ MI's of the
two-loop QED vertices introduced  in Eq.~(\ref{Eq:MIfinVertex}). In Section~\ref{Sec:vertex}, we have obtained them starting from the integrals $\top{i}$ depicted
in Fig.~\ref{Fig:Masters2}, which are normalized according to the integration measure (Minkowskian metric is understood)
\begin{displaymath}
\left ( \frac{m^{2 \epsilon}}{\Gamma(1+\epsilon)} \right )^2 \, \int \,
\frac{d^D k_1}{\pi^{D/2}} \, \int \,   \frac{d^D k_2 }{\pi^{D/2}} \, .
\end{displaymath}
The MI's exhibit uniform transcendentality. In the following  we present the expression of the coefficients
of their  expansion around $\epsilon=0$ up to $\mathcal{O}(\epsilon^4)$.
The coefficients  $g_{i}^{(a)}$ are defined as follows:
\begin{displaymath}
g_i = \sum_{a=0}^{4} \epsilon^a \ g_{i}^{(a)} \ , \qquad i=1,\ldots,17  \, .
\end{displaymath}

\begin{subequations}
\begin{align}
                   g_{1}^{(0)}   = {} &
                         -1\,,\\ \allowdisplaybreaks[1]
                   g_{1}^{(1)}   = {} &
                          0 \,,\\ \allowdisplaybreaks[1]
                   g_{1}^{(2)}   = {} &
                          0 \,,\\ \allowdisplaybreaks[1]
                   g_{1}^{(3)}   = {} &
                          0 \,,\\ \allowdisplaybreaks[1]
                   g_{1}^{(4)}   = {} &
                          0 \,,
\end{align}
\end{subequations}

\begin{subequations}
\begin{align}
                   g_{2}^{(0)}   = {} &
                          0 \,,\\ \allowdisplaybreaks[1]
                   g_{2}^{(1)}   = {} &
                         -\,\text{H}(0; x)\,,\\ \allowdisplaybreaks[1]
                   g_{2}^{(2)}   = {} &
                         2 \,\text{H}(-1,0; x)-\,\text{H}(0,0; x) + \zeta_2\,,\\ \allowdisplaybreaks[1]
                   g_{2}^{(3)}   = {} &
                         -4 \,\text{H}(-1,-1,0; x)+2 \,\text{H}(-1,0,0; x)+2 \,\text{H}(0,-1,0; x)\nonumber \\ \allowdisplaybreaks[1]
                         &-\,\text{H}(0,0,0; x) + \,\zeta_2 (\,\text{H}(0; x)-2 \,\text{H}(-1; x)) + 2 \,\zeta_3\,,\\ \allowdisplaybreaks[1]
                   g_{2}^{(4)}   = {} &
                         8 \,\text{H}(-1,-1,-1,0; x)-4 \,\text{H}(-1,-1,0,0; x)-4 \,\text{H}(-1,0,-1,0; x)\nonumber \\ \allowdisplaybreaks[1]
                         &  +2 \,\text{H}(-1,0,0,0; x)-4 \,\text{H}(0,-1,-1,0; x)+2 \,\text{H}(0,-1,0,0; x)\nonumber \\ \allowdisplaybreaks[1]
                         & +2 \,\text{H}(0,0,-1,0; x)-\,\text{H}(0,0,0,0; x)
                         + \zeta_2 (4 \,\text{H}(-1,-1; x) \nonumber \\ \allowdisplaybreaks[1]
                         & -2 \,\text{H}(-1,0; x)-2 \,\text{H}(0,-1; x)+\,\text{H}(0,0; x))  \nonumber \\ \allowdisplaybreaks[1]
                         & -2 \,\zeta_3 (2 \,\text{H}(-1; x)-\,\text{H}(0; x)) + \frac{9 \,\zeta _4}{4}\,,
\end{align}
\end{subequations}

\begin{subequations}
\begin{align}
                   g_{3}^{(0)}   = {} &
                          0 \,,\\ \allowdisplaybreaks[1]
                   g_{3}^{(1)}   = {} &
                          0 \,,\\ \allowdisplaybreaks[1]
                   g_{3}^{(2)}   = {} &
                         -2 \,\text{H}(0,0; x)\,,\\ \allowdisplaybreaks[1]
                   g_{3}^{(3)}   = {} &
                         8 \,\text{H}(-1,0,0; x)+4 \,\text{H}(0,-1,0; x)-6 \,\text{H}(0,0,0; x) + 2 \,\zeta_2 \,\text{H}(0; x)\,,\\ \allowdisplaybreaks[1]
                   g_{3}^{(4)}   = {} &
                         -32 \,\text{H}(-1,-1,0,0; x)-16 \,\text{H}(-1,0,-1,0; x)+24 \,\text{H}(-1,0,0,0; x) \nonumber \\ \allowdisplaybreaks[1]
                         &-8 \,\text{H}(0,-1,-1,0; x)+20 \,\text{H}(0,-1,0,0; x)+12 \,\text{H}(0,0,-1,0; x) \nonumber \\ \allowdisplaybreaks[1]
                         &-14 \,\text{H}(0,0,0,0; x)
                          - 2 \,\zeta_2 (4 \,\text{H}(-1,0; x)+2 \,\text{H}(0,-1; x)-3 \,\text{H}(0,0; x))  \nonumber \\ \allowdisplaybreaks[1]
                         &+4 \,\zeta_3 \,\text{H}(0; x) - \frac{5 \,\zeta _4}{2}\,,
\end{align}
\end{subequations}

\begin{subequations}
\begin{align}
                   g_{4}^{(0)}   = {} &
                         \frac{1}{4}\,,\\ \allowdisplaybreaks[1]
                   g_{4}^{(1)}   = {} &
                          0 \,,\\ \allowdisplaybreaks[1]
                   g_{4}^{(2)}   = {} &
                         \zeta_2\,,\\ \allowdisplaybreaks[1]
                   g_{4}^{(3)}   = {} &
                         2 \,\zeta_3\,,\\ \allowdisplaybreaks[1]
                   g_{4}^{(4)}   = {} &
                         16 \,\zeta _4\,,
\end{align}
\end{subequations}

\begin{subequations}
\begin{align}
                   g_{5}^{(0)}   = {} &
                          0 \,,\\ \allowdisplaybreaks[1]
                   g_{5}^{(1)}   = {} &
                          \text{H}(0; x)\,,\\ \allowdisplaybreaks[1]
                   g_{5}^{(2)}   = {} &
                         -6 \,\text{H}(-1,0; x)+5 \,\text{H}(0,0; x)+2 \,\text{H}(1,0; x) - \zeta_2\,,\\ \allowdisplaybreaks[1]
                   g_{5}^{(3)}   = {} &
                          36 \,\text{H}(-1,-1,0; x)-24 \,\text{H}(-1,0,0; x)-12 \,\text{H}(-1,1,0; x) \nonumber \\ \allowdisplaybreaks[1]
                         & -30 \,\text{H}(0,-1,0; x)+13 \,\text{H}(0,0,0; x)+10 \,\text{H}(0,1,0; x) \nonumber \\ \allowdisplaybreaks[1]
                         & -12 \,\text{H}(1,-1,0; x)+6 \,\text{H}(1,0,0; x)+4 \,\text{H}(1,1,0; x) \nonumber \\ \allowdisplaybreaks[1]
                         & +\zeta_2 (6 \,\text{H}(-1; x)-5 \,\text{H}(0; x)-2 \,\text{H}(1; x)) - 14 \,\zeta_3\,,\\ \allowdisplaybreaks[1]
                   g_{5}^{(4)}   = {} &
                          -216 \,\text{H}(-1,-1,-1,0; x)+144 \,\text{H}(-1,-1,0,0; x) \nonumber \\ \allowdisplaybreaks[1]
                         &+72 \,\text{H}(-1,-1,1,0; x)+144 \,\text{H}(-1,0,-1,0; x)-60 \,\text{H}(-1,0,0,0; x) \nonumber \\ \allowdisplaybreaks[1]
                         &-48 \,\text{H}(-1,0,1,0; x)+72 \,\text{H}(-1,1,-1,0; x)-48 \,\text{H}(-1,1,0,0; x) \nonumber \\ \allowdisplaybreaks[1]
                         &-24 \,\text{H}(-1,1,1,0; x)+180 \,\text{H}(0,-1,-1,0; x)-120 \,\text{H}(0,-1,0,0; x) \nonumber \\ \allowdisplaybreaks[1]
                         &-60 \,\text{H}(0,-1,1,0; x)-78 \,\text{H}(0,0,-1,0; x)+29 \,\text{H}(0,0,0,0; x) \nonumber \\ \allowdisplaybreaks[1]
                         &+26 \,\text{H}(0,0,1,0; x)-60 \,\text{H}(0,1,-1,0; x)+54 \,\text{H}(0,1,0,0; x) \nonumber \\ \allowdisplaybreaks[1]
                         &+20 \,\text{H}(0,1,1,0; x)+72 \,\text{H}(1,-1,-1,0; x)-48 \,\text{H}(1,-1,0,0; x) \nonumber \\ \allowdisplaybreaks[1]
                         &-24 \,\text{H}(1,-1,1,0; x)-36 \,\text{H}(1,0,-1,0; x)+14 \,\text{H}(1,0,0,0; x) \nonumber \\ \allowdisplaybreaks[1]
                         &+12 \,\text{H}(1,0,1,0; x)-24 \,\text{H}(1,1,-1,0; x)+20 \,\text{H}(1,1,0,0; x) \nonumber \\ \allowdisplaybreaks[1]
                         &+8 \,\text{H}(1,1,1,0; x) + \zeta_2 (-36 \,\text{H}(-1,-1; x)+24 \,\text{H}(-1,0; x) \nonumber \\ \allowdisplaybreaks[1]
                         &+12 \,\text{H}(-1,1; x)+30 \,\text{H}(0,-1; x)-13 \,\text{H}(0,0; x)-10 \,\text{H}(0,1; x) \nonumber \\ \allowdisplaybreaks[1]
                         &+12 \,\text{H}(1,-1; x)-6 \,\text{H}(1,0; x)-4 \,\text{H}(1,1; x))
                         + 2 \,\zeta_3 (33 \,\text{H}(-1; x) \nonumber \\ \allowdisplaybreaks[1]
                         &-17 \,\text{H}(0; x)-8 \,\text{H}(1; x)) - \frac{61 \,\zeta _4}{4}\,,
\end{align}
\end{subequations}

\begin{subequations}
\begin{align}
                   g_{6}^{(0)}   = {} &
                          0 \,,\\ \allowdisplaybreaks[1]
                   g_{6}^{(1)}   = {} &
                          0 \,,\\ \allowdisplaybreaks[1]
                   g_{6}^{(2)}   = {} &
                         2 \,\text{H}(0,0; x)\,,\\ \allowdisplaybreaks[1]
                   g_{6}^{(3)}   = {} &
                         -12 \,\text{H}(0,-1,0; x)+6 \,\text{H}(0,0,0; x)+4 \,\text{H}(0,1,0; x)-4 \,\text{H}(1,0,0; x) \nonumber \\ \allowdisplaybreaks[1]
                         &-2 \,\zeta_2 \,\text{H}(0; x) + \nonumber \\ \allowdisplaybreaks[1]
                         &-6 \,\zeta_3\,,\\ \allowdisplaybreaks[1]
                   g_{6}^{(4)}   = {} &
                         72 \,\text{H}(0,-1,-1,0; x)-48 \,\text{H}(0,-1,0,0; x)-24 \,\text{H}(0,-1,1,0; x) \nonumber \\ \allowdisplaybreaks[1]
                      & -36 \,\text{H}(0,0,-1,0; x)+14 \,\text{H}(0,0,0,0; x)+12 \,\text{H}(0,0,1,0; x) \nonumber \\ \allowdisplaybreaks[1]
                      & -24 \,\text{H}(0,1,-1,0; x)+20 \,\text{H}(0,1,0,0; x)+8 \,\text{H}(0,1,1,0; x) \nonumber \\ \allowdisplaybreaks[1]
                      & +24 \,\text{H}(1,0,-1,0; x)-12 \,\text{H}(1,0,0,0; x)-8 \,\text{H}(1,0,1,0; x) \nonumber \\ \allowdisplaybreaks[1]
                      & +8 \,\text{H}(1,1,0,0; x) + 2 \,\zeta_2 (6 \,\text{H}(0,-1; x)-3 \,\text{H}(0,0; x) \nonumber \\ \allowdisplaybreaks[1]
                      & -2 \,\text{H}(0,1; x)+2 \,\text{H}(1,0; x)) - 4 \,\zeta_3 (4 \,\text{H}(0; x)-3 \,\text{H}(1; x))
                         -\frac{13 \,\zeta _4}{2}\,,
\end{align}
\end{subequations}

\begin{subequations}
\begin{align}
                   g_{7}^{(0)}   = {} &
                          0 \,,\\ \allowdisplaybreaks[1]
                   g_{7}^{(1)}   = {} &
                          0 \,,\\ \allowdisplaybreaks[1]
                   g_{7}^{(2)}   = {} &
                         \frac{\,\zeta_2}{2}\,,\\ \allowdisplaybreaks[1]
                   g_{7}^{(3)}   = {} &
                         -3 \,\zeta_2 \log{2} + \frac{7 \,\zeta_3}{4}\,,\\ \allowdisplaybreaks[1]
                   g_{7}^{(4)}   = {} &
                         \frac{1}{2} \left(24 \,\text{Li}_4{\frac12}+\log^4{2}\right)
                         + 6 \,\zeta_2 \log^2{2} - \frac{31 \,\zeta _4}{4}\,,
\end{align}
\end{subequations}

\begin{subequations}
\begin{align}
                   g_{8}^{(0)}   = {} &
                          0 \,,\\ \allowdisplaybreaks[1]
                   g_{8}^{(1)}   = {} &
                          0 \,,\\ \allowdisplaybreaks[1]
                   g_{8}^{(2)}   = {} &
                          0 \,,\\ \allowdisplaybreaks[1]
                   g_{8}^{(3)}   = {} &
                         -4 \,\text{H}(0,0,0; x) - 4 \,\zeta_2 \,\text{H}(0; x)\,,\\ \allowdisplaybreaks[1]
                   g_{8}^{(4)}   = {} &
                         -8 \,\text{H}(-1,0,0,0; x)+24 \,\text{H}(0,0,-1,0; x)-4 \,\text{H}(0,0,0,0; x) \nonumber \\ \allowdisplaybreaks[1]
                         &-8 \,\text{H}(0,0,1,0; x)+8 \,\text{H}(0,1,0,0; x)+8 \,\text{H}(1,0,0,0; x) \nonumber \\ \allowdisplaybreaks[1]
                         &-4 \,\zeta_2 (2 \,\text{H}(-1,0; x)-3 \,\text{H}(0,0; x)-2 \,\text{H}(1,0; x)) \nonumber \\ \allowdisplaybreaks[1]
                         &+4 \,\zeta_3 \,\text{H}(0; x)  + 26 \,\zeta _4\,,
\end{align}
\end{subequations}

\begin{subequations}
\begin{align}
                   g_{9}^{(0)}   = {} &
                          0 \,,\\ \allowdisplaybreaks[1]
                   g_{9}^{(1)}   = {} &
                         -\frac{1}{2} \,\text{H}(0; x)\,,\\ \allowdisplaybreaks[1]
                   g_{9}^{(2)}   = {} &
                         2 \,\text{H}(-1,0; x)-\,\text{H}(0,0; x) + \zeta_2\,,\\ \allowdisplaybreaks[1]
                   g_{9}^{(3)}   = {} &
                         -8 \,\text{H}(-1,-1,0; x)+4 \,\text{H}(-1,0,0; x)+4 \,\text{H}(0,-1,0; x) \nonumber \\ \allowdisplaybreaks[1]
                         &-2 \,\text{H}(0,0,0; x) - 4 \,\zeta_2 \,\text{H}(-1; x) + 4 \,\zeta_3\,,\\ \allowdisplaybreaks[1]
                   g_{9}^{(4)}   = {} &
                         32 \,\text{H}(-1,-1,-1,0; x)-16 \,\text{H}(-1,-1,0,0; x)-16 \,\text{H}(-1,0,-1,0; x) \nonumber \\ \allowdisplaybreaks[1]
                         & +8 \,\text{H}(-1,0,0,0; x)-16 \,\text{H}(0,-1,-1,0; x)+8 \,\text{H}(0,-1,0,0; x) \nonumber \\ \allowdisplaybreaks[1]
                         & +8 \,\text{H}(0,0,-1,0; x)-4 \,\text{H}(0,0,0,0; x) + 8 \,\zeta_2 (2 \,\text{H}(-1,-1; x)\nonumber \\ \allowdisplaybreaks[1]
                         & -\,\text{H}(0,-1; x)) - 4 \,\zeta_3 (4 \,\text{H}(-1; x)-\,\text{H}(0; x)) + 19 \,\zeta _4\,,
\end{align}
\end{subequations}

\begin{subequations}
\begin{align}
                   g_{10}^{(0)}   = {} &
                          0 \,,\\ \allowdisplaybreaks[1]
                   g_{10}^{(1)}   = {} &
                         \frac{1}{2} \,\text{H}(0; x)\,,\\ \allowdisplaybreaks[1]
                   g_{10}^{(2)}   = {} &
                         -3 \,\text{H}(-1,0; x)+\frac{5}{2} \,\text{H}(0,0; x)+\,\text{H}(1,0; x) +\zeta_2\,,\\ \allowdisplaybreaks[1]
                   g_{10}^{(3)}   = {} &
                         18 \,\text{H}(-1,-1,0; x)-14 \,\text{H}(-1,0,0; x)-6 \,\text{H}(-1,1,0; x)\nonumber \\ \allowdisplaybreaks[1]
                         & -15 \,\text{H}(0,-1,0; x)+\frac{17}{2} \,\text{H}(0,0,0; x)+5 \,\text{H}(0,1,0; x)\nonumber \\ \allowdisplaybreaks[1]
                         & -6 \,\text{H}(1,-1,0; x)+5 \,\text{H}(1,0,0; x)+2 \,\text{H}(1,1,0; x) \nonumber \\ \allowdisplaybreaks[1]
                         & +\frac{1}{2} \,\zeta_2 (-6 \,\text{H}(-1; x)+\,\text{H}(0; x)-2 \,\text{H}(1; x)-6 \log{2})
                         -\frac{9 \,\zeta_3}{4}\,,\\ \allowdisplaybreaks[1]
                   g_{10}^{(4)}   = {} &
                         -108 \,\text{H}(-1,-1,-1,0; x)+80 \,\text{H}(-1,-1,0,0; x)\nonumber \\ \allowdisplaybreaks[1]
                         &+36 \,\text{H}(-1,-1,1,0; x)+84 \,\text{H}(-1,0,-1,0; x)-44 \,\text{H}(-1,0,0,0; x)\nonumber \\ \allowdisplaybreaks[1]
                         &-28 \,\text{H}(-1,0,1,0; x)+36 \,\text{H}(-1,1,-1,0; x)-28 \,\text{H}(-1,1,0,0; x)\nonumber \\ \allowdisplaybreaks[1]
                         &-12 \,\text{H}(-1,1,1,0; x)+90 \,\text{H}(0,-1,-1,0; x)-66 \,\text{H}(0,-1,0,0; x)\nonumber \\ \allowdisplaybreaks[1]
                         &-30 \,\text{H}(0,-1,1,0; x)-51 \,\text{H}(0,0,-1,0; x)+\frac{41}{2} \,\text{H}(0,0,0,0; x)\nonumber \\ \allowdisplaybreaks[1]
                         &+17 \,\text{H}(0,0,1,0; x)-30 \,\text{H}(0,1,-1,0; x)+29 \,\text{H}(0,1,0,0; x)\nonumber \\ \allowdisplaybreaks[1]
                         &+10 \,\text{H}(0,1,1,0; x)+36 \,\text{H}(1,-1,-1,0; x)-28 \,\text{H}(1,-1,0,0; x)\nonumber \\ \allowdisplaybreaks[1]
                         &-12 \,\text{H}(1,-1,1,0; x)-30 \,\text{H}(1,0,-1,0; x)+17 \,\text{H}(1,0,0,0; x)\nonumber \\ \allowdisplaybreaks[1]
                         &+10 \,\text{H}(1,0,1,0; x)-12 \,\text{H}(1,1,-1,0; x)+10 \,\text{H}(1,1,0,0; x)\nonumber \\ \allowdisplaybreaks[1]
                         &+4 \,\text{H}(1,1,1,0; x)+12 \,\text{Li}_4{\frac12}+\frac{\log^4{2}}{2}
                         +\frac{1}{2} \,\zeta_2 \left( 24 \log{2} \,\text{H}(-1; x) \right. \nonumber \\ \allowdisplaybreaks[1]
                         &+24 \log{2} \,\text{H}(1; x) +12 \,\text{H}(-1,-1; x)+4 \,\text{H}(-1,0; x)
                         +12 \,\text{H}(-1,1; x) \nonumber \\ \allowdisplaybreaks[1]
                         &-6 \,\text{H}(0,-1; x)-11 \,\text{H}(0,0; x)-10 \,\text{H}(0,1; x)-12 \,\text{H}(1,-1; x)\nonumber \\ \allowdisplaybreaks[1]
                         & + 2 \left.\text{H}(1,0; x)-4 \,\text{H}(1,1; x)+12 \log^2{2}\right)
                         + \zeta_3 (20 \,\text{H}(-1; x)\nonumber \\ \allowdisplaybreaks[1]
                         &-14 \,\text{H}(0; x)-15 \,\text{H}(1; x)) - \frac{95 \,\zeta _4}{8}\,,
\end{align}
\end{subequations}

\begin{subequations}
\begin{align}
                   g_{11}^{(0)}   = {} &
                          0 \,,\\ \allowdisplaybreaks[1]
                   g_{11}^{(1)}   = {} &
                          0 \,,\\ \allowdisplaybreaks[1]
                   g_{11}^{(2)}   = {} &
                          0 \,,\\ \allowdisplaybreaks[1]
                   g_{11}^{(3)}   = {} &
                         -2 \,\text{H}(0,0,0; x) -2 \,\zeta_2 \,\text{H}(0; x)\,,\\ \allowdisplaybreaks[1]
                   g_{11}^{(4)}   = {} &
                         -4 \,\text{H}(-1,0,0,0; x)+4 \,\text{H}(0,-1,0,0; x)+12 \,\text{H}(0,0,-1,0; x)\nonumber \\ \allowdisplaybreaks[1]
                         &-6 \,\text{H}(0,0,0,0; x)-4 \,\text{H}(0,0,1,0; x)+4 \,\text{H}(1,0,0,0; x)  \nonumber \\ \allowdisplaybreaks[1]
                         &-4 \,\zeta_2 (\,\text{H}(-1,0; x)-3 \,\text{H}(0,-1; x)-\,\text{H}(1,0; x)) - \frac{\,\zeta _4}{2}\,,
\end{align}
\end{subequations}

\begin{subequations}
\begin{align}
                   g_{12}^{(0)}   = {} &
                          0 \,,\\ \allowdisplaybreaks[1]
                   g_{12}^{(1)}   = {} &
                          0 \,,\\ \allowdisplaybreaks[1]
                   g_{12}^{(2)}   = {} &
                          0 \,,\\ \allowdisplaybreaks[1]
                   g_{12}^{(3)}   = {} &
                         -\,\text{H}(0,0,0; x) - \zeta_2 \text{H}(0; x) \,,\\ \allowdisplaybreaks[1]
                   g_{12}^{(4)}   = {} &
                         -2 \,\text{H}(-1,0,0,0; x)+2 \,\text{H}(0,-1,0,0; x)+2 \,\text{H}(0,0,-1,0; x)\nonumber \\ \allowdisplaybreaks[1]
                         &-3 \,\text{H}(0,0,0,0; x)-4 \,\text{H}(0,1,0,0; x) + \zeta_2 (-2 \,\text{H}(-1,0; x)\nonumber \\ \allowdisplaybreaks[1]
                         &+6 \,\text{H}(0,-1; x)-\,\text{H}(0,0; x)) + 2 \,\zeta_3 \,\text{H}(0; x) + \frac{\,\zeta _4}{4}\,,
\end{align}
\end{subequations}

\begin{subequations}
\begin{align}
                   g_{13}^{(0)}   = {} &
                          0 \,,\\ \allowdisplaybreaks[1]
                   g_{13}^{(1)}   = {} &
                          0 \,,\\ \allowdisplaybreaks[1]
                   g_{13}^{(2)}   = {} &
                         \text{H}(0,0; x) + \frac{3 \,\zeta_2}{2}\,,\\ \allowdisplaybreaks[1]
                   g_{13}^{(3)}   = {} &
                         -2 \,\text{H}(-1,0,0; x)-2 \,\text{H}(0,-1,0; x)+4 \,\text{H}(0,0,0; x)+4 \,\text{H}(1,0,0; x) \nonumber \\ \allowdisplaybreaks[1]
                         &+\,\zeta_2 (-6 \,\text{H}(-1; x)+2 \,\text{H}(0; x)-3 \log{2}) - \frac{\,\zeta_3}{4}\,,\\ \allowdisplaybreaks[1]
                   g_{13}^{(4)}   = {} &
                         4 \,\text{H}(-1,-1,0,0; x)+4 \,\text{H}(-1,0,-1,0; x)-8 \,\text{H}(-1,0,0,0; x)\nonumber \\ \allowdisplaybreaks[1]
                         & -8 \,\text{H}(-1,1,0,0; x)+4 \,\text{H}(0,-1,-1,0; x)-8 \,\text{H}(0,-1,0,0; x)\nonumber \\ \allowdisplaybreaks[1]
                         & -8 \,\text{H}(0,0,-1,0; x)+10 \,\text{H}(0,0,0,0; x)+12 \,\text{H}(0,1,0,0; x)\nonumber \\ \allowdisplaybreaks[1]
                         & -8 \,\text{H}(1,-1,0,0; x)-8 \,\text{H}(1,0,-1,0; x)+16 \,\text{H}(1,0,0,0; x)\nonumber \\ \allowdisplaybreaks[1]
                         & +16 \,\text{H}(1,1,0,0; x)+12 \,\text{Li}_4{\frac12}+\frac{\log^4{2}}{2}
                         + 2 \,\zeta_2 \left(12 \log{2} \,\text{H}(-1; x) \right. \nonumber \\ \allowdisplaybreaks[1]
                         & +12 \log{2} \,\text{H}(1; x)+6 \,\text{H}(-1,-1; x)-2 \,\text{H}(-1,0; x)-8 \,\text{H}(0,-1; x)\nonumber \\ \allowdisplaybreaks[1]
                         & \left. +\,\text{H}(0,0; x)-12 \,\text{H}(1,-1; x)+4 \,\text{H}(1,0; x)+3 \log^2{2}\right)  \nonumber \\ \allowdisplaybreaks[1]
                         & -2 \,\zeta_3 (5 \,\text{H}(-1; x)+4 \,\text{H}(0; x)+11 \,\text{H}(1; x)) -\frac{47 \,\zeta _4}{4}\,,
\end{align}
\end{subequations}

\begin{subequations}
\begin{align}
                   g_{14}^{(0)}   = {} &
                          0 \,,\\ \allowdisplaybreaks[1]
                   g_{14}^{(1)}   = {} &
                          0 \,,\\ \allowdisplaybreaks[1]
                   g_{14}^{(2)}   = {} &
                         \text{H}(0,0; x)\,,\\ \allowdisplaybreaks[1]
                   g_{14}^{(3)}   = {} &
                         -4 \,\text{H}(-1,0,0; x)-4 \,\text{H}(0,-1,0; x)+5 \,\text{H}(0,0,0; x)\nonumber \\ \allowdisplaybreaks[1]
                         &+2 \,\text{H}(0,1,0; x) + \zeta_3\,,\\ \allowdisplaybreaks[1]
                   g_{14}^{(4)}   = {} &
                         16 \,\text{H}(-1,-1,0,0; x)+16 \,\text{H}(-1,0,-1,0; x)-20 \,\text{H}(-1,0,0,0; x)\nonumber \\ \allowdisplaybreaks[1]
                         & -8 \,\text{H}(-1,0,1,0; x)+24 \,\text{H}(0,-1,-1,0; x)-26 \,\text{H}(0,-1,0,0; x)\nonumber \\ \allowdisplaybreaks[1]
                         & -12 \,\text{H}(0,-1,1,0; x)-26 \,\text{H}(0,0,-1,0; x)+9 \,\text{H}(0,0,0,0; x)\nonumber \\ \allowdisplaybreaks[1]
                         & +12 \,\text{H}(0,0,1,0; x)-12 \,\text{H}(0,1,-1,0; x)+8 \,\text{H}(0,1,0,0; x)\nonumber \\ \allowdisplaybreaks[1]
                         & +4 \,\text{H}(0,1,1,0; x) - \,\zeta_2 (13 \,\text{H}(0,0; x)+2 \,\text{H}(0,1; x)) \nonumber \\ \allowdisplaybreaks[1]
                         & -\zeta_3 (4 \,\text{H}(-1; x)+3 \,\text{H}(0; x)) - \frac{7 \,\zeta _4}{2}\,,
\end{align}
\end{subequations}

\begin{subequations}
\begin{align}
                   g_{15}^{(0)}   = {} &
                          0 \,,\\ \allowdisplaybreaks[1]
                   g_{15}^{(1)}   = {} &
                          0 \,,\\ \allowdisplaybreaks[1]
                   g_{15}^{(2)}   = {} &
                          0 \,,\\ \allowdisplaybreaks[1]
                   g_{15}^{(3)}   = {} &
                          0 \,,\\ \allowdisplaybreaks[1]
                   g_{15}^{(4)}   = {} &
                         4 \,\text{H}(0,-1,0,0; x)-2 \,\text{H}(0,0,-1,0; x)-2 \,\text{H}(0,1,0,0; x)\nonumber \\ \allowdisplaybreaks[1]
                         & +4 \,\text{H}(1,0,0,0; x) + \zeta_2 (\,\text{H}(0,0; x)+4 \,\text{H}(1,0; x)) \nonumber \\ \allowdisplaybreaks[1]
                         & -4 \,\zeta_3 \,\text{H}(0; x) + \frac{17 \,\zeta _4}{4}\,,
\end{align}
\end{subequations}

\begin{subequations}
\begin{align}
                   g_{16}^{(0)}   = {} &
                          0 \,,\\ \allowdisplaybreaks[1]
                   g_{16}^{(1)}   = {} &
                          0 \,,\\ \allowdisplaybreaks[1]
                   g_{16}^{(2)}   = {} &
                          0 \,,\\ \allowdisplaybreaks[1]
                   g_{16}^{(3)}   = {} &
                          0 \,,\\ \allowdisplaybreaks[1]
                   g_{16}^{(4)}   = {} &
                         -4 \,\text{H}(0,-1,0,0; x)+4 \,\text{H}(0,0,-1,0; x)-2 \,\text{H}(0,0,0,0; x)\nonumber \\ \allowdisplaybreaks[1]
                         &-4 \,\text{H}(0,0,1,0; x)+4 \,\text{H}(0,1,0,0; x)-4 \,\text{H}(1,0,0,0; x) \nonumber \\ \allowdisplaybreaks[1]
                         &-2 \,\zeta_2 (6 \,\text{H}(0,-1; x)-\,\text{H}(0,0; x)+2 \,\text{H}(1,0; x)) - 2 \,\zeta _4\,,
\end{align}
\end{subequations}

\begin{subequations}
\begin{align}
                   g_{17}^{(0)}   = {} &
                          0 \,,\\ \allowdisplaybreaks[1]
                   g_{17}^{(1)}   = {} &
                          0 \,,\\ \allowdisplaybreaks[1]
                   g_{17}^{(2)}   = {} &
                          0 \,,\\ \allowdisplaybreaks[1]
                   g_{17}^{(3)}   = {} &
                         2 \,(\text{H}(0,-1,0; x)-\,\text{H}(0,0,0; x)-\,\text{H}(0,1,0; x))
                         - \zeta_2 \,\text{H}(0; x) -\,\zeta_3\,,\\ \allowdisplaybreaks[1]
                   g_{17}^{(4)}   = {} &
                         -8 \,\text{H}(-1,0,-1,0; x)+8 \,\text{H}(-1,0,0,0; x)+8 \,\text{H}(-1,0,1,0; x)\nonumber \\ \allowdisplaybreaks[1]
                         &-20 \,\text{H}(0,-1,-1,0; x)+16 \,\text{H}(0,-1,0,0; x)+12 \,\text{H}(0,-1,1,0; x)\nonumber \\ \allowdisplaybreaks[1]
                         &+24 \,\text{H}(0,0,-1,0; x)-12 \,\text{H}(0,0,0,0; x)-16 \,\text{H}(0,0,1,0; x)\nonumber \\ \allowdisplaybreaks[1]
                         &+12 \,\text{H}(0,1,-1,0; x)-8 \,\text{H}(0,1,0,0; x)-4 \,\text{H}(0,1,1,0; x)\nonumber \\ \allowdisplaybreaks[1]
                         &+8 \,\text{H}(1,0,-1,0; x)-8 \,\text{H}(1,0,0,0; x)-8 \,\text{H}(1,0,1,0; x)  \nonumber \\ \allowdisplaybreaks[1]
                         &+2 \,\zeta_2 (2 \,\text{H}(-1,0; x)+\,\text{H}(0,-1; x)+\,\text{H}(0,0; x)+\,\text{H}(0,1; x)\nonumber \\ \allowdisplaybreaks[1]
                         &-2 \,\text{H}(1,0; x)) + \zeta_3 (4 \,\text{H}(-1; x)-\,\text{H}(0; x)-4 \,\text{H}(1; x))
                         -\frac{37 \,\zeta _4}{4}\,,
\end{align}
\end{subequations}

\section{Master Integrals for the two-loop non-planar Box}
\label{results:xbox2L}
In this Appendix  we present the expression of the $12$
MI's of the  two-loop non-planar Box in Eq.~(\ref{Eq:MIfinBox}).  They are obtained according to the procedure
described in Section~\ref{Sec:xbox}  starting from the integrals $\top{i}$  in Fig.~\ref{Fig:Masters3}.  The latter
are normalized according to the integration measure (Minkowskian metric is understood)
\begin{displaymath}
\left ( \frac{s^{\epsilon}\,  \Gamma(1-2\epsilon)}{\Gamma(1+\epsilon) \Gamma(1-\epsilon)^2} \right )^2 \, \int \, \frac{d^D k_1}{\pi^{D/2}} \, \int \,   \frac{d^D k_2 }{\pi^{D/2}}\, .
\end{displaymath}
In the following we collect
the coefficients $g_{i}^{(a)}$ of the expansion of the MI's around $\epsilon=0$,
\begin{displaymath}
g_i = \sum_{a=0}^{4} \epsilon^a \ g_{i}^{(a)} \ , \qquad i=1,\ldots,12  \, .
\end{displaymath}
The MI's have uniform transcendentality, as can be explicitly checked by using the expressions of the coefficients $g_{i}^{(a)}$.
\begin{subequations}
\begin{align}
                   g_{1}^{(0)} = {}
                         &-1\,,\\ \allowdisplaybreaks[1]
                   g_{1}^{(1)} = {}
                         &-2 \,i\,\pi\,,\\ \allowdisplaybreaks[1]
                   g_{1}^{(2)} = {}
                         &12 \,\zeta_2\,,\\ \allowdisplaybreaks[1]
                   g_{1}^{(3)} = {}
                         &6 \,\zeta_3 + 8 \,i\,\pi \,\zeta_2\,,\\ \allowdisplaybreaks[1]
                   g_{1}^{(4)} = {}
                         &-51 \,\zeta_4 + 12 \,i\,\pi \,\zeta_3\,,
\end{align}
\end{subequations}

\begin{subequations}
\begin{align}
                   g_{2}^{(0)} = {}
                         &-1\,,\\ \allowdisplaybreaks[1]
                   g_{2}^{(1)} = {}
                         &2 \,\text{H}(0; x)\,,\\ \allowdisplaybreaks[1]
                   g_{2}^{(2)} = {}
                         &-4 \,\text{H}(0,0; x)\,,\\ \allowdisplaybreaks[1]
                   g_{2}^{(3)} = {}
                         &8 \,\text{H}(0,0,0; x) + 6 \,\zeta_3\,,\\ \allowdisplaybreaks[1]
                   g_{2}^{(4)} = {}
                         &-16 \,\text{H}(0,0,0,0; x) -12 \,\zeta_3 \,\text{H}(0; x) + 9 \,\zeta_4\,,
\end{align}
\end{subequations}

\begin{subequations}
\begin{align}
                   g_{3}^{(0)} = {}
                         &-1\,,\\ \allowdisplaybreaks[1]
                   g_{3}^{(1)} = {}
                         &-2 \,\text{H}(1; x)\,,\\ \allowdisplaybreaks[1]
                   g_{3}^{(2)} = {}
                         &-4 \,\text{H}(1,1; x)\,,\\ \allowdisplaybreaks[1]
                   g_{3}^{(3)} = {}
                         &-8 \,\text{H}(1,1,1; x) + 6 \,\zeta_3\,,\\ \allowdisplaybreaks[1]
                   g_{3}^{(4)} = {}
                         &-16 \,\text{H}(1,1,1,1; x) +12 \,\zeta_3 \,\text{H}(1; x) + 9 \,\zeta_4\,,
\end{align}
\end{subequations}

\begin{subequations}
\begin{align}
                   g_{4}^{(0)} = {}
                         &\frac{1}{4}\,,\\ \allowdisplaybreaks[1]
                   g_{4}^{(1)} = {}
                         &\frac{\,i\,\pi}{2}\,,\\ \allowdisplaybreaks[1]
                   g_{4}^{(2)} = {}
                         &-\frac{5 \,\zeta_2}{2}\,,\\ \allowdisplaybreaks[1]
                   g_{4}^{(3)} = {}
                         &-\,\zeta_3 - \,i\,\pi \,\zeta_2\,,\\ \allowdisplaybreaks[1]
                   g_{4}^{(4)} = {}
                         &-2 \,i\,\pi \,\zeta_3\,,
\end{align}
\end{subequations}

\begin{subequations}
\begin{align}
                   g_{5}^{(0)} = {}
                         &\frac{9}{4}\,,\\ \allowdisplaybreaks[1]
                   g_{5}^{(1)} = {}
                         &-3 \,\text{H}(0; x) + \frac{3 \,i\,\pi}{2}\,,\\ \allowdisplaybreaks[1]
                   g_{5}^{(2)} = {}
                         &3 \,\text{H}(0,0; x) - \frac{15 \,\zeta_2}{2} - 3 \,i\,\pi \,\text{H}(0; x)\,,\\ \allowdisplaybreaks[1]
                   g_{5}^{(3)} = {}
                         &3 \,\text{H}(1,0,0; x) + 15 \,\zeta_2 \,\text{H}(0; x) - 15 \,\zeta_3
                          \nonumber \\ \allowdisplaybreaks[1]
                         &+ \,i\,\pi (6 \,\text{H}(0,0; x)+3 \,\text{H}(1,0; x))\,,\\ \allowdisplaybreaks[1]
                   g_{5}^{(4)} = {}
                         &-12 \,\text{H}(0,0,0,0; x)-6 \,\text{H}(0,1,0,0; x)-12 \,\text{H}(1,0,0,0; x) \nonumber \\ \allowdisplaybreaks[1]
                         &-3 \,\text{H}(1,1,0,0; x) -15 \,\zeta_2 (2 \,\text{H}(0,0; x)+\,\text{H}(1,0; x))  \nonumber \\ \allowdisplaybreaks[1]
                         &+3 \,\zeta_3 (7 \,\text{H}(0; x)+\,\text{H}(1; x)) - \frac{183 \,\zeta_4}{4}  \nonumber \\ \allowdisplaybreaks[1]
                         &+i\,\pi \left( -12 \,\text{H}(0,0,0; x)-6 \,\text{H}(0,1,0; x)-6 \,\text{H}(1,0,0; x)\right. \nonumber \\ \allowdisplaybreaks[1]
                         &-3 \,\text{H}(1,1,0; x) -3 \,\zeta_2 \,\text{H}(1; x) - 9 \,\zeta_3\left . \right) \,,
\end{align}
\end{subequations}

\begin{subequations}
\begin{align}
                   g_{6}^{(0)} = {}
                         &\frac{9}{4}\,,\\ \allowdisplaybreaks[1]
                   g_{6}^{(1)} = {}
                         &3 \,\text{H}(1; x) +\frac{3 \,i\,\pi}{2}\,,\\ \allowdisplaybreaks[1]
                   g_{6}^{(2)} = {}
                         &3 \,\text{H}(1,1; x) - \frac{15 \,\zeta_2}{2} + 3 \,i\,\pi \,\text{H}(1; x)\,,\\ \allowdisplaybreaks[1]
                   g_{6}^{(3)} = {}
                         &-3 \,\text{H}(0,1,1; x) - 15 \,\zeta_2 \,\text{H}(1; x) - 12 \,\zeta_3  \nonumber \\ \allowdisplaybreaks[1]
                         &+i\,\pi \left(
                         3 \,\text{H}(0,1; x)+6 \,\text{H}(1,1; x) - 3 \,\zeta_2  \right) \,,\\ \allowdisplaybreaks[1]
                   g_{6}^{(4)} = {}
                         &-3 \,\text{H}(0,0,1,1; x)-12 \,\text{H}(0,1,1,1; x)-6 \,\text{H}(1,0,1,1; x)\nonumber \\ \allowdisplaybreaks[1]
                         &-12 \,\text{H}(1,1,1,1; x) - 15 \,\zeta_2 (\,\text{H}(0,1; x)+2 \,\text{H}(1,1; x)) \nonumber \\ \allowdisplaybreaks[1]
                         &-15 \,\zeta_3 \,\text{H}(1; x) -\frac{27 \,\zeta_4}{2}
                         + i\,\pi \left( 3 \,\text{H}(0,0,1; x)+6 \,\text{H}(0,1,1; x) \right. \nonumber \\ \allowdisplaybreaks[1]
                         &+6 \,\text{H}(1,0,1; x)+12 \,\text{H}(1,1,1; x) - 6 \,\zeta_2 \,\text{H}(1; x) - 6 \,\zeta_3 \left. \right) \,,
\end{align}
\end{subequations}

\begin{subequations}
\begin{align}
                   g_{7}^{(0)} = {}
                         & 0 \,,\\ \allowdisplaybreaks[1]
                   g_{7}^{(1)} = {}
                         & 0 \,,\\ \allowdisplaybreaks[1]
                   g_{7}^{(2)} = {}
                         &\,\text{H}(0,0; x) +\,i\,\pi \,\text{H}(0; x)\,,\\ \allowdisplaybreaks[1]
                   g_{7}^{(3)} = {}
                         &-4 \,\text{H}(0,0,0; x)-2 \,\text{H}(1,0,0; x) - 6 \,\zeta_2 \,\text{H}(0; x)  +2 \,\zeta_3  \nonumber \\ \allowdisplaybreaks[1]
                         &-i\,\pi \left( 2 (\,\text{H}(0,0; x)+\,\text{H}(1,0; x)) + 2 \,\zeta_2  \right) \,,\\ \allowdisplaybreaks[1]
                   g_{7}^{(4)} = {}
                         &12 \,\text{H}(0,0,0,0; x)+4 \,\text{H}(0,1,0,0; x)+8 \,\text{H}(1,0,0,0; x) \nonumber \\ \allowdisplaybreaks[1]
                         &+4 \,\text{H}(1,1,0,0; x)
                         +12 \,\zeta_2 (\,\text{H}(0,0; x)+\,\text{H}(1,0; x)) \nonumber \\ \allowdisplaybreaks[1]
                         &-4 \,\zeta_3 (\,\text{H}(0; x)+\,\text{H}(1; x)) +27 \,\zeta_4 \nonumber \\ \allowdisplaybreaks[1]
                         &+i\,\pi \left( 4 (\,\text{H}(0,0,0; x)+\,\text{H}(0,1,0; x)+\,\text{H}(1,0,0; x)\right.\nonumber \\ \allowdisplaybreaks[1]
                         &+\,\text{H}(1,1,0; x)) + 4 \,\zeta_2 \,\text{H}(1; x)\left . \right) \,,
\end{align}
\end{subequations}

\begin{subequations}
\begin{align}
                   g_{8}^{(0)} = {}
                         & 0 \,,\\ \allowdisplaybreaks[1]
                   g_{8}^{(1)} = {}
                         & 0 \,,\\ \allowdisplaybreaks[1]
                   g_{8}^{(2)} = {}
                         &\,\text{H}(0,0; x)+\,\text{H}(0,1; x)+\,\text{H}(1,0; x)+\,\text{H}(1,1; x) + 3 \,\zeta_2\,,\\ \allowdisplaybreaks[1]
                   g_{8}^{(3)} = {}
                         &-4 \,\text{H}(0,0,0; x)-2 \,\text{H}(0,0,1; x)-2 \,\text{H}(0,1,0; x)+2 \,\text{H}(1,0,1; x) \nonumber \\ \allowdisplaybreaks[1]
                         &+2 \,\text{H}(1,1,0; x)+4 \,\text{H}(1,1,1; x) -6 \,\zeta_2 (\,\text{H}(0; x)-\,\text{H}(1; x)) + 2 \,\zeta_3\,,\\ \allowdisplaybreaks[1]
                   g_{8}^{(4)} = {}
                         &12 \,\text{H}(0,0,0,0; x)+4 \,\text{H}(0,0,0,1; x)+4 \,\text{H}(0,0,1,0; x)\nonumber \\ \allowdisplaybreaks[1]
                         &-4 \,\text{H}(0,1,0,1; x)-4 \,\text{H}(0,1,1,0; x)-4 \,\text{H}(0,1,1,1; x)\nonumber \\ \allowdisplaybreaks[1]
                         &-4 \,\text{H}(1,0,0,0; x)-4 \,\text{H}(1,0,0,1; x)-4 \,\text{H}(1,0,1,0; x)\nonumber \\ \allowdisplaybreaks[1]
                         &+4 \,\text{H}(1,1,0,1; x)+4 \,\text{H}(1,1,1,0; x)+12 \,\text{H}(1,1,1,1; x) \nonumber \\ \allowdisplaybreaks[1]
                         &+12 \,\zeta_2 (\,\text{H}(0,0; x)-\,\text{H}(0,1; x)-\,\text{H}(1,0; x)+\,\text{H}(1,1; x))  \nonumber \\ \allowdisplaybreaks[1]
                         &-4 \,\zeta_3 (\,\text{H}(0; x)-\,\text{H}(1; x)) +12 \,\zeta_4\,,
\end{align}
\end{subequations}

\begin{subequations}
\begin{align}
                   g_{9}^{(0)} = {}
                         & 0 \,,\\ \allowdisplaybreaks[1]
                   g_{9}^{(1)} = {}
                         & 0 \,,\\ \allowdisplaybreaks[1]
                   g_{9}^{(2)} = {}
                         &\,\text{H}(1,1; x) -\,i\,\pi \,\text{H}(1; x)\,,\\ \allowdisplaybreaks[1]
                   g_{9}^{(3)} = {}
                         &2 \,\text{H}(0,1,1; x)+4 \,\text{H}(1,1,1; x) + 6 \,\zeta_2 \,\text{H}(1; x)  \nonumber \\ \allowdisplaybreaks[1]
                         &-2 \,i\,\pi (\,\text{H}(0,1; x)+\,\text{H}(1,1; x))\,,\\ \allowdisplaybreaks[1]
                   g_{9}^{(4)} = {}
                         &4 \,\text{H}(0,0,1,1; x)+8 \,\text{H}(0,1,1,1; x)+4 \,\text{H}(1,0,1,1; x)\nonumber \\ \allowdisplaybreaks[1]
                         &+12 \,\text{H}(1,1,1,1; x) +12 \,\zeta_2 (\,\text{H}(0,1; x)+\,\text{H}(1,1; x)) \nonumber \\ \allowdisplaybreaks[1]
                         &+i\,\pi \left( -4 (\,\text{H}(0,0,1; x)+\,\text{H}(0,1,1; x)+\,\text{H}(1,0,1; x)\right.\nonumber \\ \allowdisplaybreaks[1]
                         &+\,\text{H}(1,1,1; x)) + 4 \,\zeta_2 \,\text{H}(1; x)\left . \right) \,,
\end{align}
\end{subequations}

\begin{subequations}
\begin{align}
                   g_{10}^{(0)} = {}
                         &-1\,,\\ \allowdisplaybreaks[1]
                   g_{10}^{(1)} = {}
                         &-2 \,i\,\pi\,,\\ \allowdisplaybreaks[1]
                   g_{10}^{(2)} = {}
                         &17 \,\zeta_2\,,\\ \allowdisplaybreaks[1]
                   g_{10}^{(3)} = {}
                         &23 \,\zeta_3 + 18 \,i\,\pi \,\zeta_2\,,\\ \allowdisplaybreaks[1]
                   g_{10}^{(4)} = {}
                         &-\frac{317 \,\zeta_4}{2} + 46 \,i\,\pi \,\zeta_3\,,
\end{align}
\end{subequations}

\begin{subequations}
\begin{align}
                   g_{11}^{(0)} = {}
                         & 0 \,,\\ \allowdisplaybreaks[1]
                   g_{11}^{(1)} = {}
                         &\frac{5}{2} (\,\text{H}(0; x)+\,\text{H}(1; x))\,,\\ \allowdisplaybreaks[1]
                   g_{11}^{(2)} = {}
                         &5 \,i\,\pi (\,\text{H}(0; x)+\,\text{H}(1; x))\,,\\ \allowdisplaybreaks[1]
                   g_{11}^{(3)} = {}
                         &-10 \,\text{H}(0,0,0; x)-4 \,\text{H}(0,0,1; x)-4 \,\text{H}(0,1,0; x)-10 \,\text{H}(0,1,1; x) \nonumber \\ \allowdisplaybreaks[1]
                         &-10 \,\text{H}(1,0,0; x)-4 \,\text{H}(1,0,1; x)-4 \,\text{H}(1,1,0; x)-10 \,\text{H}(1,1,1; x)  \nonumber \\ \allowdisplaybreaks[1]
                         &-50 \,\zeta_2 (\,\text{H}(0; x)+\text{H}(1; x)) + 6 \,\zeta_3 + i\,\pi \left(
                         -6 (\,\text{H}(0,0; x)-\,\text{H}(0,1; x) \right. \nonumber \\ \allowdisplaybreaks[1]
                         &+\,\text{H}(1,0; x)-\,\text{H}(1,1; x)) - 6 \,\zeta_2\left . \right) \,,\\ \allowdisplaybreaks[1]
                   g_{11}^{(4)} = {}
                         &40 \,\text{H}(0,0,0,0; x)+16 \,\text{H}(0,0,0,1; x)+16 \,\text{H}(0,0,1,0; x)\nonumber \\ \allowdisplaybreaks[1]
                         &-8 \,\text{H}(0,0,1,1; x)+8 \,\text{H}(0,1,0,0; x)-16 \,\text{H}(0,1,0,1; x)\nonumber \\ \allowdisplaybreaks[1]
                         &-16 \,\text{H}(0,1,1,0; x)-40 \,\text{H}(0,1,1,1; x)+40 \,\text{H}(1,0,0,0; x)\nonumber \\ \allowdisplaybreaks[1]
                         &+16 \,\text{H}(1,0,0,1; x)+16 \,\text{H}(1,0,1,0; x)-8 \,\text{H}(1,0,1,1; x)\nonumber \\ \allowdisplaybreaks[1]
                         &+8 \,\text{H}(1,1,0,0; x)-16 \,\text{H}(1,1,0,1; x)-16 \,\text{H}(1,1,1,0; x)\nonumber \\ \allowdisplaybreaks[1]
                         &-40 \,\text{H}(1,1,1,1; x)
                         + 96 \,\zeta_2 (\,\text{H}(0,0; x)-\,\text{H}(0,1; x)+\,\text{H}(1,0; x)\nonumber \\ \allowdisplaybreaks[1]
                         &-\,\text{H}(1,1; x)) - 55 \,\zeta_3 (\,\text{H}(0; x)+\,\text{H}(1; x)) + 282 \,\zeta_4 \nonumber \\ \allowdisplaybreaks[1]
                         &+i\,\pi \left(
                         4 \,\text{H}(0,0,0; x)+16 \,\text{H}(0,0,1; x)+16 \,\text{H}(0,1,0; x) \right.  \nonumber \\ \allowdisplaybreaks[1]
                         &+4 \,\text{H}(0,1,1; x)
                         +4 \,\text{H}(1,0,0; x)+16 \,\text{H}(1,0,1; x)+16 \,\text{H}(1,1,0; x) \nonumber \\ \allowdisplaybreaks[1]
                         &+4 \,\text{H}(1,1,1; x) -36 \,\zeta_2 (\,\text{H}(0; x)+\,\text{H}(1; x)) + 12 \,\zeta_3\left . \right) \,,
\end{align}
\end{subequations}

\begin{subequations}
\begin{align}
                   g_{12}^{(0)} = {}
                         &-\frac{1}{4}\,,\\ \allowdisplaybreaks[1]
                   g_{12}^{(1)} = {}
                         &\frac{5}{4} \,\text{H}(0; x)+\frac{11}{4} \,\text{H}(1; x) - 2 \,i\,\pi\,,\\ \allowdisplaybreaks[1]
                   g_{12}^{(2)} = {}
                         &-4 \,\text{H}(0,0; x)-4 \,\text{H}(0,1; x)-4 \,\text{H}(1,0; x)-\,\text{H}(1,1; x)  +2 \,\zeta_2  \nonumber \\ \allowdisplaybreaks[1]
                         &+\frac{5}{2} \,i\,\pi (\,\text{H}(0; x)+\,\text{H}(1; x))\,,\\ \allowdisplaybreaks[1]
                   g_{12}^{(3)} = {}
                         &11 \,\text{H}(0,0,0; x)+2 \,\text{H}(0,0,1; x)+2 \,\text{H}(0,1,0; x)-10 \,\text{H}(0,1,1; x)\nonumber \\ \allowdisplaybreaks[1]
                         &+3 \,\text{H}(1,0,0; x)-6 \,\text{H}(1,0,1; x)-6 \,\text{H}(1,1,0; x)-15 \,\text{H}(1,1,1; x)  \nonumber \\ \allowdisplaybreaks[1]
                         &+ \zeta_2 (-13 \,\text{H}(0; x)-34 \,\text{H}(1; x)) +\frac{9 \,\zeta_3}{2} +i\,\pi \left(
                         \,\text{H}(0,0; x)+4 \,\text{H}(0,1; x)\right. \nonumber \\ \allowdisplaybreaks[1]
                         &+\,\text{H}(1,0; x)+7 \,\text{H}(1,1; x) + 21 \,\zeta_2\left. \right) \,,\\ \allowdisplaybreaks[1]
                   g_{12}^{(4)} = {}
                         &-28 \,\text{H}(0,0,0,0; x)+8 \,\text{H}(0,0,0,1; x)+8 \,\text{H}(0,0,1,0; x)\nonumber \\ \allowdisplaybreaks[1]
                         &-7 \,\text{H}(0,0,1,1; x)-24 \,\text{H}(0,1,0,0; x)-24 \,\text{H}(0,1,1,1; x)\nonumber \\ \allowdisplaybreaks[1]
                         &+4 \,\text{H}(1,0,0,0; x)+16 \,\text{H}(1,0,0,1; x)+16 \,\text{H}(1,0,1,0; x)\nonumber \\ \allowdisplaybreaks[1]
                         &-26 \,\text{H}(1,0,1,1; x)-8 \,\text{H}(1,1,0,0; x)-8 \,\text{H}(1,1,0,1; x)\nonumber \\ \allowdisplaybreaks[1]
                         &-8 \,\text{H}(1,1,1,0; x)-56 \,\text{H}(1,1,1,1; x) + \zeta_2 (20 \,\text{H}(0,0; x)\nonumber \\ \allowdisplaybreaks[1]
                         &-31 \,\text{H}(0,1; x)+44 \,\text{H}(1,0; x)-70 \,\text{H}(1,1; x)) \nonumber \\ \allowdisplaybreaks[1]
                         &+ \frac{1}{2} \,\zeta_3 (-15 \,\text{H}(0; x)-11 \,\text{H}(1; x)) + \frac{125 \,\zeta_4}{4} \nonumber \\ \allowdisplaybreaks[1]
                         &+i\,\pi \left( -14 \,\text{H}(0,0,0; x)+19 \,\text{H}(0,0,1; x)-20 \,\text{H}(0,1,0; x)\right. \nonumber \\ \allowdisplaybreaks[1]
                         &+4 \,\text{H}(0,1,1; x)-6 \,\text{H}(1,0,0; x)+30 \,\text{H}(1,0,1; x)-12 \,\text{H}(1,1,0; x)\nonumber \\ \allowdisplaybreaks[1]
                         &+18 \,\text{H}(1,1,1; x) - 6 \,\zeta_2 (9 \,\text{H}(0; x) +4 \,\text{H}(1; x)) + 78 \,\zeta_3\left . \right) \,,
\end{align}
\end{subequations}

\section{Magnus Theorem}
\label{App:Magnus}

We closely follow the discussion of Ref.\cite{PolySdeMany}.
Given  an operator, $\Omega$,
we define the
derivative of $\Omega^k$ w.r.t. $\Omega$ by its
action on a generic operator $H$:
\begin{eqnarray}
\bigg(
{d \over d\Omega} \Omega^k
\bigg) H \equiv H \Omega^{k-1} + \Omega H \Omega^{k-2} + \ldots +
\Omega^{k-1} H \, .
\label{Eq:DefDer}
\end{eqnarray}
This definition guarantees that,  when $\Omega=\Omega(x)$ and  $H=\partial_x\Omega$,
 \begin{eqnarray}
\partial_x \Omega^k    = \bigg( {d \over d\Omega} \Omega^k   \bigg) \,  \partial_x\Omega.
\label{Eq:ProDer}
\end{eqnarray}
The definition~(\ref{Eq:DefDer}) reduces to  $k H \Omega^{k-1}$ when $[\Omega, H] = 0$, therefore  it is
natural to write it as $k H \Omega^{k-1}$ plus correction terms involving (iterated) commutators. Using the
adjoint operator
\begin{equation}
\mbox{ad}_\Omega (H) \equiv [\Omega, H],
\end{equation}
and its iterated application  $\mbox{ad}^{\; i}_\Omega$ we obtain
\begin{eqnarray}
\bigg(
{d \over d\Omega} \Omega^2
\bigg) H  &=&  H\Omega + \Omega H = 2 H \Omega +\mbox{ad}_\Omega (H) \nonumber \\
\bigg(
{d \over d\Omega} \Omega^3
\bigg) H  &=&  H\Omega^2 + \Omega H \Omega  + \Omega^2 H = 3 H \Omega^2 + 3 [\Omega,H]\Omega + \mbox{ad}^{\, 2}_\Omega (H)
 \nonumber \\ [0.5ex]
\vdots \qquad    && \qquad \vdots \nonumber \\  [0.5ex]
\bigg(
{d \over d\Omega} \Omega^k
\bigg) H  &=& \sum_{i=0}^{k-1}
\begin{pmatrix}
k \\
i+1
\end{pmatrix}
 \mbox{ad}^{\, i}_\Omega (H) \ \Omega^{k-i-1}\, .
\end{eqnarray}
The last equation can be obtained by  induction using the relation
\begin{eqnarray}
\Omega \; \mbox{ad}^{\, i}_\Omega (H)    =   \mbox{ad}^{\, i}_\Omega (H) \;  \Omega +  \mbox{ad}^{\, i+1}_\Omega (H)
\end{eqnarray}
\bigskip

\noindent
The exponential of a matrix $\Omega$ is defined via a series expansion:
\begin{eqnarray}
e^\Omega \equiv \sum_{k \ge 0} {1 \over k!} \Omega^k \ .
\label{Eq:DefExp}
\end{eqnarray}
The derivative and the inverse of the exponential of a matrix can be straightforwardly obtained
by using the previous results:

\begin{lemma}[Derivative of the exponential]
The derivative of the matrix exponential can be derived from its
action on a generic operator $H$ and reads as follows
\label{lemma1}
\end{lemma}
\begin{eqnarray}
\bigg(
{d \over d\Omega} e^\Omega\bigg) H
= d \exp_\Omega(H)  \, e^\Omega \ ,  \qquad
d \exp_\Omega(H)
\equiv
\sum_{k \ge 0} {1 \over (k+1)!} \; \mbox{ad}^{\, i}_\Omega (H)
\end{eqnarray}

\begin{lemma}[Inverse of the exponential]
If the eigenvalues of $\mbox{ad}_\Omega$
are different from $2 \ell \pi i $ with $\ell \in \{ \pm 1, \pm 2, \ldots \}$, then $d \exp_\Omega$
is invertible, and
\label{lemma2}
\end{lemma}
\begin{eqnarray}
d \exp_\Omega^{-1}(H)
=
\sum_{k \ge 0} {\beta_k \over k!}  \; \mbox{ad}^{\, i}_\Omega (H)   \ ,
\end{eqnarray}
where $\beta_k$ are the Bernoulli numbers,  whose generating function is
\begin{eqnarray}
{t \over e^t-1} = \sum_{k=0}^\infty \ {\beta_k \over k!} \ t^k \ .
\end{eqnarray}
\bigskip

\noindent
We have now all the ingredients to prove the  following~\cite{Magnus}
\begin{theo}[Magnus]
The solution of  a generic linear matrix differential equation
\begin{eqnarray}
\partial_x Y = A(x) Y \ , \quad Y(x_0) = Y_0 \,
\end{eqnarray}
can be written as
\begin{eqnarray}
Y(x) = e^{\Omega(x,x_0)} \ Y(x_0) \equiv  e^{\Omega(x)} \ Y_0  \,
\label{Eq:solMagnus}
\end{eqnarray}
where $\Omega(x)$ can be computed by solving the differential equation,
\begin{eqnarray}
\partial_x \Omega = d \exp_\Omega^{-1} \Big( A(x) \Big)  \ ,
\qquad \Omega(x_0) = 0 \ .
\label{Eq:Thesis}
\end{eqnarray}
\label{theoMagnus}
\end{theo}
{\bf Proof} \quad  Let us
consider the derivative of~(\ref{Eq:solMagnus}). Using the definition~(\ref{Eq:DefExp})
and the property~(\ref{Eq:ProDer}) we have
\begin{eqnarray}
\partial_x Y
=
\bigg({d \over d \Omega} e^\Omega \bigg) \
\partial_x \Omega \ Y_0
=
d \exp_\Omega ( \partial_x  \Omega)  \ e^\Omega \ Y_0
=  d \exp_\Omega ( \partial_x  \Omega )  Y(x) \, . \nonumber
\end{eqnarray}
The {\it r.h.s.} is equal to $A(x) Y(x)$ when
\begin{eqnarray}
 d \exp_\Omega (  \partial_x   \Omega ) = A(x) \ .
\label{eq:dexpOmega}
\end{eqnarray}
The relation~(\ref{Eq:Thesis}) is thus proven by applying the operator $d
\exp_\Omega^{-1} $ to both sides of~(\ref{eq:dexpOmega}). \; \hfill $\Box$
\bigskip

\noindent
The differential equation for $\Omega$ explicitly reads,
\begin{eqnarray}
\partial_x \Omega = A(x) - {1\over 2} [\Omega, A(x)]
+ {1 \over 12} [\Omega,[\Omega, A(x)] ]
+ \ldots \ ,
\end{eqnarray}
and the solution can be written as a series, called  {\it Magnus expansion},
\begin{eqnarray}
\Omega=\sum_{n=1}^\infty \Omega_n(x) \ ,   \qquad
\Omega_n(x)  &=& \sum_{j=1}^{n-1} {\beta_j \over j !} \int_{x_0}^x
S_n^{(j)}(\tau) d\tau  \, .
\end{eqnarray}
The coefficients $\beta_j$ are the  Bernoulli numbers while the  integrands $S_n^{(j)}$ can be computed recursively,
\begin{eqnarray}
S_n^{(1)} &=&  [\Omega_{n-1},A] \ , \nonumber \\
S_n^{(j)} &=& \sum_{m=j-1}^{n-j} \,  \left [ \Omega_m, S_{n-m}^{(j-1)} \right ] \
\qquad 2 \le j \le n-2 \ , \nonumber \\
S_n^{(n-1)} &=& [\Omega_1, A] \ .
\end{eqnarray}

\section{One-Loop Bhabha scattering}
The calculation of the one-loop Bhabha scattering within  the DE's
method was discussed in \cite{BoboThesis,Bonciani:2003cj}, and more recently in
Ref.~\cite{Henn:2013woa}.  A selection of the Feynman diagrams contributing to this process
is depicted in Fig.~\ref{Fig:Diagrams1}.  In this section,  we  compute
a set of MI's with a slightly different definition from the ones in~\cite{Henn:2013woa},
which will be also adopted for the the one-loop $\times$ one-loop
subtopologies
of the QED vertices in the next section.
\FIGURE[t]{
\input{PICTURES/Bhabha/Diagrams}
\caption{Selection of Feynman diagrams entering the Bhabha scattering at one loop.}
\label{Fig:Diagrams1}
}

The diagrams  depend on the invariants
$s=(p_1+p_2)^2$, $t=(p_1+p_3)^2$,  $u=(p_2+p_3)^2$  and on the fermion
mass $m$.  Momentum  conservation and the on-shellness
of the external legs render these variables not independent as  they are related by the
condition $s+t+u = 4m^2$.   The integrals can be expressed in terms  of   the Landau
auxiliary variables $x$ and $y$,    defined as follows
\begin{equation}
s= -\frac{m^2 (1-x)^2}{x} \ , \qquad
t=-\frac{m^2 (1-y)^2}{y} \ .
\end{equation}
We identify the following basis $f$ of scalar integrals,
\begin{align}
f_1 &= \epsilon \top{1}\ ,  &
f_2 &= \epsilon \top{2}(t) \ ,   &
f_3 &= \epsilon \top{3}(s)  \ ,  \nonumber \\
f_4 &= \epsilon^2 \top{4}(t) \ ,  &
f_5 &= \epsilon^2 \top{5}(s,t)  \ ,
\end{align}
in terms of the integrals $\mathcal{T}$  in Fig.~\ref{Fig:Masters1}.
The basis $f$  fulfills the following systems of differential equations $(\sigma=x,y)$
\begin{eqnarray}
\partial_\sigma f(\epsilon,x,y) = A_\sigma(\epsilon, x,y) f(\epsilon,x,y) \ , \quad A_\sigma(\epsilon,x,y) = D_{\sigma,0}(x,y) + \epsilon A_{\sigma,1}(x,y)  \ .
\end{eqnarray}
Both systems are linear in $\epsilon$ and in both cases the $\mathcal{O}(\epsilon^0)$ term, $D_{\sigma,0}$,
is diagonal.  The systems can be brought in the canonical form by performing the
transformation
\begin{equation}
f(\epsilon,x,y) = B_0(x,y) g(\epsilon,x,y) \qquad
B_0(x,y)= e^{\int_{x_0}^x d \tau D_{ x,0}(\tau,y)}
 e^{\int_{y_0}^y d \tau D_{y,0}(x, \tau)} \ .
\end{equation}
The new basis $g$,
\begin{align}
g_1 &=  f_1  \ ,  &
g_2 &=  t \,f_2 \ , &
g_3 &=  \sqrt{\left(-s \right) \left(4 m^2-s\right)} \, f_3  \nonumber \\
g_4 &=  \sqrt{\left(-t \right) \left(4 m^2-t\right)} \, f_4  \ ,  &
g_5 &=  \sqrt{\left(-s \right) \left(4 m^2-s\right)} \, t \, f_5 \ . &
\end{align}
fulfills  the canonical systems
\begin{equation}
\partial_x g(\epsilon,x,y)=\epsilon \hat{A}_{x,1}(x,y)  \ g(\epsilon,x,y) \ , \qquad
\partial_y g(\epsilon, x,y)=\epsilon \hat{A}_{y,1}(x,y) \ g(\epsilon,x,y) \ ,
\label{eq:bhabha:cansys}
\end{equation}
with
\begin{align}
\hat{A}_{x,1}(x,y) &= \left(
\begin{array}{ccccc}
 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 \\
 \frac{1}{x} & 0 & \frac{1-x}{x (1+x)} & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 \\
 0 & \frac{2}{x} & \frac{2 (1-x) (1-y)^2}{(1+x) (x+y) (1+x y)} &
 -\frac{2 (1-y) (1+y)}{(x+y) (1+x y)} &
 \frac{(1-x) (1-y)^2}{(1+x) (x+y) (1+x y)} \\
\end{array}
\right)  \, ,  \nonumber \\[0.5ex] \allowdisplaybreaks[1]
 \hat{A}_{y,1}(x,y) &=\left(
\begin{array}{ccccc}
 0 & 0 & 0 & 0 & 0 \\
 0 & \frac{1+y}{(1-y) y} & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 \\
 \frac{1}{y} & \frac{1}{y} & 0 & \frac{4}{(1-y) (y+1)} & 0 \\
 0 & 0 & -\frac{2 x (1-y) (1+y)}{y (x+y) (1+x y)} &
-\frac{2 (1-x) (1+x)}{(x+y) (1+x y)} & \frac{(1+x)^2 (1+y)}{(1-y) (x+y) (1+x y)} \\
\end{array}
\right)  \, . \allowdisplaybreaks[1]
\end{align}
\FIGURE[t]{
\input{PICTURES/Bhabha/Masters}
\caption{MI's for the one-loop corrections to the Bhabha scattering. All the external momenta are incoming. A dot denotes a squared propagator.}
\label{Fig:Masters1}
}
The two systems of DE's in Eq.(\ref{eq:bhabha:cansys})
can be combined in a full differential form, along the  lines of Ref.~\cite{Henn:2013woa},
\begin{eqnarray}
d g(\epsilon, x,y)  = \epsilon \ d {\hat {\cal A}}_1(x,y) \ g(\epsilon, x,y) \ ,
\label{Eq:SystemBhabha}
\end{eqnarray}
where the matrix ${\hat{\cal A}}_1$ fulfills the relations,
\begin{eqnarray}
 {\partial_x {\hat {\cal A}}_1(x,y)} = {\hat A}_{x,1}(x,y) \ , \qquad
  {\partial_y {\hat {\cal A}}_1(x,y)} = {\hat A}_{y,1}(x,y) \ .
 \end{eqnarray}
 and the integrability condition
 \begin{eqnarray}
\epsilon \left (  \partial_x  \partial_y \hat {\cal A}_1(x,y)  -  \partial_y  \partial_x \hat {\cal A}_1(x,y) \right ) + \epsilon^2 \,
 \left [  \partial_x \hat {\cal A}_1(x,y) ,  \partial_y \hat {\cal A}_1(x,y)  \right ] =0\ .
 \end{eqnarray}
The matrix  $\hat{{\cal A}}_1$ is  logarithmic in the variables $x$ and  $y$,
\begin{eqnarray}
\hat{{\cal A}}_1(x,y)&=&
   M_1 \, \text{log}(x)
+M_2 \, \text{log}(1+x)
+ M_3 \, \text{log}(y)
+ M_4 \, \text{log}(1+y) + \nonumber \\
&&
+ M_5 \, \text{log}(1-y)
+ M_6 \, \text{log}(x+y)
+ M_7 \, \text{log}(1+xy) \ ,
\end{eqnarray}
with
\begin{align}
M_1 &=\left(
\begin{array}{ccccc}
 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 \\
 1 & 0 & 1 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 \\
 0 & 2 & 0 & 0 & 0 \\
\end{array}
\right) \ , &
M_2 &=\left(
\begin{array}{ccccc}
 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & -2 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & -4 & 0 & -2 \\
\end{array}
\right) \ ,  &
M_3 &=\left(
\begin{array}{ccccc}
 0 & 0 & 0 & 0 & 0 \\
 0 & 1 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 \\
 1 & 1 & 0 & 0 & 0 \\
 0 & 0 & -2 & 0 & 0 \\
\end{array}
\right) \ ,\nonumber \\[0.5ex] \allowdisplaybreaks[1]
M_4&=\left(
\begin{array}{ccccc}
 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 2 & 0 \\
 0 & 0 & 0 & 0 & 0 \\
\end{array}
\right)\ , &
M_5&=\left(
\begin{array}{ccccc}
 0 & 0 & 0 & 0 & 0 \\
 0 & -2 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & -2 & 0 \\
 0 & 0 & 0 & 0 & 2 \\
\end{array}
\right)\ , &
M_6&=\left(
\begin{array}{ccccc}
 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & -2 & -2 & -1 \\
\end{array}
\right)\ , \nonumber \\[0.5ex] \allowdisplaybreaks[1]
M_7&=\left(
\begin{array}{ccccc}
 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 2 & 2 & 1 \\
\end{array}
\right) \ .
\end{align}
The position of the non-zero  entries of the sparse matrices $M_i$ agrees with the result obtained in Ref.~\cite{Henn:2013woa}.
The actual value of the non-zero entries, however, are different, owing to the different normalization of the elements
of the basis of MI's. The solution of the system~(\ref{Eq:SystemBhabha}) can be computed along the lines
of Ref.~\cite{Henn:2013woa}. In particular, the solution is computed in the Euclidean region $0<x, y <1$
by using the analytic structures of the $g_i$ and then extended
in the physical region by analytic continuation~\cite{Bonciani:2003cj}.

\section{Conclusions}
In this article we elaborated on the method of differential
equations for Feynman integrals within the $D$-dimensional
regularization scheme.

The freedom in the choice of the MI's allowed us to
analyze the paradigmatic case of systems of
differential equations whose matrix is linear
in the dimensional parameter, $\epsilon = (4 - D)/2$.
We showed that these systems admit a {\it canonical} form,
where the dependence on $\epsilon$ is factorized from the
kinematic variables, as recently suggested by Henn.

We used  Magnus series  to  obtain the matrix
implementing the transformation from the linear to the canonical
form. The solution of the canonical system is obtained by using   either Dyson
series or Magnus series.  Both series require multiple integrations
which allow one to naturally express the MI's  in
terms of polylogarithms and of  their generalization.

We demonstrated that
the one-loop Bhabha scattering,
the two-loop electron form factors in QED and
the two-loop $2 \to 2$ massless scattering
exhibit a basis of MI's leading to linear systems of DE's.
We then obtained the corresponding canonical bases, in terms of
uniform transcendentality functions.

Finally, we have shown that our procedure can be extended to the more
general case of systems of DE's that are polynomial in $\epsilon$.

The  range of applicability  of the algorithm is rather  wide and  can be used to compute generic sets of MI's,
provided that the matrix associated to the system of DE's can be Taylor expanded  around $\epsilon =0$.

\section{Differential equations for Master Integrals}
\label{Sec:diffeq}
We consider  a linear system of first order differential equations
\begin{eqnarray}
\partial_x f(\epsilon,x) =  A(\epsilon, x) \ f(\epsilon,x)\ ,
\label{eq:MIdiffeq}
\end{eqnarray}
where  $f$ is a vector of MI's, while $x$ is a variable depending on kinematic invariants and masses.
We suppose that $A$ depends linearly on $\epsilon$,
\begin{eqnarray}
A(\epsilon,x) = A_0(x) + \epsilon A_1(x) \ ,
\end{eqnarray}
and we change the basis of MI's via the Magnus series  obtained by using   $A_0$
as kernel,
\begin{eqnarray}
f(\epsilon,x) = B_0(x) \  g(\epsilon,x) \ , \qquad  B_0(x) \equiv e^{ \magnus{A_0} (x,x_0) } \,  .
\label{Eq:change0}
\label{eq:MagnusA0}
\end{eqnarray}
Using Eq.~(\ref{eq:dexpOmega}), one can show that $B_0$ obeys the equation,
\begin{eqnarray}
\partial_x B_0(x) = A_0(x) B_0(x) \ ,
\label{eq:dotB}
\end{eqnarray}
which, analogously to the quantum-mechanical case, Eq.~(\ref{Eq:eqB}),
implies that the new basis $g$ of MI's fulfills a
system of differential equations  in the {\it canonical} factorized form,
\begin{eqnarray}
\partial_x g(\epsilon,x) = \epsilon {\hat A}_1(x) g(\epsilon,x) \,  .
\label{eq:cansys}
\end{eqnarray}
The matrix ${\hat A}_1$ is related to  $A_1$ by a similarity map,
\begin{eqnarray}
{\hat A}_1(x) = B_0^{-1}(x) A_1(x) B_0(x) \ ,
\label{eq:cantfm}
\end{eqnarray}
and does not depend on $\epsilon$.
The solution of Eq.~(\ref{eq:cansys}) can be found by using the Magnus theorem
with $\epsilon \hat A_1$ as kernel
\begin{eqnarray}
g(\epsilon,x) = B_1(\epsilon,x)  g_0(\epsilon) \ ,  \qquad  B_1(\epsilon,x)=e^{\magnus{\epsilon \hat A_1}(x,x_0)}  \, ,
\label{eq:MagnusA1}
\end{eqnarray}
where the vector $g_0$ corresponds to the boundary values of the MI's.
Therefore, the solution of the original system Eq.~(\ref{eq:MIdiffeq})
finally reads,
\begin{eqnarray}
f (\epsilon,x)   =B_0(x)   B_1(\epsilon,x)  g_0(\epsilon) \, .
\end{eqnarray}

\noindent
It is worth to notice
that $\magnus{\epsilon \hat A_1}$ in Eq.~(\ref{eq:MagnusA1}) depends on
$\epsilon$, while $\magnus{A_0}$ in Eq.~(\ref{eq:MagnusA0})
does not.

Let us remark that the previously described two-step procedure
is equivalent to solving, first, the {\it homogeneous} system
 \begin{equation}
 \partial_x f_{\mbox{\tiny H}}(\epsilon,x) = A_0(x) f_{\mbox{\tiny H}}(\epsilon,x) \, ,
\end{equation}
whose solution reads,
\begin{equation}
  f_{\mbox{\tiny H}}(\epsilon,x) = B_0(x) g(\epsilon) \ ,
 \end{equation}
 and,  then, to find the solution of the full system by Euler
 constants' variation. In fact, by promoting $g$ to be function of $x$,
 \begin{equation}
f_{\mbox{\tiny H}}(\epsilon,x) \rightarrow  f(\epsilon,x) =  B_0(x) g(\epsilon,x)  \ ,
\end{equation}
and by requiring $f$ to be solution of Eq.~(\ref{eq:MIdiffeq}), one
finds that $g(\epsilon,x)$ obeys the differential equation~(\ref{eq:cansys}).
\medskip

\noindent
The matrix $B_0$,  implementing the transformation from the
linear to the canonical form, is simply given as the
product of two matrix exponentials. Indeed one can split $A_0$
into a diagonal  term, $D_0$, and a matrix with vanishing diagonal entries $N_0$,
\begin{eqnarray}
A_0(x) = D_0(x) + N_0(x) \ .
\end{eqnarray}
The transformation $B$ is then  obtained by the composition of two   transformations
\begin{eqnarray}
B(x) =  e^{\magnus{D_0}(x,x_0)} e^{\magnus{\hat N_0}(x,x_0)}    =e^{\int_{x_0}^x d \tau  \ D_0(\tau)} e^{\magnus{\hat N_0}(x,x_0)} \ ,
\label{Eq:Bversion2}
\end{eqnarray}
where $\hat{N}_0$ is given by
\begin{eqnarray}
\hat N_0(x) =    e^{-  \int_{x_0}^x d \tau \ D_0(\tau)}  \ N_0(x) \ e^{\int_{x_0}^x d \tau \ D_0(\tau)}
\end{eqnarray}
In the last step of Eq.~(\ref{Eq:Bversion2}) we have used the  commutativity of the diagonal matrix $D_0$
with its own integral.  The leftmost expansion performs a transformation  that ``rotates" away $D_0$, while the
second expansion gets rid of the $\mathcal{O}(\epsilon^0)$ contribution coming from $\hat N_0$, i.e. coming from
the image of $N_0$ under the first transformation.
 \medskip

\noindent
In the examples hereby discussed it is possible, by trials and errors, to find a set of MI's obeying a system of DE's linear in $\epsilon$.
Moreover in these cases  one finds that $\magnus{\hat N_0}$
contains just the first term of the series, except for the non-planar box,
where also the second order is non vanishing.

\section{Introduction}
\label{sec:intro}

The {\it method of differential equations} (DE's), developed by
Kotikov, Remiddi and Gehrmann ~\cite{Kotikov:1990kg,Remiddi:1997ny,Gehrmann:1999as} and reviewed in Ref.~\cite{Argeri:2007up,Smirnov:2012gma}, is one of the most effective
techniques for computing dimensionally regulated multi-loop integrals and has led to significant achievements in the
context of multi-loop corrections. Within the continuous dimensional
regularization scheme, Feynman integrals can be related by using integration-by-parts identities
(IBP-id's)~\cite{Chetyrkin:1980pr,Tkachov:1981wb}, Lorentz invariance
identities~\cite{Gehrmann:1999as}, Gram
identities~\cite{Gluza:2010ws}, and quasi-Shouten
identities~\cite{Remiddi:2013joa}.   These relations can be exploited
in order to {\it identify} a set of independent integrals, dubbed {\it
  master integrals} (MI's),  that can be used as a basis of functions for the
virtual contributions to scattering amplitudes.

The MI's are functions of the kinematic invariants
constructed from the external momenta, of
the masses of the external  particles and  of the
particles running  in the loops, as well as of the number of spacetime dimensions.
Remarkably,  the existence of the aforementioned  relations forces the MI's
to obey linear systems of first-order differential equations in
the kinematic invariants,
which can be used for the determination of their expression.
When possible, these systems are solved exactly for generic values of the space-time dimension  $D$. Alternatively, they can be Laurent-expanded around suitable values of the dimensional parameter up to the required order,
obtaining a system of chained differential equations for the
coefficients of the expansions.
In the most general case, the latter are finally  integrated by using the method of  Euler's variation of
constants.

The nested structure of the Laurent expansion of the
linear system  leads to an iterative structure for the solution that, order-by-order in $\epsilon = (4-D)/2$, is
written in terms of repeated integrals, starting from the  kernels dictated by the homogeneous solution. The
transcendentality of the solution is  associated to the number of repeated integrations and increases by one unit as
the order of the $\epsilon$-expansion increases.
{\it The} solution of the system, namely the MI's, is finally determined by imposing the {\it
  boundary conditions} at special values of the kinematic
variables,  properly chosen either in correspondence of configurations that reduce the MI's to
 simpler integrals or in correspondence of pseudo-thresholds. In this
 latter case, the boundary conditions
 are obtained  by imposing  the {\it regularity} of the MI's  around  unphysical
singularities, ruling out divergent behavior of the general solution of the systems.

For any given scattering process the set of MI's is not unique, and, in practice,
their choice is rather arbitrary. Usually MI's are identified
after applying the Laporta reduction algorithm~\cite{Laporta:2001dd}.
Afterward, convenient manipulations of the basis of MI's may be performed.

Proper choices of MI's can simplify the form of the
systems of differential equations and, hence, of their solution, although general criteria for
determining such optimal sets are not available.
An important step in this direction has been recently taken in
Ref.~\cite{Henn:2013pwa}, where Henn proposes to solve the systems of DE's
for MI's with algebraic methods. The key observation
 is that a {\it good} choice of MI's allows one to cast  the system of DE's in a {\it canonical
form}, where the dependence on $\epsilon$, is factorized from the
kinematic.
The integration of a system in canonical form trivializes and the analytic properties of its general  solution are manifestly
inherited  from the matrix associated to the
system, which is  the kernel of the
representation of the solutions in terms of repeated integrations.

This novel idea has  been applied in a number of cases by Henn,
Smirnov, and Smirnov \cite{Henn:2013tua,Henn:2013woa,Henn:2013nsa},
showing the effectiveness of this approach.
As pointed out in \cite{Henn:2013pwa}, finding an algorithmic  procedure which,
starting from a generic set of MI's, leads to a set MI's fulfilling a
canonical system of DE's is a formidable task. In practice, the
quest for the suitable basis of MI's is determined by qualitative
properties required for the solution, such as finiteness in the $\epsilon \to 0$
limit, and homogeneous transcendentality, which turn into quantitative
tools like the unit leading singularity criterion and the {\it dlog}
representation in terms of Feynman parameters \cite{Gehrmann:2011xn}.

\bigskip
In this article, we suggest a convenient form  for the
initial system of MI's, and we  propose an algorithm to find
the transformation matrix yielding  a canonical system.
In particular, we choose a set of MI's obeying to a system of
DE's which has a {\it linear} $\epsilon$-dependence,
and we find a transformation which absorbs the $\mathcal{O}(\epsilon^0)$
term and leads to a new system of DE's where the $\epsilon$-dependence
is factorized. This transformation, as well as the integration of the
canonical system, are obtained by using  Magnus and Dyson series
expansions \cite{Magnus,Dyson:1949bp,Blanes:arXiv0810.5488}.
The  procedure we propose can be generalized to the case of
systems that are {\it polynomial} in $\epsilon$. Nevertheless, for the
cases at hand, we have succeeded to begin from a set of
MI's obeying a system that is linear  in $\epsilon$.
We show the effectiveness of our method by applying it to
non-trivial integrals. In particular,  we apply our procedure to  determine the MI's of
 the two-loop vertex diagrams contributing to the massive fermion
form-factors in QED~\cite{Bonciani:2003te,Bonciani:2003ai} and of  the non-planar two-loop diagrams
contributing to the $2 \to 2$ scattering of massless particles~\cite{Tausk:1999vh,Anastasiou:2000mf}.
Together with the ones in Ref.~\cite{Henn:2013woa}, the set of MI's for the two-loop QED vertices hereby presented
constitute a transcendentally-homogeneous subset for tackling the
analytic calculation of the still unknown non-planar two-loop box diagrams contributing
to the  massive Bhabha scattering in QED~\cite{Czakon:2004wm,Bonciani:2005im,Czakon:2006pa}.  It may enter  as well in  more general classes
of scattering processes involving massive particles.

Let us finally remark, that, while the canonical form of the system guaranties an easy integration procedure,
it alone does not directly imply the homogeneous transcendentality of the solution.
Indeed, this property may be affected by the analytic properties of additional
inputs such as the boundary conditions and the integrals which appear in the system of
DE's but cannot be determined from it. The latter are   integrals whose differential equation is
homogeneous and carries only the scaling information.

The  paper is structured as follows.
In section 2, we  present the Quantum Mechanical example that inspired this
study. The definition of Magnus series and its connection to the Dyson series are presented  in Section 3.
In Section 4 we show how to derive the canonical system starting from a
linear $\epsilon$-dependent system.
In Section 5, 6 and 7,  we apply our procedure to
the one-loop massive Bhabha scattering in QED,
to the two-loop vertex diagrams contributing to the massive electron
form factors in QED,
and to the two-loop non-planar box diagram respectively.
In Section 8, we show how our method generalizes to the case of a system
of differential equations which is polynomial in $\epsilon$.
The properties of the matrix exponential and
the proof of Magnus theorem are shown in Appendix~A,
while the Appendices~B, and~C collect the expressions of the MI's
of  the two-loop QED vertex diagrams and of  the two-loop non-planar
box diagram respectively.

We used the computer code {\sc Reduze2}~\cite{Studerus:2009ye,vonManteuffel:2012np} for the generation of the
systems of differential equations.

This manuscript is accompanied by two ancillary files, containing the
results of the canonical MI's for the two-loop QED vertices, and for
the two-loop non-planar box, respectively.

\section{Magnus series expansion}

Consider a generic linear matrix differential equation~\cite{Blanes:arXiv0810.5488}
\begin{eqnarray}
\partial_x Y(x) = A(x) Y(x) \ , \quad Y(x_0) = Y_0 \ .
\label{Eq:mainSY}
\end{eqnarray}
 If  $A(x)$ commutes with its integral $\int_{x_0}^x d\tau A(\tau)$, {\it e.g.}  in the scalar case,
 the solution can be written as
\begin{eqnarray}
Y(x) = e^{\int_{x_0}^x d\tau A(\tau)} \ Y_0 \ .
\label{eq:solMagn0}
\end{eqnarray}
In the general non-commutative case, one can use the {\it Magnus theorem}
\cite{Magnus} to write the solution as,
\begin{eqnarray}
Y(x) = e^{\Omega(x,x_0)} \ Y(x_0) \equiv  e^{\Omega(x)} \ Y_0 \,  ,
\label{eq:solMagn1}
\end{eqnarray}
where $\Omega(x)$ is written as a  series expansion,  called {\it Magnus expansion},
\begin{eqnarray}
\Omega(x)=\sum_{n=1}^\infty \Omega_n(x) \ .
\label{eq:Magnus}
\end{eqnarray}
The proof of the Magnus theorem is presented in the Appendix~\ref{App:Magnus}, together with
the actual expression of the  terms $\Omega_n$.
The first three terms of the expansion~(\ref{eq:Magnus}) read as follows:
\begin{eqnarray}
\Omega_1(x) &=& \int_{x_0}^x d\tau_1 A(\tau_1)\  ,\nonumber  \\
\Omega_2(x) &=&
{1 \over 2}
\int_{x_0}^x d\tau_1
\int_{x_0}^{\tau_1} d\tau_2  \
[A(\tau_1), A(\tau_2)] \ ,\nonumber  \\
\Omega_3(x) &=&
{1 \over 6}
\int_{x_0}^t d\tau_1
\int_{x_0}^{\tau_1} d\tau_2
\int_{x_0}^{\tau_2} d\tau_3  \
[A(\tau_1),[A(\tau_2),A(\tau_3)]] +
[A(\tau_3),[A(\tau_2),A(\tau_1)]]
\label{eq:Magnus:exp} \ . \qquad
\end{eqnarray}
We remark that  if  $A$  and its integral commute, the series~(\ref{eq:Magnus})
is truncated at the first order,  $\Omega = \Omega_1$, and we recover
the solution~(\ref{eq:solMagn0}).
As a notational aside, in the following we
will use the symbol
$\magnus{A}(x)$ to denote the Magnus expansion obtained using  $A$ as kernel.

\subsection{Magnus and Dyson series expansion}

Magnus series is related to the Dyson series~\cite{Blanes:arXiv0810.5488},
and their connection  can be obtained starting from the Dyson
expansion of the solution of the system~(\ref{Eq:mainSY}),
\begin{eqnarray}
Y(x) = Y_0 + \sum_{n=1}^\infty Y_n(x) \ , \quad Y_n(x) \equiv \int_{x_0}^x d\tau_1 \ldots \int_{x_0}^{\tau_{n-1}} d\tau_n \
A(\tau_1) A(\tau_2) \cdots A(\tau_n) \ ,
\label{eq:Dyson}
\label{eq:Dyson:exp}
\end{eqnarray}
in terms of the  {\it time-ordered} integrals $Y_n$.
Comparing Eq.~(\ref{eq:solMagn1}) and~(\ref{eq:Dyson}) we have
\begin{eqnarray}
\sum_{j=1}^\infty \Omega_j(x) = {\rm log}
\left( Y_0 + \sum_{n=1}^\infty Y_n(x) \right) \ ,
\end{eqnarray}
and the following relations
\begin{eqnarray}
Y_1 &=& \Omega_1 \ , \nonumber \\
Y_2 &=& \Omega_2 + {1 \over 2 !} \Omega_1^2 \ , \nonumber  \\
Y_3 &=& \Omega_3 + {1 \over 2 !}
(\Omega_1 \Omega_2 +
\Omega_2 \Omega_1)
+ {1 \over 3 !} \Omega_1^3 \ ,  \nonumber \\ [0.5ex]
\vdots & & \qquad \vdots \nonumber \\[0.5ex]
Y_n &=& \Omega_n + \sum_{j=2}^n {1 \over j} Q_n^{(j)} \, .
\end{eqnarray}
The matrices $Q_n^{(j)}$ are defined as
\begin{eqnarray}
Q_n^{(j)} = \sum_{m=1}^{n-j+1} Q_m^{(1)} Q_{n-m}^{(j-1)} \ ,
\quad
Q_n^{(1)} \equiv \Omega_n \ , \quad
Q_n^{(n)} \equiv \Omega_1^n \ .
\end{eqnarray}
In the following, we will use both Magnus and Dyson series.
The former allows us to easily demonstrate how a system of DE's,
whose matrix is linear in $\epsilon$, can be cast in the canonical form.
The latter can be more conveniently used for the explicit
representation of the solution.
\section*{Acknowledgments}
We wish to thank Roberto Bonciani for interesting discussions and comments on the manuscript.
The work of  P.M. and U.S. was supported by the Alexander von Humboldt
Foundation, in the framework of the Sofja Kovalevskaja Award Project ``Mathematical Methods for Particle Physics'', endowed by the German Federal Ministry of Education and Research.
The work of L.T. was supported in part by the Swiss National Science
Foundation (SNF) under  contract PDFMP2-135101. M.A.  wishes to acknowledge the kind hospitality of the Max Planck Institut f\"ur Physik
in Munich  during the completion of  this project. The Feynman diagrams depicted in this
paper are drawn using {\sc FeynArts}~\cite{Hahn:2000kx}.

\section{Polynomial $\epsilon$ dependence}
The cases discussed above admitted an initial choice of MI's $f$
obeying a system of differential equations linear in $\epsilon$.
We cannot be sure that this feature is general, and holds for any
scattering process in dimensional regularization.
Nevertheless, the use of Magnus series enables us to generalize our
algorithm to the case of systems of DE's whose matrix is a
{\it polynomial} in
$\epsilon$.
In fact, let us consider a system of equations where $A$ is of degree
$\kappa$ in $\epsilon$,
\begin{eqnarray}
\partial_x f(\epsilon,x) = A(\epsilon,x) \,  f(\epsilon,x)  \ , \qquad
A(\epsilon, x ) \equiv \sum_{k=0}^{\kappa} \epsilon^k { A}_k(x) \ .
\label{Eq:DiffEqPol}
\end{eqnarray}
By iterating the algorithm described in Section~\ref{Sec:diffeq},
the solution of the differential equation~(\ref{Eq:DiffEqPol}) can be expressed in terms of
a chain of products of Magnus exponentials,
\begin{eqnarray}
f(\epsilon,x) =
B_0(x)
B_1(\epsilon, x)
\cdots
B_\kappa(\epsilon, x)
\ef{\kappa}(\epsilon) \ ,  \qquad
B_k(\epsilon,x) \equiv e^{\magnus{\epsilon^k \hat A_k}(x,x_0)} \ ,
\label{eq:polycase:sol}
\end{eqnarray}
where the kernel  $\hat A_k$ is defined as
\begin{eqnarray}
 {\hat { A}}_k(\epsilon, x)  &=&  {\hat { A}}_k^{(k)}(\epsilon, x)  \, , \nonumber \\
 {\hat { A}}_k^{(j)}(\epsilon, x) &=&
B_{j-1}^{-1}(\epsilon,x)  \cdots
B_{1}^{-1}(\epsilon,x)
B_{0}^{-1}(x)
\; { A}_k(x)  \; B_{0}(x) B_{1}(\epsilon,x)
 \cdots
 B_{j-1}(\epsilon,x) \ .  \qquad
\label{eq:polycase:map}
\end{eqnarray}
It is worth to observe that, within our construction, the solution $f$ is given by
repeated transformations. Starting from
\begin{eqnarray}
f(\epsilon,x) =  B_0(x)\ef{0}(\epsilon,x) \, ,
\end{eqnarray}
we iteratively write
$\ef{k}$ as,
\begin{eqnarray}
f_{k}(\epsilon,x) = B_{k+1}(\epsilon,x)f_{k+1}(\epsilon,x) \ , \qquad (0 \le k \le \kappa-1) \, ,
\end{eqnarray}
which obeys the system
\begin{eqnarray}
\partial_x  \ef{k}(\epsilon, x) &=&
\epsilon^{k}
\bigg ( \,
\sum_{j=1}^{\kappa - k}
\epsilon^{j} {\hat { A}}_{k+j}^{(k+1)}(\epsilon,x)  \,
\bigg  )
\ \ef{k}(\epsilon, x)
 \ .
\end{eqnarray}
The generalization of the canonical system Eq.~(\ref{eq:cansys}) is obtained at the last
step of the iteration, when $k=\kappa-1$,
\begin{eqnarray}
\ef{\kappa-1}(\epsilon,x) = B_\kappa(\epsilon,x)\ef{\kappa}(\epsilon) \ , \qquad
\partial_x  \ef{\kappa-1}(\epsilon, x) &=& \epsilon^{\kappa}
{\hat { A}}_{\kappa}(\epsilon,x)
\ \ef{\kappa-1}(\epsilon, x) \ .
\end{eqnarray}

It is important to remark that the complete factorization of $\epsilon$ is
achieved only if $\kappa=1$, i.e. if the system is linear in $\epsilon$, because, although
${\hat { A}}_1$ is independent of $\epsilon$, ${\hat { A}}_k$ acquires a dependence on
$\epsilon$  for $k> 1$, {\it cfr.}   Eq.~(\ref{eq:polycase:map}). \\

The algorithm here described has a wide range of applicability and  can be used to compute generic sets of MI's,
provided that the matrix associated to the system of DE's can be Taylor expanded  around $\epsilon =0$. In this case,
the MI's are obtained perturbatively by  truncating the $\epsilon$
expansion of the  matrices associated to the systems of DE's.

\section{On time-dependent perturbation theory}
Given an Hamiltonian operator $H$, we consider the Schr\"odinger equation ($\partial_t \equiv \partial / \partial t$)
\begin{eqnarray}
i\hbar \ \partial_t | \Psi(t) \rangle =   H(t)   | \Psi(t) \rangle  \, .
\end{eqnarray}
Let us assume that ${H}$ can be split in two terms as
\begin{eqnarray}
{ H(t)}  = { H} _0(t) + \epsilon { H}_1(t) \ ,
\label{eq:Ham}
\end{eqnarray}
where  ${H} _0$ is a solvable Hamiltonian and $\epsilon \ll 1$ is a
small perturbation parameter. We may move to the  {\it interaction picture}
by performing a  transformation via a unitary  operator  $B$. In this
representation  any operator $A$ transforms according to
\begin{eqnarray}
A(t) = B(t) A_I(t) B^\dagger(t) \; .
\end{eqnarray}
In the interaction picture one imposes that only $H_1$ ($H_0$)  enters the time evolution of the states (of the operators),  thus $B$ is obtained by
imposing
\begin{eqnarray}
i \hbar \ \partial_t U_I(t)= \epsilon \, {H}_{1,I}(t) U_I(t)   +\left ( H_{0,I}(t) -i \hbar \,  B^\dagger(t) \ \partial_t B(t) \right ) U_I(t) \stackrel{!}{=}\epsilon \, {H}_{1,I}(t) U_I(t),  \;
\end{eqnarray}
so that  $B$ fulfills
\begin{eqnarray}
 i \hbar \,  \partial_t B(t) =H_{0}(t) B(t)\ .
 \label{Eq:eqB}
\end{eqnarray}
In the interaction picture the  Schr\"odinger equation  can be cast in a {\it canonical form},
\begin{eqnarray}
i\hbar \ \partial_t | \Psi_I(t) \rangle = \epsilon \,   H_{1,I}(t)   | \Psi_I(t) \rangle \, ,
\end{eqnarray}
where the $\epsilon$-dependence is factorized.
If the Hamiltonian $H_0$ at different times commute, the solution of Eq.~(\ref{Eq:eqB}) is
\begin{eqnarray}
B(t) = e^{- {i \over \hbar} \int_{t_0}^t d\tau H_0(\tau)}\ .
\end{eqnarray}

The important remark in this derivation is that, as a consequence of
the linear $\epsilon$-dependence of the original Hamiltonian Eq.~(\ref{eq:Ham}),
the states fulfill  an equation in a canonical form
by means of a transformation matrix $B$ that obeys the differential
equation~(\ref{Eq:eqB}). This simple quantum mechanical example contains the
two main guiding principles for building  canonical systems of differential equations for Feynman integrals:
\begin{itemize}
\item choose a set of Master Integrals obeying a system of
  differential equations linear in $\epsilon$;
\item find the transformation matrix by solving a differential
  equation governed by the constant term.
\end{itemize}
In this example $H_0(t)$ and $B(t)$ commute.
In the case of Feynman integrals, no assumption can be made on the
properties of the matrix associated to the systems of DE's built out of IBP-id's. Therefore,
in the following, we need to consider the generic case of non-commutative operators.

\section{Two-Loop QED Vertices}
\label{Sec:vertex}
A basis of MI's for the electron form factor at two loops in QED \cite{Bonciani:2003ai}
was computed in Ref.~\cite{Bonciani:2003te}, for arbitrary kinematics and finite electron mass.
The diagrams contributing to such corrections are depicted in Fig.~\ref{Fig:Diagrams2} and
depend on $s=(p_1+p_2)^2$ and $p_1^2=p_2^2=m^2$. In this example we
start from an alternative set of MI's,
\begin{align}
f_1 &= \epsilon^2 \top{1} \ , &
f_2 &= \epsilon^2   \top{2}  \ , &
f_3 &= \epsilon^2  \top{3}  \ , &
f_4 &= \epsilon^2 \top{4}  \ , &
f_5 &= \epsilon^2   \top{5}  \ , \nonumber \\ \displaybreak[1]
f_6 &= \epsilon^2  \top{6}   \ , &
f_7 &= \epsilon^2  \top{7}  \ , &
f_8 &= \epsilon^3  \top{8}  \ , &
f_9 &= \epsilon^3   \top{9}  \ , &
f_{10} &= \epsilon^2  \top{10}  \ ,  \nonumber \\\displaybreak[1]
f_{11} &= \epsilon^3  \top{11}  \ ,  &
f_{12} &= \epsilon^3  \top{12}  \ , &
f_{13} &= \epsilon^2    \top{13}  \ ,  &
f_{14} &= \epsilon^3  \top{14}   \ ,  &
f_{15} &= \epsilon^4   \top{15}  \ , \nonumber \\ \displaybreak[1]
f_{16} &= \epsilon^4  \top{16}   \ ,  &
f_{17} &=  \epsilon^4  \top{17}  \ ,  &
\end{align}
where the integrals $\top{i}$ are collected in Fig.~\ref{Fig:Masters2}.
The system of differential equation for $f$, in the auxiliary
variable $x$, defined through
 \begin{equation}
s= -\frac{m^2 (1-x)^2}{x} \, ,
\end{equation}
is linear in $\epsilon$,
\begin{eqnarray}
\partial_x f(\epsilon,x) = A(\epsilon, x) \  f(\epsilon,x) \ , \qquad A(\epsilon, x) = A_{0}(x) + \epsilon A_{1}(x)  \ .
\end{eqnarray}
The canonical form can be obtained performing the transformation described in Section~\ref{Sec:diffeq},
\begin{eqnarray}
f(\epsilon, x)  = B_0(x) \ g(\epsilon,x), \qquad B_0(x) = e^{\magnus{A_0}(x)  }  \ .
\end{eqnarray}
The new basis $g$ is given by
\begin{align}
g_1 &= f_1 \ , &   g_2 &=\lambda_1 f_2\ ,  \nonumber \\
 g_3 &=  (-s)\lambda_2 f_3 \ , & g_4 &=  m^2 f_4\ , \nonumber \\
 g_{5} &=\lambda_1  \left ( f_5 +\frac{f_6}{2} \right ) - \frac{s}{2} f_6  \ ,& g_6 &=   (-s) f_6 \ ,\nonumber \\
g_7 &=  m^2 f_7 \ ,&   g_8 &=  \lambda_1 f_8 \ ,\nonumber \\
g_9 &=   \lambda_1 f_9  \ ,&
g_{10} &= \lambda_3  \left ( 2  f_5 +f_6\right ) +  m^2 \lambda_2 f_{10} \ , \nonumber \\
g_{11} &=  \lambda_1 f_{11} \ , &  g_{12} &=  \lambda_1 f_{12} \ ,\nonumber \\
g_{13} &=  3\, \left( m^2-\frac{ s}{2}\right) f_7 -  s \lambda_2 f_{13} \ , & g_{14} &=   (-s) \lambda_2 f_{14}\ , \nonumber \\
 g_{15}&=  \lambda_1 f_{15} \ ,& g_{16} &=  \lambda_1 f_{16} \ ,\nonumber \\
  g_{17} &= (-s) \lambda_2 f_{17} \ ,
  \label{Eq:MIfinVertex}
\end{align}
where
\begin{eqnarray}
\lambda_1 =   \sqrt{-s}\sqrt{4m^2-s} \ ,   \qquad \lambda_2 =  (4m^2-s) \ ,   \qquad \lambda_3  = \frac{\lambda_1+\lambda_2}{4} \ .
\end{eqnarray}
The new basis of MI's obeys a system of DE's in the canonical form,
\FIGURE[t]{
\input{PICTURES/Vertex2/Masters}
\caption{MI's for the two-loop corrections to the QED vertex. All the external momenta depicted are incoming. In the integral  $\top{16}$
the loop momenta $k_1$, $k_2$ are fixed according to the first diagram of Fig.~\ref{Fig:Diagrams2} and a term
$(k_1+k_2)^2$ has to be included in  the numerator of the integrand. A dot indicates a squared propagator.}
\label{Fig:Masters2}
}
\begin{equation}
\partial_x g(\epsilon, x) =\epsilon \hat{A}_1(x) g(\epsilon, x) \ , \qquad \hat{A}_1(x)= \frac{M_1}{x}+ \frac{M_2}{1+x} + \frac{M_3}{1-x} \, ,
\label{Eq:SystemVertex}
\end{equation}
with
\begin{align}
M_1 &= \scriptsize \left(
\begin{array}{ccccccccccccccccc}
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
   \\
 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
   \\
 0 & 2 & 2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
   \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
   \\
 -1 & 0 & 0 & 0 & 5 & -6 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 &
   0 \\
 0 & 0 & 0 & 0 & 2 & -2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
   \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
   \\
 -1 & 0 & 0 & -4 & 0 & -2 & 0 & -2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
   & 0 \\
 0 & 0 & 0 & -2 & 0 & 0 & 0 & 0 & 2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
   \\
 -\frac{1}{2} & 0 & 0 & 0 & 1 & -2 & -3 & 0 & 0 & 3 & 3 & 0 & 0 & 0
   & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 1 & -1 & 2 & 0 & 0 & -2 & -2 & 0 & 0 & 0 & 0 & 0 &
   0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & -1 & -1 & 0 & 0 & 0 &
   0 \\
 0 & -1 & 0 & 0 & 0 & 0 & -3 & 0 & 0 & 0 & 0 & 3 & 3 & 0 & 0 & 0 &
   0 \\
 0 & -1 & 0 & 0 & 1 & -\frac{1}{2} & 0 & 2 & 2 & 0 & 0 & 0 & 0 & 2
   & 2 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & \frac{1}{2} & 0 & -\frac{1}{2} & 0 & 0 & 0 & 0
   & 0 & -1 & -1 & 0 & 0 \\
 -\frac{1}{2} & 0 & 0 & -2 & -1 & 0 & -2 & 1 & 0 & 2 & 0 & -2 & 0 &
   0 & -2 & -2 & 2 \\
 0 & 0 & 0 & 0 & -1 & \frac{1}{2} & 0 & 3 & -2 & 0 & -6 & -2 & 0 &
   0 & -4 & -4 & 4 \\
\end{array}
\right)  \normalsize \ , \nonumber \\[0.5ex] \displaybreak[1]
M_2 &=  \scriptsize \left(
\begin{array}{ccccccccccccccccc}
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
   \\
 0 & -2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
   \\
 0 & 0 & -4 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
   \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
   \\
 0 & 0 & 0 & 0 & -6 & 3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
   \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
   \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
   \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
   \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -4 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
   \\
 0 & 0 & 0 & 0 & -1 & \frac{1}{2} & 0 & 0 & 0 & -4 & 0 & 0 & 0 & 0
   & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 2 & 0 & 0 & 0 & 0 & 0 & 0
   \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 2 & 0 & 0 & 0 & 0 & 0
   \\
 0 & 0 & 0 & 0 & 0 & 0 & -6 & 0 & 0 & 0 & 0 & 0 & -2 & 0 & 0 & 0 &
   0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -4 & 0 & 0 & 0
   \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 2 & 0 & 0
   \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 2 & 0
   \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -4
   \\
\end{array}
\right) \ ,\normalsize  \nonumber \\  [0.5ex] \displaybreak[1]
M_3&=  \scriptsize \left(
\begin{array}{ccccccccccccccccc}
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
   \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
   \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
   \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
   \\
 0 & 0 & 0 & 0 & 2 & -2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
   \\
 0 & 0 & 0 & 0 & 0 & -2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
   \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
   \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & -2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
   \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
   \\
 0 & 0 & 0 & 0 & 0 & 0 & -6 & 0 & 0 & 2 & 0 & 0 & 0 & 0 & 0 & 0 & 0
   \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -2 & 0 & 0 & 0 & 0 & 0 & 0
   \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
   \\
 0 & 0 & 0 & 0 & 0 & 0 & -12 & 0 & 0 & 0 & 0 & 0 & 4 & 0 & 0 & 0 &
   0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
   \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
   \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 2 & 0 & 0 & 0 & -4 & 0 & 0 & -4 & -2 & 0
   \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 4
   \\
\end{array}
\right) \ . \normalsize \displaybreak[1]
\end{align}
The solution of the system can be expressed as Dyson series, as well
as Magnus series, in terms
of one-dimensional  Harmonic Polylogarithms  (HPL's)~\cite{Remiddi:1999ew}. The requirements
that the MI's are real-valued in the Euclidean region and regular
in $x=1$ ($s=0$), or simply the matching against the known integrals
at $x=1$, fix all but three boundary conditions,
corresponding to the {\it constant} MI's
$g_1$, $g_4$ and  $g_7$ (that do not depend on $x$).   The integrals $g_1$ and $g_4$ can be easily
computed by direct integration, while
$g_7$ can be determined from the results of Ref.~\cite{Argeri:2002wz}. Our results were checked analytically,
using the code {\sc HPL}~\cite{Maitre:2005uu, Maitre:2007kp},  against
the results available in the literature~\cite{Bonciani:2003te}.
The expressions of the transcendentally homogenous MI's $g$ are shown
in Appendix~\ref{results:vertex2L}, and collected in the ancillary file {\tt <vertex2L.m>}.

\section{Two-Loop non-planar Box}
\label{Sec:xbox}
The evaluation of the two-loop non-planar box diagram in Fig.~\ref{Fig:Diagrams3}, contributing
to the $2 \to 2$ scattering among massless particles, has already been
considered in the literature \cite{Tausk:1999vh,Anastasiou:2000mf}. Recently, for its planar
partner, a set of MI's with homogeneous transcendentality was presented
in Ref.~\cite{Henn:2013pwa}.
Our method can be easily applied to it, but instead of showing the case
of the ladder-box diagram, in this section, we compute the additional MI's
required for determining the non-planar contribution, having
expressions with manifest homogeneous transcendentality as well.

The integrals, in this case, are functions of the invariants
$s = (p_1+p_2)^2$,
$t = (p_1+p_3)^2$, and
$u = (p_2+p_3)^2$, with $p_i^2=0$, and $s+t+u=0$.

\noindent
We adopt the following initial choice of MI's,
\begin{align}
&
\begin{aligned}
f_1 &= \epsilon^2\, s\, \top{a}(s)\,,   &
f_2 &=  \epsilon^2\, t\, \top{a}(t)\,, &
f_3 &=  \epsilon^2\, u\, \top{a}(u)\,,   \nonumber \\\allowdisplaybreaks[1]
f_4 &=  \epsilon^3\, s\, \top{b}(s)\,,  &
f_5 &=  \epsilon^3\, s\, t\, \top{c}(s,t) \,, &
f_6 &=  \epsilon^3\, s\, u\, \top{c}(s,u) \,,   \nonumber \\\allowdisplaybreaks[1]
f_7 &=  \epsilon^4\, u\, \top{d}(s,t) \,, &
f_8 &= \epsilon^4\, s\, \top{d}(t,u) \,,  &
f_9 &=  \epsilon^4\, t\, \top{d}(u,s)  \,,  \nonumber \\\allowdisplaybreaks[1]
f_{10} &= \epsilon^4\, s^2\, \top{e}(s)\,, &&&&  \nonumber
\end{aligned}  \nonumber \\\allowdisplaybreaks[1]
f_{11} = {} &         \epsilon^4\, s\, t\, u\, \top{f}(s,t)
                -\frac{3}{4\, s\, (4 \epsilon+1)} \left[    \epsilon^2 \left(s^2\, \top{a}(s) + t^2\, \top{a}(t) + u^2\, \top{a}(u) \,\right) \right. \nonumber \\\allowdisplaybreaks[1]
      &          \left. -\,4 \epsilon^4 \left(u^2\, \top{d}(s,t) + s^2\, \top{d}(t,u) + t^2\, \top{d}(u,s) \right)    \right]       \, ,          \nonumber \\\allowdisplaybreaks[1]
 f_{12} = {} &      \epsilon^4\, s\, t\, \top{g}(s,t)
                -\frac{3}{8\, u\, (4 \epsilon+1)} \left[    \epsilon^2 \left(s^2\, \top{a}(s) + t^2\, \top{a}(t) + u^2\, \top{a}(u) \,\right) \right. \nonumber \\\allowdisplaybreaks[1]
      &          \left. -\,4 \epsilon^4 \left(u^2\, \top{d}(s,t) + s^2\, \top{d}(t,u) + t^2\, \top{d}(u,s) \right)    \right] \, ,
\end{align}
where the integrals $\top{i}$ correspond to the diagrams in
Fig.~\ref{Fig:Masters3}.
We notice that the integrals $f_1, \ldots, f_9$  are common to the two-loop planar box
diagram~\cite{Henn:2013pwa}.
The set $f$ of MI's obeys a system of differential equations the variable $x$, defined as,
\begin{eqnarray}
x = - {t\over s} \ ,
\end{eqnarray}
which is linear in $\epsilon$.
According to the procedure in Section~\ref{Sec:diffeq}, we can
build the matrix $B_0(x)$ ruling the change of basis $f(\epsilon,x) =
B_0(x) g(\epsilon,x)$, so that the new MI's,
\begin{align}
g_i = {} & f_i \ ,  \qquad  i=1,\ldots, 10\,  ,  \nonumber \\\allowdisplaybreaks[1]
g_{11}  = {} & \frac{s}{8\, t\, u}   \,  \big[ 3 f_1 (3\, t-5\, u)-3 f_2 (t+4\, u)+3 f_3 (2\,t+u)-16 f_5 \,u+8 f_6 \,t \nonumber \\\allowdisplaybreaks[1]
       &  -60 f_ 7 \,u-12 f_ 8 (t-u)+36 f_ 9 \,t-8 f_{11}\, u-8 f_{12}\, u \big] \ , \nonumber \\\allowdisplaybreaks[1]
g_{12} = {} & \frac{s}{8\, u} \,  \left( 9 f_1-3 f_2+6 f_3+8 f_6-12 f_8+36 f_9 \right) +f_{12} \, ,
\label{Eq:MIfinBox}
\end{align}
obey the canonical system,
\begin{equation}
\partial_x g(\epsilon,x)=\epsilon \hat{A}_1(x) g(\epsilon,x) \ , \qquad \hat{A}(x)= \frac{M_1}{x}+ \frac{M_2}{1-x} \, ,
\end{equation}
with
\FIGURE[t]{
\input{PICTURES/Box2C/Masters}
\caption{MI's for the two-loop diagram in Fig.~\ref{Fig:Diagrams3}. All the external
momenta depicted are incoming. In the last integral the loop momenta have to be
fixed according to Fig.~\ref{Fig:Diagrams3} and a term
$(k_2+p_1)^2$ enters the numerator of its integrand. A dot indicates a squared propagator.}
\label{Fig:Masters3}
}
\begin{align}
M_1 &=  \small
\left(
\begin{array}{cccccccccccc}
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 0 & -2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 0 & -\frac{3}{2} & 0 & 0 & -2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & \frac{3}{2} & -3 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\
 -\frac{1}{2} & \frac{1}{2} & 0 & 0 & 0 & 0 & -2 & 0 & 0 & 0 & 0 & 0 \\
 0 & \frac{1}{2} & -\frac{1}{2} & 0 & 0 & 0 & 0 & -2 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 2 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 -6 & -6 & -\frac{9}{2} & 0 & -4 & -2 & -18 & -12 & -12 & 1 & 1 & -2 \\
 \frac{3}{4} & \frac{9}{4} & -\frac{21}{4} & 3 & 2 & -3 & 12 & -6 & -18 & 0 & 0 & -2 \\
\end{array}
\right)
\ ,  \normalsize\nonumber \\  \allowdisplaybreaks[1]
M_2  &=  \small  \left(
\begin{array}{cccccccccccc}
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 2 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 0 & -\frac{3}{2} & 0 & 3 & -1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & \frac{3}{2} & 0 & 0 & 2 & 0 & 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & -2 & 0 & 0 & 0 & 0 & 0 \\
 0 & \frac{1}{2} & -\frac{1}{2} & 0 & 0 & 0 & 0 & 2 & 0 & 0 & 0 & 0 \\
 \frac{1}{2} & 0 & -\frac{1}{2} & 0 & 0 & 0 & 0 & 0 & 2 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 -6 & -6 & -\frac{9}{2} & 0 & -4 & -2 & -18 & -12 & -12 & 1 & 1 & -2 \\
 -\frac{21}{4} & \frac{9}{4} & -\frac{27}{4} & -6 & 2 & -4 & 12 & -6 & -24 & 1 & -1 & 0 \\
\end{array}
\right) \, \normalsize.%
\end{align}
The solution of the system can be expressed as Dyson series, as well as
Magnus series, in terms of one-dimensional HPL's~\cite{Remiddi:1999ew}. All MI's have been
computed in the scattering kinematics, i.e. $s>0$, $t<0$, $u<0$ with $|s|>|t|$, which gives $0 < x < 1$. As long as the planar sub
topologies are concerned, one can fix the boundary conditions using the
regularity properties of the integrals in some special kinematical points.
On the other hand, the analyticity structure of the crossed box is more complicated, since it involves at the same time
cuts in all three Mandelstam variables $s$, $t$, $u$.
 Nevertheless, in this particular case, the boundaries can be fixed
by direct comparison with the results presented in~\cite{Tausk:1999vh,Anastasiou:2000mf}.
The expressions of the transcendentally homogeneous MI's $g$ are shown
in Appendix~\ref{results:xbox2L}, and collected in the ancillary file {\tt <xbox2L.m>}.

redpajama_set_name: RedPajamaArXiv


September 2015 – ADVOCATES FOR FINE ART ENAMELING INC.
I have always been fascinated by colors. I almost became a chemist because I was fascinated by a solution of copper nitrate. The blue was hypnotizing. My first attempt at a dissertation in psychology was a study of the emotionally stimulating capacities of colors. I was first attracted to enameling by the deep, limpid quality of transparent enamels. If I have tried to make myself an artist, it was, in part, an attempt to capture and own these colors which have fascinated me for so long.
Thus, I was a bit disappointed by the last Enamelist Society exhibition, Alchemy 3. There was no doubt that the art work was superior and original, in some cases fiercely original. But where was the color? Certainly, there was color here and there. I found the color I so desired to look at in the bowls of Judy Stone and Fay Rooke, as well as the jewelry done jointly by Sharon Kree and Diane Buettner, and that of Toni Strassler. The wall art of Gabrielle Castonguay used color in an inventive way; Averill Shepps is doing interesting color experiments, and Cynthia Miller's abstraction and Terry Gay Puckett's portrait caught my color sensitive eye.
The color, however, was muted in so many other pieces. Herb Friedson, known for his carnival color extravaganzas, was represented by a muted piece. There was a lot of gray around. There was certainly not much celebration of the particular qualities of enamel colors. I have catalogs from recent , similar enamel exhibitions in Europe. They show a plethora of enamel pieces with the vibrant and pellucid colors which we expect from enamel art, so this dimming and graying of enamel art is not a world wide phenomena.
Hal Nelson put this wonderful enamel on facebook. This is enamel color! This is what I like to see.


text:
Man pleads guilty to being a felon in possession of gun
A New Hampshire man has pleaded guilty in federal court to being a felon in possession of a firearm and ammunition
CONCORD, N.H. (AP) � A New Hampshire man has pleaded guilty in federal court to being a felon in possession of a firearm and ammunition, prosecutors said Wednesday.
Court documents said David Weekly, 23, of Concord, was in a vehicle stopped by police in Manchester in January.
Weekly admitted to an officer he had a gun in his waistband and a loaded extended 31-round magazine. The firearm was reported stolen in 2018.
Weekly is scheduled to be sentenced on March 26.
The case was investigated by the Bureau of Alcohol, Tobacco, Firearms and Explosives and the Manchester Police Department.

redpajama_set_name: RedPajamaCommonCrawl


text:
Q: Using HTML in CSS After Tag for Hover State I am currently trying to pass in some html into the the :after pseudo class in css to add an arrow after a link once hovered on.  I have seen this done quite some time ago in the stock WordPress theme that used to ship with the default install.  I feel like I must me missing something simple.
a:hover:after{content: "»"}

thank in advance!

A: You need to use the hex value \00BB.
http://jsfiddle.net/durilai/AGDbH/

A: You're looking for, I think:
a:hover:after
{
content: "\00bb";
}

JS Fiddle demo
Link to reference: http://css-tricks.com/snippets/html/glyphs/

redpajama_set_name: RedPajamaStackExchange